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 29 : void Ensemble::registerKeywords(Keywords& keys) {
75 29 : Function::registerKeywords(keys);
76 29 : keys.addFlag("REWEIGHT",false,"simple REWEIGHT using the latest ARG as energy");
77 29 : keys.addFlag("CENTRAL",false,"calculate a central moment instead of a standard moment");
78 29 : keys.add("optional","TEMP","the system temperature - this is only needed if you are reweighting");
79 29 : keys.add("optional","MOMENT","the moment you want to calculate in alternative to the mean or the variance");
80 29 : keys.add("optional","POWER","the power of the mean (and moment)");
81 29 : ActionWithValue::useCustomisableComponents(keys);
82 29 : }
83 :
84 27 : Ensemble::Ensemble(const ActionOptions&ao):
85 : Action(ao),
86 : Function(ao),
87 27 : do_reweight(false),
88 27 : do_moments(false),
89 27 : do_central(false),
90 27 : do_powers(false),
91 27 : kbt(-1.0),
92 27 : moment(0),
93 27 : power(0) {
94 27 : parseFlag("REWEIGHT", do_reweight);
95 27 : if(do_reweight) {
96 12 : kbt=getkBT();
97 12 : 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 15 : double temp=0.0;
102 30 : parse("TEMP",temp);
103 : }
104 :
105 27 : parse("MOMENT",moment);
106 27 : if(moment==1) {
107 0 : error("MOMENT can be any number but for 0 and 1");
108 : }
109 27 : if(moment!=0) {
110 0 : do_moments=true;
111 : }
112 27 : parseFlag("CENTRAL", do_central);
113 27 : if(!do_moments&&do_central) {
114 0 : error("To calculate a CENTRAL moment you need to define for which MOMENT");
115 : }
116 :
117 27 : parse("POWER",power);
118 27 : if(power==1) {
119 0 : error("POWER can be any number but for 0 and 1");
120 : }
121 27 : if(power!=0) {
122 2 : do_powers=true;
123 : }
124 :
125 27 : checkRead();
126 :
127 27 : master = (comm.Get_rank()==0);
128 27 : ens_dim=0;
129 27 : my_repl=0;
130 27 : if(master) {
131 17 : ens_dim=multi_sim_comm.Get_size();
132 17 : my_repl=multi_sim_comm.Get_rank();
133 : }
134 27 : comm.Bcast(ens_dim,0);
135 27 : comm.Bcast(my_repl,0);
136 27 : if(ens_dim<2) {
137 1 : log.printf("WARNING: ENSEMBLE with one replica is not doing any averaging!\n");
138 : }
139 :
140 : // prepare output components, the number depending on reweighing or not
141 27 : narg = getNumberOfArguments();
142 27 : if(do_reweight) {
143 12 : narg--;
144 : }
145 :
146 : // these are the averages
147 3044 : for(unsigned i=0; i<narg; i++) {
148 3017 : std::string s=getPntrToArgument(i)->getName();
149 3017 : addComponentWithDerivatives(s);
150 3017 : getPntrToComponent(i)->setNotPeriodic();
151 : }
152 : // these are the moments
153 27 : if(do_moments) {
154 0 : for(unsigned i=0; i<narg; i++) {
155 0 : std::string s=getPntrToArgument(i)->getName()+"_m";
156 0 : addComponentWithDerivatives(s);
157 0 : getPntrToComponent(i+narg)->setNotPeriodic();
158 : }
159 : }
160 :
161 27 : log.printf(" averaging over %u replicas.\n", ens_dim);
162 27 : if(do_reweight) {
163 12 : log.printf(" doing simple REWEIGHT using the latest ARGUMENT as energy.\n");
164 : }
165 27 : if(do_moments&&!do_central) {
166 0 : log.printf(" calculating also the %lf standard moment\n", moment);
167 : }
168 27 : if(do_moments&&do_central) {
169 0 : log.printf(" calculating also the %lf central moment\n", moment);
170 : }
171 27 : if(do_powers) {
172 2 : log.printf(" calculating the %lf power of the mean (and moment)\n", power);
173 : }
174 27 : }
175 :
176 0 : std::string Ensemble::getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const {
177 0 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
178 0 : if( cname==getPntrToArgument(i)->getName() ) {
179 0 : return "the average for argument " + cname;
180 : }
181 0 : if( cname==getPntrToArgument(i)->getName() + "_m" ) {
182 0 : return "the moment for argument " + cname;
183 : }
184 : }
185 0 : plumed_error();
186 : return "";
187 : }
188 :
189 :
190 125 : void Ensemble::calculate() {
191 : double norm = 0.0;
192 125 : double fact = 0.0;
193 :
194 : // calculate the weights either from BIAS
195 125 : if(do_reweight) {
196 : std::vector<double> bias;
197 0 : bias.resize(ens_dim);
198 0 : if(master) {
199 0 : bias[my_repl] = getArgument(narg);
200 0 : if(ens_dim>1) {
201 0 : multi_sim_comm.Sum(&bias[0], ens_dim);
202 : }
203 : }
204 0 : comm.Sum(&bias[0], ens_dim);
205 0 : const double maxbias = *(std::max_element(bias.begin(), bias.end()));
206 0 : for(unsigned i=0; i<ens_dim; ++i) {
207 0 : bias[i] = exp((bias[i]-maxbias)/kbt);
208 0 : norm += bias[i];
209 : }
210 0 : fact = bias[my_repl]/norm;
211 : // or arithmetic ones
212 : } else {
213 125 : norm = static_cast<double>(ens_dim);
214 125 : fact = 1.0/norm;
215 : }
216 :
217 125 : const double fact_kbt = fact/kbt;
218 :
219 125 : std::vector<double> mean(narg);
220 125 : std::vector<double> dmean(narg,fact);
221 : // calculate the mean
222 125 : if(master) {
223 2106 : for(unsigned i=0; i<narg; ++i) {
224 2007 : mean[i] = fact*getArgument(i);
225 : }
226 99 : if(ens_dim>1) {
227 98 : multi_sim_comm.Sum(&mean[0], narg);
228 : }
229 : }
230 125 : comm.Sum(&mean[0], narg);
231 :
232 : std::vector<double> v_moment, dv_moment;
233 : // calculate other moments
234 125 : if(do_moments) {
235 0 : v_moment.resize(narg);
236 0 : dv_moment.resize(narg);
237 : // standard moment
238 0 : if(!do_central) {
239 0 : if(master) {
240 0 : for(unsigned i=0; i<narg; ++i) {
241 0 : const double tmp = fact*std::pow(getArgument(i),moment-1);
242 0 : v_moment[i] = tmp*getArgument(i);
243 0 : dv_moment[i] = moment*tmp;
244 : }
245 0 : if(ens_dim>1) {
246 0 : multi_sim_comm.Sum(&v_moment[0], narg);
247 : }
248 : } else {
249 0 : for(unsigned i=0; i<narg; ++i) {
250 0 : const double tmp = fact*std::pow(getArgument(i),moment-1);
251 0 : dv_moment[i] = moment*tmp;
252 : }
253 : }
254 : // central moment
255 : } else {
256 0 : if(master) {
257 0 : for(unsigned i=0; i<narg; ++i) {
258 0 : const double tmp = std::pow(getArgument(i)-mean[i],moment-1);
259 0 : v_moment[i] = fact*tmp*(getArgument(i)-mean[i]);
260 0 : dv_moment[i] = moment*tmp*(fact-fact/norm);
261 : }
262 0 : if(ens_dim>1) {
263 0 : multi_sim_comm.Sum(&v_moment[0], narg);
264 : }
265 : } else {
266 0 : for(unsigned i=0; i<narg; ++i) {
267 0 : const double tmp = std::pow(getArgument(i)-mean[i],moment-1);
268 0 : dv_moment[i] = moment*tmp*(fact-fact/norm);
269 : }
270 : }
271 : }
272 0 : comm.Sum(&v_moment[0], narg);
273 : }
274 :
275 : // calculate powers of moments
276 125 : if(do_powers) {
277 72 : for(unsigned i=0; i<narg; ++i) {
278 48 : const double tmp1 = std::pow(mean[i],power-1);
279 48 : mean[i] *= tmp1;
280 48 : dmean[i] *= power*tmp1;
281 48 : if(do_moments) {
282 0 : const double tmp2 = std::pow(v_moment[i],power-1);
283 0 : v_moment[i] *= tmp2;
284 0 : dv_moment[i] *= power*tmp2;
285 : }
286 : }
287 : }
288 :
289 : // set components
290 3358 : for(unsigned i=0; i<narg; ++i) {
291 : // set mean
292 3233 : Value* v=getPntrToComponent(i);
293 3233 : v->set(mean[i]);
294 3233 : setDerivative(v, i, dmean[i]);
295 3233 : if(do_reweight) {
296 0 : const double w_tmp = fact_kbt*(getArgument(i) - mean[i]);
297 0 : setDerivative(v, narg, w_tmp);
298 : }
299 3233 : if(do_moments) {
300 : // set moments
301 0 : Value* u=getPntrToComponent(i+narg);
302 0 : u->set(v_moment[i]);
303 0 : setDerivative(u, i, dv_moment[i]);
304 0 : if(do_reweight) {
305 0 : const double w_tmp = fact_kbt*(pow(getArgument(i),moment) - v_moment[i]);
306 0 : setDerivative(u, narg, w_tmp);
307 : }
308 : }
309 : }
310 125 : }
311 :
312 : }
313 : }
|
