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 [ECV_MULTITHERMAL](ECV_MULTITHERMAL.md).
31 :
32 : The difference in Hamiltonian $\Delta U$ is expected as ARG.
33 :
34 : $$
35 : \Delta u_\lambda=\beta \lambda \Delta U\, .
36 : $$
37 :
38 : Use the DIMENSIONLESS flag to avoid multiplying for the inverse temperature $\beta$.
39 :
40 : By default the needed steps in lambda are automatically guessed from few initial unbiased MD steps, as described in the paper cited below.
41 : Otherwise one can set this number with LAMBDA_STEPS.
42 : In both cases the steps will be uniformly distriuted.
43 : Finally, one can use instead the keyword LAMBDA_SET_ALL and explicitly provide each lambda value.
44 :
45 : ## Examples
46 :
47 : Typical multibaric simulation:
48 :
49 : ```plumed
50 : vol: VOLUME
51 : ecv: ECV_LINEAR ...
52 : ARG=vol
53 : TEMP=300
54 : LAMBDA=0.06022140857*2000 #2 kbar
55 : LAMBDA_MIN=0.06022140857 #1 bar
56 : LAMBDA_MAX=0.06022140857*4000 #4 kbar
57 : ...
58 : opes: OPES_EXPANDED ARG=ecv.vol PACE=500
59 : ```
60 :
61 : Typical thermodynamic integration:
62 :
63 : ```plumed
64 : DeltaU: EXTRACV NAME=energy_difference
65 : ecv: ECV_LINEAR ARG=DeltaU TEMP=300
66 : opes: OPES_EXPANDED ARG=ecv.* PACE=100
67 : ```
68 :
69 : Notice that by default LAMBDA=0, LAMBDA_MIN=0 and LAMBDA_MAX=1, which is the typical case for thermodynamic integration.
70 :
71 : */
72 : //+ENDPLUMEDOC
73 :
74 : class ECVlinear :
75 : public ExpansionCVs {
76 : private:
77 : bool todoAutomatic_;
78 : bool geom_spacing_;
79 : double beta0_;
80 : double lambda0_;
81 : std::vector<double> ECVs_;
82 : std::vector<double> derECVs_; //beta0*(lambda_k-lambda0)
83 : void initECVs();
84 :
85 : public:
86 : explicit ECVlinear(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(ECVlinear,"ECV_LINEAR")
97 :
98 8 : void ECVlinear::registerKeywords(Keywords& keys) {
99 8 : ExpansionCVs::registerKeywords(keys);
100 16 : keys.addInputKeyword("compulsory","ARG","scalar","the label of the Hamiltonian difference. Delta U");
101 8 : keys.add("compulsory","LAMBDA","0","the lambda at which the underlying simulation runs");
102 8 : keys.add("optional","LAMBDA_MIN","( default=0 ) the minimum of the lambda range");
103 8 : keys.add("optional","LAMBDA_MAX","( default=1 ) the maximum of the lambda range");
104 8 : keys.add("optional","LAMBDA_STEPS","uniformly place the lambda values, for a total of LAMBDA_STEPS");
105 8 : keys.add("optional","LAMBDA_SET_ALL","manually set all the lamdbas");
106 8 : keys.addFlag("DIMENSIONLESS",false,"ARG is considered dimensionless rather than an energy, thus is not multiplied by beta");
107 8 : keys.addFlag("GEOM_SPACING",false,"use geometrical spacing in lambda instead of linear spacing");
108 8 : keys.addDOI("10.1103/PhysRevX.10.041034");
109 8 : }
110 :
111 4 : ECVlinear::ECVlinear(const ActionOptions&ao)
112 : : Action(ao)
113 : , ExpansionCVs(ao)
114 4 : , todoAutomatic_(false)
115 4 : , beta0_(1./kbt_) {
116 4 : plumed_massert(getNumberOfArguments()==1,"only DeltaU should be given as ARG");
117 :
118 : //set beta0_
119 : bool dimensionless;
120 4 : parseFlag("DIMENSIONLESS",dimensionless);
121 4 : if(dimensionless) {
122 1 : beta0_=1;
123 : }
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 4 : parse("LAMBDA_STEPS",lambda_steps);
134 4 : 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 1 : plumed_massert(lambda_steps==0,"cannot set both LAMBDA_STEPS and LAMBDA_SET_ALL");
143 1 : plumed_massert(lambda_min==myNone && lambda_max==myNone,"cannot set both LAMBDA_SET_ALL and LAMBDA_MIN/MAX");
144 1 : plumed_massert(lambdas.size()>=2,"set at least 2 lambdas with LAMBDA_SET_ALL");
145 4 : for(unsigned k=0; k<lambdas.size()-1; k++) {
146 3 : plumed_massert(lambdas[k]<=lambdas[k+1],"LAMBDA_SET_ALL must be properly ordered");
147 : }
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 0 : lambda_min=0;
158 0 : log.printf(" no LAMBDA_MIN provided, using LAMBDA_MIN = %g\n",lambda_min);
159 : }
160 3 : if(lambda_max==myNone) {
161 1 : lambda_max=1;
162 1 : log.printf(" no LAMBDA_MAX provided, using LAMBDA_MAX = %g\n",lambda_max);
163 : }
164 3 : plumed_massert(lambda_max>=lambda_min,"LAMBDA_MAX should be bigger than LAMBDA_MIN");
165 3 : derECVs_.resize(2);
166 3 : derECVs_[0]=beta0_*(lambda_min-lambda0_);
167 3 : derECVs_[1]=beta0_*(lambda_max-lambda0_);
168 3 : if(lambda_min==lambda_max && lambda_steps==0) {
169 0 : lambda_steps=1;
170 : }
171 3 : if(lambda_steps>0) {
172 3 : derECVs_=getSteps(derECVs_[0],derECVs_[1],lambda_steps,"LAMBDA",geom_spacing_,beta0_*lambda0_);
173 : } else {
174 2 : todoAutomatic_=true;
175 : }
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 :
181 : //print some info
182 4 : log.printf(" running at LAMBDA=%g\n",lambda0_);
183 4 : log.printf(" targeting a lambda range from LAMBDA_MIN=%g to LAMBDA_MAX=%g\n",lambda_min,lambda_max);
184 4 : if(dimensionless) {
185 1 : log.printf(" -- DIMENSIONLESS: the ARG is not multiplied by beta\n");
186 : }
187 4 : if(geom_spacing_) {
188 1 : log.printf(" -- GEOM_SPACING: lambdas will be geometrically spaced\n");
189 : }
190 4 : }
191 :
192 182 : void ECVlinear::calculateECVs(const double * DeltaU) {
193 1587 : for(unsigned k=0; k<derECVs_.size(); k++) {
194 1405 : ECVs_[k]=derECVs_[k]*DeltaU[0];
195 : }
196 : // derivatives never change: derECVs_k=beta0*(lambda_k-lambda0)
197 182 : }
198 :
199 4 : const double * ECVlinear::getPntrToECVs(unsigned j) {
200 4 : plumed_massert(isReady_,"cannot access ECVs before initialization");
201 4 : plumed_massert(j==0,getName()+" has only one CV, the DeltaU");
202 4 : return &ECVs_[0];
203 : }
204 :
205 4 : const double * ECVlinear::getPntrToDerECVs(unsigned j) {
206 4 : plumed_massert(isReady_,"cannot access ECVs before initialization");
207 4 : plumed_massert(j==0,getName()+" has only one CV, the DeltaU");
208 4 : return &derECVs_[0];
209 : }
210 :
211 4 : std::vector<std::string> ECVlinear::getLambdas() const {
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 31 : std::ostringstream subs;
216 31 : subs<<derECVs_[k]/beta0_+lambda0_; //lambda_k
217 31 : lambdas[k]=subs.str();
218 31 : }
219 4 : return lambdas;
220 0 : }
221 :
222 4 : void ECVlinear::initECVs() {
223 4 : plumed_massert(!isReady_,"initialization should not be called twice");
224 4 : plumed_massert(!todoAutomatic_,"this should not happen");
225 4 : totNumECVs_=derECVs_.size();
226 4 : ECVs_.resize(derECVs_.size());
227 4 : isReady_=true;
228 4 : log.printf(" *%4lu lambdas for %s\n",derECVs_.size(),getName().c_str());
229 4 : }
230 :
231 3 : void ECVlinear::initECVs_observ(const std::vector<double>& all_obs_cvs,const unsigned ncv,const unsigned index_j) {
232 3 : if(todoAutomatic_) { //estimate the steps in lambda from observations
233 1 : plumed_massert(all_obs_cvs.size()%ncv==0 && index_j<ncv,"initECVs_observ parameters are inconsistent");
234 1 : std::vector<double> obs_cv(all_obs_cvs.size()/ncv); //copy only useful observation (would be better not to copy...)
235 11 : for(unsigned t=0; t<obs_cv.size(); t++) {
236 10 : obs_cv[t]=all_obs_cvs[t*ncv+index_j];
237 : }
238 1 : const unsigned lambda_steps=estimateNumSteps(derECVs_[0],derECVs_[1],obs_cv,"LAMBDA");
239 1 : if(beta0_!=1) {
240 0 : log.printf(" (spacing is in beta0 units)\n");
241 : }
242 1 : derECVs_=getSteps(derECVs_[0],derECVs_[1],lambda_steps,"LAMBDA",geom_spacing_,beta0_*lambda0_);
243 1 : todoAutomatic_=false;
244 1 : }
245 3 : initECVs();
246 3 : calculateECVs(&all_obs_cvs[index_j]);
247 3 : }
248 :
249 1 : void ECVlinear::initECVs_restart(const std::vector<std::string>& lambdas) {
250 1 : std::size_t pos=lambdas[0].find("_");
251 1 : plumed_massert(pos==std::string::npos,"this should not happen, only one CV is used in "+getName());
252 1 : if(todoAutomatic_) {
253 1 : derECVs_=getSteps(derECVs_[0],derECVs_[1],lambdas.size(),"LAMBDA",geom_spacing_,beta0_*lambda0_);
254 1 : todoAutomatic_=false;
255 : }
256 1 : std::vector<std::string> myLambdas=getLambdas();
257 1 : plumed_massert(myLambdas.size()==lambdas.size(),"RESTART - mismatch in number of "+getName());
258 1 : plumed_massert(std::equal(myLambdas.begin(),myLambdas.end(),lambdas.begin()),"RESTART - mismatch in lambda values of "+getName());
259 :
260 1 : initECVs();
261 1 : }
262 :
263 : }
264 : }
|
