Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2014-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 "Function.h"
23 : #include "tools/Communicator.h"
24 : #include "core/ActionRegister.h"
25 : #include "core/PlumedMain.h"
26 :
27 : namespace PLMD {
28 : namespace function {
29 :
30 : //+PLUMEDOC FUNCTION ENSEMBLE
31 : /*
32 : Calculates the replica averaging of a collective variable over multiple replicas.
33 :
34 : Each collective variable is averaged separately and stored in a component labelled <em>label</em>.cvlabel.
35 :
36 : \par Examples
37 :
38 : The following input tells plumed to calculate the distance between atoms 3 and 5
39 : and the average it over the available replicas.
40 : \plumedfile
41 : dist: DISTANCE ATOMS=3,5
42 : ens: ENSEMBLE ARG=dist
43 : PRINT ARG=dist,ens.dist
44 : \endplumedfile
45 :
46 : */
47 : //+ENDPLUMEDOC
48 :
49 :
50 : class Ensemble :
51 : public Function {
52 : unsigned ens_dim;
53 : unsigned my_repl;
54 : unsigned narg;
55 : bool master;
56 : bool do_reweight;
57 : bool do_moments;
58 : bool do_central;
59 : bool do_powers;
60 : double kbt;
61 : double moment;
62 : double power;
63 : public:
64 : explicit Ensemble(const ActionOptions&);
65 : std::string getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const override ;
66 : void calculate() override;
67 : static void registerKeywords(Keywords& keys);
68 : };
69 :
70 :
71 : PLUMED_REGISTER_ACTION(Ensemble,"ENSEMBLE")
72 :
73 164 : void Ensemble::registerKeywords(Keywords& keys) {
74 164 : Function::registerKeywords(keys);
75 164 : keys.use("ARG");
76 328 : keys.addFlag("REWEIGHT",false,"simple REWEIGHT using the latest ARG as energy");
77 328 : keys.addFlag("CENTRAL",false,"calculate a central moment instead of a standard moment");
78 328 : keys.add("optional","TEMP","the system temperature - this is only needed if you are reweighting");
79 328 : keys.add("optional","MOMENT","the moment you want to calculate in alternative to the mean or the variance");
80 328 : keys.add("optional","POWER","the power of the mean (and moment)");
81 164 : ActionWithValue::useCustomisableComponents(keys);
82 164 : }
83 :
84 160 : Ensemble::Ensemble(const ActionOptions&ao):
85 : Action(ao),
86 : Function(ao),
87 160 : do_reweight(false),
88 160 : do_moments(false),
89 160 : do_central(false),
90 160 : do_powers(false),
91 160 : kbt(-1.0),
92 160 : moment(0),
93 160 : power(0) {
94 160 : parseFlag("REWEIGHT", do_reweight);
95 160 : if(do_reweight) {
96 77 : kbt=getkBT();
97 77 : if(kbt==0.0) {
98 0 : error("Unless the MD engine passes the temperature to plumed, with REWEIGHT you must specify TEMP");
99 : }
100 : } else {
101 83 : double temp=0.0;
102 166 : parse("TEMP",temp);
103 : }
104 :
105 160 : parse("MOMENT",moment);
106 160 : if(moment==1) {
107 0 : error("MOMENT can be any number but for 0 and 1");
108 : }
109 160 : if(moment!=0) {
110 132 : do_moments=true;
111 : }
112 160 : parseFlag("CENTRAL", do_central);
113 160 : if(!do_moments&&do_central) {
114 0 : error("To calculate a CENTRAL moment you need to define for which MOMENT");
115 : }
116 :
117 160 : parse("POWER",power);
118 160 : if(power==1) {
119 0 : error("POWER can be any number but for 0 and 1");
120 : }
121 160 : if(power!=0) {
122 67 : do_powers=true;
123 : }
124 :
125 160 : master = (comm.Get_rank()==0);
126 160 : ens_dim=0;
127 160 : my_repl=0;
128 160 : if(master) {
129 86 : ens_dim=multi_sim_comm.Get_size();
130 86 : my_repl=multi_sim_comm.Get_rank();
131 : }
132 160 : comm.Bcast(ens_dim,0);
133 160 : comm.Bcast(my_repl,0);
134 160 : if(ens_dim<2) {
135 6 : log.printf("WARNING: ENSEMBLE with one replica is not doing any averaging!\n");
136 : }
137 :
138 : // prepare output components, the number depending on reweighing or not
139 160 : narg = getNumberOfArguments();
140 160 : if(do_reweight) {
141 77 : narg--;
142 : }
143 :
144 : // these are the averages
145 3310 : for(unsigned i=0; i<narg; i++) {
146 3150 : std::string s=getPntrToArgument(i)->getName();
147 3150 : addComponentWithDerivatives(s);
148 3150 : getPntrToComponent(i)->setNotPeriodic();
149 : }
150 : // these are the moments
151 160 : if(do_moments) {
152 264 : for(unsigned i=0; i<narg; i++) {
153 132 : std::string s=getPntrToArgument(i)->getName()+"_m";
154 132 : addComponentWithDerivatives(s);
155 132 : getPntrToComponent(i+narg)->setNotPeriodic();
156 : }
157 : }
158 :
159 : // checkRead() may calculate constant arguments immediately.
160 : // All state and output components must be ready first.
161 160 : checkRead();
162 :
163 160 : log.printf(" averaging over %u replicas.\n", ens_dim);
164 160 : if(do_reweight) {
165 77 : log.printf(" doing simple REWEIGHT using the latest ARGUMENT as energy.\n");
166 : }
167 160 : if(do_moments&&!do_central) {
168 67 : log.printf(" calculating also the %lf standard moment\n", moment);
169 : }
170 160 : if(do_moments&&do_central) {
171 65 : log.printf(" calculating also the %lf central moment\n", moment);
172 : }
173 160 : if(do_powers) {
174 67 : log.printf(" calculating the %lf power of the mean (and moment)\n", power);
175 : }
176 160 : }
177 :
178 0 : std::string Ensemble::getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const {
179 0 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
180 0 : if( cname==getPntrToArgument(i)->getName() ) {
181 0 : return "the average for argument " + cname;
182 : }
183 0 : if( cname==getPntrToArgument(i)->getName() + "_m" ) {
184 0 : return "the moment for argument " + cname;
185 : }
186 : }
187 0 : plumed_error();
188 : return "";
189 : }
190 :
191 :
192 1031 : void Ensemble::calculate() {
193 : double norm = 0.0;
194 1031 : double fact = 0.0;
195 :
196 : // calculate the weights either from BIAS
197 1031 : if(do_reweight) {
198 : std::vector<double> bias;
199 578 : bias.resize(ens_dim);
200 578 : if(master) {
201 290 : bias[my_repl] = getArgument(narg);
202 290 : if(ens_dim>1) {
203 288 : multi_sim_comm.Sum(&bias[0], ens_dim);
204 : }
205 : }
206 578 : comm.Sum(&bias[0], ens_dim);
207 578 : const double maxbias = *(std::max_element(bias.begin(), bias.end()));
208 1732 : for(unsigned i=0; i<ens_dim; ++i) {
209 1154 : bias[i] = exp((bias[i]-maxbias)/kbt);
210 1154 : norm += bias[i];
211 : }
212 578 : fact = bias[my_repl]/norm;
213 : // or arithmetic ones
214 : } else {
215 453 : norm = static_cast<double>(ens_dim);
216 453 : fact = 1.0/norm;
217 : }
218 :
219 1031 : const double fact_kbt = fact/kbt;
220 :
221 1031 : std::vector<double> mean(narg);
222 1031 : std::vector<double> dmean(narg,fact);
223 : // calculate the mean
224 1031 : if(master) {
225 3022 : for(unsigned i=0; i<narg; ++i) {
226 2465 : mean[i] = fact*getArgument(i);
227 : }
228 557 : if(ens_dim>1) {
229 546 : multi_sim_comm.Sum(&mean[0], narg);
230 : }
231 : }
232 1031 : comm.Sum(&mean[0], narg);
233 :
234 : std::vector<double> v_moment, dv_moment;
235 1031 : std::vector<double> central_lower_moment(do_central ? narg : 0,0.0);
236 : // calculate other moments
237 1031 : if(do_moments) {
238 904 : v_moment.resize(narg);
239 904 : dv_moment.resize(narg);
240 : // standard moment
241 904 : if(!do_central) {
242 454 : if(master) {
243 460 : for(unsigned i=0; i<narg; ++i) {
244 230 : const double tmp = fact*std::pow(getArgument(i),moment-1);
245 230 : v_moment[i] = tmp*getArgument(i);
246 230 : dv_moment[i] = moment*tmp;
247 : }
248 230 : if(ens_dim>1) {
249 224 : multi_sim_comm.Sum(&v_moment[0], narg);
250 : }
251 : } else {
252 448 : for(unsigned i=0; i<narg; ++i) {
253 224 : const double tmp = fact*std::pow(getArgument(i),moment-1);
254 224 : dv_moment[i] = moment*tmp;
255 : }
256 : }
257 : // A central moment also depends on the replica-averaged mean.
258 : } else {
259 450 : if(master) {
260 452 : for(unsigned i=0; i<narg; ++i) {
261 226 : const double delta=getArgument(i)-mean[i];
262 226 : central_lower_moment[i]=fact*std::pow(delta,moment-1);
263 226 : v_moment[i]=central_lower_moment[i]*delta;
264 : }
265 226 : if(ens_dim>1) {
266 224 : multi_sim_comm.Sum(&v_moment[0], narg);
267 224 : multi_sim_comm.Sum(¢ral_lower_moment[0], narg);
268 : }
269 : }
270 450 : comm.Sum(¢ral_lower_moment[0], narg);
271 900 : for(unsigned i=0; i<narg; ++i) {
272 450 : dv_moment[i]=moment*fact*(std::pow(getArgument(i)-mean[i],moment-1)-central_lower_moment[i]);
273 : }
274 : }
275 904 : comm.Sum(&v_moment[0], narg);
276 : }
277 :
278 : // Differentiate normalized weights before applying an optional power.
279 1031 : std::vector<double> dmean_bias(do_reweight ? narg : 0,0.0);
280 1031 : std::vector<double> dmoment_bias(do_reweight && do_moments ? narg : 0,0.0);
281 1031 : if(do_reweight) {
282 1156 : for(unsigned i=0; i<narg; ++i) {
283 578 : const double delta=getArgument(i)-mean[i];
284 578 : dmean_bias[i]=fact_kbt*delta;
285 578 : if(do_moments) {
286 578 : if(do_central) {
287 288 : dmoment_bias[i]=fact_kbt*(std::pow(delta,moment)-v_moment[i]
288 288 : -moment*delta*central_lower_moment[i]);
289 : } else {
290 290 : dmoment_bias[i]=fact_kbt*(std::pow(getArgument(i),moment)-v_moment[i]);
291 : }
292 : }
293 : }
294 : }
295 :
296 : // calculate powers of moments
297 1031 : if(do_powers) {
298 972 : for(unsigned i=0; i<narg; ++i) {
299 498 : const double tmp1 = std::pow(mean[i],power-1);
300 498 : mean[i] *= tmp1;
301 498 : dmean[i] *= power*tmp1;
302 498 : if(do_reweight) {
303 288 : dmean_bias[i]*=power*tmp1;
304 : }
305 498 : if(do_moments) {
306 450 : const double tmp2 = std::pow(v_moment[i],power-1);
307 450 : v_moment[i] *= tmp2;
308 450 : dv_moment[i] *= power*tmp2;
309 450 : if(do_reweight) {
310 288 : dmoment_bias[i]*=power*tmp2;
311 : }
312 : }
313 : }
314 : }
315 :
316 : // set components
317 5170 : for(unsigned i=0; i<narg; ++i) {
318 : // set mean
319 4139 : Value* v=getPntrToComponent(i);
320 4139 : v->set(mean[i]);
321 4139 : setDerivative(v, i, dmean[i]);
322 4139 : if(do_reweight) {
323 578 : setDerivative(v, narg, dmean_bias[i]);
324 : }
325 4139 : if(do_moments) {
326 : // set moments
327 904 : Value* u=getPntrToComponent(i+narg);
328 904 : u->set(v_moment[i]);
329 904 : setDerivative(u, i, dv_moment[i]);
330 904 : if(do_reweight) {
331 578 : setDerivative(u, narg, dmoment_bias[i]);
332 : }
333 : }
334 : }
335 1031 : }
336 :
337 : }
338 : }
|
