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_MULTIBARIC
28 : /*
29 : Expand a simulation to sample multiple temperatures and pressures.
30 :
31 : The potential \ref ENERGY, \f$E\f$, and the \ref VOLUME, \f$V\f$, of the system should be used as ARG.
32 : \f[
33 : \Delta u_{\beta',p'}=(\beta'-\beta) E + (\beta' p' -\beta p) V\, ,
34 : \f]
35 : where \f$\beta', p'\f$ are the temperatures and pressures to be sampled, while \f$\beta, p\f$ is the temperature and pressure at which the simulation is conducted.
36 :
37 : If instead you wish to sample multiple temperatures and a single pressure, you should use \ref ECV_MULTITHERMAL with as ARG the internal energy \f$U=E+pV\f$.
38 :
39 : The TEMP_STEPS and PRESSURE_STEPS are automatically guessed from the initial unbiased steps (see OBSERVATION_STEPS in \ref OPES_EXPANDED), unless explicitly set.
40 : The algorithm for this guess is described in \cite Invernizzi2020unified should provide a rough estimate useful for most applications.
41 : The pressures are uniformely spaced, while the temperatures steps are geometrically spaced.
42 : Use instead the keyword NO_GEOM_SPACING for a linear spacing in inverse temperature (beta).
43 : For more detailed control you can use instead TEMP_SET_ALL and/or PRESSURE_SET_ALL to explicitly set all of them.
44 : The temperatures and pressures are then combined in a 2D grid.
45 :
46 : You can use CUT_CORNER to avoid a high-temperature/low-pressure region.
47 : This can be useful e.g. to increase the temperature for greater ergodicity, while avoiding water to vaporize, as in Ref.\cite Invernizzi2020unified.
48 :
49 : You can reweight the resulting simulation at any temperature and pressure in chosen target, using e.g. \ref REWEIGHT_TEMP_PRESS.
50 : A similar target distribution can be sampled using \ref TD_MULTITHERMAL_MULTIBARIC.
51 :
52 : \par Examples
53 :
54 : \plumedfile
55 : ene: ENERGY
56 : vol: VOLUME
57 : ecv: ECV_MULTITHERMAL_MULTIBARIC ...
58 : ARG=ene,vol
59 : TEMP=500
60 : TEMP_MIN=270
61 : TEMP_MAX=800
62 : PRESSURE=0.06022140857*2000 #2 kbar
63 : PRESSURE_MIN=0.06022140857 #1 bar
64 : PRESSURE_MAX=0.06022140857*4000 #4 kbar
65 : CUT_CORNER=500,0.06022140857,800,0.06022140857*1000
66 : ...
67 : opes: OPES_EXPANDED ARG=ecv.* FILE=DeltaF.data PACE=500 WALKERS_MPI
68 : \endplumedfile
69 :
70 : Notice that \f$p=0.06022140857\f$ corresponds to 1 bar only when using the default PLUMED units.
71 : If you modify them via the \ref UNITS command, then the pressure has to be rescaled accordingly.
72 :
73 : */
74 : //+ENDPLUMEDOC
75 :
76 : class ECVmultiThermalBaric :
77 : public ExpansionCVs
78 : {
79 : private:
80 : bool todoAutomatic_beta_;
81 : bool todoAutomatic_pres_;
82 : bool geom_spacing_;
83 : double pres0_;
84 : std::vector<double> pres_;
85 : std::vector<double> ECVs_beta_;
86 : std::vector<double> ECVs_pres_;
87 : std::vector<double> derECVs_beta_; //(beta_k-beta0) or (temp0/temp_k-1)/kbt
88 : std::vector<double> derECVs_pres_; //(beta_k*pres_kk-beta0*pres0) or (temp0/temp_k*pres_kk-pres0)/kbt
89 : void initECVs();
90 :
91 : //CUT_CORNER stuff
92 : double coeff_;
93 : double pres_low_;
94 : double kB_temp_low_;
95 : //SET_ALL_TEMP_PRESSURE stuff
96 : std::vector<std::string> custom_lambdas_;
97 :
98 : public:
99 : explicit ECVmultiThermalBaric(const ActionOptions&);
100 : static void registerKeywords(Keywords& keys);
101 : void calculateECVs(const double *) override;
102 : const double * getPntrToECVs(unsigned) override;
103 : const double * getPntrToDerECVs(unsigned) override;
104 : std::vector< std::vector<unsigned> > getIndex_k() const override;
105 : std::vector<std::string> getLambdas() const override;
106 : void initECVs_observ(const std::vector<double>&,const unsigned,const unsigned) override;
107 : void initECVs_restart(const std::vector<std::string>&) override;
108 : };
109 :
110 10437 : PLUMED_REGISTER_ACTION(ECVmultiThermalBaric,"ECV_MULTITHERMAL_MULTIBARIC")
111 :
112 10 : void ECVmultiThermalBaric::registerKeywords(Keywords& keys)
113 : {
114 10 : ExpansionCVs::registerKeywords(keys);
115 10 : keys.remove("ARG");
116 20 : keys.add("compulsory","ARG","the labels of the potential energy and of the volume of the system. You can calculate them with \\ref ENERGY and \\ref VOLUME respectively");
117 : //temperature
118 20 : keys.add("optional","TEMP_MIN","the minimum of the temperature range");
119 20 : keys.add("optional","TEMP_MAX","the maximum of the temperature range");
120 20 : keys.add("optional","TEMP_STEPS","the number of steps in temperature");
121 20 : keys.add("optional","TEMP_SET_ALL","manually set all the temperatures");
122 20 : keys.addFlag("NO_GEOM_SPACING",false,"do not use geometrical spacing in temperature, but instead linear spacing in inverse temperature");
123 : //pressure
124 20 : keys.add("compulsory","PRESSURE","pressure. Use the proper units");
125 20 : keys.add("optional","PRESSURE_MIN","the minimum of the pressure range");
126 20 : keys.add("optional","PRESSURE_MAX","the maximum of the pressure range");
127 20 : keys.add("optional","PRESSURE_STEPS","the number of steps in pressure");
128 30 : keys.add("optional","PRESSURE_SET_ALL","manually set all the pressures");
129 : //other
130 30 : keys.add("optional","SET_ALL_TEMP_PRESSURE","manually set all the target temperature_pressure pairs. An underscore separates temperature and pressure, while different points are comma-separated, e.g.: temp1_pres1,temp1_pres2,...");
131 20 : keys.add("optional","CUT_CORNER","avoid region of high temperature and low pressure. Exclude all points below a line in the temperature-pressure plane, defined by two points: \\f$T_{\\text{low}},P_{\\text{low}},T_{\\text{high}},P_{\\text{high}}\\f$");
132 10 : }
133 :
134 9 : ECVmultiThermalBaric::ECVmultiThermalBaric(const ActionOptions&ao)
135 : : Action(ao)
136 : , ExpansionCVs(ao)
137 9 : , todoAutomatic_beta_(false)
138 9 : , todoAutomatic_pres_(false)
139 9 : , coeff_(0)
140 9 : , pres_low_(0)
141 9 : , kB_temp_low_(0)
142 : {
143 9 : plumed_massert(getNumberOfArguments()==2,"ENERGY and VOLUME should be given as ARG");
144 :
145 : //set temp0
146 9 : const double kB=plumed.getAtoms().getKBoltzmann();
147 9 : const double temp0=kbt_/kB;
148 :
149 : //parse temp range
150 9 : double temp_min=-1;
151 9 : double temp_max=-1;
152 9 : parse("TEMP_MIN",temp_min);
153 9 : parse("TEMP_MAX",temp_max);
154 9 : unsigned temp_steps=0;
155 18 : parse("TEMP_STEPS",temp_steps);
156 : std::vector<double> temps;
157 9 : parseVector("TEMP_SET_ALL",temps);
158 9 : parseFlag("NO_GEOM_SPACING",geom_spacing_);
159 9 : geom_spacing_=!geom_spacing_;
160 : //parse pressures
161 9 : parse("PRESSURE",pres0_);
162 : const double myNone=std::numeric_limits<double>::lowest(); //quiet_NaN is not supported by some intel compiler
163 9 : double pres_min=myNone; //-1 might be a meaningful pressure
164 9 : double pres_max=myNone;
165 9 : parse("PRESSURE_MIN",pres_min);
166 9 : parse("PRESSURE_MAX",pres_max);
167 9 : unsigned pres_steps=0;
168 9 : parse("PRESSURE_STEPS",pres_steps);
169 18 : parseVector("PRESSURE_SET_ALL",pres_);
170 : //other
171 : std::vector<double> cut_corner;
172 9 : parseVector("CUT_CORNER",cut_corner);
173 9 : parseVector("SET_ALL_TEMP_PRESSURE",custom_lambdas_);
174 :
175 9 : checkRead();
176 :
177 9 : if(custom_lambdas_.size()>0)
178 : {
179 : //make sure no incompatible options are used
180 2 : plumed_massert(temps.size()==0,"cannot set both SET_ALL_TEMP_PRESSURE and TEMP_SET_ALL");
181 2 : plumed_massert(pres_.size()==0,"cannot set both SET_ALL_TEMP_PRESSURE and PRESSURE_SET_ALL");
182 2 : plumed_massert(temp_steps==0,"cannot set both SET_ALL_TEMP_PRESSURE and TEMP_STEPS");
183 2 : plumed_massert(pres_steps==0,"cannot set both SET_ALL_TEMP_PRESSURE and PRESSURE_STEPS");
184 2 : plumed_massert(temp_min==-1 && temp_max==-1,"cannot set both SET_ALL_TEMP_PRESSURE and TEMP_MIN/MAX");
185 2 : plumed_massert(pres_min==myNone && pres_max==myNone,"cannot set both SET_ALL_TEMP_PRESSURE and PRESSURE_MIN/MAX");
186 2 : plumed_massert(cut_corner.size()==0,"cannot set both SET_ALL_TEMP_PRESSURE and CUT_CORNER");
187 : //setup the target temperature-pressure grid
188 2 : derECVs_beta_.resize(custom_lambdas_.size());
189 2 : derECVs_pres_.resize(custom_lambdas_.size());
190 2 : const std::string error_msg="SET_ALL_TEMP_PRESSURE: two underscore-separated values are expected for each comma-separated point, cannot understand: ";
191 22 : for(unsigned i=0; i<custom_lambdas_.size(); i++)
192 : {
193 : try
194 : {
195 : std::size_t pos1;
196 : const double temp_i=std::stod(custom_lambdas_[i],&pos1);
197 20 : plumed_massert(pos1+1<custom_lambdas_[i].size(),error_msg+custom_lambdas_[i]);
198 20 : plumed_massert(custom_lambdas_[i][pos1]=='_',error_msg+custom_lambdas_[i]);
199 : std::size_t pos2;
200 20 : const double pres_i=std::stod(custom_lambdas_[i].substr(pos1+1),&pos2);
201 20 : plumed_massert(pos1+1+pos2==custom_lambdas_[i].size(),error_msg+custom_lambdas_[i]);
202 :
203 20 : derECVs_beta_[i]=(temp0/temp_i-1.)/kbt_;
204 20 : derECVs_pres_[i]=(temp0/temp_i*pres_i-pres0_)/kbt_;
205 : }
206 0 : catch (std::exception &ex)
207 : {
208 0 : plumed_merror(error_msg+custom_lambdas_[i]);
209 0 : }
210 : }
211 : }
212 : else
213 : {
214 : //set the intermediate temperatures
215 7 : if(temps.size()>0)
216 : {
217 1 : plumed_massert(temp_steps==0,"cannot set both TEMP_STEPS and TEMP_SET_ALL");
218 1 : plumed_massert(temp_min==-1 && temp_max==-1,"cannot set both TEMP_SET_ALL and TEMP_MIN/MAX");
219 1 : plumed_massert(temps.size()>=2,"set at least 2 temperatures");
220 1 : temp_min=temps[0];
221 1 : temp_max=temps[temps.size()-1];
222 1 : derECVs_beta_.resize(temps.size());
223 5 : for(unsigned k=0; k<derECVs_beta_.size(); k++)
224 : {
225 4 : derECVs_beta_[k]=(temp0/temps[k]-1.)/kbt_;
226 4 : if(k<derECVs_beta_.size()-1)
227 3 : plumed_massert(temps[k]<=temps[k+1],"TEMP_SET_ALL must be properly ordered");
228 : }
229 : }
230 : else
231 : { //get TEMP_MIN and TEMP_MAX
232 6 : if(temp_min==-1)
233 : {
234 0 : temp_min=temp0;
235 0 : log.printf(" no TEMP_MIN provided, using TEMP_MIN=TEMP\n");
236 : }
237 6 : if(temp_max==-1)
238 : {
239 1 : temp_max=temp0;
240 1 : log.printf(" no TEMP_MAX provided, using TEMP_MAX=TEMP\n");
241 : }
242 6 : plumed_massert(temp_max>=temp_min,"TEMP_MAX should be bigger than TEMP_MIN");
243 6 : derECVs_beta_.resize(2);
244 6 : derECVs_beta_[0]=(temp0/temp_min-1.)/kbt_;
245 6 : derECVs_beta_[1]=(temp0/temp_max-1.)/kbt_;
246 6 : if(temp_min==temp_max && temp_steps==0)
247 0 : temp_steps=1;
248 6 : if(temp_steps>0)
249 4 : derECVs_beta_=getSteps(derECVs_beta_[0],derECVs_beta_[1],temp_steps,"TEMP",geom_spacing_,1./kbt_);
250 : else
251 4 : todoAutomatic_beta_=true;
252 : }
253 : const double tol=1e-3; //if temp is taken from MD engine it might be numerically slightly different
254 7 : if(temp0<(1-tol)*temp_min || temp0>(1+tol)*temp_max)
255 1 : log.printf(" +++ WARNING +++ running at TEMP=%g which is outside the chosen temperature range\n",temp0);
256 :
257 : //set the intermediate pressures
258 7 : if(pres_.size()>0)
259 : {
260 1 : plumed_massert(pres_steps==0,"cannot set both PRESSURE_STEPS and PRESSURE_SET_ALL");
261 1 : plumed_massert(pres_min==myNone && pres_max==myNone,"cannot set both PRESSURE_SET_ALL and PRESSURE_MIN/MAX");
262 1 : plumed_massert(pres_.size()>=2,"set at least 2 pressures");
263 6 : for(unsigned kk=0; kk<pres_.size()-1; kk++)
264 5 : plumed_massert(pres_[kk]<=pres_[kk+1],"PRESSURE_SET_ALL must be properly ordered");
265 1 : pres_min=pres_[0];
266 1 : pres_max=pres_[pres_.size()-1];
267 : }
268 : else
269 : { //get PRESSURE_MIN and PRESSURE_MAX
270 6 : if(pres_min==myNone)
271 : {
272 3 : pres_min=pres0_;
273 3 : log.printf(" no PRESSURE_MIN provided, using PRESSURE_MIN=PRESSURE\n");
274 : }
275 6 : if(pres_max==myNone)
276 : {
277 2 : pres_max=pres0_;
278 2 : log.printf(" no PRESSURE_MAX provided, using PRESSURE_MAX=PRESSURE\n");
279 : }
280 6 : plumed_massert(pres_max>=pres_min,"PRESSURE_MAX should be bigger than PRESSURE_MIN");
281 6 : if(pres_min==pres_max && pres_steps==0)
282 0 : pres_steps=1;
283 6 : if(pres_steps>0)
284 4 : pres_=getSteps(pres_min,pres_max,pres_steps,"PRESSURE",false,0);
285 : else
286 : {
287 4 : pres_.resize(2);
288 4 : pres_[0]=pres_min;
289 4 : pres_[1]=pres_max;
290 4 : todoAutomatic_pres_=true;
291 : }
292 : }
293 7 : if(pres0_<pres_min || pres0_>pres_max)
294 0 : log.printf(" +++ WARNING +++ running at PRESSURE=%g which is outside the chosen pressure range\n",pres0_);
295 :
296 : //set CUT_CORNER
297 7 : std::string cc_usage("CUT_CORNER=temp_low,pres_low,temp_high,pres_high");
298 7 : if(cut_corner.size()==4)
299 : {
300 6 : const double temp_low=cut_corner[0];
301 6 : const double pres_low=cut_corner[1];
302 6 : const double temp_high=cut_corner[2];
303 6 : const double pres_high=cut_corner[3];
304 6 : plumed_massert(temp_low<temp_high,"temp_low="+std::to_string(temp_low)+" should be smaller than temp_high="+std::to_string(temp_high)+", "+cc_usage);
305 6 : plumed_massert(temp_low>=temp_min && temp_low<=temp_max,"temp_low="+std::to_string(temp_low)+" is out of temperature range. "+cc_usage);
306 6 : plumed_massert(temp_high>=temp_min && temp_high<=temp_max,"temp_high="+std::to_string(temp_high)+" is out of temperature range. "+cc_usage);
307 6 : plumed_massert(pres_low<pres_high,"pres_low="+std::to_string(pres_low)+" should be smaller than pres_high="+std::to_string(pres_high)+", "+cc_usage);
308 6 : plumed_massert(pres_low>=pres_min && pres_low<=pres_max,"pres_low="+std::to_string(pres_low)+" is out of pressure range. "+cc_usage);
309 6 : plumed_massert(pres_high>=pres_min && pres_high<=pres_max,"pres_high="+std::to_string(pres_high)+" is out of pressure range. "+cc_usage);
310 6 : kB_temp_low_=kB*temp_low;
311 6 : coeff_=(pres_high-pres_low)/(temp_high-temp_low)/kB;
312 6 : plumed_massert(coeff_!=0,"this should not be possible");
313 6 : const double small_value=(temp_high-pres_low)/1e4;
314 6 : pres_low_=pres_low-small_value; //make sure pres_max is included
315 6 : plumed_massert(pres_max>=coeff_*(kB*temp_max-kB_temp_low_)+pres_low_,"please chose a pres_high slightly smaller than PRESSURE_MAX in "+cc_usage);
316 : }
317 : else
318 : {
319 1 : plumed_massert(cut_corner.size()==0,"expected 4 values: "+cc_usage);
320 : }
321 : }
322 :
323 : //print some info
324 9 : log.printf(" running at TEMP=%g and PRESSURE=%g\n",temp0,pres0_);
325 9 : log.printf(" targeting a temperature range from TEMP_MIN=%g to TEMP_MAX=%g\n",temp_min,temp_max);
326 9 : if(temp_min==temp_max)
327 2 : log.printf(" +++ WARNING +++ if you only need a multibaric simulation it is more efficient to set it up with ECV_LINEAR\n");
328 9 : log.printf(" and a pressure range from PRESSURE_MIN=%g to PRESSURE_MAX=%g\n",pres_min,pres_max);
329 9 : if(pres_min==pres_max)
330 2 : log.printf(" +++ WARNING +++ if you only need a multithermal simulation it is more efficient to set it up with ECV_MULTITHERMAL\n");
331 9 : if(geom_spacing_)
332 8 : log.printf(" -- NO_GEOM_SPACING: inverse temperatures will be linearly spaced\n");
333 9 : if(coeff_!=0)
334 6 : log.printf(" -- CUT_CORNER: ignoring some high temperature and low pressure values\n");
335 9 : }
336 :
337 463 : void ECVmultiThermalBaric::calculateECVs(const double * ene_vol)
338 : {
339 5925 : for(unsigned k=0; k<derECVs_beta_.size(); k++)
340 5462 : ECVs_beta_[k]=derECVs_beta_[k]*ene_vol[0];
341 50075 : for(unsigned i=0; i<derECVs_pres_.size(); i++)
342 49612 : ECVs_pres_[i]=derECVs_pres_[i]*ene_vol[1];
343 : // derivatives are constant, as usual in linear expansions
344 463 : }
345 :
346 18 : const double * ECVmultiThermalBaric::getPntrToECVs(unsigned j)
347 : {
348 18 : plumed_massert(isReady_,"cannot access ECVs before initialization");
349 18 : plumed_massert(j==0 || j==1,getName()+" has only two CVs, the ENERGY and the VOLUME");
350 18 : if(j==0)
351 9 : return &ECVs_beta_[0];
352 : else //if (j==1)
353 9 : return &ECVs_pres_[0];
354 : }
355 :
356 18 : const double * ECVmultiThermalBaric::getPntrToDerECVs(unsigned j)
357 : {
358 18 : plumed_massert(isReady_,"cannot access ECVs before initialization");
359 18 : plumed_massert(j==0 || j==1,getName()+" has only two CVs, the ENERGY and the VOLUME");
360 18 : if(j==0)
361 9 : return &derECVs_beta_[0];
362 : else //if (j==1)
363 9 : return &derECVs_pres_[0];
364 : }
365 :
366 9 : std::vector< std::vector<unsigned> > ECVmultiThermalBaric::getIndex_k() const
367 : {
368 9 : plumed_massert(isReady_ && totNumECVs_>0,"cannot access getIndex_k() of ECV before initialization");
369 : std::vector< std::vector<unsigned> > index_k;
370 9 : if(custom_lambdas_.size()>0)
371 : { //same as default getIndex_k() function
372 2 : plumed_massert(totNumECVs_==custom_lambdas_.size(),"this should not happen");
373 22 : for(unsigned i=0; i<totNumECVs_; i++)
374 40 : index_k.emplace_back(std::vector<unsigned> {i,i});
375 : }
376 : else
377 : {
378 : unsigned i=0;
379 146 : for(unsigned k=0; k<derECVs_beta_.size(); k++)
380 : {
381 139 : const double kB_temp_k=kbt_/(derECVs_beta_[k]*kbt_+1);
382 139 : const double line_k=coeff_*(kB_temp_k-kB_temp_low_)+pres_low_;
383 2594 : for(unsigned kk=0; kk<pres_.size(); kk++)
384 : {
385 2455 : if(coeff_==0 || pres_[kk]>=line_k) //important to be inclusive, thus >=, not just >
386 : {
387 2254 : index_k.emplace_back(std::vector<unsigned> {k,i});
388 2254 : i++;
389 : }
390 : }
391 : }
392 7 : plumed_massert(totNumECVs_==index_k.size(),"this should not happen, is something wrong with CUT_CORNER ?");
393 : }
394 9 : return index_k;
395 0 : }
396 :
397 18 : std::vector<std::string> ECVmultiThermalBaric::getLambdas() const
398 : {
399 18 : if(custom_lambdas_.size()>0)
400 4 : return custom_lambdas_;
401 :
402 14 : plumed_massert(!todoAutomatic_beta_ && !todoAutomatic_pres_,"cannot access lambdas before initializing them");
403 : std::vector<std::string> lambdas;
404 14 : const double kB=plumed.getAtoms().getKBoltzmann();
405 292 : for(unsigned k=0; k<derECVs_beta_.size(); k++)
406 : {
407 278 : const double kB_temp_k=kbt_/(derECVs_beta_[k]*kbt_+1);
408 278 : const double line_k=coeff_*(kB_temp_k-kB_temp_low_)+pres_low_;
409 5188 : for(unsigned kk=0; kk<pres_.size(); kk++)
410 : {
411 4910 : if(coeff_==0 || pres_[kk]>=line_k)
412 : {
413 4508 : std::ostringstream subs;
414 4508 : subs<<kB_temp_k/kB<<"_"<<pres_[kk];
415 4508 : lambdas.emplace_back(subs.str());
416 4508 : }
417 : }
418 : }
419 : return lambdas;
420 14 : }
421 :
422 9 : void ECVmultiThermalBaric::initECVs()
423 : {
424 9 : plumed_massert(!isReady_,"initialization should not be called twice");
425 9 : plumed_massert(!todoAutomatic_beta_ && !todoAutomatic_pres_,"this should not happen");
426 9 : totNumECVs_=getLambdas().size(); //slow, but runs only once
427 9 : if(custom_lambdas_.size()>0)
428 : {
429 2 : log.printf(" *%4lu temperatures for %s\n",derECVs_beta_.size(),getName().c_str());
430 2 : log.printf(" *%4lu beta-pressures for %s\n",derECVs_pres_.size(),getName().c_str());
431 2 : log.printf(" -- SET_ALL_TEMP_PRESSURE: total number of temp-pres points is %u\n",totNumECVs_);
432 : }
433 : else
434 : {
435 7 : plumed_massert(derECVs_beta_.size()*pres_.size()>=totNumECVs_,"this should not happen, is something wrong with CUT_CORNER ?");
436 7 : derECVs_pres_.resize(totNumECVs_); //pres is mixed with temp (beta*p*V), thus we need to store all possible
437 : //initialize the derECVs.
438 : //this could be done before and one could avoid storing also beta0, beta_k, etc. but this way the code should be more readable
439 : unsigned i=0;
440 146 : for(unsigned k=0; k<derECVs_beta_.size(); k++)
441 : {
442 139 : const double kB_temp_k=kbt_/(derECVs_beta_[k]*kbt_+1);
443 139 : const double line_k=coeff_*(kB_temp_k-kB_temp_low_)+pres_low_;
444 2594 : for(unsigned kk=0; kk<pres_.size(); kk++)
445 : {
446 2455 : if(coeff_==0 || pres_[kk]>=line_k)
447 : {
448 2254 : derECVs_pres_[i]=(pres_[kk]/kB_temp_k-pres0_/kbt_);
449 2254 : i++;
450 : }
451 : }
452 : }
453 7 : log.printf(" *%4lu temperatures for %s\n",derECVs_beta_.size(),getName().c_str());
454 7 : log.printf(" *%4lu pressures for %s\n",pres_.size(),getName().c_str());
455 7 : if(coeff_!=0)
456 6 : log.printf(" -- CUT_CORNER: %lu temp-pres points were excluded, thus total is %u\n",derECVs_beta_.size()*pres_.size()-totNumECVs_,totNumECVs_);
457 : }
458 9 : ECVs_beta_.resize(derECVs_beta_.size());
459 9 : ECVs_pres_.resize(derECVs_pres_.size());
460 9 : isReady_=true;
461 9 : }
462 :
463 6 : void ECVmultiThermalBaric::initECVs_observ(const std::vector<double>& all_obs_cvs,const unsigned ncv,const unsigned index_j)
464 : {
465 6 : if(todoAutomatic_beta_) //estimate the steps in beta from observations
466 : {
467 2 : plumed_massert(all_obs_cvs.size()%ncv==0 && index_j<ncv,"initECVs_observ parameters are inconsistent");
468 2 : std::vector<double> obs_ene(all_obs_cvs.size()/ncv); //copy only useful observations
469 17 : for(unsigned t=0; t<obs_ene.size(); t++)
470 15 : obs_ene[t]=all_obs_cvs[t*ncv+index_j]+pres0_*all_obs_cvs[t*ncv+index_j+1]; //U=E+pV
471 2 : const unsigned temp_steps=estimateNumSteps(derECVs_beta_[0],derECVs_beta_[1],obs_ene,"TEMP");
472 2 : log.printf(" (spacing is on beta, not on temperature)\n");
473 4 : derECVs_beta_=getSteps(derECVs_beta_[0],derECVs_beta_[1],temp_steps,"TEMP",geom_spacing_,1./kbt_);
474 2 : todoAutomatic_beta_=false;
475 : }
476 6 : if(todoAutomatic_pres_) //estimate the steps in pres from observations
477 : {
478 2 : plumed_massert(all_obs_cvs.size()%ncv==0 && index_j+1<ncv,"initECVs_observ parameters are inconsistent");
479 2 : std::vector<double> obs_vol(all_obs_cvs.size()/ncv); //copy only useful observations
480 17 : for(unsigned t=0; t<obs_vol.size(); t++)
481 15 : obs_vol[t]=all_obs_cvs[t*ncv+index_j+1];
482 2 : const unsigned pres_steps=estimateNumSteps((pres_[0]-pres0_)/kbt_,(pres_[1]-pres0_)/kbt_,obs_vol,"PRESSURE");
483 2 : log.printf(" (spacing is in beta0 units)\n");
484 4 : pres_=getSteps(pres_[0],pres_[1],pres_steps,"PRESSURE",false,0);
485 2 : todoAutomatic_pres_=false;
486 : }
487 6 : initECVs();
488 6 : calculateECVs(&all_obs_cvs[index_j]);
489 6 : }
490 :
491 3 : void ECVmultiThermalBaric::initECVs_restart(const std::vector<std::string>& lambdas)
492 : {
493 3 : std::size_t pos=lambdas[0].find("_");
494 3 : plumed_massert(pos!=std::string::npos,"this should not happen, two CVs are used in "+getName()+", not less");
495 3 : pos=lambdas[0].find("_",pos+1);
496 3 : plumed_massert(pos==std::string::npos,"this should not happen, two CVs are used in "+getName()+", not more");
497 :
498 230 : auto getPres=[&lambdas](const unsigned i) {return lambdas[i].substr(lambdas[i].find("_")+1);};
499 3 : if(todoAutomatic_pres_)
500 : {
501 : unsigned pres_steps=1;
502 2 : std::string pres_min=getPres(0);
503 20 : for(unsigned i=1; i<lambdas.size(); i++) //pres is second, thus increas by 1
504 : {
505 20 : if(getPres(i)==pres_min)
506 : break;
507 18 : pres_steps++;
508 : }
509 4 : pres_=getSteps(pres_[0],pres_[1],pres_steps,"PRESSURE",false,0);
510 2 : todoAutomatic_pres_=false;
511 : }
512 3 : if(todoAutomatic_beta_)
513 : {
514 : unsigned temp_steps=1;
515 2 : std::string pres_max=getPres(pres_.size()-1);
516 208 : for(unsigned i=pres_.size(); i<lambdas.size(); i++)
517 : { //even if CUT_CORNER, the max pressures are all present, for each temp
518 206 : if(getPres(i)==pres_max)
519 24 : temp_steps++;
520 : }
521 4 : derECVs_beta_=getSteps(derECVs_beta_[0],derECVs_beta_[1],temp_steps,"TEMP",geom_spacing_,1./kbt_);
522 2 : todoAutomatic_beta_=false;
523 : }
524 3 : std::vector<std::string> myLambdas=getLambdas();
525 3 : plumed_assert(myLambdas.size()==lambdas.size())<<"RESTART - mismatch in number of "<<getName()<<".\nFrom "<<lambdas.size()<<" labels "<<derECVs_beta_.size()<<" temperatures and "<<pres_.size()<<" pressures were found, for a total of "<<myLambdas.size()<<" estimated steps.\nCheck if the CUT_CORNER or the SET_ALL_TEMP_PRESSURE options are consistent\n";
526 3 : plumed_massert(std::equal(myLambdas.begin(),myLambdas.end(),lambdas.begin()),"RESTART - mismatch in lambda values of "+getName());
527 :
528 3 : initECVs();
529 3 : }
530 :
531 : }
532 : }
|
