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 5 : keys.use("ARG");
102 10 : keys.add("compulsory","FILE","the name of the file containing the umbrellas");
103 10 : keys.add("optional","BARRIER","a guess of the free energy barrier to be overcome (better to stay higher than lower)");
104 10 : keys.addFlag("ADD_P0",false,"add the unbiased Boltzmann distribution to the target distribution, to make sure to sample it");
105 10 : keys.addFlag("LOWER_HALF_ONLY",false,"use only the lower half of each umbrella potentials");
106 5 : }
107 :
108 3 : ECVumbrellasFile::ECVumbrellasFile(const ActionOptions&ao):
109 : Action(ao),
110 3 : ExpansionCVs(ao)
111 : {
112 : //get number of CVs
113 : const unsigned ncv=getNumberOfArguments();
114 3 : centers_.resize(ncv);
115 3 : sigmas_.resize(ncv);
116 :
117 : //set P0_contribution_
118 3 : bool add_P0=false;
119 3 : parseFlag("ADD_P0",add_P0);
120 3 : if(add_P0)
121 2 : P0_contribution_=1;
122 : else
123 1 : P0_contribution_=0;
124 :
125 : //set barrier_
126 3 : barrier_=std::numeric_limits<double>::infinity();
127 3 : parse("BARRIER",barrier_);
128 6 : parseFlag("LOWER_HALF_ONLY",lower_only_);
129 :
130 : //set umbrellas
131 : std::string umbrellasFileName;
132 3 : parse("FILE",umbrellasFileName);
133 3 : IFile ifile;
134 3 : ifile.link(*this);
135 3 : if(ifile.FileExist(umbrellasFileName))
136 : {
137 3 : log.printf(" reading from FILE '%s'\n",umbrellasFileName.c_str());
138 3 : ifile.open(umbrellasFileName);
139 3 : ifile.allowIgnoredFields();
140 : double time; //first field is ignored
141 1332 : while(ifile.scanField("time",time))
142 : {
143 1989 : for(unsigned j=0; j<ncv; j++)
144 : {
145 : double centers_j;
146 1326 : ifile.scanField(getPntrToArgument(j)->getName(),centers_j);
147 1326 : centers_[j].push_back(centers_j); //this might be slow
148 : }
149 1989 : for(unsigned j=0; j<ncv; j++)
150 : {
151 : double sigmas_j;
152 2652 : ifile.scanField("sigma_"+getPntrToArgument(j)->getName(),sigmas_j);
153 1326 : sigmas_[j].push_back(sigmas_j);
154 : }
155 663 : ifile.scanField();
156 : }
157 : }
158 : else
159 0 : plumed_merror("Umbrellas FILE '"+umbrellasFileName+"' not found");
160 :
161 3 : checkRead();
162 :
163 : //extra consistency checks
164 3 : const unsigned sizeUmbrellas=centers_[0].size();
165 9 : for(unsigned j=0; j<ncv; j++)
166 : {
167 6 : plumed_massert(centers_[j].size()==sizeUmbrellas,"mismatch in the number of centers read from file");
168 6 : plumed_massert(sigmas_[j].size()==sizeUmbrellas,"mismatch in the number of sigmas read from file");
169 : }
170 :
171 : //set ECVs stuff
172 3 : totNumECVs_=sizeUmbrellas+P0_contribution_;
173 3 : ECVs_.resize(ncv,std::vector<double>(totNumECVs_));
174 3 : derECVs_.resize(ncv,std::vector<double>(totNumECVs_));
175 :
176 : //printing some info
177 3 : log.printf(" total number of umbrellas = %u\n",sizeUmbrellas);
178 3 : if(barrier_!=std::numeric_limits<double>::infinity())
179 1 : log.printf(" guess for free energy BARRIER = %g\n",barrier_);
180 3 : if(P0_contribution_==1)
181 2 : log.printf(" -- ADD_P0: the target includes also the unbiased probability itself\n");
182 3 : if(lower_only_)
183 1 : log.printf(" -- LOWER_HALF_ONLY: the ECVs are set to zero for values of the CV above the respective center\n");
184 6 : }
185 :
186 131 : void ECVumbrellasFile::calculateECVs(const double * cv)
187 : {
188 131 : if(lower_only_)
189 : {
190 153 : for(unsigned j=0; j<getNumberOfArguments(); j++)
191 : {
192 22644 : for(unsigned k=P0_contribution_; k<totNumECVs_; k++) //if ADD_P0, the first ECVs=0
193 : {
194 22542 : const unsigned kk=k-P0_contribution_;
195 22542 : const double dist_jk=difference(j,centers_[j][kk],cv[j])/sigmas_[j][kk]; //PBC might be present
196 22542 : if(dist_jk>=0)
197 : {
198 11483 : ECVs_[j][k]=0;
199 11483 : derECVs_[j][k]=0;
200 : }
201 : else
202 : {
203 11059 : ECVs_[j][k]=0.5*std::pow(dist_jk,2);
204 11059 : derECVs_[j][k]=dist_jk/sigmas_[j][kk];
205 : }
206 : }
207 : }
208 : }
209 : else
210 : {
211 240 : for(unsigned j=0; j<getNumberOfArguments(); j++)
212 : {
213 35520 : for(unsigned k=P0_contribution_; k<totNumECVs_; k++) //if ADD_P0, the first ECVs=0
214 : {
215 35360 : const unsigned kk=k-P0_contribution_;
216 35360 : const double dist_jk=difference(j,centers_[j][kk],cv[j])/sigmas_[j][kk]; //PBC might be present
217 35360 : ECVs_[j][k]=0.5*std::pow(dist_jk,2);
218 35360 : derECVs_[j][k]=dist_jk/sigmas_[j][kk];
219 : }
220 : }
221 : }
222 131 : }
223 :
224 6 : const double * ECVumbrellasFile::getPntrToECVs(unsigned j)
225 : {
226 6 : plumed_massert(isReady_,"cannot access ECVs before initialization");
227 6 : plumed_massert(j<getNumberOfArguments(),getName()+" has fewer CVs");
228 6 : return &ECVs_[j][0];
229 : }
230 :
231 6 : const double * ECVumbrellasFile::getPntrToDerECVs(unsigned j)
232 : {
233 6 : plumed_massert(isReady_,"cannot access ECVs before initialization");
234 6 : plumed_massert(j<getNumberOfArguments(),getName()+" has fewer CVs");
235 6 : return &derECVs_[j][0];
236 : }
237 :
238 3 : std::vector<std::string> ECVumbrellasFile::getLambdas() const
239 : { //notice that sigmas are not considered!
240 3 : std::vector<std::string> lambdas(totNumECVs_);
241 3 : if(P0_contribution_==1)
242 : {
243 2 : std::ostringstream subs;
244 2 : subs<<"P0";
245 4 : for(unsigned j=1; j<getNumberOfArguments(); j++)
246 2 : subs<<"_P0";
247 2 : lambdas[0]=subs.str();
248 2 : }
249 666 : for(unsigned k=P0_contribution_; k<totNumECVs_; k++)
250 : {
251 663 : const unsigned kk=k-P0_contribution_;
252 663 : std::ostringstream subs;
253 663 : subs<<centers_[0][kk];
254 1326 : for(unsigned j=1; j<getNumberOfArguments(); j++)
255 663 : subs<<"_"<<centers_[j][kk];
256 663 : lambdas[k]=subs.str();
257 663 : }
258 3 : return lambdas;
259 0 : }
260 :
261 3 : void ECVumbrellasFile::initECVs()
262 : {
263 3 : plumed_massert(!isReady_,"initialization should not be called twice");
264 3 : isReady_=true;
265 3 : log.printf(" *%4u windows for %s\n",totNumECVs_,getName().c_str());
266 3 : }
267 :
268 2 : void ECVumbrellasFile::initECVs_observ(const std::vector<double>& all_obs_cvs,const unsigned ncv,const unsigned index_j)
269 : {
270 : //this non-linear exansion never uses automatic initialization
271 2 : initECVs();
272 2 : calculateECVs(&all_obs_cvs[index_j]); //use only first obs point
273 6 : for(unsigned j=0; j<getNumberOfArguments(); j++)
274 888 : for(unsigned k=P0_contribution_; k<totNumECVs_; k++)
275 1680 : ECVs_[j][k]=std::min(barrier_/kbt_,ECVs_[j][k]);
276 2 : }
277 :
278 1 : void ECVumbrellasFile::initECVs_restart(const std::vector<std::string>& lambdas)
279 : {
280 : std::size_t pos=0;
281 2 : for(unsigned j=0; j<getNumberOfArguments()-1; j++)
282 1 : pos=lambdas[0].find("_",pos+1); //checking only lambdas[0] is hopefully enough
283 1 : plumed_massert(pos<lambdas[0].length(),"this should not happen, fewer '_' than expected in "+getName());
284 1 : pos=lambdas[0].find("_",pos+1);
285 1 : plumed_massert(pos>lambdas[0].length(),"this should not happen, more '_' than expected in "+getName());
286 :
287 1 : std::vector<std::string> myLambdas=getLambdas();
288 1 : plumed_massert(myLambdas.size()==lambdas.size(),"RESTART - mismatch in number of "+getName());
289 1 : plumed_massert(std::equal(myLambdas.begin(),myLambdas.end(),lambdas.begin()),"RESTART - mismatch in lambda values of "+getName());
290 :
291 1 : initECVs();
292 1 : }
293 :
294 : }
295 : }
|
