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 "core/PlumedMain.h"
22 : #include "core/Atoms.h"
23 :
24 : namespace PLMD {
25 : namespace opes {
26 :
27 : //+PLUMEDOC OPES_EXPANSION_CV ECV_MULTITHERMAL
28 : /*
29 : Expand a simulation to sample multiple temperatures simultaneously.
30 :
31 : The internal energy \f$U\f$ of of the system should be used as ARG.
32 : \f[
33 : \Delta u_{\beta'}=(\beta'-\beta) U\, ,
34 : \f]
35 : where \f$\beta'\f$ are the temperatures to be sampled and \f$\beta\f$ is the temperature at which the simulation is conducted.
36 : In case of fixed volume, the internal energy is simply the potential energy given by the \ref ENERGY colvar\f$U=E\f$, and you will run a multicanonical simulation.
37 : If instead the simulation is at fixed pressure \f$p\f$, the contribution of the volume must be added \f$U=E+pV\f$ (see example below).
38 :
39 : By defauly the needed steps in temperatures are automatically guessed from few initial unbiased MD steps, as descibed in \cite Invernizzi2020unified.
40 : Otherwise you can manually set this number with TEMP_STEPS.
41 : In both cases the steps will be geometrically spaced in temperature.
42 : Use instead the keyword NO_GEOM_SPACING for a linear spacing in the inverse temperature (beta), that typically increases the focus on lower temperatures.
43 : Finally, you can use instead the keyword TEMP_SET_ALL and explicitly provide each temperature.
44 :
45 : You can reweight the resulting simulation at any temperature in the chosen range, using e.g. \ref REWEIGHT_TEMP_PRESS.
46 : A similar target distribution can be sampled using \ref TD_MULTICANONICAL.
47 :
48 :
49 : \par Examples
50 :
51 : Fixed volume, multicanonical simulation:
52 :
53 : \plumedfile
54 : ene: ENERGY
55 : ecv: ECV_MULTITHERMAL ARG=ene TEMP=300 TEMP_MIN=300 TEMP_MAX=800
56 : opes: OPES_EXPANDED ARG=ecv.ene PACE=500
57 : \endplumedfile
58 :
59 : which, if your MD code passes the temperature to PLUMED, is equivalent to:
60 :
61 : \plumedfile
62 : ene: ENERGY
63 : ecv: ECV_MULTITHERMAL ARG=ene TEMP_MAX=800
64 : opes: OPES_EXPANDED ARG=ecv.ene PACE=500
65 : \endplumedfile
66 :
67 : If instead the pressure is fixed and the volume changes, you shuld calculate the internal energy first, \f$U=E+pV\f$
68 :
69 : \plumedfile
70 : ene: ENERGY
71 : vol: VOLUME
72 : intEne: CUSTOM PERIODIC=NO ARG=ene,vol FUNC=x+0.06022140857*y
73 : ecv: ECV_MULTITHERMAL ARG=intEne TEMP_MAX=800
74 : opes: OPES_EXPANDED ARG=ecv.intEne PACE=500
75 : \endplumedfile
76 :
77 : Notice that \f$p=0.06022140857\f$ corresponds to 1 bar when using the default PLUMED units.
78 :
79 : */
80 : //+ENDPLUMEDOC
81 :
82 : class ECVmultiThermal :
83 : public ExpansionCVs
84 : {
85 : private:
86 : bool todoAutomatic_;
87 : bool geom_spacing_;
88 : std::vector<double> ECVs_;
89 : std::vector<double> derECVs_; //(beta_k-beta0) or (temp0/temp_k-1)/kbt
90 : void initECVs();
91 :
92 : public:
93 : explicit ECVmultiThermal(const ActionOptions&);
94 : static void registerKeywords(Keywords& keys);
95 : void calculateECVs(const double *) override;
96 : const double * getPntrToECVs(unsigned) override;
97 : const double * getPntrToDerECVs(unsigned) override;
98 : std::vector<std::string> getLambdas() const override;
99 : void initECVs_observ(const std::vector<double>&,const unsigned,const unsigned) override;
100 : void initECVs_restart(const std::vector<std::string>&) override;
101 : };
102 :
103 10441 : PLUMED_REGISTER_ACTION(ECVmultiThermal,"ECV_MULTITHERMAL")
104 :
105 12 : void ECVmultiThermal::registerKeywords(Keywords& keys)
106 : {
107 12 : ExpansionCVs::registerKeywords(keys);
108 12 : keys.remove("ARG");
109 24 : keys.add("compulsory","ARG","the label of the internal energy of the system. If volume is fixed it is calculated by the \\ref ENERGY colvar");
110 24 : keys.add("optional","TEMP_MIN","the minimum of the temperature range");
111 24 : keys.add("optional","TEMP_MAX","the maximum of the temperature range");
112 24 : keys.add("optional","TEMP_STEPS","the number of steps in temperature");
113 24 : keys.add("optional","TEMP_SET_ALL","manually set all the temperatures");
114 24 : keys.addFlag("NO_GEOM_SPACING",false,"do not use geometrical spacing in temperature, but instead linear spacing in inverse temperature");
115 12 : }
116 :
117 11 : ECVmultiThermal::ECVmultiThermal(const ActionOptions&ao)
118 : : Action(ao)
119 : , ExpansionCVs(ao)
120 11 : , todoAutomatic_(false)
121 : {
122 11 : plumed_massert(getNumberOfArguments()==1,"only the internal energy should be given as ARG");
123 :
124 : //set temp0
125 11 : const double temp0=kbt_/plumed.getAtoms().getKBoltzmann();
126 :
127 : //parse temp range
128 11 : double temp_min=-1;
129 11 : double temp_max=-1;
130 11 : parse("TEMP_MIN",temp_min);
131 11 : parse("TEMP_MAX",temp_max);
132 11 : unsigned temp_steps=0;
133 22 : parse("TEMP_STEPS",temp_steps);
134 : std::vector<double> temps;
135 11 : parseVector("TEMP_SET_ALL",temps);
136 11 : parseFlag("NO_GEOM_SPACING",geom_spacing_);
137 11 : geom_spacing_=!geom_spacing_;
138 :
139 11 : checkRead();
140 :
141 : //set the intermediate temperatures
142 11 : if(temps.size()>0)
143 : {
144 2 : plumed_massert(temp_steps==0,"cannot set both TEMP_STEPS and TEMP_SET_ALL");
145 2 : plumed_massert(temp_min==-1 && temp_max==-1,"cannot set both TEMP_SET_ALL and TEMP_MIN/MAX");
146 2 : plumed_massert(temps.size()>=2,"set at least 2 temperatures");
147 2 : temp_min=temps[0];
148 2 : temp_max=temps[temps.size()-1];
149 2 : derECVs_.resize(temps.size());
150 10 : for(unsigned k=0; k<derECVs_.size(); k++)
151 : {
152 8 : derECVs_[k]=(temp0/temps[k]-1.)/kbt_;
153 8 : if(k<derECVs_.size()-1)
154 6 : plumed_massert(temps[k]<=temps[k+1],"TEMP_SET_ALL must be properly ordered");
155 : }
156 : }
157 : else
158 : { //get TEMP_MIN and TEMP_MAX
159 9 : plumed_massert(temp_min!=-1 || temp_max!=-1,"TEMP_MIN, TEMP_MAX or both, should be set");
160 9 : if(temp_min==-1)
161 : {
162 2 : temp_min=temp0;
163 2 : log.printf(" no TEMP_MIN provided, using TEMP_MIN=TEMP\n");
164 : }
165 9 : if(temp_max==-1)
166 : {
167 0 : temp_max=temp0;
168 0 : log.printf(" no TEMP_MAX provided, using TEMP_MAX=TEMP\n");
169 : }
170 9 : plumed_massert(temp_max>=temp_min,"TEMP_MAX should be bigger than TEMP_MIN");
171 9 : derECVs_.resize(2);
172 9 : derECVs_[0]=(temp0/temp_min-1.)/kbt_;
173 9 : derECVs_[1]=(temp0/temp_max-1.)/kbt_;
174 9 : if(temp_min==temp_max && temp_steps==0)
175 0 : temp_steps=1;
176 9 : if(temp_steps>0)
177 14 : derECVs_=getSteps(derECVs_[0],derECVs_[1],temp_steps,"TEMP",geom_spacing_,1./kbt_);
178 : else
179 2 : todoAutomatic_=true;
180 : }
181 : const double tol=1e-3; //if temp is taken from MD engine it might be numerically slightly different
182 11 : if(temp0<(1-tol)*temp_min || temp0>(1+tol)*temp_max)
183 0 : log.printf(" +++ WARNING +++ running at TEMP=%g which is outside the chosen temperature range\n",temp0);
184 :
185 : //print some info
186 11 : log.printf(" targeting a temperature range from TEMP_MIN=%g to TEMP_MAX=%g\n",temp_min,temp_max);
187 11 : if(!geom_spacing_)
188 1 : log.printf(" -- NO_GEOM_SPACING: inverse temperatures will be linearly spaced\n");
189 11 : }
190 :
191 586 : void ECVmultiThermal::calculateECVs(const double * ene)
192 : {
193 3618 : for(unsigned k=0; k<derECVs_.size(); k++)
194 3032 : ECVs_[k]=derECVs_[k]*ene[0];
195 : // derivatives never change: derECVs_k=(beta_k-beta0)
196 586 : }
197 :
198 11 : const double * ECVmultiThermal::getPntrToECVs(unsigned j)
199 : {
200 11 : plumed_massert(isReady_,"cannot access ECVs before initialization");
201 11 : plumed_massert(j==0,getName()+" has only one CV, the ENERGY");
202 11 : return &ECVs_[0];
203 : }
204 :
205 11 : const double * ECVmultiThermal::getPntrToDerECVs(unsigned j)
206 : {
207 11 : plumed_massert(isReady_,"cannot access ECVs before initialization");
208 11 : plumed_massert(j==0,getName()+" has only one CV, the ENERGY");
209 11 : return &derECVs_[0];
210 : }
211 :
212 11 : std::vector<std::string> ECVmultiThermal::getLambdas() const
213 : {
214 11 : plumed_massert(!todoAutomatic_,"cannot access lambdas before initializing them");
215 11 : const double temp0=kbt_/plumed.getAtoms().getKBoltzmann();
216 11 : std::vector<std::string> lambdas(derECVs_.size());
217 70 : for(unsigned k=0; k<derECVs_.size(); k++)
218 : {
219 59 : std::ostringstream subs;
220 59 : subs<<temp0/(derECVs_[k]*kbt_+1); //temp_k
221 59 : lambdas[k]=subs.str();
222 59 : }
223 11 : return lambdas;
224 0 : }
225 :
226 11 : void ECVmultiThermal::initECVs()
227 : {
228 11 : plumed_massert(!isReady_,"initialization should not be called twice");
229 11 : plumed_massert(!todoAutomatic_,"this should not happen");
230 11 : totNumECVs_=derECVs_.size();
231 11 : ECVs_.resize(derECVs_.size());
232 11 : isReady_=true;
233 11 : log.printf(" *%4lu temperatures for %s\n",derECVs_.size(),getName().c_str());
234 11 : }
235 :
236 8 : void ECVmultiThermal::initECVs_observ(const std::vector<double>& all_obs_cvs,const unsigned ncv,const unsigned index_j)
237 : {
238 8 : if(todoAutomatic_) //estimate the steps in beta from observations
239 : {
240 1 : plumed_massert(all_obs_cvs.size()%ncv==0 && index_j<ncv,"initECVs_observ parameters are inconsistent");
241 1 : std::vector<double> obs_ene(all_obs_cvs.size()/ncv); //copy only useful observation (would be better not to copy...)
242 11 : for(unsigned t=0; t<obs_ene.size(); t++)
243 10 : obs_ene[t]=all_obs_cvs[t*ncv+index_j];
244 1 : const unsigned temp_steps=estimateNumSteps(derECVs_[0],derECVs_[1],obs_ene,"TEMP");
245 1 : log.printf(" (spacing is in beta, not in temperature)\n");
246 2 : derECVs_=getSteps(derECVs_[0],derECVs_[1],temp_steps,"TEMP",geom_spacing_,1./kbt_);
247 1 : todoAutomatic_=false;
248 : }
249 8 : initECVs();
250 8 : calculateECVs(&all_obs_cvs[index_j]);
251 8 : }
252 :
253 3 : void ECVmultiThermal::initECVs_restart(const std::vector<std::string>& lambdas)
254 : {
255 3 : std::size_t pos=lambdas[0].find("_");
256 3 : plumed_massert(pos==std::string::npos,"this should not happen, only one CV is used in "+getName());
257 3 : if(todoAutomatic_)
258 : {
259 2 : derECVs_=getSteps(derECVs_[0],derECVs_[1],lambdas.size(),"TEMP",geom_spacing_,1./kbt_);
260 1 : todoAutomatic_=false;
261 : }
262 3 : std::vector<std::string> myLambdas=getLambdas();
263 3 : plumed_massert(myLambdas.size()==lambdas.size(),"RESTART - mismatch in number of "+getName());
264 3 : plumed_massert(std::equal(myLambdas.begin(),myLambdas.end(),lambdas.begin()),"RESTART - mismatch in lambda values of "+getName());
265 :
266 3 : initECVs();
267 3 : }
268 :
269 : }
270 : }
|
