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 "ReweightBase.h"
24 :
25 : //+PLUMEDOC REWEIGHTING REWEIGHT_TEMP_PRESS
26 : /*
27 : Calculate weights for ensemble averages at temperatures and/or pressures different than those used in your original simulation.
28 :
29 : 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
30 : contained in trajectories and obtain ensemble averages at different temperatures and/or pressures.
31 :
32 : 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$.
33 : This observable is in practice any collective variable (CV) calculated by Plumed.
34 : 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:
35 : \f[
36 : \langle O(\mathbf{R},\mathcal{V}) \rangle_{\xi'} = \frac{\langle O(\mathbf{R},\mathcal{V}) w(\mathbf{R},\mathcal{V}) \rangle_{\xi}}
37 : {\langle w(\mathbf{R},\mathcal{V}) \rangle_{\xi}}
38 : \f]
39 : 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$.
40 : This action calculates the weights \f$ w (\mathbf{R},\mathcal{V}) \f$ and handles 4 different cases:
41 : 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.
42 : 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$.
43 : 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$.
44 : 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$.
45 :
46 : These weights can be used in any action that computes ensemble averages.
47 : For example this action can be used in tandem with \ref HISTOGRAM or \ref AVERAGE.
48 :
49 :
50 : 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.
51 : For this reason an unbiased simulation is of little use to reweight at different temperatures and/or pressures.
52 : A successful approach has been altering the probability of observing a configuration in order to increase this overlap \cite wanglandau.
53 : 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).
54 : In order to calculate ensemble averages, also the effect of this bias must be taken into account.
55 : 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:
56 : \f[
57 : \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}}
58 : {\langle w (\mathbf{R},\mathcal{V}) e^{\beta V(\mathbf{s})} \rangle_{\xi,V}}
59 : \f]
60 : where \f$\langle \cdot \rangle_{\xi,V}\f$ is a mean value in the biased ensemble with static bias \f$ V(\mathbf{s}) \f$.
61 : 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).
62 :
63 : The bias potential \f$ V(\mathbf{s}) \f$ can be constructed with \ref METAD using \ref ENERGY as a CV \cite mich+04prl.
64 : 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.
65 : In the latter algorithms the interval of temperatures and pressures in which the trajectory can be reweighted is chosen explicitly.
66 :
67 : \par Examples
68 :
69 : We consider the 4 cases described above.
70 :
71 : 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.
72 :
73 : \plumedfile
74 : energy: READ FILE=COLVAR VALUES=energy IGNORE_TIME
75 : distance: READ FILE=COLVAR VALUES=distance IGNORE_TIME
76 : mybias: READ FILE=COLVAR VALUES=mybias.bias IGNORE_TIME
77 :
78 : # Shift energy (to avoid numerical issues)
79 : renergy: COMBINE ARG=energy PARAMETERS=-13250 PERIODIC=NO
80 :
81 : # Weights
82 : bias_weights: REWEIGHT_BIAS TEMP=500 ARG=mybias.bias
83 : temp_press_weights: REWEIGHT_TEMP_PRESS TEMP=500 REWEIGHT_TEMP=300 ENERGY=renergy
84 :
85 : # Ensemble average of the distance at 300 K
86 : avg_dist: AVERAGE ARG=distance LOGWEIGHTS=bias_weights,temp_press_weights
87 :
88 : PRINT ARG=avg_dist FILE=COLVAR_REWEIGHT STRIDE=1
89 : \endplumedfile
90 :
91 : 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
92 : to calculate the ensemble average of the distance at 300 K.
93 :
94 : \auxfile{COLVAR}
95 : #! FIELDS time energy volume mybias.bias distance
96 : 10000.000000 -13133.769283 7.488921 63.740530 0.10293
97 : 10001.000000 -13200.239722 7.116548 36.691988 0.16253
98 : 10002.000000 -13165.108850 7.202273 44.408815 0.17625
99 : \endauxfile
100 :
101 : 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.
102 :
103 : 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).
104 :
105 : \plumedfile
106 : energy: READ FILE=COLVAR VALUES=energy IGNORE_TIME
107 : volume: READ FILE=COLVAR VALUES=volume IGNORE_TIME
108 : mybias: READ FILE=COLVAR VALUES=mybias.bias IGNORE_TIME
109 :
110 : # Shift energy and volume (to avoid numerical issues)
111 : rvol: COMBINE ARG=volume PARAMETERS=7.8 PERIODIC=NO
112 : renergy: COMBINE ARG=energy PARAMETERS=-13250 PERIODIC=NO
113 :
114 : # Weights
115 : bias_weights: REWEIGHT_BIAS TEMP=500 ARG=mybias.bias
116 : temp_press_weights: REWEIGHT_TEMP_PRESS TEMP=500 REWEIGHT_TEMP=300 PRESSURE=0.06022140857 ENERGY=renergy VOLUME=rvol
117 :
118 : # Ensemble average of the volume at 300 K
119 : avg_vol: AVERAGE ARG=volume LOGWEIGHTS=bias_weights,temp_press_weights
120 : # Ensemble average of the density at 300 K
121 : avg_density: CUSTOM ARG=avg_vol FUNC=1000/x PERIODIC=NO
122 :
123 : PRINT ARG=avg_density FILE=COLVAR_REWEIGHT STRIDE=1
124 : \endplumedfile
125 :
126 : 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).
127 :
128 : \plumedfile
129 : volume: READ FILE=COLVAR VALUES=volume IGNORE_TIME
130 : mybias: READ FILE=COLVAR VALUES=mybias.bias IGNORE_TIME
131 :
132 : # Shift volume (to avoid numerical issues)
133 : rvol: COMBINE ARG=volume PARAMETERS=7.8 PERIODIC=NO
134 :
135 : # Weights
136 : bias_weights: REWEIGHT_BIAS TEMP=500 ARG=mybias.bias
137 : temp_press_weights: REWEIGHT_TEMP_PRESS TEMP=500 PRESSURE=0.06022140857 REWEIGHT_PRESSURE=180.66422571 VOLUME=volume
138 :
139 : # Ensemble average of the volume at 300 K and 300 MPa
140 : avg_vol: AVERAGE ARG=volume LOGWEIGHTS=bias_weights,temp_press_weights
141 : # Ensemble average of the density at 300 K and 300 MPa
142 : avg_density: CUSTOM ARG=avg_vol FUNC=1000/x PERIODIC=NO
143 :
144 : PRINT ARG=avg_density FILE=COLVAR_REWEIGHT STRIDE=1
145 : \endplumedfile
146 :
147 :
148 : 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).
149 :
150 : \plumedfile
151 : energy: READ FILE=COLVAR VALUES=energy IGNORE_TIME
152 : volume: READ FILE=COLVAR VALUES=volume IGNORE_TIME
153 : mybias: READ FILE=COLVAR VALUES=mybias.bias IGNORE_TIME
154 :
155 : # Shift energy and volume (to avoid numerical issues)
156 : rvol: COMBINE ARG=volume PARAMETERS=7.8 PERIODIC=NO
157 : renergy: COMBINE ARG=energy PARAMETERS=-13250 PERIODIC=NO
158 :
159 : # Weights
160 : bias_weights: REWEIGHT_BIAS TEMP=500 ARG=mybias.bias
161 : temp_press_weights: REWEIGHT_TEMP_PRESS TEMP=500 REWEIGHT_TEMP=300 PRESSURE=0.06022140857 REWEIGHT_PRESSURE=180.66422571 ENERGY=renergy VOLUME=rvol
162 :
163 : # Ensemble average of the volume at 300 K and 300 MPa
164 : avg_vol: AVERAGE ARG=volume LOGWEIGHTS=bias_weights,temp_press_weights
165 : # Ensemble average of the density at 300 K and 300 MPa
166 : avg_density: CUSTOM ARG=avg_vol FUNC=1000/x PERIODIC=NO
167 :
168 : PRINT ARG=avg_density FILE=COLVAR_REWEIGHT STRIDE=1
169 : \endplumedfile
170 :
171 : */
172 : //+ENDPLUMEDOC
173 :
174 : namespace PLMD {
175 : namespace bias {
176 :
177 : class ReweightTemperaturePressure : public ReweightBase {
178 : private:
179 : ///
180 : double rpress_, press_, rtemp_;
181 : std::vector<Value*> myenergy, myvol;
182 : public:
183 : static void registerKeywords(Keywords&);
184 : explicit ReweightTemperaturePressure(const ActionOptions&ao);
185 : double getLogWeight() override;
186 : };
187 :
188 : PLUMED_REGISTER_ACTION(ReweightTemperaturePressure,"REWEIGHT_TEMP_PRESS")
189 :
190 6 : void ReweightTemperaturePressure::registerKeywords(Keywords& keys ) {
191 6 : ReweightBase::registerKeywords( keys );
192 6 : keys.remove("ARG");
193 12 : keys.add("optional","ENERGY","Energy");
194 12 : keys.add("optional","VOLUME","Volume");
195 12 : keys.add("optional","REWEIGHT_PRESSURE","Reweighting pressure");
196 12 : keys.add("optional","PRESSURE","The system pressure");
197 12 : keys.add("optional","REWEIGHT_TEMP","Reweighting temperature");
198 6 : keys.setValueDescription("the weight to use for this frame to determine its contribution at a different temperature/pressure");
199 6 : }
200 :
201 4 : ReweightTemperaturePressure::ReweightTemperaturePressure(const ActionOptions&ao):
202 : Action(ao),
203 4 : ReweightBase(ao)
204 : {
205 : // Initialize to not defined (negative)
206 4 : rpress_=-1;
207 4 : press_=-1;
208 4 : rtemp_=-1;
209 4 : parse("REWEIGHT_PRESSURE",rpress_);
210 4 : parse("PRESSURE",press_);
211 4 : parse("REWEIGHT_TEMP",rtemp_);
212 4 : rtemp_*=getKBoltzmann();
213 :
214 8 : parseArgumentList("ENERGY",myenergy);
215 4 : if(!myenergy.empty()) {
216 3 : log.printf(" with energies: ");
217 6 : for(unsigned i=0; i<myenergy.size(); i++) log.printf(" %s",myenergy[i]->getName().c_str());
218 3 : log.printf("\n");
219 : }
220 : //requestArguments(myenergy);
221 :
222 8 : parseArgumentList("VOLUME",myvol);
223 4 : if(!myvol.empty()) {
224 3 : log.printf(" with volumes: ");
225 6 : for(unsigned i=0; i<myvol.size(); i++) log.printf(" %s",myvol[i]->getName().c_str());
226 3 : log.printf("\n");
227 : }
228 :
229 : std::vector<Value*> conc;
230 4 : conc.insert(conc.begin(), myenergy.begin(), myenergy.end());
231 4 : conc.insert(conc.end(), myvol.begin(), myvol.end());
232 4 : requestArguments(conc);
233 :
234 : // 4 possible cases
235 : // Case 1) Reweight from T to T' with V=const (canonical)
236 4 : if (rtemp_>=0 && press_<0 && rpress_<0 && !myenergy.empty() && myvol.empty() ) {
237 1 : log.printf(" reweighting simulation from temperature %f to temperature %f at constant volume \n",simtemp/getKBoltzmann(),rtemp_/getKBoltzmann() );
238 1 : log.printf(" WARNING: If the simulation is performed at constant pressure add the keywords PRESSURE and VOLUME \n" );
239 : }
240 : // Case 2) Reweight from T to T' with P=const (isothermal-isobaric)
241 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/getKBoltzmann(),rtemp_/getKBoltzmann(), press_ );
242 : // Case 3) Reweight from P to P' with T=const (isothermal-isobaric)
243 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/getKBoltzmann() );
244 : // Case 4) Reweight from T,P to T',P' (isothermal-isobaric)
245 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/getKBoltzmann(), press_, rtemp_/getKBoltzmann(), rpress_);
246 0 : else error("Combination of ENERGY, VOLUME, REWEIGHT_PRESSURE, PRESSURE and REWEIGHT_TEMP not supported. Please refer to the manual for supported combinations.");
247 4 : }
248 :
249 1001 : double ReweightTemperaturePressure::getLogWeight() {
250 2002 : double energy=0.0; for(unsigned i=0; i<myenergy.size(); ++i) energy+=getArgument(i);
251 2002 : double volume=0.0; for(unsigned i=0; i<myvol.size(); ++i) volume+=getArgument(myenergy.size()+i);
252 : // 4 possible cases
253 : // Case 1) Reweight from T to T' with V=const (canonical)
254 1001 : if (rtemp_>=0 && press_<0 && rpress_<0) return ((1.0/simtemp)- (1.0/rtemp_) )*energy;
255 : // Case 2) Reweight from T to T' with P=const (isothermal-isobaric)
256 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;
257 : // Case 3) Reweight from P to P' with T=const (isothermal-isobaric)
258 1001 : else if (rtemp_<0 && press_>=0 && rpress_>=0) return (1.0/simtemp)*(press_ - rpress_)*volume;
259 : // Case 4) Reweight from T,P to T',P' (isothermal-isobaric)
260 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;
261 : else return 0;
262 : }
263 :
264 : }
265 : }
|
