Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2019-2023 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "core/ActionRegister.h"
23 : #include "core/PlumedMain.h"
24 : #include "core/Atoms.h"
25 : #include "ReweightBase.h"
26 :
27 : //+PLUMEDOC REWEIGHTING REWEIGHT_TEMP_PRESS
28 : /*
29 : Calculate weights for ensemble averages at temperatures and/or pressures different than those used in your original simulation.
30 :
31 : We can use our knowledge of the probability distribution in the canonical (N\f$\mathcal{V}\f$T) or the isothermal-isobaric ensemble (NPT) to reweight the data
32 : contained in trajectories and obtain ensemble averages at different temperatures and/or pressures.
33 :
34 : Consider the ensemble average of an observable \f$O(\mathbf{R},\mathcal{V})\f$ that depends on the atomic coordinates \f$\mathbf{R}\f$ and the volume \f$\mathcal{V}\f$.
35 : This observable is in practice any collective variable (CV) calculated by Plumed.
36 : The ensemble average of the observable in an ensemble \f$ \xi' \f$ can be calculated from a simulation performed in an ensemble \f$ \xi \f$ using:
37 : \f[
38 : \langle O(\mathbf{R},\mathcal{V}) \rangle_{\xi'} = \frac{\langle O(\mathbf{R},\mathcal{V}) w(\mathbf{R},\mathcal{V}) \rangle_{\xi}}
39 : {\langle w(\mathbf{R},\mathcal{V}) \rangle_{\xi}}
40 : \f]
41 : where \f$\langle \cdot \rangle_{\xi}\f$ and \f$\langle \cdot \rangle_{\xi'}\f$ are mean values in the simulated and targeted ensemble, respectively, \f$ E(\mathbf{R}) \f$ is the potential energy of the system, and \f$ w (\mathbf{R},\mathcal{V}) \f$ are the appropriate weights to take from \f$ \xi \f$ to \f$ \xi' \f$.
42 : This action calculates the weights \f$ w (\mathbf{R},\mathcal{V}) \f$ and handles 4 different cases:
43 : 1. Change of temperature from T to T' at constant volume. That is to say, from a simulation performed in the N\f$\mathcal{V}\f$T (canonical) ensemble, obtain an ensemble average in the N\f$\mathcal{V}\f$T' ensemble. The weights in this case are \f$ w(\mathbf{R},\mathcal{V}) = e^{(\beta-\beta')E(\mathbf{R})} \f$ with \f$ \beta \f$ and \f$ \beta' \f$ the inverse temperatures.
44 : 2. Change of temperature from T to T' at constant pressure. That is to say, from a simulation performed in the NPT (isothermal-isobaric) ensemble, obtain an ensemble average in the NPT' ensemble. The weights in this case are \f$ w(\mathbf{R},\mathcal{V}) = e^{(\beta-\beta')(E(\mathbf{R}) + P\mathcal{V}) } \f$.
45 : 3. Change of pressure from P to P' at constant temperature. That is to say, from a simulation performed in the NPT (isothermal-isobaric) ensemble, obtain an ensemble average in the NP'T ensemble. The weights in this case are \f$ w(\mathbf{R},\mathcal{V}) = e^{\beta (P - P') \mathcal{V}} \f$.
46 : 4. Change of temperature and pressure from T,P to T',P'. That is to say, from a simulation performed in the NPT (isothermal-isobaric) ensemble, obtain an ensemble average in the NP'T' ensemble. The weights in this case are \f$ w(\mathbf{R},\mathcal{V}) = e^{(\beta-\beta')E(\mathbf{R}) + (\beta P - \beta' P') \mathcal{V}} \f$.
47 :
48 : These weights can be used in any action that computes ensemble averages.
49 : For example this action can be used in tandem with \ref HISTOGRAM or \ref AVERAGE.
50 :
51 :
52 : The above equation is often impractical since the overlap between the distributions of energy and volume at different temperatures and pressures is only significant for neighboring temperatures and pressures.
53 : For this reason an unbiased simulation is of little use to reweight at different temperatures and/or pressures.
54 : A successful approach has been altering the probability of observing a configuration in order to increase this overlap \cite wanglandau.
55 : This is done through a bias potential \f$ V(\mathbf{s}) \f$ where \f$ \mathbf{s} \f$ is a set of CVs, that often is the energy (and possibly the volume).
56 : In order to calculate ensemble averages, also the effect of this bias must be taken into account.
57 : The ensemble average of the observable in the ensemble \f$ \xi' \f$ can be calculated from a biased simulation performed in the ensemble \f$\xi\f$ with bias \f$ V(\mathbf{s}) \f$ using:
58 : \f[
59 : \langle O(\mathbf{R},\mathcal{V}) \rangle_{\xi'} = \frac{\langle O(\mathbf{R},\mathcal{V}) w (\mathbf{R},\mathcal{V}) e^{\beta V(\mathbf{s})} \rangle_{\xi,V}}
60 : {\langle w (\mathbf{R},\mathcal{V}) e^{\beta V(\mathbf{s})} \rangle_{\xi,V}}
61 : \f]
62 : where \f$\langle \cdot \rangle_{\xi,V}\f$ is a mean value in the biased ensemble with static bias \f$ V(\mathbf{s}) \f$.
63 : Therefore in order to reweight the trajectory at different temperatures and/or pressures one must use the weights calculated by this action \f$ w (\mathbf{R},\mathcal{V}) \f$ together with the weights of \ref REWEIGHT_BIAS (see the examples below).
64 :
65 : The bias potential \f$ V(\mathbf{s}) \f$ can be constructed with \ref METAD using \ref ENERGY as a CV \cite mich+04prl.
66 : More specialized tools are available, for instance using bespoke target distributions such as \ref TD_MULTICANONICAL and \ref TD_MULTITHERMAL_MULTIBARIC \cite Piaggi-PRL-2019 \cite Piaggi-JCP-2019 within \ref VES.
67 : In the latter algorithms the interval of temperatures and pressures in which the trajectory can be reweighted is chosen explicitly.
68 :
69 : \par Examples
70 :
71 : We consider the 4 cases described above.
72 :
73 : The following input can be used to postprocess a molecular dynamics trajectory of a system of 1000 particles run at 500 K and constant volume using a static bias potential.
74 :
75 : \plumedfile
76 : energy: READ FILE=COLVAR VALUES=energy IGNORE_TIME
77 : distance: READ FILE=COLVAR VALUES=distance IGNORE_TIME
78 : mybias: READ FILE=COLVAR VALUES=mybias.bias IGNORE_TIME
79 :
80 : # Shift energy (to avoid numerical issues)
81 : renergy: COMBINE ARG=energy PARAMETERS=-13250 PERIODIC=NO
82 :
83 : # Weights
84 : bias_weights: REWEIGHT_BIAS TEMP=500 ARG=mybias.bias
85 : temp_press_weights: REWEIGHT_TEMP_PRESS TEMP=500 REWEIGHT_TEMP=300 ENERGY=renergy
86 :
87 : # Ensemble average of the distance at 300 K
88 : avg_dist: AVERAGE ARG=distance LOGWEIGHTS=bias_weights,temp_press_weights
89 :
90 : PRINT ARG=avg_dist FILE=COLVAR_REWEIGHT STRIDE=1
91 : \endplumedfile
92 :
93 : Clearly, in performing the analysis above we would read from the potential energy, a distance, and the value of the bias potential from a COLVAR file like the one shown below. We would then be able
94 : to calculate the ensemble average of the distance at 300 K.
95 :
96 : \auxfile{COLVAR}
97 : #! FIELDS time energy volume mybias.bias distance
98 : 10000.000000 -13133.769283 7.488921 63.740530 0.10293
99 : 10001.000000 -13200.239722 7.116548 36.691988 0.16253
100 : 10002.000000 -13165.108850 7.202273 44.408815 0.17625
101 : \endauxfile
102 :
103 : The next three inputs can be used to postprocess a molecular dynamics trajectory of a system of 1000 particles run at 500 K and 1 bar using a static bias potential.
104 :
105 : We read from a file COLVAR the potential energy, the volume, and the value of the bias potential and calculate the ensemble average of the (particle) density at 300 K and 1 bar (the simulation temperature was 500 K).
106 :
107 : \plumedfile
108 : energy: READ FILE=COLVAR VALUES=energy IGNORE_TIME
109 : volume: READ FILE=COLVAR VALUES=volume IGNORE_TIME
110 : mybias: READ FILE=COLVAR VALUES=mybias.bias IGNORE_TIME
111 :
112 : # Shift energy and volume (to avoid numerical issues)
113 : rvol: COMBINE ARG=volume PARAMETERS=7.8 PERIODIC=NO
114 : renergy: COMBINE ARG=energy PARAMETERS=-13250 PERIODIC=NO
115 :
116 : # Weights
117 : bias_weights: REWEIGHT_BIAS TEMP=500 ARG=mybias.bias
118 : temp_press_weights: REWEIGHT_TEMP_PRESS TEMP=500 REWEIGHT_TEMP=300 PRESSURE=0.06022140857 ENERGY=renergy VOLUME=rvol
119 :
120 : # Ensemble average of the volume at 300 K
121 : avg_vol: AVERAGE ARG=volume LOGWEIGHTS=bias_weights,temp_press_weights
122 : # Ensemble average of the density at 300 K
123 : avg_density: CUSTOM ARG=avg_vol FUNC=1000/x PERIODIC=NO
124 :
125 : PRINT ARG=avg_density FILE=COLVAR_REWEIGHT STRIDE=1
126 : \endplumedfile
127 :
128 : In the next example we calculate the ensemble average of the (particle) density at 500 K and 300 MPa (the simulation pressure was 1 bar).
129 :
130 : \plumedfile
131 : volume: READ FILE=COLVAR VALUES=volume IGNORE_TIME
132 : mybias: READ FILE=COLVAR VALUES=mybias.bias IGNORE_TIME
133 :
134 : # Shift volume (to avoid numerical issues)
135 : rvol: COMBINE ARG=volume PARAMETERS=7.8 PERIODIC=NO
136 :
137 : # Weights
138 : bias_weights: REWEIGHT_BIAS TEMP=500 ARG=mybias.bias
139 : temp_press_weights: REWEIGHT_TEMP_PRESS TEMP=500 PRESSURE=0.06022140857 REWEIGHT_PRESSURE=180.66422571 VOLUME=volume
140 :
141 : # Ensemble average of the volume at 300 K and 300 MPa
142 : avg_vol: AVERAGE ARG=volume LOGWEIGHTS=bias_weights,temp_press_weights
143 : # Ensemble average of the density at 300 K and 300 MPa
144 : avg_density: CUSTOM ARG=avg_vol FUNC=1000/x PERIODIC=NO
145 :
146 : PRINT ARG=avg_density FILE=COLVAR_REWEIGHT STRIDE=1
147 : \endplumedfile
148 :
149 :
150 : In this final example we calculate the ensemble average of the (particle) density at 300 K and 300 MPa (the simulation temperature and pressure were 500 K and 1 bar).
151 :
152 : \plumedfile
153 : energy: READ FILE=COLVAR VALUES=energy IGNORE_TIME
154 : volume: READ FILE=COLVAR VALUES=volume IGNORE_TIME
155 : mybias: READ FILE=COLVAR VALUES=mybias.bias IGNORE_TIME
156 :
157 : # Shift energy and volume (to avoid numerical issues)
158 : rvol: COMBINE ARG=volume PARAMETERS=7.8 PERIODIC=NO
159 : renergy: COMBINE ARG=energy PARAMETERS=-13250 PERIODIC=NO
160 :
161 : # Weights
162 : bias_weights: REWEIGHT_BIAS TEMP=500 ARG=mybias.bias
163 : temp_press_weights: REWEIGHT_TEMP_PRESS TEMP=500 REWEIGHT_TEMP=300 PRESSURE=0.06022140857 REWEIGHT_PRESSURE=180.66422571 ENERGY=renergy VOLUME=rvol
164 :
165 : # Ensemble average of the volume at 300 K and 300 MPa
166 : avg_vol: AVERAGE ARG=volume LOGWEIGHTS=bias_weights,temp_press_weights
167 : # Ensemble average of the density at 300 K and 300 MPa
168 : avg_density: CUSTOM ARG=avg_vol FUNC=1000/x PERIODIC=NO
169 :
170 : PRINT ARG=avg_density FILE=COLVAR_REWEIGHT STRIDE=1
171 : \endplumedfile
172 :
173 : */
174 : //+ENDPLUMEDOC
175 :
176 : namespace PLMD {
177 : namespace bias {
178 :
179 : class ReweightTemperaturePressure : public ReweightBase {
180 : private:
181 : ///
182 : double rpress_, press_, rtemp_;
183 : std::vector<Value*> myenergy, myvol;
184 : public:
185 : static void registerKeywords(Keywords&);
186 : explicit ReweightTemperaturePressure(const ActionOptions&ao);
187 : double getLogWeight() override;
188 : };
189 :
190 10427 : PLUMED_REGISTER_ACTION(ReweightTemperaturePressure,"REWEIGHT_TEMP_PRESS")
191 :
192 5 : void ReweightTemperaturePressure::registerKeywords(Keywords& keys ) {
193 5 : ReweightBase::registerKeywords( keys );
194 5 : keys.remove("ARG");
195 10 : keys.add("optional","ENERGY","Energy");
196 10 : keys.add("optional","VOLUME","Volume");
197 15 : keys.add("optional","REWEIGHT_PRESSURE","Reweighting pressure");
198 10 : keys.add("optional","PRESSURE","The system pressure");
199 10 : keys.add("optional","REWEIGHT_TEMP","Reweighting temperature");
200 5 : }
201 :
202 4 : ReweightTemperaturePressure::ReweightTemperaturePressure(const ActionOptions&ao):
203 : Action(ao),
204 4 : ReweightBase(ao)
205 : {
206 : // Initialize to not defined (negative)
207 4 : rpress_=-1;
208 4 : press_=-1;
209 4 : rtemp_=-1;
210 4 : parse("REWEIGHT_PRESSURE",rpress_);
211 4 : parse("PRESSURE",press_);
212 4 : parse("REWEIGHT_TEMP",rtemp_);
213 4 : rtemp_*=plumed.getAtoms().getKBoltzmann();
214 :
215 8 : parseArgumentList("ENERGY",myenergy);
216 4 : if(!myenergy.empty()) {
217 3 : log.printf(" with energies: ");
218 6 : for(unsigned i=0; i<myenergy.size(); i++) log.printf(" %s",myenergy[i]->getName().c_str());
219 3 : log.printf("\n");
220 : }
221 : //requestArguments(myenergy);
222 :
223 8 : parseArgumentList("VOLUME",myvol);
224 4 : if(!myvol.empty()) {
225 3 : log.printf(" with volumes: ");
226 6 : for(unsigned i=0; i<myvol.size(); i++) log.printf(" %s",myvol[i]->getName().c_str());
227 3 : log.printf("\n");
228 : }
229 :
230 : std::vector<Value*> conc;
231 4 : conc.insert(conc.begin(), myenergy.begin(), myenergy.end());
232 4 : conc.insert(conc.end(), myvol.begin(), myvol.end());
233 4 : requestArguments(conc);
234 :
235 : // 4 possible cases
236 : // Case 1) Reweight from T to T' with V=const (canonical)
237 4 : if (rtemp_>=0 && press_<0 && rpress_<0 && !myenergy.empty() && myvol.empty() ) {
238 1 : log.printf(" reweighting simulation from temperature %f to temperature %f at constant volume \n",simtemp/plumed.getAtoms().getKBoltzmann(),rtemp_/plumed.getAtoms().getKBoltzmann() );
239 1 : log.printf(" WARNING: If the simulation is performed at constant pressure add the keywords PRESSURE and VOLUME \n" );
240 : }
241 : // Case 2) Reweight from T to T' with P=const (isothermal-isobaric)
242 3 : else if (rtemp_>=0 && press_>=0 && rpress_<0 && !myenergy.empty() && !myvol.empty() ) log.printf(" reweighting simulation from temperature %f to temperature %f at constant pressure %f \n",simtemp/plumed.getAtoms().getKBoltzmann(),rtemp_/plumed.getAtoms().getKBoltzmann(), press_ );
243 : // Case 3) Reweight from P to P' with T=const (isothermal-isobaric)
244 2 : else if (rtemp_<0 && press_>=0 && rpress_>=0 && myenergy.empty() && !myvol.empty() ) log.printf(" reweighting simulation from pressure %f to pressure %f at constant temperature %f\n",press_,rpress_,simtemp/plumed.getAtoms().getKBoltzmann() );
245 : // Case 4) Reweight from T,P to T',P' (isothermal-isobaric)
246 1 : else if (rtemp_>0 && press_>=0 && rpress_>=0 && !myenergy.empty() && !myvol.empty() ) log.printf(" reweighting simulation from temperature %f and pressure %f to temperature %f and pressure %f \n",simtemp/plumed.getAtoms().getKBoltzmann(), press_, rtemp_/plumed.getAtoms().getKBoltzmann(), rpress_);
247 0 : else error("Combination of ENERGY, VOLUME, REWEIGHT_PRESSURE, PRESSURE and REWEIGHT_TEMP not supported. Please refer to the manual for supported combinations.");
248 4 : }
249 :
250 1001 : double ReweightTemperaturePressure::getLogWeight() {
251 2002 : double energy=0.0; for(unsigned i=0; i<myenergy.size(); ++i) energy+=getArgument(i);
252 2002 : double volume=0.0; for(unsigned i=0; i<myvol.size(); ++i) volume+=getArgument(myenergy.size()+i);
253 : // 4 possible cases
254 : // Case 1) Reweight from T to T' with V=const (canonical)
255 1001 : if (rtemp_>=0 && press_<0 && rpress_<0) return ((1.0/simtemp)- (1.0/rtemp_) )*energy;
256 : // Case 2) Reweight from T to T' with P=const (isothermal-isobaric)
257 1001 : else if (rtemp_>=0 && press_>=0 && rpress_<0) return ((1.0/simtemp)- (1.0/rtemp_) )*energy + ((1.0/simtemp) - (1.0/rtemp_))*press_*volume;
258 : // Case 3) Reweight from P to P' with T=const (isothermal-isobaric)
259 1001 : else if (rtemp_<0 && press_>=0 && rpress_>=0) return (1.0/simtemp)*(press_ - rpress_)*volume;
260 : // Case 4) Reweight from T,P to T',P' (isothermal-isobaric)
261 1001 : else if (rtemp_>0 && press_>=0 && rpress_>=0) return ((1.0/simtemp)- (1.0/rtemp_) )*energy + ((1.0/simtemp)*press_ - (1.0/rtemp_)*rpress_ )*volume;
262 : else return 0;
263 : }
264 :
265 : }
266 : }
|
