Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2019 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-arXiv-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 : We read from a file COLVAR the potential energy, a distance, and the value of the bias potential and calculate the ensemble average of the distance at 300 K.
75 :
76 : \plumedfile
77 : energy: READ FILE=COLVAR VALUES=energy IGNORE_TIME
78 : distance: READ FILE=COLVAR VALUES=distance IGNORE_TIME
79 : mybias: READ FILE=COLVAR VALUES=mybias.bias IGNORE_TIME
80 :
81 : # Shift energy (to avoid numerical issues)
82 : renergy: COMBINE ARG=energy PARAMETERS=-13250 PERIODIC=NO
83 :
84 : # Weights
85 : bias_weights: REWEIGHT_BIAS TEMP=500 ARG=mybias.bias
86 : temp_press_weights: REWEIGHT_TEMP_PRESS TEMP=500 REWEIGHT_TEMP=300 ENERGY=renergy
87 :
88 : # Ensemble average of the distance at 300 K
89 : avg_dist: AVERAGE ARG=distance LOGWEIGHTS=bias_weights,temp_press_weights
90 :
91 : PRINT ARG=avg_dist FILE=COLVAR_REWEIGHT STRIDE=1
92 : \endplumedfile
93 :
94 : 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.
95 :
96 : 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).
97 :
98 : \plumedfile
99 : energy: READ FILE=COLVAR VALUES=energy IGNORE_TIME
100 : volume: READ FILE=COLVAR VALUES=volume IGNORE_TIME
101 : mybias: READ FILE=COLVAR VALUES=mybias.bias IGNORE_TIME
102 :
103 : # Shift energy and volume (to avoid numerical issues)
104 : rvol: COMBINE ARG=vol PARAMETERS=7.8 PERIODIC=NO
105 : renergy: COMBINE ARG=energy PARAMETERS=-13250 PERIODIC=NO
106 :
107 : # Weights
108 : bias_weights: REWEIGHT_BIAS TEMP=500 ARG=mybias.bias
109 : temp_press_weights: REWEIGHT_TEMP_PRESS TEMP=500 REWEIGHT_TEMP=300 PRESSURE=0.06022140857 ENERGY=renergy VOLUME=rvolume
110 :
111 : # Ensemble average of the volume at 300 K
112 : avg_vol: AVERAGE ARG=volume LOGWEIGHTS=bias_weights,temp_press_weights
113 : # Ensemble average of the density at 300 K
114 : avg_density: CUSTOM ARG=avg_vol FUNC=1000/x PERIODIC=NO
115 :
116 : PRINT ARG=avg_density FILE=COLVAR_REWEIGHT STRIDE=1
117 : \endplumedfile
118 :
119 : 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).
120 :
121 : \plumedfile
122 : volume: READ FILE=COLVAR VALUES=volume IGNORE_TIME
123 : mybias: READ FILE=COLVAR VALUES=mybias.bias IGNORE_TIME
124 :
125 : # Shift volume (to avoid numerical issues)
126 : rvol: COMBINE ARG=vol PARAMETERS=7.8 PERIODIC=NO
127 :
128 : # Weights
129 : bias_weights: REWEIGHT_BIAS TEMP=500 ARG=mybias.bias
130 : temp_press_weights: REWEIGHT_TEMP_PRESS TEMP=500 PRESSURE=0.06022140857 REWEIGHT_PRESSURE=180.66422571 VOLUME=volume
131 :
132 : # Ensemble average of the volume at 300 K and 300 MPa
133 : avg_vol: AVERAGE ARG=volume LOGWEIGHTS=bias_weights,temp_press_weights
134 : # Ensemble average of the density at 300 K and 300 MPa
135 : avg_density: CUSTOM ARG=avg_vol FUNC=1000/x PERIODIC=NO
136 :
137 : PRINT ARG=avg_density FILE=COLVAR_REWEIGHT STRIDE=1
138 : \endplumedfile
139 :
140 :
141 : 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).
142 :
143 : \plumedfile
144 : energy: READ FILE=COLVAR VALUES=energy IGNORE_TIME
145 : volume: READ FILE=COLVAR VALUES=volume IGNORE_TIME
146 : mybias: READ FILE=COLVAR VALUES=mybias.bias IGNORE_TIME
147 :
148 : # Shift energy and volume (to avoid numerical issues)
149 : rvol: COMBINE ARG=vol PARAMETERS=7.8 PERIODIC=NO
150 : renergy: COMBINE ARG=energy PARAMETERS=-13250 PERIODIC=NO
151 :
152 : # Weights
153 : bias_weights: REWEIGHT_BIAS TEMP=500 ARG=mybias.bias
154 : temp_press_weights: REWEIGHT_TEMP_PRESS TEMP=500 REWEIGHT_TEMP=300 PRESSURE=0.06022140857 REWEIGHT_PRESSURE=180.66422571 ENERGY=renergy VOLUME=rvolume
155 :
156 : # Ensemble average of the volume at 300 K and 300 MPa
157 : avg_vol: AVERAGE ARG=volume LOGWEIGHTS=bias_weights,temp_press_weights
158 : # Ensemble average of the density at 300 K and 300 MPa
159 : avg_density: CUSTOM ARG=avg_vol FUNC=1000/x PERIODIC=NO
160 :
161 : PRINT ARG=avg_density FILE=COLVAR_REWEIGHT STRIDE=1
162 : \endplumedfile
163 :
164 : */
165 : //+ENDPLUMEDOC
166 :
167 : namespace PLMD {
168 : namespace bias {
169 :
170 12 : class ReweightTemperaturePressure : public ReweightBase {
171 : private:
172 : ///
173 : double rpress_, press_, rtemp_;
174 : std::vector<Value*> myenergy, myvol;
175 : public:
176 : static void registerKeywords(Keywords&);
177 : explicit ReweightTemperaturePressure(const ActionOptions&ao);
178 : double getLogWeight() override;
179 : };
180 :
181 7840 : PLUMED_REGISTER_ACTION(ReweightTemperaturePressure,"REWEIGHT_TEMP_PRESS")
182 :
183 5 : void ReweightTemperaturePressure::registerKeywords(Keywords& keys ) {
184 5 : ReweightBase::registerKeywords( keys );
185 10 : keys.remove("ARG");
186 20 : keys.add("optional","ENERGY","Energy");
187 20 : keys.add("optional","VOLUME","Volume");
188 20 : keys.add("optional","REWEIGHT_PRESSURE","Reweighting pressure");
189 20 : keys.add("optional","PRESSURE","The system pressure");
190 20 : keys.add("optional","REWEIGHT_TEMP","Reweighting temperature");
191 5 : }
192 :
193 4 : ReweightTemperaturePressure::ReweightTemperaturePressure(const ActionOptions&ao):
194 : Action(ao),
195 4 : ReweightBase(ao)
196 : {
197 : // Initialize to not defined (negative)
198 4 : rpress_=-1;
199 4 : press_=-1;
200 4 : rtemp_=-1;
201 8 : parse("REWEIGHT_PRESSURE",rpress_);
202 8 : parse("PRESSURE",press_);
203 8 : parse("REWEIGHT_TEMP",rtemp_);
204 8 : rtemp_*=plumed.getAtoms().getKBoltzmann();
205 :
206 8 : parseArgumentList("ENERGY",myenergy);
207 4 : if(!myenergy.empty()) {
208 3 : log.printf(" with energies: ");
209 12 : for(unsigned i=0; i<myenergy.size(); i++) log.printf(" %s",myenergy[i]->getName().c_str());
210 3 : log.printf("\n");
211 : }
212 : //requestArguments(myenergy);
213 :
214 8 : parseArgumentList("VOLUME",myvol);
215 4 : if(!myvol.empty()) {
216 3 : log.printf(" with volumes: ");
217 12 : for(unsigned i=0; i<myvol.size(); i++) log.printf(" %s",myvol[i]->getName().c_str());
218 3 : log.printf("\n");
219 : }
220 :
221 : std::vector<Value*> conc;
222 4 : conc.insert(conc.begin(), myenergy.begin(), myenergy.end());
223 4 : conc.insert(conc.end(), myvol.begin(), myvol.end());
224 4 : requestArguments(conc);
225 :
226 : // 4 possible cases
227 : // Case 1) Reweight from T to T' with V=const (canonical)
228 6 : if (rtemp_>=0 && press_<0 && rpress_<0 && !myenergy.empty() && myvol.empty() ) {
229 3 : log.printf(" reweighting simulation from temperature %f to temperature %f at constant volume \n",simtemp/plumed.getAtoms().getKBoltzmann(),rtemp_/plumed.getAtoms().getKBoltzmann() );
230 1 : log.printf(" WARNING: If the simulation is performed at constant pressure add the keywords PRESSURE and VOLUME \n" );
231 : }
232 : // Case 2) Reweight from T to T' with P=const (isothermal-isobaric)
233 7 : 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_ );
234 : // Case 3) Reweight from P to P' with T=const (isothermal-isobaric)
235 5 : 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() );
236 : // Case 4) Reweight from T,P to T',P' (isothermal-isobaric)
237 5 : 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_);
238 0 : else error("Combination of ENERGY, VOLUME, REWEIGHT_PRESSURE, PRESSURE and REWEIGHT_TEMP not supported. Please refer to the manual for supported combinations.");
239 4 : }
240 :
241 1001 : double ReweightTemperaturePressure::getLogWeight() {
242 4004 : double energy=0.0; for(unsigned i=0; i<myenergy.size(); ++i) energy+=getArgument(i);
243 4004 : double volume=0.0; for(unsigned i=0; i<myvol.size(); ++i) volume+=getArgument(myenergy.size()+i);
244 : // 4 possible cases
245 : // Case 1) Reweight from T to T' with V=const (canonical)
246 1001 : if (rtemp_>=0 && press_<0 && rpress_<0) return ((1.0/simtemp)- (1.0/rtemp_) )*energy;
247 : // Case 2) Reweight from T to T' with P=const (isothermal-isobaric)
248 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;
249 : // Case 3) Reweight from P to P' with T=const (isothermal-isobaric)
250 1001 : else if (rtemp_<0 && press_>=0 && rpress_>=0) return (1.0/simtemp)*(press_ - rpress_)*volume;
251 : // Case 4) Reweight from T,P to T',P' (isothermal-isobaric)
252 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;
253 : else return 0;
254 : }
255 :
256 : }
257 5874 : }
|
