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 :
22 : namespace PLMD {
23 : namespace opes {
24 :
25 : //+PLUMEDOC OPES_EXPANSION_CV ECV_LINEAR
26 : /*
27 : Linear expansion, according to a parameter lambda.
28 :
29 : This can be used e.g. for thermodynamic integration, or for multibaric simulations, in which case lambda=pressure.
30 : It can also be used for multithermal simulations, but for simplicity it is more convenient to use \ref ECV_MULTITHERMAL.
31 :
32 : The difference in Hamiltonian \f$\Delta U\f$ is expected as ARG.
33 : \f[
34 : \Delta u_\lambda=\beta \lambda \Delta U\, .
35 : \f]
36 : Use the DIMENSIONLESS flag to avoid multiplying for the inverse temperature \f$\beta\f$.
37 :
38 : By defauly the needed steps in lambda are automatically guessed from few initial unbiased MD steps, as descibed in \cite Invernizzi2020unified.
39 : Otherwise one can set this number with LAMBDA_STEPS.
40 : In both cases the steps will be uniformly distriuted.
41 : Finally, one can use instead the keyword LAMBDA_SET_ALL and explicitly provide each lambda value.
42 :
43 : \par Examples
44 :
45 : Typical multibaric simulation:
46 :
47 : \plumedfile
48 : vol: VOLUME
49 : ecv: ECV_LINEAR ...
50 : ARG=vol
51 : TEMP=300
52 : LAMBDA=0.06022140857*2000 #2 kbar
53 : LAMBDA_MIN=0.06022140857 #1 bar
54 : LAMBDA_MAX=0.06022140857*4000 #4 kbar
55 : ...
56 : opes: OPES_EXPANDED ARG=ecv.vol PACE=500
57 : \endplumedfile
58 :
59 : Typical thermodynamic integration:
60 :
61 : \plumedfile
62 : DeltaU: EXTRACV NAME=energy_difference
63 : ecv: ECV_LINEAR ARG=DeltaU TEMP=300
64 : opes: OPES_EXPANDED ARG=ecv.* PACE=100
65 : \endplumedfile
66 :
67 : Notice that by defauly LAMBDA=0, LAMBDA_MIN=0 and LAMBDA_MAX=1, which is the typical case for thermodynamic integration.
68 :
69 : */
70 : //+ENDPLUMEDOC
71 :
72 : class ECVlinear :
73 : public ExpansionCVs
74 : {
75 : private:
76 : bool todoAutomatic_;
77 : bool geom_spacing_;
78 : double beta0_;
79 : double lambda0_;
80 : std::vector<double> ECVs_;
81 : std::vector<double> derECVs_; //beta0*(lambda_k-lambda0)
82 : void initECVs();
83 :
84 : public:
85 : explicit ECVlinear(const ActionOptions&);
86 : static void registerKeywords(Keywords& keys);
87 : void calculateECVs(const double *) override;
88 : const double * getPntrToECVs(unsigned) override;
89 : const double * getPntrToDerECVs(unsigned) override;
90 : std::vector<std::string> getLambdas() const override;
91 : void initECVs_observ(const std::vector<double>&,const unsigned,const unsigned) override;
92 : void initECVs_restart(const std::vector<std::string>&) override;
93 : };
94 :
95 : PLUMED_REGISTER_ACTION(ECVlinear,"ECV_LINEAR")
96 :
97 6 : void ECVlinear::registerKeywords(Keywords& keys)
98 : {
99 6 : ExpansionCVs::registerKeywords(keys);
100 6 : keys.remove("ARG");
101 12 : keys.add("compulsory","ARG","the label of the Hamiltonian difference. Delta U");
102 12 : keys.add("compulsory","LAMBDA","0","the lambda at which the underlying simulation runs");
103 12 : keys.add("optional","LAMBDA_MIN","( default=0 ) the minimum of the lambda range");
104 12 : keys.add("optional","LAMBDA_MAX","( default=1 ) the maximum of the lambda range");
105 12 : keys.add("optional","LAMBDA_STEPS","uniformly place the lambda values, for a total of LAMBDA_STEPS");
106 12 : keys.add("optional","LAMBDA_SET_ALL","manually set all the lamdbas");
107 12 : keys.addFlag("DIMENSIONLESS",false,"ARG is considered dimensionless rather than an energy, thus is not multiplied by beta");
108 12 : keys.addFlag("GEOM_SPACING",false,"use geometrical spacing in lambda instead of linear spacing");
109 6 : }
110 :
111 4 : ECVlinear::ECVlinear(const ActionOptions&ao)
112 : : Action(ao)
113 : , ExpansionCVs(ao)
114 4 : , todoAutomatic_(false)
115 4 : , beta0_(1./kbt_)
116 : {
117 4 : plumed_massert(getNumberOfArguments()==1,"only DeltaU should be given as ARG");
118 :
119 : //set beta0_
120 : bool dimensionless;
121 4 : parseFlag("DIMENSIONLESS",dimensionless);
122 4 : if(dimensionless)
123 1 : beta0_=1;
124 :
125 : //parse lambda info
126 4 : parse("LAMBDA",lambda0_);
127 : const double myNone=std::numeric_limits<double>::lowest(); //quiet_NaN is not supported by some intel compiler
128 4 : double lambda_min=myNone;
129 4 : double lambda_max=myNone;
130 4 : parse("LAMBDA_MIN",lambda_min);
131 4 : parse("LAMBDA_MAX",lambda_max);
132 4 : unsigned lambda_steps=0;
133 8 : parse("LAMBDA_STEPS",lambda_steps);
134 : std::vector<double> lambdas;
135 4 : parseVector("LAMBDA_SET_ALL",lambdas);
136 4 : parseFlag("GEOM_SPACING",geom_spacing_);
137 :
138 4 : checkRead();
139 :
140 : //set the diff vector using lambdas
141 4 : if(lambdas.size()>0)
142 : {
143 1 : plumed_massert(lambda_steps==0,"cannot set both LAMBDA_STEPS and LAMBDA_SET_ALL");
144 1 : plumed_massert(lambda_min==myNone && lambda_max==myNone,"cannot set both LAMBDA_SET_ALL and LAMBDA_MIN/MAX");
145 1 : plumed_massert(lambdas.size()>=2,"set at least 2 lambdas with LAMBDA_SET_ALL");
146 4 : for(unsigned k=0; k<lambdas.size()-1; k++)
147 3 : plumed_massert(lambdas[k]<=lambdas[k+1],"LAMBDA_SET_ALL must be properly ordered");
148 1 : lambda_min=lambdas[0];
149 1 : lambda_max=lambdas[lambdas.size()-1];
150 1 : derECVs_.resize(lambdas.size());
151 5 : for(unsigned k=0; k<derECVs_.size(); k++)
152 4 : derECVs_[k]=beta0_*(lambdas[k]-lambda0_);
153 : }
154 : else
155 : { //get LAMBDA_MIN and LAMBDA_MAX
156 3 : if(lambda_min==myNone)
157 : {
158 0 : lambda_min=0;
159 0 : log.printf(" no LAMBDA_MIN provided, using LAMBDA_MIN = %g\n",lambda_min);
160 : }
161 3 : if(lambda_max==myNone)
162 : {
163 1 : lambda_max=1;
164 1 : log.printf(" no LAMBDA_MAX provided, using LAMBDA_MAX = %g\n",lambda_max);
165 : }
166 3 : plumed_massert(lambda_max>=lambda_min,"LAMBDA_MAX should be bigger than LAMBDA_MIN");
167 3 : derECVs_.resize(2);
168 3 : derECVs_[0]=beta0_*(lambda_min-lambda0_);
169 3 : derECVs_[1]=beta0_*(lambda_max-lambda0_);
170 3 : if(lambda_min==lambda_max && lambda_steps==0)
171 0 : lambda_steps=1;
172 3 : if(lambda_steps>0)
173 2 : derECVs_=getSteps(derECVs_[0],derECVs_[1],lambda_steps,"LAMBDA",geom_spacing_,beta0_*lambda0_);
174 : else
175 2 : todoAutomatic_=true;
176 : }
177 4 : if(lambda0_<lambda_min || lambda0_>lambda_max)
178 1 : log.printf(" +++ WARNING +++ running at LAMBDA=%g which is outside the chosen lambda range\n",lambda0_);
179 :
180 : //print some info
181 4 : log.printf(" running at LAMBDA=%g\n",lambda0_);
182 4 : log.printf(" targeting a lambda range from LAMBDA_MIN=%g to LAMBDA_MAX=%g\n",lambda_min,lambda_max);
183 4 : if(dimensionless)
184 1 : log.printf(" -- DIMENSIONLESS: the ARG is not multiplied by beta\n");
185 4 : if(geom_spacing_)
186 1 : log.printf(" -- GEOM_SPACING: lambdas will be geometrically spaced\n");
187 4 : }
188 :
189 182 : void ECVlinear::calculateECVs(const double * DeltaU)
190 : {
191 1587 : for(unsigned k=0; k<derECVs_.size(); k++)
192 1405 : ECVs_[k]=derECVs_[k]*DeltaU[0];
193 : // derivatives never change: derECVs_k=beta0*(lambda_k-lambda0)
194 182 : }
195 :
196 4 : const double * ECVlinear::getPntrToECVs(unsigned j)
197 : {
198 4 : plumed_massert(isReady_,"cannot access ECVs before initialization");
199 4 : plumed_massert(j==0,getName()+" has only one CV, the DeltaU");
200 4 : return &ECVs_[0];
201 : }
202 :
203 4 : const double * ECVlinear::getPntrToDerECVs(unsigned j)
204 : {
205 4 : plumed_massert(isReady_,"cannot access ECVs before initialization");
206 4 : plumed_massert(j==0,getName()+" has only one CV, the DeltaU");
207 4 : return &derECVs_[0];
208 : }
209 :
210 4 : std::vector<std::string> ECVlinear::getLambdas() const
211 : {
212 4 : plumed_massert(!todoAutomatic_,"cannot access lambdas before initializing them");
213 4 : std::vector<std::string> lambdas(derECVs_.size());
214 35 : for(unsigned k=0; k<derECVs_.size(); k++)
215 : {
216 31 : std::ostringstream subs;
217 31 : subs<<derECVs_[k]/beta0_+lambda0_; //lambda_k
218 31 : lambdas[k]=subs.str();
219 31 : }
220 4 : return lambdas;
221 0 : }
222 :
223 4 : void ECVlinear::initECVs()
224 : {
225 4 : plumed_massert(!isReady_,"initialization should not be called twice");
226 4 : plumed_massert(!todoAutomatic_,"this should not happen");
227 4 : totNumECVs_=derECVs_.size();
228 4 : ECVs_.resize(derECVs_.size());
229 4 : isReady_=true;
230 4 : log.printf(" *%4lu lambdas for %s\n",derECVs_.size(),getName().c_str());
231 4 : }
232 :
233 3 : void ECVlinear::initECVs_observ(const std::vector<double>& all_obs_cvs,const unsigned ncv,const unsigned index_j)
234 : {
235 3 : if(todoAutomatic_) //estimate the steps in lambda from observations
236 : {
237 1 : plumed_massert(all_obs_cvs.size()%ncv==0 && index_j<ncv,"initECVs_observ parameters are inconsistent");
238 1 : std::vector<double> obs_cv(all_obs_cvs.size()/ncv); //copy only useful observation (would be better not to copy...)
239 11 : for(unsigned t=0; t<obs_cv.size(); t++)
240 10 : obs_cv[t]=all_obs_cvs[t*ncv+index_j];
241 1 : const unsigned lambda_steps=estimateNumSteps(derECVs_[0],derECVs_[1],obs_cv,"LAMBDA");
242 1 : if(beta0_!=1)
243 0 : log.printf(" (spacing is in beta0 units)\n");
244 2 : derECVs_=getSteps(derECVs_[0],derECVs_[1],lambda_steps,"LAMBDA",geom_spacing_,beta0_*lambda0_);
245 1 : todoAutomatic_=false;
246 : }
247 3 : initECVs();
248 3 : calculateECVs(&all_obs_cvs[index_j]);
249 3 : }
250 :
251 1 : void ECVlinear::initECVs_restart(const std::vector<std::string>& lambdas)
252 : {
253 1 : std::size_t pos=lambdas[0].find("_");
254 1 : plumed_massert(pos==std::string::npos,"this should not happen, only one CV is used in "+getName());
255 1 : if(todoAutomatic_)
256 : {
257 2 : derECVs_=getSteps(derECVs_[0],derECVs_[1],lambdas.size(),"LAMBDA",geom_spacing_,beta0_*lambda0_);
258 1 : todoAutomatic_=false;
259 : }
260 1 : std::vector<std::string> myLambdas=getLambdas();
261 1 : plumed_massert(myLambdas.size()==lambdas.size(),"RESTART - mismatch in number of "+getName());
262 1 : plumed_massert(std::equal(myLambdas.begin(),myLambdas.end(),lambdas.begin()),"RESTART - mismatch in lambda values of "+getName());
263 :
264 1 : initECVs();
265 1 : }
266 :
267 : }
268 : }
|
