1 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 Copyright (c) 2016-2021 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 "Bias.h"
23 #include "core/PlumedMain.h"
24 #include "core/Atoms.h"
25 #include <string>
26 #include <cstring>
27 //#include "ActionRegister.h"
28 #include "core/ActionRegister.h"
29 #include "core/ActionWithValue.h"
30 #include "tools/Communicator.h"
31 #include "tools/File.h"
32 #include <iostream>
33
34 //#include "Analysis.h"
35 //#include "core/PlumedMain.h"
36 //#include "core/ActionRegister.h"
37 //#include "tools/Grid.h"
38 //#include "tools/KernelFunctions.h"
39 //#include "tools/IFile.h"
40 //#include "tools/OFile.h"
41
42 // The original implementation of this method was contributed
43 // by Andrea Cesari (andreacesari90@gmail.com).
44 // Copyright has been then transferred to PLUMED developers
45 // (see https://github.com/plumed/plumed2/blob/master/.github/CONTRIBUTING.md)
46
47 using namespace std;
48
49
50 namespace PLMD {
51 namespace bias {
52
53 //+PLUMEDOC BIAS MAXENT
54 /*
55 Add a linear biasing potential on one or more variables \f$f_{i}\left(\boldsymbol{x}\right)\f$ satisfying the maximum entropy principle as proposed in Ref. \cite cesari2016maxent .
56
57 \warning
58 Notice that syntax is still under revision and might change
59
60 The resulting biasing potential is given by:
61 \f[
62 V_{BIAS}(\boldsymbol{x},t)=K_{B}T\sum_{i=1}^{\#arguments}f_{i}(\boldsymbol{x},t)\lambda_{i}(t)
63 \f]
64 Lagrangian multipliers \f$ \lambda_{i}\f$ are updated, every PACE steps, according to the following update rule:
65 \f[
66 \lambda_{i}=\lambda_{i}+\frac{k_{i}}{1+\frac{t}{\tau_{i}}}\left(f_{exp,i}+\xi_{i}\lambda_{i}-f_{i}(\boldsymbol{x})\right)
67 \f]
68 \f$k\f$ set the initial value of the learning rate and its units are \f$[observable]^{-2}ps^{-1}\f$. This can be set with the keyword KAPPA.
69 The number of components for any KAPPA vector must be equal to the number of arguments of the action.
70
71 Variable \f$ \xi_{i}(\lambda) \f$ is related to the chosen prior to model experimental errors. If a GAUSSIAN prior is used then:
72 \f[
73 \xi_{i}(\lambda)=-\lambda_{i}\sigma^{2}
74 \f]
75 where \f$ \sigma \f$ is the typical expected error on the observable \f$ f_i\f$.
76 For a LAPLACE prior:
77 \f[
78 \xi_{i}(\lambda)=-\frac{\lambda_{i}\sigma^{2}}{1-\frac{\lambda^{2}\sigma^{2}}{2}}
79
80 \f]
81 The value of \f$ \xi(\lambda,t)\f$ is written in output as a component named: argument name followed by the string _error.
82 Setting \f$ \sigma =0\f$ is equivalent to enforce a pure Maximum Entropy restraint without any noise modelling.
83 This method can be also used to enforce inequality restraint as shown in following examples.
84
85 Notice that a similar method is available as \ref EDS, although with different features and using a different optimization algorithm.
86
87 \par Examples
88
89 The following input tells plumed to restrain the distance between atoms 7 and 15
90 and the distance between atoms 2 and 19, at different equilibrium
91 values, and to print the energy of the restraint.
92 Lagrangian multiplier will be printed on a file called restraint.LAGMULT with a stride set by the variable PACE to 200ps.
93 Moreover plumed will compute the average of each Lagrangian multiplier in the window [TSTART,TEND] and use that to continue the simulations with fixed Lagrangian multipliers.
94 \plumedfile
95 DISTANCE ATOMS=7,15 LABEL=d1
96 DISTANCE ATOMS=2,19 LABEL=d2
97 MAXENT ...
98 ARG=d1,d2
99 TYPE=EQUAL
100 AT=0.2,0.5
101 KAPPA=35000.0,35000.0
102 TAU=0.02,0.02
103 PACE=200
104 TSTART=100
105 TEND=500
106 LABEL=restraint
107 ... MAXENT
108 PRINT ARG=restraint.bias
109 \endplumedfile
110 Lagrangian multipliers will be printed on a file called restraint.bias
111 The following input tells plumed to restrain the distance between atoms 7 and 15
112 to be greater than 0.2 and to print the energy of the restraint
113 \plumedfile
114 DISTANCE ATOMS=7,15 LABEL=d
115 MAXENT ARG=d TYPE=INEQUAL> AT=0.02 KAPPA=35000.0 TAU=3 LABEL=restraint
116 PRINT ARG=restraint.bias
117 \endplumedfile
118
119 (See also \ref DISTANCE and \ref PRINT).
120
121 */
122 //+ENDPLUMEDOC
123
124 class MaxEnt : public Bias {
125 std::vector<double> at;
126 std::vector<double> kappa;
127 std::vector<double> lambda;
128 std::vector<double> avgx;
129 std::vector<double> work;
130 std::vector<double> oldlambda;
131 std::vector<double> tau;
132 std::vector<double> avglambda;
133 std::vector<double> avglambda_restart;
134 std::vector<double> expected_eps;
135 std::vector<double> apply_weights;
136 double sigma;
137 double tstart;
138 double tend;
139 double avgstep; //current number of samples over which to compute the average. Check if could be replaced bu getStep()
140 long int pace_;
141 long int stride_;
142 double totalWork;
143 double BetaReweightBias;
144 double simtemp;
145 vector<ActionWithValue*> biases;
146 std::string type;
147 std::string error_type;
148 double alpha;
149 double avg_counter;
150 int learn_replica;
151 Value* valueForce2;
152 Value* valueWork;
153 OFile lagmultOfile_;
154 IFile ifile;
155 string lagmultfname;
156 string ifilesnames;
157 string fmt;
158 bool isFirstStep;
159 bool reweight;
160 bool no_broadcast;
161 bool printFirstStep;
162 std::vector<bool> done_average;
163 int myrep,nrep;
164 public:
165 explicit MaxEnt(const ActionOptions&);
166 void calculate() override;
167 void update() override;
168 void update_lambda();
169 static void registerKeywords(Keywords& keys);
170 void ReadLagrangians(IFile &ifile);
171 void WriteLagrangians(vector<double> &lagmult,OFile &file);
172 double compute_error(string &err_type,double &l);
173 double convert_lambda(string &type,double lold);
174 void check_lambda_boundaries(string &err_type,double &l);
175 };
176 PLUMED_REGISTER_ACTION(MaxEnt,"MAXENT")
177
registerKeywords(Keywords & keys)178 void MaxEnt::registerKeywords(Keywords& keys) {
179 Bias::registerKeywords(keys);
180 componentsAreNotOptional(keys);
181 keys.use("ARG");
182 keys.add("compulsory","KAPPA","0.0","specifies the initial value for the learning rate");
183 keys.add("compulsory","TAU","Specify the dumping time for the learning rate.");
184 keys.add("compulsory","TYPE","specify the restraint type. "
185 "EQUAL to restrain the variable at a given equilibrium value "
186 "INEQUAL< to restrain the variable to be smaller than a given value "
187 "INEQUAL> to restrain the variable to be greater than a given value");
188 keys.add("optional","ERROR_TYPE","specify the prior on the error to use."
189 "GAUSSIAN: use a Gaussian prior "
190 "LAPLACE: use a Laplace prior");
191 keys.add("optional","TSTART","time from where to start averaging the Lagrangian multiplier. By default no average is computed, hence lambda is updated every PACE steps");
192 keys.add("optional","TEND","time in ps where to stop to compute the average of Lagrangian multiplier. From this time until the end of the simulation Lagrangian multipliers are kept fix to the average computed between TSTART and TEND;");
193 keys.add("optional","ALPHA","default=1.0; To be used with LAPLACE KEYWORD, allows to choose a prior function proportional to a Gaussian times an exponential function. ALPHA=1 correspond to the LAPLACE prior.");
194 keys.add("compulsory","AT","the position of the restraint");
195 keys.add("optional","SIGMA","The typical errors expected on observable");
196 keys.add("optional","FILE","Lagrangian multipliers output file. The default name is: label name followed by the string .LAGMULT ");
197 keys.add("optional","LEARN_REPLICA","In a multiple replica environment specify which is the reference replica. By default replica 0 will be used.");
198 keys.add("optional","APPLY_WEIGHTS","Vector of weights containing 1 in correspondence of each replica that will receive the Lagrangian multiplier from the current one.");
199 keys.add("optional","PACE","the frequency for Lagrangian multipliers update");
200 keys.add("optional","PRINT_STRIDE","stride of Lagrangian multipliers output file. If no STRIDE is passed they are written every time they are updated (PACE).");
201 keys.add("optional","FMT","specify format for Lagrangian multipliers files (useful to decrease the number of digits in regtests)");
202 keys.addFlag("REWEIGHT",false,"to be used with plumed driver in order to reweight a trajectory a posteriori");
203 keys.addFlag("NO_BROADCAST",false,"If active will avoid Lagrangian multipliers to be communicated to other replicas.");
204 keys.add("optional","TEMP","the system temperature. This is required if you are reweighting.");
205 keys.addOutputComponent("force2","default","the instantaneous value of the squared force due to this bias potential");
206 keys.addOutputComponent("work","default","the instantaneous value of the work done by the biasing force");
207 keys.addOutputComponent("_work","default","the instantaneous value of the work done by the biasing force for each argument. "
208 "These quantities will named with the arguments of the bias followed by "
209 "the character string _work.");
210 keys.addOutputComponent("_error","default","Instantaneous values of the discrepancy between the observable and the restraint center");
211 keys.addOutputComponent("_coupling","default","Instantaneous values of Lagrangian multipliers. They are also written by default in a separate output file.");
212 keys.use("RESTART");
213 }
MaxEnt(const ActionOptions & ao)214 MaxEnt::MaxEnt(const ActionOptions&ao):
215 PLUMED_BIAS_INIT(ao),
216 at(getNumberOfArguments()),
217 kappa(getNumberOfArguments(),0.0),
218 lambda(getNumberOfArguments(),0.0),
219 avgx(getNumberOfArguments(),0.0),
220 oldlambda(getNumberOfArguments(),0.0),
221 tau(getNumberOfArguments(),getTimeStep()),
222 avglambda(getNumberOfArguments(),0.0),
223 avglambda_restart(getNumberOfArguments(),0.0),
224 expected_eps(getNumberOfArguments(),0.0),
225 sigma(0.0),
226 pace_(100),
227 stride_(100),
228 alpha(1.0),
229 avg_counter(0.0),
230 isFirstStep(true),
231 reweight(false),
232 no_broadcast(false),
233 printFirstStep(true),
234 done_average(getNumberOfArguments(),false)
235 {
236 if(comm.Get_rank()==0) nrep=multi_sim_comm.Get_size();
237 if(comm.Get_rank()==0) myrep=multi_sim_comm.Get_rank();
238 comm.Bcast(nrep,0);
239 comm.Bcast(myrep,0);
240 parseFlag("NO_BROADCAST",no_broadcast);
241 //if(no_broadcast){
242 //for(int irep=0;irep<nrep;irep++){
243 // if(irep!=myrep)
244 // apply_weights[irep]=0.0;}
245 //}
246 avgstep=1.0;
247 tstart=-1.0;
248 tend=-1.0;
249 totalWork=0.0;
250 learn_replica=0;
251
252 parse("LEARN_REPLICA",learn_replica);
253 parseVector("APPLY_WEIGHTS",apply_weights);
254 if(apply_weights.size()==0) apply_weights.resize(nrep,1.0);
255 parseVector("KAPPA",kappa);
256 parseVector("AT",at);
257 parseVector("TAU",tau);
258 parse("TYPE",type);
259 error_type="GAUSSIAN";
260 parse("ERROR_TYPE",error_type);
261 parse("ALPHA",alpha);
262 parse("SIGMA",sigma);
263 parse("TSTART",tstart);
264 if(tstart <0 && tstart != -1.0) error("TSTART should be a positive number");
265 parse("TEND",tend);
266 if(tend<0 && tend != -1.0) error("TSTART should be a positive number");
267 if(tend<tstart) error("TEND should be >= TSTART");
268 lagmultfname=getLabel()+".LAGMULT";
269 parse("FILE",lagmultfname);
270 parse("FMT",fmt);
271 parse("PACE",pace_);
272 if(pace_<=0 ) error("frequency for Lagrangian multipliers update (PACE) is nonsensical");
273 stride_=pace_; //if no STRIDE is passed, then Lagrangian multipliers willbe printed at each update
274 parse("PRINT_STRIDE",stride_);
275 if(stride_<=0 ) error("frequency for Lagrangian multipliers printing (STRIDE) is nonsensical");
276 simtemp=0.;
277 parse("TEMP",simtemp);
278 if(simtemp>0) simtemp*=plumed.getAtoms().getKBoltzmann();
279 else simtemp=plumed.getAtoms().getKbT();
280 parseFlag("REWEIGHT",reweight);
281 if(simtemp<=0 && reweight) error("Set the temperature (TEMP) if you want to do reweighting.");
282
283 checkRead();
284
285 log.printf(" at");
286 for(unsigned i=0; i<at.size(); i++) log.printf(" %f",at[i]);
287 log.printf("\n");
288 log.printf(" with initial learning rate for optimization of");
289 for(unsigned i=0; i<kappa.size(); i++) log.printf(" %f",kappa[i]);
290 log.printf("\n");
291 log.printf("Dumping times for the learning rates are (ps): ");
292 for(unsigned i=0; i<tau.size(); i++) log.printf(" %f",tau[i]);
293 log.printf("\n");
294 log.printf("Lagrangian multipliers are updated every %ld steps (PACE)\n",pace_);
295 log.printf("Lagrangian multipliers output file %s\n",lagmultfname.c_str());
296 log.printf("Lagrangian multipliers are written every %ld steps (PRINT_STRIDE)\n",stride_);
297 if(fmt.length()>0)
298 log.printf("The format for real number in Lagrangian multipliers file is %s\n",fmt.c_str());
299 if(tstart >-1.0 && tend>-1.0)
300 log.printf("Lagrangian multipliers are averaged from %lf ps to %lf ps\n",tstart,tend);
301 if(no_broadcast)
302 log.printf("Using NO_BROADCAST options. Lagrangian multipliers will not be comunicated to other replicas.\n");
303 //for(int irep=0;irep<nrep;irep++){
304 // if(apply_weights[irep]!=0)
305 // log.printf("%d",irep);
306 // }
307 addComponent("force2"); componentIsNotPeriodic("force2");
308 addComponent("work"); componentIsNotPeriodic("work");
309 valueForce2=getPntrToComponent("force2");
310 valueWork=getPntrToComponent("work");
311
312 std::string comp;
313 for(unsigned i=0; i< getNumberOfArguments() ; i++) {
314 comp=getPntrToArgument(i)->getName()+"_coupling";
315 addComponent(comp); componentIsNotPeriodic(comp);
316 comp=getPntrToArgument(i)->getName()+"_work";
317 addComponent(comp); componentIsNotPeriodic(comp);
318 work.push_back(0.); // initialize the work value
319 comp=getPntrToArgument(i)->getName()+"_error";
320 addComponent(comp); componentIsNotPeriodic(comp);
321 }
322 string fname;
323 fname=lagmultfname;
324 ifile.link(*this);
325 if(ifile.FileExist(fname)) {
326 ifile.open(fname);
327 if(getRestart()) {
328 log.printf(" Restarting from: %s\n",fname.c_str());
329 ReadLagrangians(ifile);
330 printFirstStep=false;
331 }
332 ifile.reset(false);
333 }
334
335 lagmultOfile_.link(*this);
336 lagmultOfile_.open(fname);
337 if(fmt.length()>0) {fmt=" "+fmt; lagmultOfile_.fmtField(fmt);}
338 }
339 ////MEMBER FUNCTIONS
ReadLagrangians(IFile & ifile)340 void MaxEnt::ReadLagrangians(IFile &ifile)
341 {
342 double dummy;
343 while(ifile.scanField("time",dummy)) {
344 for(unsigned j=0; j<getNumberOfArguments(); ++j) {
345 ifile.scanField(getPntrToArgument(j)->getName()+"_coupling",lambda[j]);
346 if(dummy>=tstart && dummy <=tend)
347 avglambda[j]+=lambda[j];
348 if(dummy>=tend) {
349 avglambda[j]=lambda[j];
350 done_average[j]=true;
351 }
352 }
353 if(dummy>=tstart && dummy <=tend)
354 avg_counter++;
355 ifile.scanField();
356 }
357 }
WriteLagrangians(vector<double> & lagmult,OFile & file)358 void MaxEnt::WriteLagrangians(vector<double> &lagmult,OFile &file) {
359 if(printFirstStep) {
360 unsigned ncv=getNumberOfArguments();
361 file.printField("time",getTimeStep()*getStep());
362 for(unsigned i=0; i<ncv; ++i)
363 file.printField(getPntrToArgument(i)->getName()+"_coupling",lagmult[i]);
364 file.printField();
365 } else {
366 if(!isFirstStep) {
367 unsigned ncv=getNumberOfArguments();
368 file.printField("time",getTimeStep()*getStep());
369 for(unsigned i=0; i<ncv; ++i)
370 file.printField(getPntrToArgument(i)->getName()+"_coupling",lagmult[i]);
371 file.printField();
372 }
373 }
374 }
compute_error(string & err_type,double & l)375 double MaxEnt::compute_error(string &err_type,double &l) {
376 double sigma2=pow(sigma,2.0);
377 double l2=convert_lambda(type,l);
378 double return_error=0;
379 if(err_type=="GAUSSIAN" && sigma!=0.0)
380 return_error=-l2*sigma2;
381 else {
382 if(err_type=="LAPLACE" && sigma!=0) {
383 return_error=-l2*sigma2/(1.0-l2*l2*sigma2/(alpha+1));
384 }
385 }
386 return return_error;
387 }
convert_lambda(string & type,double lold)388 double MaxEnt::convert_lambda(string &type,double lold) {
389 double return_lambda=0;
390 if(type=="EQUAL")
391 return_lambda=lold;
392 else {
393 if(type=="INEQUAL>") {
394 if(lold>0.0)
395 return_lambda=0.0;
396 else
397 return_lambda=lold;
398 }
399 else {
400 if(type=="INEQUAL<") {
401 if(lold<0.0)
402 return_lambda=0.0;
403 else
404 return_lambda=lold;
405 }
406 }
407 }
408 return return_lambda;
409 }
check_lambda_boundaries(string & err_type,double & l)410 void MaxEnt::check_lambda_boundaries(string &err_type,double &l) {
411 if(err_type=="LAPLACE" && sigma !=0 ) {
412 double l2=convert_lambda(err_type,l);
413 if(l2 <-(sqrt(alpha+1)/sigma-0.01)) {
414 l=-(sqrt(alpha+1)/sigma-0.01);
415 log.printf("Lambda exceeded the allowed range\n");
416 }
417 if(l2>(sqrt(alpha+1)/sigma-0.01)) {
418 l=sqrt(alpha+1)/sigma-0.01;
419 log.printf("Lambda exceeded the allowed range\n");
420 }
421 }
422 }
423
update_lambda()424 void MaxEnt::update_lambda() {
425
426 double totalWork_=0.0;
427 const double time=getTime();
428 const double step=getStep();
429 double KbT=simtemp;
430 double learning_rate;
431 if(reweight)
432 BetaReweightBias=plumed.getBias()/KbT;
433 else
434 BetaReweightBias=0.0;
435
436 for(unsigned i=0; i<getNumberOfArguments(); ++i) {
437 const double k=kappa[i];
438 double cv=(getArgument(i)+compute_error(error_type,lambda[i])-at[i]);
439 if(reweight)
440 learning_rate=1.0*k/(1+step/tau[i]);
441 else
442 learning_rate=1.0*k/(1+time/tau[i]);
443 lambda[i]+=learning_rate*cv*exp(-BetaReweightBias); //update Lagrangian multipliers and reweight them if REWEIGHT is set
444 check_lambda_boundaries(error_type,lambda[i]); //check that Lagrangians multipliers not exceed the allowed range
445 if(time>=tstart && time <=tend && !done_average[i]) {
446 avglambda[i]+=convert_lambda(type,lambda[i]); //compute the average of Lagrangian multipliers over the required time window
447 }
448 if(time>=tend && tend >=0) { //setting tend<0 will disable this feature
449 if(!done_average[i]) {
450 avglambda[i]=avglambda[i]/avg_counter;
451 done_average[i]=true;
452 lambda[i]=avglambda[i];
453 }
454 else
455 lambda[i]=avglambda[i]; //keep Lagrangian multipliers fixed to the previously computed average.
456 }
457 work[i]+=(convert_lambda(type,lambda[i])-oldlambda[i])*getArgument(i); //compute the work performed in updating lambda
458 totalWork_+=work[i];
459 totalWork=totalWork_;
460 oldlambda[i]=convert_lambda(type,lambda[i]);
461 };
462 if(time>=tstart && time <=tend)
463 avg_counter++;
464 }
465
calculate()466 void MaxEnt::calculate() {
467 double totf2=0.0;
468 double ene=0.0;
469 double KbT=simtemp;
470 for(unsigned i=0; i<getNumberOfArguments(); ++i) {
471 getPntrToComponent(getPntrToArgument(i)->getName()+"_error")->set(expected_eps[i]);
472 getPntrToComponent(getPntrToArgument(i)->getName()+"_work")->set(work[i]);
473 valueWork->set(totalWork);
474 getPntrToComponent(getPntrToArgument(i)->getName()+"_coupling")->set(lambda[i]);
475 const double f=-KbT*convert_lambda(type,lambda[i])*apply_weights[myrep];
476 totf2+=f*f;
477 ene+=KbT*convert_lambda(type,lambda[i])*getArgument(i)*apply_weights[myrep];
478 setOutputForce(i,f);
479 }
480 setBias(ene);
481 valueForce2->set(totf2);
482 }
483
update()484 void MaxEnt::update() {
485
486 if(getStep()%stride_ == 0)
487 WriteLagrangians(lambda,lagmultOfile_);
488 if(getStep()%pace_ == 0) {
489 update_lambda();
490 if(!no_broadcast) {
491 if(comm.Get_rank()==0) //Comunicate Lagrangian multipliers from reference replica to higher ones
492 multi_sim_comm.Bcast(lambda,learn_replica);
493 }
494 comm.Bcast(lambda,0);
495 }
496 isFirstStep=false;
497 }
498
499 }
500
501 }
502