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$\mathcal{V}$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 $O(\mathbf{R},\mathcal{V})$ that depends on the atomic coordinates $\mathbf{R}$ and the volume $\mathcal{V}$.
33 : This observable is in practice any collective variable (CV) calculated by PLUMED.
34 : The ensemble average of the observable in an ensemble $\xi'$ can be calculated from a simulation performed in an ensemble $\xi$ using:
35 :
36 : $$
37 : \langle O(\mathbf{R},\mathcal{V}) \rangle_{\xi'} = \frac{\langle O(\mathbf{R},\mathcal{V}) w(\mathbf{R},\mathcal{V}) \rangle_{\xi}}
38 : {\langle w(\mathbf{R},\mathcal{V}) \rangle_{\xi}}
39 : $$
40 :
41 : where $\langle \cdot \rangle_{\xi}$ and $\langle \cdot \rangle_{\xi'}$ are mean values in the simulated and targeted ensemble, respectively, $E(\mathbf{R})$ is the potential energy of the system, and $w (\mathbf{R},\mathcal{V})$ are the appropriate weights to take from $\xi$ to $\xi'$.
42 : This action calculates the weights $ w (\mathbf{R},\mathcal{V})$ and handles 4 different cases:
43 :
44 : 1. Change of temperature from T to T' at constant volume. That is to say, from a simulation performed in the N$\mathcal{V}$T (canonical) ensemble, obtain an ensemble average in the N$\mathcal{V}$T' ensemble. The weights in this case are $ w(\mathbf{R},\mathcal{V}) = e^{(\beta-\beta')E(\mathbf{R})}$ with $\beta$ and $\beta'$ the inverse temperatures.
45 : 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 $w(\mathbf{R},\mathcal{V}) = e^{(\beta-\beta')(E(\mathbf{R}) + P\mathcal{V}) }$.
46 : 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 $w(\mathbf{R},\mathcal{V}) = e^{\beta (P - P') \mathcal{V}}$.
47 : 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 $w(\mathbf{R},\mathcal{V}) = e^{(\beta-\beta')E(\mathbf{R}) + (\beta P - \beta' P') \mathcal{V}}$.
48 :
49 : These weights can be used in any action that computes ensemble averages. For example this action can be used in tandem with [HISTOGRAM](HISTOGRAM.md) or [AVERAGE](AVERAGE.md).
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 that is discussed in the first paper cited below is altering the probability of observing a configuration in order to increase this overlap.
55 : This is done through a bias potential $V(\mathbf{s})$ where $\mathbf{s}$ 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 $\xi'$ can be calculated from a biased simulation performed in the ensemble $\xi$ with bias $V(\mathbf{s})$ using:
58 :
59 : $$
60 : \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}}
61 : {\langle w (\mathbf{R},\mathcal{V}) e^{\beta V(\mathbf{s})} \rangle_{\xi,V}}
62 : $$
63 :
64 : where $\langle \cdot \rangle_{\xi,V}$ is a mean value in the biased ensemble with static bias $V(\mathbf{s})$.
65 : Therefore in order to reweight the trajectory at different temperatures and/or pressures one must use the weights calculated by this action $w (\mathbf{R},\mathcal{V})$ together with the weights of [REWEIGHT_BIAS](REWEIGHT_BIAS.md) (see the examples below).
66 :
67 : The bias potential $V(\mathbf{s})$ can be constructed with [METAD](METAD.md) using [ENERGY](ENERGY.md) as a CV as discussed in the second paper cited below.
68 : The remaining papers cited below discuss more specialized tools that available, that, for instance, use bespoke target distributions such as [TD_MULTICANONICAL](TD_MULTICANONICAL.md) and [TD_MULTITHERMAL_MULTIBARIC](TD_MULTITHERMAL_MULTIBARIC.md) that are implemented within the [VES](module_ves.md).
69 : In the latter algorithms the interval of temperatures and pressures in which the trajectory can be reweighted is chosen explicitly.
70 :
71 : ## Examples
72 :
73 : We consider the 4 cases described above in the following examples
74 :
75 : 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.
76 :
77 : ```plumed
78 : #SETTINGS INPUTFILES=regtest/landmarks/rt-reweight-temp-press/COLVAR-READ
79 :
80 : energy: READ FILE=regtest/landmarks/rt-reweight-temp-press/COLVAR-READ VALUES=energy IGNORE_TIME
81 : distance: READ FILE=regtest/landmarks/rt-reweight-temp-press/COLVAR-READ VALUES=distance IGNORE_TIME
82 : mybias: READ FILE=regtest/landmarks/rt-reweight-temp-press/COLVAR-READ VALUES=mybias.bias IGNORE_TIME
83 :
84 : # Shift energy (to avoid numerical issues)
85 : renergy: COMBINE ARG=energy PARAMETERS=-13250 PERIODIC=NO
86 :
87 : # Weights
88 : bias_weights: REWEIGHT_BIAS TEMP=500 ARG=mybias.bias
89 : temp_press_weights: REWEIGHT_TEMP_PRESS TEMP=500 REWEIGHT_TEMP=300 ENERGY=renergy
90 :
91 : # Ensemble average of the distance at 300 K
92 : avg_dist: AVERAGE ARG=distance LOGWEIGHTS=bias_weights,temp_press_weights
93 :
94 : PRINT ARG=avg_dist FILE=COLVAR_REWEIGHT STRIDE=1
95 : ```
96 :
97 : Clearly, in performing the analysis above we read from the potential energy, a distance, and the value of the bias potential from a COLVAR file. Having done the reweighting
98 : we can then calculate the ensemble average of the distance at 300 K.
99 :
100 : 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.
101 :
102 : 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).
103 :
104 : ```plumed
105 : #SETTINGS INPUTFILES=regtest/landmarks/rt-reweight-temp-press/COLVAR-READ
106 :
107 : energy: READ FILE=regtest/landmarks/rt-reweight-temp-press/COLVAR-READ VALUES=energy IGNORE_TIME
108 : volume: READ FILE=regtest/landmarks/rt-reweight-temp-press/COLVAR-READ VALUES=volume IGNORE_TIME
109 : mybias: READ FILE=regtest/landmarks/rt-reweight-temp-press/COLVAR-READ VALUES=mybias.bias IGNORE_TIME
110 :
111 : # Shift energy and volume (to avoid numerical issues)
112 : rvol: COMBINE ARG=volume PARAMETERS=7.8 PERIODIC=NO
113 : renergy: COMBINE ARG=energy PARAMETERS=-13250 PERIODIC=NO
114 :
115 : # Weights
116 : bias_weights: REWEIGHT_BIAS TEMP=500 ARG=mybias.bias
117 : temp_press_weights: REWEIGHT_TEMP_PRESS TEMP=500 REWEIGHT_TEMP=300 PRESSURE=0.06022140857 ENERGY=renergy VOLUME=rvol
118 :
119 : # Ensemble average of the volume at 300 K
120 : avg_vol: AVERAGE ARG=volume LOGWEIGHTS=bias_weights,temp_press_weights
121 : # Ensemble average of the density at 300 K
122 : avg_density: CUSTOM ARG=avg_vol FUNC=1000/x PERIODIC=NO
123 :
124 : PRINT ARG=avg_density FILE=COLVAR_REWEIGHT STRIDE=1
125 : ```
126 :
127 : 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).
128 :
129 : ```plumed
130 : #SETTINGS INPUTFILES=regtest/landmarks/rt-reweight-temp-press/COLVAR-READ
131 :
132 : volume: READ FILE=regtest/landmarks/rt-reweight-temp-press/COLVAR-READ VALUES=volume IGNORE_TIME
133 : mybias: READ FILE=regtest/landmarks/rt-reweight-temp-press/COLVAR-READ VALUES=mybias.bias IGNORE_TIME
134 :
135 : # Shift volume (to avoid numerical issues)
136 : rvol: COMBINE ARG=volume PARAMETERS=7.8 PERIODIC=NO
137 :
138 : # Weights
139 : bias_weights: REWEIGHT_BIAS TEMP=500 ARG=mybias.bias
140 : temp_press_weights: REWEIGHT_TEMP_PRESS TEMP=500 PRESSURE=0.06022140857 REWEIGHT_PRESSURE=180.66422571 VOLUME=volume
141 :
142 : # Ensemble average of the volume at 300 K and 300 MPa
143 : avg_vol: AVERAGE ARG=volume LOGWEIGHTS=bias_weights,temp_press_weights
144 : # Ensemble average of the density at 300 K and 300 MPa
145 : avg_density: CUSTOM ARG=avg_vol FUNC=1000/x PERIODIC=NO
146 :
147 : PRINT ARG=avg_density FILE=COLVAR_REWEIGHT STRIDE=1
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 : ```plumed
153 : #SETTINGS INPUTFILES=regtest/landmarks/rt-reweight-temp-press/COLVAR-READ
154 :
155 : energy: READ FILE=regtest/landmarks/rt-reweight-temp-press/COLVAR-READ VALUES=energy IGNORE_TIME
156 : volume: READ FILE=regtest/landmarks/rt-reweight-temp-press/COLVAR-READ VALUES=volume IGNORE_TIME
157 : mybias: READ FILE=regtest/landmarks/rt-reweight-temp-press/COLVAR-READ VALUES=mybias.bias IGNORE_TIME
158 :
159 : # Shift energy and volume (to avoid numerical issues)
160 : rvol: COMBINE ARG=volume PARAMETERS=7.8 PERIODIC=NO
161 : renergy: COMBINE ARG=energy PARAMETERS=-13250 PERIODIC=NO
162 :
163 : # Weights
164 : bias_weights: REWEIGHT_BIAS TEMP=500 ARG=mybias.bias
165 : temp_press_weights: REWEIGHT_TEMP_PRESS TEMP=500 REWEIGHT_TEMP=300 PRESSURE=0.06022140857 REWEIGHT_PRESSURE=180.66422571 ENERGY=renergy VOLUME=rvol
166 :
167 : # Ensemble average of the volume at 300 K and 300 MPa
168 : avg_vol: AVERAGE ARG=volume LOGWEIGHTS=bias_weights,temp_press_weights
169 : # Ensemble average of the density at 300 K and 300 MPa
170 : avg_density: CUSTOM ARG=avg_vol FUNC=1000/x PERIODIC=NO
171 :
172 : PRINT ARG=avg_density FILE=COLVAR_REWEIGHT STRIDE=1
173 : ```
174 :
175 : */
176 : //+ENDPLUMEDOC
177 :
178 : namespace PLMD {
179 : namespace bias {
180 :
181 : class ReweightTemperaturePressure : public ReweightBase {
182 : private:
183 : ///
184 : double rpress_, press_, rtemp_;
185 : std::vector<Value*> myenergy, myvol;
186 : public:
187 : static void registerKeywords(Keywords&);
188 : explicit ReweightTemperaturePressure(const ActionOptions&ao);
189 : double getLogWeight() override;
190 : };
191 :
192 : PLUMED_REGISTER_ACTION(ReweightTemperaturePressure,"REWEIGHT_TEMP_PRESS")
193 :
194 6 : void ReweightTemperaturePressure::registerKeywords(Keywords& keys ) {
195 6 : ReweightBase::registerKeywords( keys );
196 12 : keys.addInputKeyword("optional","ENERGY","scalar","Energy");
197 12 : keys.addInputKeyword("optional","VOLUME","scalar","Volume");
198 6 : keys.add("optional","REWEIGHT_PRESSURE","Reweighting pressure");
199 6 : keys.add("optional","PRESSURE","The system pressure");
200 6 : keys.add("optional","REWEIGHT_TEMP","Reweighting temperature");
201 12 : keys.setValueDescription("scalar","the weight to use for this frame to determine its contribution at a different temperature/pressure");
202 6 : keys.addDOI("10.1103/PhysRevLett.86.2050");
203 6 : keys.addDOI("10.1103/PhysRevLett.92.170601");
204 6 : keys.addDOI("10.1103/PhysRevLett.122.050601");
205 6 : keys.addDOI("10.1063/1.5102104");
206 6 : }
207 :
208 4 : ReweightTemperaturePressure::ReweightTemperaturePressure(const ActionOptions&ao):
209 : Action(ao),
210 4 : ReweightBase(ao) {
211 : // Initialize to not defined (negative)
212 4 : rpress_=-1;
213 4 : press_=-1;
214 4 : rtemp_=-1;
215 4 : parse("REWEIGHT_PRESSURE",rpress_);
216 4 : parse("PRESSURE",press_);
217 4 : parse("REWEIGHT_TEMP",rtemp_);
218 4 : rtemp_*=getKBoltzmann();
219 :
220 8 : parseArgumentList("ENERGY",myenergy);
221 4 : if(!myenergy.empty()) {
222 3 : log.printf(" with energies: ");
223 6 : for(unsigned i=0; i<myenergy.size(); i++) {
224 3 : log.printf(" %s",myenergy[i]->getName().c_str());
225 : }
226 3 : log.printf("\n");
227 : }
228 : //requestArguments(myenergy);
229 :
230 8 : parseArgumentList("VOLUME",myvol);
231 4 : if(!myvol.empty()) {
232 3 : log.printf(" with volumes: ");
233 6 : for(unsigned i=0; i<myvol.size(); i++) {
234 3 : log.printf(" %s",myvol[i]->getName().c_str());
235 : }
236 3 : log.printf("\n");
237 : }
238 :
239 : std::vector<Value*> conc;
240 4 : conc.insert(conc.begin(), myenergy.begin(), myenergy.end());
241 4 : conc.insert(conc.end(), myvol.begin(), myvol.end());
242 4 : requestArguments(conc);
243 :
244 : // 4 possible cases
245 : // Case 1) Reweight from T to T' with V=const (canonical)
246 4 : if (rtemp_>=0 && press_<0 && rpress_<0 && !myenergy.empty() && myvol.empty() ) {
247 1 : log.printf(" reweighting simulation from temperature %f to temperature %f at constant volume \n",simtemp/getKBoltzmann(),rtemp_/getKBoltzmann() );
248 1 : log.printf(" WARNING: If the simulation is performed at constant pressure add the keywords PRESSURE and VOLUME \n" );
249 : }
250 : // Case 2) Reweight from T to T' with P=const (isothermal-isobaric)
251 3 : else if (rtemp_>=0 && press_>=0 && rpress_<0 && !myenergy.empty() && !myvol.empty() ) {
252 1 : log.printf(" reweighting simulation from temperature %f to temperature %f at constant pressure %f \n",simtemp/getKBoltzmann(),rtemp_/getKBoltzmann(), press_ );
253 : }
254 : // Case 3) Reweight from P to P' with T=const (isothermal-isobaric)
255 2 : else if (rtemp_<0 && press_>=0 && rpress_>=0 && myenergy.empty() && !myvol.empty() ) {
256 1 : log.printf(" reweighting simulation from pressure %f to pressure %f at constant temperature %f\n",press_,rpress_,simtemp/getKBoltzmann() );
257 : }
258 : // Case 4) Reweight from T,P to T',P' (isothermal-isobaric)
259 1 : else if (rtemp_>0 && press_>=0 && rpress_>=0 && !myenergy.empty() && !myvol.empty() ) {
260 1 : log.printf(" reweighting simulation from temperature %f and pressure %f to temperature %f and pressure %f \n",simtemp/getKBoltzmann(), press_, rtemp_/getKBoltzmann(), rpress_);
261 : } else {
262 0 : error("Combination of ENERGY, VOLUME, REWEIGHT_PRESSURE, PRESSURE and REWEIGHT_TEMP not supported. Please refer to the manual for supported combinations.");
263 : }
264 4 : }
265 :
266 1001 : double ReweightTemperaturePressure::getLogWeight() {
267 : double energy=0.0;
268 2002 : for(unsigned i=0; i<myenergy.size(); ++i) {
269 1001 : energy+=getArgument(i);
270 : }
271 : double volume=0.0;
272 2002 : for(unsigned i=0; i<myvol.size(); ++i) {
273 1001 : volume+=getArgument(myenergy.size()+i);
274 : }
275 : // 4 possible cases
276 : // Case 1) Reweight from T to T' with V=const (canonical)
277 1001 : if (rtemp_>=0 && press_<0 && rpress_<0) {
278 0 : return ((1.0/simtemp)- (1.0/rtemp_) )*energy;
279 : }
280 : // Case 2) Reweight from T to T' with P=const (isothermal-isobaric)
281 1001 : else if (rtemp_>=0 && press_>=0 && rpress_<0) {
282 0 : return ((1.0/simtemp)- (1.0/rtemp_) )*energy + ((1.0/simtemp) - (1.0/rtemp_))*press_*volume;
283 : }
284 : // Case 3) Reweight from P to P' with T=const (isothermal-isobaric)
285 1001 : else if (rtemp_<0 && press_>=0 && rpress_>=0) {
286 0 : return (1.0/simtemp)*(press_ - rpress_)*volume;
287 : }
288 : // Case 4) Reweight from T,P to T',P' (isothermal-isobaric)
289 1001 : else if (rtemp_>0 && press_>=0 && rpress_>=0) {
290 1001 : return ((1.0/simtemp)- (1.0/rtemp_) )*energy + ((1.0/simtemp)*press_ - (1.0/rtemp_)*rpress_ )*volume;
291 : } else {
292 : return 0;
293 : }
294 : }
295 :
296 : }
297 : }
|
