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 _label_.cvlabel.
35 :
36 : ## 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 :
41 : ```plumed
42 : dist: DISTANCE ATOMS=3,5
43 : ens: ENSEMBLE ARG=dist
44 : PRINT ARG=dist,ens.dist
45 : ```
46 :
47 : */
48 : //+ENDPLUMEDOC
49 :
50 :
51 : class Ensemble :
52 : public Function {
53 : unsigned ens_dim;
54 : unsigned my_repl;
55 : unsigned narg;
56 : bool master;
57 : bool do_reweight;
58 : bool do_moments;
59 : bool do_central;
60 : bool do_powers;
61 : double kbt;
62 : double moment;
63 : double power;
64 : public:
65 : explicit Ensemble(const ActionOptions&);
66 : std::string getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const override ;
67 : void calculate() override;
68 : static void registerKeywords(Keywords& keys);
69 : };
70 :
71 :
72 : PLUMED_REGISTER_ACTION(Ensemble,"ENSEMBLE")
73 :
74 407 : void Ensemble::registerKeywords(Keywords& keys) {
75 407 : Function::registerKeywords(keys);
76 407 : keys.addFlag("REWEIGHT",false,"simple REWEIGHT using the latest ARG as energy");
77 407 : keys.addFlag("CENTRAL",false,"calculate a central moment instead of a standard moment");
78 407 : keys.add("optional","TEMP","the system temperature - this is only needed if you are reweighting");
79 407 : keys.add("optional","MOMENT","the moment you want to calculate in alternative to the mean or the variance");
80 407 : keys.add("optional","POWER","the power of the mean (and moment)");
81 407 : ActionWithValue::useCustomisableComponents(keys);
82 407 : }
83 :
84 403 : Ensemble::Ensemble(const ActionOptions&ao):
85 : Action(ao),
86 : Function(ao),
87 403 : do_reweight(false),
88 403 : do_moments(false),
89 403 : do_central(false),
90 403 : do_powers(false),
91 403 : kbt(-1.0),
92 403 : moment(0),
93 403 : power(0) {
94 403 : parseFlag("REWEIGHT", do_reweight);
95 403 : if(do_reweight) {
96 76 : kbt=getkBT();
97 76 : 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 327 : double temp=0.0;
102 654 : parse("TEMP",temp);
103 : }
104 :
105 403 : parse("MOMENT",moment);
106 403 : if(moment==1) {
107 0 : error("MOMENT can be any number but for 0 and 1");
108 : }
109 403 : if(moment!=0) {
110 128 : do_moments=true;
111 : }
112 403 : parseFlag("CENTRAL", do_central);
113 403 : if(!do_moments&&do_central) {
114 0 : error("To calculate a CENTRAL moment you need to define for which MOMENT");
115 : }
116 :
117 403 : parse("POWER",power);
118 403 : if(power==1) {
119 0 : error("POWER can be any number but for 0 and 1");
120 : }
121 403 : if(power!=0) {
122 66 : do_powers=true;
123 : }
124 :
125 403 : checkRead();
126 :
127 403 : master = (comm.Get_rank()==0);
128 403 : ens_dim=getNumberOfReplicas();
129 403 : my_repl=getReplicaIndex();
130 403 : if(ens_dim<2) {
131 25 : log.printf("WARNING: ENSEMBLE with one replica is not doing any averaging!\n");
132 : }
133 :
134 : // prepare output components, the number depending on reweighing or not
135 403 : narg = getNumberOfArguments();
136 403 : if(do_reweight) {
137 76 : narg--;
138 : }
139 :
140 : // these are the averages
141 3796 : for(unsigned i=0; i<narg; i++) {
142 3393 : std::string s=getPntrToArgument(i)->getName();
143 3393 : addComponentWithDerivatives(s);
144 3393 : getPntrToComponent(i)->setNotPeriodic();
145 : }
146 : // these are the moments
147 403 : if(do_moments) {
148 256 : for(unsigned i=0; i<narg; i++) {
149 128 : std::string s=getPntrToArgument(i)->getName()+"_m";
150 128 : addComponentWithDerivatives(s);
151 128 : getPntrToComponent(i+narg)->setNotPeriodic();
152 : }
153 : }
154 :
155 403 : log.printf(" averaging over %u replicas.\n", ens_dim);
156 403 : if(do_reweight) {
157 76 : log.printf(" doing simple REWEIGHT using the latest ARGUMENT as energy.\n");
158 : }
159 403 : if(do_moments&&!do_central) {
160 64 : log.printf(" calculating also the %lf standard moment\n", moment);
161 : }
162 403 : if(do_moments&&do_central) {
163 64 : log.printf(" calculating also the %lf central moment\n", moment);
164 : }
165 403 : if(do_powers) {
166 66 : log.printf(" calculating the %lf power of the mean (and moment)\n", power);
167 : }
168 403 : }
169 :
170 0 : std::string Ensemble::getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const {
171 0 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
172 0 : if( cname==getPntrToArgument(i)->getName() ) {
173 0 : return "the average for argument " + cname;
174 : }
175 0 : if( cname==getPntrToArgument(i)->getName() + "_m" ) {
176 0 : return "the moment for argument " + cname;
177 : }
178 : }
179 0 : plumed_error();
180 : return "";
181 : }
182 :
183 :
184 2113 : void Ensemble::calculate() {
185 : double norm = 0.0;
186 2113 : double fact = 0.0;
187 :
188 : // calculate the weights either from BIAS
189 2113 : if(do_reweight) {
190 576 : std::vector<double> bias;
191 576 : bias.resize(ens_dim);
192 576 : if(master) {
193 288 : bias[my_repl] = getArgument(narg);
194 288 : if(ens_dim>1) {
195 288 : multi_sim_comm.Sum(&bias[0], ens_dim);
196 : }
197 : }
198 576 : comm.Sum(&bias[0], ens_dim);
199 576 : const double maxbias = *(std::max_element(bias.begin(), bias.end()));
200 1728 : for(unsigned i=0; i<ens_dim; ++i) {
201 1152 : bias[i] = exp((bias[i]-maxbias)/kbt);
202 1152 : norm += bias[i];
203 : }
204 576 : fact = bias[my_repl]/norm;
205 : // or arithmetic ones
206 576 : } else {
207 1537 : norm = static_cast<double>(ens_dim);
208 1537 : fact = 1.0/norm;
209 : }
210 :
211 2113 : const double fact_kbt = fact/kbt;
212 :
213 2113 : std::vector<double> mean(narg);
214 2113 : std::vector<double> dmean(narg,fact);
215 : // calculate the mean
216 2113 : if(master) {
217 4652 : for(unsigned i=0; i<narg; ++i) {
218 3280 : mean[i] = fact*getArgument(i);
219 : }
220 1372 : if(ens_dim>1) {
221 1344 : multi_sim_comm.Sum(&mean[0], narg);
222 : }
223 : }
224 2113 : comm.Sum(&mean[0], narg);
225 :
226 2113 : std::vector<double> v_moment, dv_moment;
227 2113 : std::vector<double> central_lower_moment(do_central ? narg : 0,0.0);
228 : // calculate other moments
229 2113 : if(do_moments) {
230 896 : v_moment.resize(narg);
231 896 : dv_moment.resize(narg);
232 : // standard moment
233 896 : if(!do_central) {
234 448 : if(master) {
235 448 : for(unsigned i=0; i<narg; ++i) {
236 224 : const double tmp = fact*std::pow(getArgument(i),moment-1);
237 224 : v_moment[i] = tmp*getArgument(i);
238 224 : dv_moment[i] = moment*tmp;
239 : }
240 224 : if(ens_dim>1) {
241 224 : multi_sim_comm.Sum(&v_moment[0], narg);
242 : }
243 : } else {
244 448 : for(unsigned i=0; i<narg; ++i) {
245 224 : const double tmp = fact*std::pow(getArgument(i),moment-1);
246 224 : dv_moment[i] = moment*tmp;
247 : }
248 : }
249 : // A central moment also depends on the replica-averaged mean.
250 : } else {
251 448 : if(master) {
252 448 : for(unsigned i=0; i<narg; ++i) {
253 224 : const double delta=getArgument(i)-mean[i];
254 224 : central_lower_moment[i]=fact*std::pow(delta,moment-1);
255 224 : v_moment[i]=central_lower_moment[i]*delta;
256 : }
257 224 : if(ens_dim>1) {
258 224 : multi_sim_comm.Sum(&v_moment[0], narg);
259 224 : multi_sim_comm.Sum(¢ral_lower_moment[0], narg);
260 : }
261 : }
262 448 : comm.Sum(¢ral_lower_moment[0], narg);
263 896 : for(unsigned i=0; i<narg; ++i) {
264 448 : dv_moment[i]=moment*fact*(std::pow(getArgument(i)-mean[i],moment-1)-central_lower_moment[i]);
265 : }
266 : }
267 896 : comm.Sum(&v_moment[0], narg);
268 : }
269 :
270 : // Differentiate normalized weights before applying an optional power.
271 2113 : std::vector<double> dmean_bias(do_reweight ? narg : 0,0.0);
272 2113 : std::vector<double> dmoment_bias(do_reweight && do_moments ? narg : 0,0.0);
273 2113 : if(do_reweight) {
274 1152 : for(unsigned i=0; i<narg; ++i) {
275 576 : const double delta=getArgument(i)-mean[i];
276 576 : dmean_bias[i]=fact_kbt*delta;
277 576 : if(do_moments) {
278 576 : if(do_central) {
279 288 : dmoment_bias[i]=fact_kbt*(std::pow(delta,moment)-v_moment[i]
280 288 : -moment*delta*central_lower_moment[i]);
281 : } else {
282 288 : dmoment_bias[i]=fact_kbt*(std::pow(getArgument(i),moment)-v_moment[i]);
283 : }
284 : }
285 : }
286 : }
287 :
288 : // calculate powers of moments
289 2113 : if(do_powers) {
290 968 : for(unsigned i=0; i<narg; ++i) {
291 496 : const double tmp1 = std::pow(mean[i],power-1);
292 496 : mean[i] *= tmp1;
293 496 : dmean[i] *= power*tmp1;
294 496 : if(do_reweight) {
295 288 : dmean_bias[i]*=power*tmp1;
296 : }
297 496 : if(do_moments) {
298 448 : const double tmp2 = std::pow(v_moment[i],power-1);
299 448 : v_moment[i] *= tmp2;
300 448 : dv_moment[i] *= power*tmp2;
301 448 : if(do_reweight) {
302 288 : dmoment_bias[i]*=power*tmp2;
303 : }
304 : }
305 : }
306 : }
307 :
308 : // set components
309 7334 : for(unsigned i=0; i<narg; ++i) {
310 : // set mean
311 5221 : Value* v=getPntrToComponent(i);
312 5221 : v->set(mean[i]);
313 5221 : setDerivative(v, i, dmean[i]);
314 5221 : if(do_reweight) {
315 576 : setDerivative(v, narg, dmean_bias[i]);
316 : }
317 5221 : if(do_moments) {
318 : // set moments
319 896 : Value* u=getPntrToComponent(i+narg);
320 896 : u->set(v_moment[i]);
321 896 : setDerivative(u, i, dv_moment[i]);
322 896 : if(do_reweight) {
323 576 : setDerivative(u, narg, dmoment_bias[i]);
324 : }
325 : }
326 : }
327 2113 : }
328 :
329 : }
330 : }
|
