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 : private:
75 : bool todoAutomatic_;
76 : bool geom_spacing_;
77 : double beta0_;
78 : double lambda0_;
79 : std::vector<double> ECVs_;
80 : std::vector<double> derECVs_; //beta0*(lambda_k-lambda0)
81 : void initECVs();
82 :
83 : public:
84 : explicit ECVlinear(const ActionOptions&);
85 : static void registerKeywords(Keywords& keys);
86 : void calculateECVs(const double *) override;
87 : const double * getPntrToECVs(unsigned) override;
88 : const double * getPntrToDerECVs(unsigned) override;
89 : std::vector<std::string> getLambdas() const override;
90 : void initECVs_observ(const std::vector<double>&,const unsigned,const unsigned) override;
91 : void initECVs_restart(const std::vector<std::string>&) override;
92 : };
93 :
94 : PLUMED_REGISTER_ACTION(ECVlinear,"ECV_LINEAR")
95 :
96 6 : void ECVlinear::registerKeywords(Keywords& keys) {
97 6 : ExpansionCVs::registerKeywords(keys);
98 12 : keys.addInputKeyword("compulsory","ARG","scalar","the label of the Hamiltonian difference. Delta U");
99 6 : keys.add("compulsory","LAMBDA","0","the lambda at which the underlying simulation runs");
100 6 : keys.add("optional","LAMBDA_MIN","( default=0 ) the minimum of the lambda range");
101 6 : keys.add("optional","LAMBDA_MAX","( default=1 ) the maximum of the lambda range");
102 6 : keys.add("optional","LAMBDA_STEPS","uniformly place the lambda values, for a total of LAMBDA_STEPS");
103 6 : keys.add("optional","LAMBDA_SET_ALL","manually set all the lamdbas");
104 6 : keys.addFlag("DIMENSIONLESS",false,"ARG is considered dimensionless rather than an energy, thus is not multiplied by beta");
105 6 : keys.addFlag("GEOM_SPACING",false,"use geometrical spacing in lambda instead of linear spacing");
106 6 : }
107 :
108 4 : ECVlinear::ECVlinear(const ActionOptions&ao)
109 : : Action(ao)
110 : , ExpansionCVs(ao)
111 4 : , todoAutomatic_(false)
112 4 : , beta0_(1./kbt_) {
113 4 : plumed_massert(getNumberOfArguments()==1,"only DeltaU should be given as ARG");
114 :
115 : //set beta0_
116 : bool dimensionless;
117 4 : parseFlag("DIMENSIONLESS",dimensionless);
118 4 : if(dimensionless) {
119 1 : beta0_=1;
120 : }
121 :
122 : //parse lambda info
123 4 : parse("LAMBDA",lambda0_);
124 : const double myNone=std::numeric_limits<double>::lowest(); //quiet_NaN is not supported by some intel compiler
125 4 : double lambda_min=myNone;
126 4 : double lambda_max=myNone;
127 4 : parse("LAMBDA_MIN",lambda_min);
128 4 : parse("LAMBDA_MAX",lambda_max);
129 4 : unsigned lambda_steps=0;
130 8 : parse("LAMBDA_STEPS",lambda_steps);
131 : std::vector<double> lambdas;
132 4 : parseVector("LAMBDA_SET_ALL",lambdas);
133 4 : parseFlag("GEOM_SPACING",geom_spacing_);
134 :
135 4 : checkRead();
136 :
137 : //set the diff vector using lambdas
138 4 : if(lambdas.size()>0) {
139 1 : plumed_massert(lambda_steps==0,"cannot set both LAMBDA_STEPS and LAMBDA_SET_ALL");
140 1 : plumed_massert(lambda_min==myNone && lambda_max==myNone,"cannot set both LAMBDA_SET_ALL and LAMBDA_MIN/MAX");
141 1 : plumed_massert(lambdas.size()>=2,"set at least 2 lambdas with LAMBDA_SET_ALL");
142 4 : for(unsigned k=0; k<lambdas.size()-1; k++) {
143 3 : plumed_massert(lambdas[k]<=lambdas[k+1],"LAMBDA_SET_ALL must be properly ordered");
144 : }
145 1 : lambda_min=lambdas[0];
146 1 : lambda_max=lambdas[lambdas.size()-1];
147 1 : derECVs_.resize(lambdas.size());
148 5 : for(unsigned k=0; k<derECVs_.size(); k++) {
149 4 : derECVs_[k]=beta0_*(lambdas[k]-lambda0_);
150 : }
151 : } else {
152 : //get LAMBDA_MIN and LAMBDA_MAX
153 3 : if(lambda_min==myNone) {
154 0 : lambda_min=0;
155 0 : log.printf(" no LAMBDA_MIN provided, using LAMBDA_MIN = %g\n",lambda_min);
156 : }
157 3 : if(lambda_max==myNone) {
158 1 : lambda_max=1;
159 1 : log.printf(" no LAMBDA_MAX provided, using LAMBDA_MAX = %g\n",lambda_max);
160 : }
161 3 : plumed_massert(lambda_max>=lambda_min,"LAMBDA_MAX should be bigger than LAMBDA_MIN");
162 3 : derECVs_.resize(2);
163 3 : derECVs_[0]=beta0_*(lambda_min-lambda0_);
164 3 : derECVs_[1]=beta0_*(lambda_max-lambda0_);
165 3 : if(lambda_min==lambda_max && lambda_steps==0) {
166 0 : lambda_steps=1;
167 : }
168 3 : if(lambda_steps>0) {
169 2 : derECVs_=getSteps(derECVs_[0],derECVs_[1],lambda_steps,"LAMBDA",geom_spacing_,beta0_*lambda0_);
170 : } else {
171 2 : todoAutomatic_=true;
172 : }
173 : }
174 4 : if(lambda0_<lambda_min || lambda0_>lambda_max) {
175 1 : log.printf(" +++ WARNING +++ running at LAMBDA=%g which is outside the chosen lambda range\n",lambda0_);
176 : }
177 :
178 : //print some info
179 4 : log.printf(" running at LAMBDA=%g\n",lambda0_);
180 4 : log.printf(" targeting a lambda range from LAMBDA_MIN=%g to LAMBDA_MAX=%g\n",lambda_min,lambda_max);
181 4 : if(dimensionless) {
182 1 : log.printf(" -- DIMENSIONLESS: the ARG is not multiplied by beta\n");
183 : }
184 4 : if(geom_spacing_) {
185 1 : log.printf(" -- GEOM_SPACING: lambdas will be geometrically spaced\n");
186 : }
187 4 : }
188 :
189 182 : void ECVlinear::calculateECVs(const double * DeltaU) {
190 1587 : for(unsigned k=0; k<derECVs_.size(); k++) {
191 1405 : ECVs_[k]=derECVs_[k]*DeltaU[0];
192 : }
193 : // derivatives never change: derECVs_k=beta0*(lambda_k-lambda0)
194 182 : }
195 :
196 4 : const double * ECVlinear::getPntrToECVs(unsigned j) {
197 4 : plumed_massert(isReady_,"cannot access ECVs before initialization");
198 4 : plumed_massert(j==0,getName()+" has only one CV, the DeltaU");
199 4 : return &ECVs_[0];
200 : }
201 :
202 4 : const double * ECVlinear::getPntrToDerECVs(unsigned j) {
203 4 : plumed_massert(isReady_,"cannot access ECVs before initialization");
204 4 : plumed_massert(j==0,getName()+" has only one CV, the DeltaU");
205 4 : return &derECVs_[0];
206 : }
207 :
208 4 : std::vector<std::string> ECVlinear::getLambdas() const {
209 4 : plumed_massert(!todoAutomatic_,"cannot access lambdas before initializing them");
210 4 : std::vector<std::string> lambdas(derECVs_.size());
211 35 : for(unsigned k=0; k<derECVs_.size(); k++) {
212 31 : std::ostringstream subs;
213 31 : subs<<derECVs_[k]/beta0_+lambda0_; //lambda_k
214 31 : lambdas[k]=subs.str();
215 31 : }
216 4 : return lambdas;
217 0 : }
218 :
219 4 : void ECVlinear::initECVs() {
220 4 : plumed_massert(!isReady_,"initialization should not be called twice");
221 4 : plumed_massert(!todoAutomatic_,"this should not happen");
222 4 : totNumECVs_=derECVs_.size();
223 4 : ECVs_.resize(derECVs_.size());
224 4 : isReady_=true;
225 4 : log.printf(" *%4lu lambdas for %s\n",derECVs_.size(),getName().c_str());
226 4 : }
227 :
228 3 : void ECVlinear::initECVs_observ(const std::vector<double>& all_obs_cvs,const unsigned ncv,const unsigned index_j) {
229 3 : if(todoAutomatic_) { //estimate the steps in lambda from observations
230 1 : plumed_massert(all_obs_cvs.size()%ncv==0 && index_j<ncv,"initECVs_observ parameters are inconsistent");
231 1 : std::vector<double> obs_cv(all_obs_cvs.size()/ncv); //copy only useful observation (would be better not to copy...)
232 11 : for(unsigned t=0; t<obs_cv.size(); t++) {
233 10 : obs_cv[t]=all_obs_cvs[t*ncv+index_j];
234 : }
235 1 : const unsigned lambda_steps=estimateNumSteps(derECVs_[0],derECVs_[1],obs_cv,"LAMBDA");
236 1 : if(beta0_!=1) {
237 0 : log.printf(" (spacing is in beta0 units)\n");
238 : }
239 2 : derECVs_=getSteps(derECVs_[0],derECVs_[1],lambda_steps,"LAMBDA",geom_spacing_,beta0_*lambda0_);
240 1 : todoAutomatic_=false;
241 : }
242 3 : initECVs();
243 3 : calculateECVs(&all_obs_cvs[index_j]);
244 3 : }
245 :
246 1 : void ECVlinear::initECVs_restart(const std::vector<std::string>& lambdas) {
247 1 : std::size_t pos=lambdas[0].find("_");
248 1 : plumed_massert(pos==std::string::npos,"this should not happen, only one CV is used in "+getName());
249 1 : if(todoAutomatic_) {
250 2 : derECVs_=getSteps(derECVs_[0],derECVs_[1],lambdas.size(),"LAMBDA",geom_spacing_,beta0_*lambda0_);
251 1 : todoAutomatic_=false;
252 : }
253 1 : std::vector<std::string> myLambdas=getLambdas();
254 1 : plumed_massert(myLambdas.size()==lambdas.size(),"RESTART - mismatch in number of "+getName());
255 1 : plumed_massert(std::equal(myLambdas.begin(),myLambdas.end(),lambdas.begin()),"RESTART - mismatch in lambda values of "+getName());
256 :
257 1 : initECVs();
258 1 : }
259 :
260 : }
261 : }
|
