Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2014-2018 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 "ActionRegister.h"
24 : #include "core/PlumedMain.h"
25 : #include "core/Atoms.h"
26 :
27 : using namespace std;
28 :
29 : namespace PLMD {
30 : namespace function {
31 :
32 : //+PLUMEDOC FUNCTION ENSEMBLE
33 : /*
34 : Calculates the replica averaging of a collective variable over multiple replicas.
35 :
36 : Each collective variable is averaged separately and stored in a component labelled <em>label</em>.cvlabel.
37 :
38 : \par Examples
39 :
40 : The following input tells plumed to calculate the distance between atoms 3 and 5
41 : and the average it over the available replicas.
42 : \verbatim
43 : dist: DISTANCE ATOMS=3,5
44 : ens: ENSEMBLE ARG=dist
45 : PRINT ARG=dist,ens.dist
46 : \endverbatim
47 : (See also \ref PRINT and \ref DISTANCE).
48 :
49 : */
50 : //+ENDPLUMEDOC
51 :
52 :
53 46 : class Ensemble :
54 : public Function
55 : {
56 : unsigned ens_dim;
57 : unsigned my_repl;
58 : unsigned narg;
59 : bool master;
60 : bool do_reweight;
61 : bool do_moments;
62 : bool do_central;
63 : bool do_powers;
64 : double kbt;
65 : double moment;
66 : double power;
67 : public:
68 : explicit Ensemble(const ActionOptions&);
69 : void calculate();
70 : static void registerKeywords(Keywords& keys);
71 : };
72 :
73 :
74 2546 : PLUMED_REGISTER_ACTION(Ensemble,"ENSEMBLE")
75 :
76 24 : void Ensemble::registerKeywords(Keywords& keys) {
77 24 : Function::registerKeywords(keys);
78 24 : keys.use("ARG");
79 24 : keys.addFlag("REWEIGHT",false,"simple REWEIGHT using the latest ARG as energy");
80 24 : keys.addFlag("CENTRAL",false,"calculate a central moment instead of a standard moment");
81 24 : keys.add("optional","TEMP","the system temperature - this is only needed if you are reweighting");
82 24 : keys.add("optional","MOMENT","the moment you want to calculate in alternative to the mean or the variance");
83 24 : keys.add("optional","POWER","the power of the mean (and moment)");
84 24 : ActionWithValue::useCustomisableComponents(keys);
85 24 : }
86 :
87 23 : Ensemble::Ensemble(const ActionOptions&ao):
88 : Action(ao),
89 : Function(ao),
90 : do_reweight(false),
91 : do_moments(false),
92 : do_central(false),
93 : do_powers(false),
94 : kbt(-1.0),
95 : moment(0),
96 23 : power(0)
97 : {
98 23 : parseFlag("REWEIGHT", do_reweight);
99 23 : double temp=0.0;
100 23 : parse("TEMP",temp);
101 23 : if(do_reweight) {
102 4 : if(temp>0.0) kbt=plumed.getAtoms().getKBoltzmann()*temp;
103 0 : else kbt=plumed.getAtoms().getKbT();
104 4 : if(kbt==0.0) error("Unless the MD engine passes the temperature to plumed, with REWEIGHT you must specify TEMP");
105 : }
106 :
107 23 : parse("MOMENT",moment);
108 23 : if(moment==1) error("MOMENT can be any number but for 0 and 1");
109 23 : if(moment!=0) do_moments=true;
110 23 : parseFlag("CENTRAL", do_central);
111 23 : if(!do_moments&&do_central) error("To calculate a CENTRAL moment you need to define for which MOMENT");
112 :
113 23 : parse("POWER",power);
114 23 : if(power==1) error("POWER can be any number but for 0 and 1");
115 23 : if(power!=0) do_powers=true;
116 :
117 23 : checkRead();
118 :
119 23 : master = (comm.Get_rank()==0);
120 23 : ens_dim=0;
121 23 : my_repl=0;
122 23 : if(master) {
123 15 : ens_dim=multi_sim_comm.Get_size();
124 15 : my_repl=multi_sim_comm.Get_rank();
125 : }
126 23 : comm.Bcast(ens_dim,0);
127 23 : comm.Bcast(my_repl,0);
128 23 : if(ens_dim<2) log.printf("WARNING: ENSEMBLE with one replica is not doing any averaging!\n");
129 :
130 : // prepare output components, the number depending on reweighing or not
131 23 : narg = getNumberOfArguments();
132 23 : if(do_reweight) narg--;
133 :
134 : // these are the averages
135 3024 : for(unsigned i=0; i<narg; i++) {
136 3001 : std::string s=getPntrToArgument(i)->getName();
137 3001 : addComponentWithDerivatives(s);
138 3001 : getPntrToComponent(i)->setNotPeriodic();
139 3001 : }
140 : // these are the moments
141 23 : if(do_moments) {
142 0 : for(unsigned i=0; i<narg; i++) {
143 0 : std::string s=getPntrToArgument(i)->getName()+"_m";
144 0 : addComponentWithDerivatives(s);
145 0 : getPntrToComponent(i+narg)->setNotPeriodic();
146 0 : }
147 : }
148 :
149 23 : log.printf(" averaging over %u replicas.\n", ens_dim);
150 23 : if(do_reweight) log.printf(" doing simple REWEIGHT using the latest ARGUMENT as energy.\n");
151 23 : if(do_moments&&!do_central) log.printf(" calculating also the %lf standard moment\n", moment);
152 23 : if(do_moments&&do_central) log.printf(" calculating also the %lf central moment\n", moment);
153 23 : if(do_powers) log.printf(" calculating the %lf power of the mean (and moment)\n", power);
154 23 : }
155 :
156 221 : void Ensemble::calculate() {
157 221 : double norm = 0.0;
158 221 : double fact = 0.0;
159 :
160 : // calculate the weights either from BIAS
161 221 : if(do_reweight) {
162 48 : vector<double> bias;
163 48 : bias.resize(ens_dim);
164 48 : if(master) {
165 24 : bias[my_repl] = getArgument(narg);
166 24 : if(ens_dim>1) multi_sim_comm.Sum(&bias[0], ens_dim);
167 : }
168 48 : comm.Sum(&bias[0], ens_dim);
169 48 : const double maxbias = *(std::max_element(bias.begin(), bias.end()));
170 144 : for(unsigned i=0; i<ens_dim; ++i) {
171 96 : bias[i] = exp((bias[i]-maxbias)/kbt);
172 96 : norm += bias[i];
173 : }
174 48 : fact = bias[my_repl]/norm;
175 : // or arithmetic ones
176 : } else {
177 173 : norm = static_cast<double>(ens_dim);
178 173 : fact = 1.0/norm;
179 : }
180 :
181 221 : const double fact_kbt = fact/kbt;
182 :
183 221 : vector<double> mean(narg);
184 442 : vector<double> dmean(narg,fact);
185 : // calculate the mean
186 221 : if(master) {
187 147 : for(unsigned i=0; i<narg; ++i) mean[i] = fact*getArgument(i);
188 147 : if(ens_dim>1) multi_sim_comm.Sum(&mean[0], narg);
189 : }
190 221 : comm.Sum(&mean[0], narg);
191 :
192 442 : vector<double> v_moment, dv_moment;
193 : // calculate other moments
194 221 : if(do_moments) {
195 0 : v_moment.resize(narg);
196 0 : dv_moment.resize(narg);
197 : // standard moment
198 0 : if(!do_central) {
199 0 : if(master) {
200 0 : for(unsigned i=0; i<narg; ++i) {
201 0 : const double tmp = fact*pow(getArgument(i),moment-1);
202 0 : v_moment[i] = tmp*getArgument(i);
203 0 : dv_moment[i] = moment*tmp;
204 : }
205 0 : if(ens_dim>1) multi_sim_comm.Sum(&v_moment[0], narg);
206 : } else {
207 0 : for(unsigned i=0; i<narg; ++i) {
208 0 : const double tmp = fact*pow(getArgument(i),moment-1);
209 0 : dv_moment[i] = moment*tmp;
210 : }
211 : }
212 : // central moment
213 : } else {
214 0 : if(master) {
215 0 : for(unsigned i=0; i<narg; ++i) {
216 0 : const double tmp = pow(getArgument(i)-mean[i],moment-1);
217 0 : v_moment[i] = fact*tmp*(getArgument(i)-mean[i]);
218 0 : dv_moment[i] = moment*tmp*(fact-fact/norm);
219 : }
220 0 : if(ens_dim>1) multi_sim_comm.Sum(&v_moment[0], narg);
221 : } else {
222 0 : for(unsigned i=0; i<narg; ++i) {
223 0 : const double tmp = pow(getArgument(i)-mean[i],moment-1);
224 0 : dv_moment[i] = moment*tmp*(fact-fact/norm);
225 : }
226 : }
227 : }
228 0 : comm.Sum(&v_moment[0], narg);
229 : }
230 :
231 : // calculate powers of moments
232 221 : if(do_powers) {
233 72 : for(unsigned i=0; i<narg; ++i) {
234 48 : const double tmp1 = pow(mean[i],power-1);
235 48 : mean[i] *= tmp1;
236 48 : dmean[i] *= power*tmp1;
237 48 : if(do_moments) {
238 0 : const double tmp2 = pow(v_moment[i],power-1);
239 0 : v_moment[i] *= tmp2;
240 0 : dv_moment[i] *= power*tmp2;
241 : }
242 : }
243 : }
244 :
245 : // set components
246 3838 : for(unsigned i=0; i<narg; ++i) {
247 : // set mean
248 3617 : Value* v=getPntrToComponent(i);
249 3617 : v->set(mean[i]);
250 3617 : setDerivative(v, i, dmean[i]);
251 3617 : if(do_reweight) {
252 192 : const double w_tmp = fact_kbt*(getArgument(i) - mean[i]);
253 192 : setDerivative(v, narg, w_tmp);
254 : }
255 3617 : if(do_moments) {
256 : // set moments
257 0 : Value* u=getPntrToComponent(i+narg);
258 0 : u->set(v_moment[i]);
259 0 : setDerivative(u, i, dv_moment[i]);
260 0 : if(do_reweight) {
261 0 : const double w_tmp = fact_kbt*(pow(getArgument(i),moment) - v_moment[i]);
262 0 : setDerivative(u, narg, w_tmp);
263 : }
264 : }
265 221 : }
266 221 : }
267 :
268 : }
269 2523 : }
|
