Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2020-2021 of Michele Invernizzi.
3 :
4 : This file is part of the OPES plumed module.
5 :
6 : The OPES plumed module is free software: you can redistribute it and/or modify
7 : it under the terms of the GNU Lesser General Public License as published by
8 : the Free Software Foundation, either version 3 of the License, or
9 : (at your option) any later version.
10 :
11 : The OPES plumed module is distributed in the hope that it will be useful,
12 : but WITHOUT ANY WARRANTY; without even the implied warranty of
13 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 : GNU Lesser General Public License for more details.
15 :
16 : You should have received a copy of the GNU Lesser General Public License
17 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
18 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
19 : #include "ExpansionCVs.h"
20 : #include "core/ActionRegister.h"
21 : #include "tools/File.h"
22 :
23 : namespace PLMD {
24 : namespace opes {
25 :
26 : //+PLUMEDOC OPES_EXPANSION_CV ECV_UMBRELLAS_FILE
27 : /*
28 : Target a multiumbrella ensemble, by combining systems each with a parabolic bias potential at a different location.
29 :
30 : Any set of collective variables \f$\mathbf{s}\f$ can be used as ARG.
31 : The positions \f$\mathbf{s}_i\f$ and dimension \f$\mathbf{\sigma}_i\f$ of the umbrellas are read from file.
32 : \f[
33 : \Delta u_{\mathbf{s}_i}(\mathbf{s})=\sum_j^{\text{dim}}\frac{([s]_j-[s_i]_j)^2}{2[\sigma_i]_j^2}\, .
34 : \f]
35 : Notice that \f$\mathbf{\sigma}_i\f$ is diagonal, thus only one SIGMA per CV has to be specified for each umbrella.
36 : You can choose the umbrellas manually, or place them on a grid, or along a path, similar to \ref PATH.
37 : They must cover all the CV space that one wishes to sample.
38 :
39 : The first column of the umbrellas file is always ignored and must be called "time".
40 : You can also use as input file a STATE file from an earlier \ref OPES_METAD run (or an \ref OPES_MEAD_EXPLORE run, if you combine it with other ECVs).
41 :
42 : Similarly to \ref ECV_UMBRELLAS_LINE, you should set the flag ADD_P0 if you think your umbrellas might not properly cover all the CV region relevant for the unbiased distribution.
43 : You can also use BARRIER to set the maximum barrier height to be explored, and avoid huge biases at the beginning of your simulation.
44 : See also Appendix B of Ref.\cite Invernizzi2020unified for more details on these last two options.
45 :
46 : \par Examples
47 :
48 : \plumedfile
49 : cv1: DISTANCE ATOMS=1,2
50 : cv2: DISTANCE ATOMS=3,4
51 : cv3: DISTANCE ATOMS=4,1
52 : ecv: ECV_UMBRELLAS_FILE ARG=cv1,cv2,cv3 FILE=Umbrellas.data ADD_P0 BARRIER=70
53 : opes: OPES_EXPANDED ARG=ecv.* PACE=500
54 : PRINT FILE=COLVAR STRIDE=500 ARG=cv1,cv2,cv3,opes.bias
55 : \endplumedfile
56 :
57 : The umbrellas file might look like this:
58 : \auxfile{Umbrellas.data}
59 : #! FIELDS time cv1 cv2 cv3 sigma_cv1 sigma_cv2 sigma_cv3
60 : 1 1.17958 2.93697 1.06109 0.19707 0.28275 0.32427
61 : 2 2.04023 2.69714 1.84770 0.22307 0.25933 0.31783
62 : 3 1.99693 1.10299 1.13351 0.19517 0.26260 0.37427
63 : 4 1.15954 1.37447 2.25975 0.20096 0.27168 0.33353
64 : 5 1.10126 2.45936 2.40260 0.19747 0.24215 0.35523
65 : \endauxfile
66 :
67 : */
68 : //+ENDPLUMEDOC
69 :
70 : class ECVumbrellasFile :
71 : public ExpansionCVs
72 : {
73 : private:
74 : double barrier_;
75 : unsigned P0_contribution_;
76 : bool lower_only_;
77 :
78 : std::vector< std::vector<double> > centers_;
79 : std::vector< std::vector<double> > sigmas_;
80 :
81 : std::vector< std::vector<double> > ECVs_;
82 : std::vector< std::vector<double> > derECVs_;
83 : void initECVs();
84 :
85 : public:
86 : explicit ECVumbrellasFile(const ActionOptions&);
87 : static void registerKeywords(Keywords& keys);
88 : void calculateECVs(const double *) override;
89 : const double * getPntrToECVs(unsigned) override;
90 : const double * getPntrToDerECVs(unsigned) override;
91 : std::vector<std::string> getLambdas() const override;
92 : void initECVs_observ(const std::vector<double>&,const unsigned,const unsigned) override;
93 : void initECVs_restart(const std::vector<std::string>&) override;
94 : };
95 :
96 : PLUMED_REGISTER_ACTION(ECVumbrellasFile,"ECV_UMBRELLAS_FILE")
97 :
98 5 : void ECVumbrellasFile::registerKeywords(Keywords& keys)
99 : {
100 5 : ExpansionCVs::registerKeywords(keys);
101 10 : keys.add("compulsory","FILE","the name of the file containing the umbrellas");
102 10 : keys.add("optional","BARRIER","a guess of the free energy barrier to be overcome (better to stay higher than lower)");
103 10 : keys.addFlag("ADD_P0",false,"add the unbiased Boltzmann distribution to the target distribution, to make sure to sample it");
104 10 : keys.addFlag("LOWER_HALF_ONLY",false,"use only the lower half of each umbrella potentials");
105 5 : }
106 :
107 3 : ECVumbrellasFile::ECVumbrellasFile(const ActionOptions&ao):
108 : Action(ao),
109 3 : ExpansionCVs(ao)
110 : {
111 : //get number of CVs
112 : const unsigned ncv=getNumberOfArguments();
113 3 : centers_.resize(ncv);
114 3 : sigmas_.resize(ncv);
115 :
116 : //set P0_contribution_
117 3 : bool add_P0=false;
118 3 : parseFlag("ADD_P0",add_P0);
119 3 : if(add_P0)
120 2 : P0_contribution_=1;
121 : else
122 1 : P0_contribution_=0;
123 :
124 : //set barrier_
125 3 : barrier_=std::numeric_limits<double>::infinity();
126 3 : parse("BARRIER",barrier_);
127 6 : parseFlag("LOWER_HALF_ONLY",lower_only_);
128 :
129 : //set umbrellas
130 : std::string umbrellasFileName;
131 3 : parse("FILE",umbrellasFileName);
132 3 : IFile ifile;
133 3 : ifile.link(*this);
134 3 : if(ifile.FileExist(umbrellasFileName))
135 : {
136 3 : log.printf(" reading from FILE '%s'\n",umbrellasFileName.c_str());
137 3 : ifile.open(umbrellasFileName);
138 3 : ifile.allowIgnoredFields();
139 : double time; //first field is ignored
140 1332 : while(ifile.scanField("time",time))
141 : {
142 1989 : for(unsigned j=0; j<ncv; j++)
143 : {
144 : double centers_j;
145 1326 : ifile.scanField(getPntrToArgument(j)->getName(),centers_j);
146 1326 : centers_[j].push_back(centers_j); //this might be slow
147 : }
148 1989 : for(unsigned j=0; j<ncv; j++)
149 : {
150 : double sigmas_j;
151 2652 : ifile.scanField("sigma_"+getPntrToArgument(j)->getName(),sigmas_j);
152 1326 : sigmas_[j].push_back(sigmas_j);
153 : }
154 663 : ifile.scanField();
155 : }
156 : }
157 : else
158 0 : plumed_merror("Umbrellas FILE '"+umbrellasFileName+"' not found");
159 :
160 3 : checkRead();
161 :
162 : //extra consistency checks
163 3 : const unsigned sizeUmbrellas=centers_[0].size();
164 9 : for(unsigned j=0; j<ncv; j++)
165 : {
166 6 : plumed_massert(centers_[j].size()==sizeUmbrellas,"mismatch in the number of centers read from file");
167 6 : plumed_massert(sigmas_[j].size()==sizeUmbrellas,"mismatch in the number of sigmas read from file");
168 : }
169 :
170 : //set ECVs stuff
171 3 : totNumECVs_=sizeUmbrellas+P0_contribution_;
172 3 : ECVs_.resize(ncv,std::vector<double>(totNumECVs_));
173 3 : derECVs_.resize(ncv,std::vector<double>(totNumECVs_));
174 :
175 : //printing some info
176 3 : log.printf(" total number of umbrellas = %u\n",sizeUmbrellas);
177 3 : if(barrier_!=std::numeric_limits<double>::infinity())
178 1 : log.printf(" guess for free energy BARRIER = %g\n",barrier_);
179 3 : if(P0_contribution_==1)
180 2 : log.printf(" -- ADD_P0: the target includes also the unbiased probability itself\n");
181 3 : if(lower_only_)
182 1 : log.printf(" -- LOWER_HALF_ONLY: the ECVs are set to zero for values of the CV above the respective center\n");
183 6 : }
184 :
185 131 : void ECVumbrellasFile::calculateECVs(const double * cv)
186 : {
187 131 : if(lower_only_)
188 : {
189 153 : for(unsigned j=0; j<getNumberOfArguments(); j++)
190 : {
191 22644 : for(unsigned k=P0_contribution_; k<totNumECVs_; k++) //if ADD_P0, the first ECVs=0
192 : {
193 22542 : const unsigned kk=k-P0_contribution_;
194 22542 : const double dist_jk=difference(j,centers_[j][kk],cv[j])/sigmas_[j][kk]; //PBC might be present
195 22542 : if(dist_jk>=0)
196 : {
197 11483 : ECVs_[j][k]=0;
198 11483 : derECVs_[j][k]=0;
199 : }
200 : else
201 : {
202 11059 : ECVs_[j][k]=0.5*std::pow(dist_jk,2);
203 11059 : derECVs_[j][k]=dist_jk/sigmas_[j][kk];
204 : }
205 : }
206 : }
207 : }
208 : else
209 : {
210 240 : for(unsigned j=0; j<getNumberOfArguments(); j++)
211 : {
212 35520 : for(unsigned k=P0_contribution_; k<totNumECVs_; k++) //if ADD_P0, the first ECVs=0
213 : {
214 35360 : const unsigned kk=k-P0_contribution_;
215 35360 : const double dist_jk=difference(j,centers_[j][kk],cv[j])/sigmas_[j][kk]; //PBC might be present
216 35360 : ECVs_[j][k]=0.5*std::pow(dist_jk,2);
217 35360 : derECVs_[j][k]=dist_jk/sigmas_[j][kk];
218 : }
219 : }
220 : }
221 131 : }
222 :
223 6 : const double * ECVumbrellasFile::getPntrToECVs(unsigned j)
224 : {
225 6 : plumed_massert(isReady_,"cannot access ECVs before initialization");
226 6 : plumed_massert(j<getNumberOfArguments(),getName()+" has fewer CVs");
227 6 : return &ECVs_[j][0];
228 : }
229 :
230 6 : const double * ECVumbrellasFile::getPntrToDerECVs(unsigned j)
231 : {
232 6 : plumed_massert(isReady_,"cannot access ECVs before initialization");
233 6 : plumed_massert(j<getNumberOfArguments(),getName()+" has fewer CVs");
234 6 : return &derECVs_[j][0];
235 : }
236 :
237 3 : std::vector<std::string> ECVumbrellasFile::getLambdas() const
238 : { //notice that sigmas are not considered!
239 3 : std::vector<std::string> lambdas(totNumECVs_);
240 3 : if(P0_contribution_==1)
241 : {
242 2 : std::ostringstream subs;
243 2 : subs<<"P0";
244 4 : for(unsigned j=1; j<getNumberOfArguments(); j++)
245 2 : subs<<"_P0";
246 2 : lambdas[0]=subs.str();
247 2 : }
248 666 : for(unsigned k=P0_contribution_; k<totNumECVs_; k++)
249 : {
250 663 : const unsigned kk=k-P0_contribution_;
251 663 : std::ostringstream subs;
252 663 : subs<<centers_[0][kk];
253 1326 : for(unsigned j=1; j<getNumberOfArguments(); j++)
254 663 : subs<<"_"<<centers_[j][kk];
255 663 : lambdas[k]=subs.str();
256 663 : }
257 3 : return lambdas;
258 0 : }
259 :
260 3 : void ECVumbrellasFile::initECVs()
261 : {
262 3 : plumed_massert(!isReady_,"initialization should not be called twice");
263 3 : isReady_=true;
264 3 : log.printf(" *%4u windows for %s\n",totNumECVs_,getName().c_str());
265 3 : }
266 :
267 2 : void ECVumbrellasFile::initECVs_observ(const std::vector<double>& all_obs_cvs,const unsigned ncv,const unsigned index_j)
268 : {
269 : //this non-linear exansion never uses automatic initialization
270 2 : initECVs();
271 2 : calculateECVs(&all_obs_cvs[index_j]); //use only first obs point
272 6 : for(unsigned j=0; j<getNumberOfArguments(); j++)
273 888 : for(unsigned k=P0_contribution_; k<totNumECVs_; k++)
274 1680 : ECVs_[j][k]=std::min(barrier_/kbt_,ECVs_[j][k]);
275 2 : }
276 :
277 1 : void ECVumbrellasFile::initECVs_restart(const std::vector<std::string>& lambdas)
278 : {
279 : std::size_t pos=0;
280 2 : for(unsigned j=0; j<getNumberOfArguments()-1; j++)
281 1 : pos=lambdas[0].find("_",pos+1); //checking only lambdas[0] is hopefully enough
282 1 : plumed_massert(pos<lambdas[0].length(),"this should not happen, fewer '_' than expected in "+getName());
283 1 : pos=lambdas[0].find("_",pos+1);
284 1 : plumed_massert(pos>lambdas[0].length(),"this should not happen, more '_' than expected in "+getName());
285 :
286 1 : std::vector<std::string> myLambdas=getLambdas();
287 1 : plumed_massert(myLambdas.size()==lambdas.size(),"RESTART - mismatch in number of "+getName());
288 1 : plumed_massert(std::equal(myLambdas.begin(),myLambdas.end(),lambdas.begin()),"RESTART - mismatch in lambda values of "+getName());
289 :
290 1 : initECVs();
291 1 : }
292 :
293 : }
294 : }
|
