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