Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2020 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 200 : 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();
167 : void update();
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 7456 : PLUMED_REGISTER_ACTION(MaxEnt,"MAXENT")
177 :
178 51 : void MaxEnt::registerKeywords(Keywords& keys) {
179 51 : Bias::registerKeywords(keys);
180 51 : componentsAreNotOptional(keys);
181 102 : keys.use("ARG");
182 255 : keys.add("compulsory","KAPPA","0.0","specifies the initial value for the learning rate");
183 204 : keys.add("compulsory","TAU","Specify the dumping time for the learning rate.");
184 204 : 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 204 : 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 204 : 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 204 : 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 204 : 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 204 : keys.add("compulsory","AT","the position of the restraint");
195 204 : keys.add("optional","SIGMA","The typical errors expected on observable");
196 204 : keys.add("optional","FILE","Lagrangian multipliers output file. The default name is: label name followed by the string .LAGMULT ");
197 204 : keys.add("optional","LEARN_REPLICA","In a multiple replica environment specify which is the reference replica. By default replica 0 will be used.");
198 204 : 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 204 : keys.add("optional","PACE","the frequency for Lagrangian multipliers update");
200 204 : 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 204 : keys.add("optional","FMT","specify format for Lagrangian multipliers files (useful to decrease the number of digits in regtests)");
202 153 : keys.addFlag("REWEIGHT",false,"to be used with plumed driver in order to reweight a trajectory a posteriori");
203 153 : keys.addFlag("NO_BROADCAST",false,"If active will avoid Lagrangian multipliers to be communicated to other replicas.");
204 204 : keys.add("optional","TEMP","the system temperature. This is required if you are reweighting.");
205 204 : keys.addOutputComponent("force2","default","the instantaneous value of the squared force due to this bias potential");
206 204 : keys.addOutputComponent("work","default","the instantaneous value of the work done by the biasing force");
207 204 : 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 204 : keys.addOutputComponent("_error","default","Instantaneous values of the discrepancy between the observable and the restraint center");
211 204 : keys.addOutputComponent("_coupling","default","Instantaneous values of Lagrangian multipliers. They are also written by default in a separate output file.");
212 102 : keys.use("RESTART");
213 51 : }
214 50 : 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 100 : 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 950 : done_average(getNumberOfArguments(),false)
235 : {
236 50 : if(comm.Get_rank()==0) nrep=multi_sim_comm.Get_size();
237 50 : if(comm.Get_rank()==0) myrep=multi_sim_comm.Get_rank();
238 50 : comm.Bcast(nrep,0);
239 50 : comm.Bcast(myrep,0);
240 100 : 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 50 : avgstep=1.0;
247 50 : tstart=-1.0;
248 50 : tend=-1.0;
249 50 : totalWork=0.0;
250 50 : learn_replica=0;
251 :
252 100 : parse("LEARN_REPLICA",learn_replica);
253 100 : parseVector("APPLY_WEIGHTS",apply_weights);
254 50 : if(apply_weights.size()==0) apply_weights.resize(nrep,1.0);
255 100 : parseVector("KAPPA",kappa);
256 100 : parseVector("AT",at);
257 100 : parseVector("TAU",tau);
258 100 : parse("TYPE",type);
259 : error_type="GAUSSIAN";
260 100 : parse("ERROR_TYPE",error_type);
261 100 : parse("ALPHA",alpha);
262 100 : parse("SIGMA",sigma);
263 100 : parse("TSTART",tstart);
264 50 : if(tstart <0 && tstart != -1.0) error("TSTART should be a positive number");
265 100 : parse("TEND",tend);
266 50 : if(tend<0 && tend != -1.0) error("TSTART should be a positive number");
267 50 : if(tend<tstart) error("TEND should be >= TSTART");
268 100 : lagmultfname=getLabel()+".LAGMULT";
269 100 : parse("FILE",lagmultfname);
270 100 : parse("FMT",fmt);
271 100 : parse("PACE",pace_);
272 50 : if(pace_<=0 ) error("frequency for Lagrangian multipliers update (PACE) is nonsensical");
273 50 : stride_=pace_; //if no STRIDE is passed, then Lagrangian multipliers willbe printed at each update
274 100 : parse("PRINT_STRIDE",stride_);
275 50 : if(stride_<=0 ) error("frequency for Lagrangian multipliers printing (STRIDE) is nonsensical");
276 50 : simtemp=0.;
277 100 : parse("TEMP",simtemp);
278 100 : if(simtemp>0) simtemp*=plumed.getAtoms().getKBoltzmann();
279 0 : else simtemp=plumed.getAtoms().getKbT();
280 100 : parseFlag("REWEIGHT",reweight);
281 50 : if(simtemp<=0 && reweight) error("Set the temperature (TEMP) if you want to do reweighting.");
282 :
283 50 : checkRead();
284 :
285 50 : log.printf(" at");
286 1546 : for(unsigned i=0; i<at.size(); i++) log.printf(" %f",at[i]);
287 50 : log.printf("\n");
288 50 : log.printf(" with initial learning rate for optimization of");
289 1546 : for(unsigned i=0; i<kappa.size(); i++) log.printf(" %f",kappa[i]);
290 50 : log.printf("\n");
291 50 : log.printf("Dumping times for the learning rates are (ps): ");
292 1546 : for(unsigned i=0; i<tau.size(); i++) log.printf(" %f",tau[i]);
293 50 : log.printf("\n");
294 50 : log.printf("Lagrangian multipliers are updated every %ld steps (PACE)\n",pace_);
295 100 : log.printf("Lagrangian multipliers output file %s\n",lagmultfname.c_str());
296 50 : log.printf("Lagrangian multipliers are written every %ld steps (PRINT_STRIDE)\n",stride_);
297 50 : if(fmt.length()>0)
298 100 : log.printf("The format for real number in Lagrangian multipliers file is %s\n",fmt.c_str());
299 50 : if(tstart >-1.0 && tend>-1.0)
300 14 : log.printf("Lagrangian multipliers are averaged from %lf ps to %lf ps\n",tstart,tend);
301 50 : if(no_broadcast)
302 0 : 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 150 : addComponent("force2"); componentIsNotPeriodic("force2");
308 150 : addComponent("work"); componentIsNotPeriodic("work");
309 100 : valueForce2=getPntrToComponent("force2");
310 100 : valueWork=getPntrToComponent("work");
311 :
312 : std::string comp;
313 1014 : for(unsigned i=0; i< getNumberOfArguments() ; i++) {
314 964 : comp=getPntrToArgument(i)->getName()+"_coupling";
315 482 : addComponent(comp); componentIsNotPeriodic(comp);
316 964 : comp=getPntrToArgument(i)->getName()+"_work";
317 482 : addComponent(comp); componentIsNotPeriodic(comp);
318 964 : work.push_back(0.); // initialize the work value
319 964 : comp=getPntrToArgument(i)->getName()+"_error";
320 482 : addComponent(comp); componentIsNotPeriodic(comp);
321 : }
322 : string fname;
323 : fname=lagmultfname;
324 50 : ifile.link(*this);
325 50 : if(ifile.FileExist(fname)) {
326 37 : ifile.open(fname);
327 37 : if(getRestart()) {
328 74 : log.printf(" Restarting from: %s\n",fname.c_str());
329 37 : ReadLagrangians(ifile);
330 37 : printFirstStep=false;
331 : }
332 37 : ifile.reset(false);
333 : }
334 :
335 50 : lagmultOfile_.link(*this);
336 50 : lagmultOfile_.open(fname);
337 100 : if(fmt.length()>0) {fmt=" "+fmt; lagmultOfile_.fmtField(fmt);}
338 50 : }
339 : ////MEMBER FUNCTIONS
340 37 : void MaxEnt::ReadLagrangians(IFile &ifile)
341 : {
342 : double dummy;
343 888 : while(ifile.scanField("time",dummy)) {
344 9009 : for(unsigned j=0; j<getNumberOfArguments(); ++j) {
345 12903 : ifile.scanField(getPntrToArgument(j)->getName()+"_coupling",lambda[j]);
346 4301 : if(dummy>=tstart && dummy <=tend)
347 84 : avglambda[j]+=lambda[j];
348 4301 : if(dummy>=tend) {
349 4231 : avglambda[j]=lambda[j];
350 : done_average[j]=true;
351 : }
352 : }
353 407 : if(dummy>=tstart && dummy <=tend)
354 6 : avg_counter++;
355 407 : ifile.scanField();
356 : }
357 37 : }
358 550 : void MaxEnt::WriteLagrangians(vector<double> &lagmult,OFile &file) {
359 550 : if(printFirstStep) {
360 143 : unsigned ncv=getNumberOfArguments();
361 286 : file.printField("time",getTimeStep()*getStep());
362 2145 : for(unsigned i=0; i<ncv; ++i)
363 4004 : file.printField(getPntrToArgument(i)->getName()+"_coupling",lagmult[i]);
364 143 : file.printField();
365 : } else {
366 407 : if(!isFirstStep) {
367 370 : unsigned ncv=getNumberOfArguments();
368 740 : file.printField("time",getTimeStep()*getStep());
369 8190 : for(unsigned i=0; i<ncv; ++i)
370 15640 : file.printField(getPntrToArgument(i)->getName()+"_coupling",lagmult[i]);
371 370 : file.printField();
372 : }
373 : }
374 550 : }
375 5302 : double MaxEnt::compute_error(string &err_type,double &l) {
376 5302 : double sigma2=pow(sigma,2.0);
377 5302 : double l2=convert_lambda(type,l);
378 : double return_error=0;
379 5302 : if(err_type=="GAUSSIAN" && sigma!=0.0)
380 0 : return_error=-l2*sigma2;
381 : else {
382 5302 : if(err_type=="LAPLACE" && sigma!=0) {
383 5302 : return_error=-l2*sigma2/(1.0-l2*l2*sigma2/(alpha+1));
384 : }
385 : }
386 5302 : return return_error;
387 : }
388 119118 : double MaxEnt::convert_lambda(string &type,double lold) {
389 : double return_lambda=0;
390 119118 : if(type=="EQUAL")
391 : return_lambda=lold;
392 : else {
393 5302 : if(type=="INEQUAL>") {
394 0 : if(lold>0.0)
395 : return_lambda=0.0;
396 : else
397 : return_lambda=lold;
398 : }
399 : else {
400 5302 : if(type=="INEQUAL<") {
401 0 : if(lold<0.0)
402 : return_lambda=0.0;
403 : else
404 : return_lambda=lold;
405 : }
406 : }
407 : }
408 119118 : return return_lambda;
409 : }
410 5302 : void MaxEnt::check_lambda_boundaries(string &err_type,double &l) {
411 5302 : if(err_type=="LAPLACE" && sigma !=0 ) {
412 5302 : double l2=convert_lambda(err_type,l);
413 5302 : if(l2 <-(sqrt(alpha+1)/sigma-0.01)) {
414 0 : l=-(sqrt(alpha+1)/sigma-0.01);
415 0 : log.printf("Lambda exceeded the allowed range\n");
416 : }
417 5302 : if(l2>(sqrt(alpha+1)/sigma-0.01)) {
418 0 : l=sqrt(alpha+1)/sigma-0.01;
419 0 : log.printf("Lambda exceeded the allowed range\n");
420 : }
421 : }
422 5302 : }
423 :
424 550 : void MaxEnt::update_lambda() {
425 :
426 : double totalWork_=0.0;
427 550 : const double time=getTime();
428 550 : const double step=getStep();
429 550 : double KbT=simtemp;
430 : double learning_rate;
431 550 : if(reweight)
432 396 : BetaReweightBias=plumed.getBias()/KbT;
433 : else
434 154 : BetaReweightBias=0.0;
435 :
436 11154 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
437 10604 : const double k=kappa[i];
438 10604 : double cv=(getArgument(i)+compute_error(error_type,lambda[i])-at[i]);
439 5302 : if(reweight)
440 4224 : learning_rate=1.0*k/(1+step/tau[i]);
441 : else
442 1078 : learning_rate=1.0*k/(1+time/tau[i]);
443 10604 : lambda[i]+=learning_rate*cv*exp(-BetaReweightBias); //update Lagrangian multipliers and reweight them if REWEIGHT is set
444 5302 : check_lambda_boundaries(error_type,lambda[i]); //check that Lagrangians multipliers not exceed the allowed range
445 5890 : if(time>=tstart && time <=tend && !done_average[i]) {
446 1092 : avglambda[i]+=convert_lambda(type,lambda[i]); //compute the average of Lagrangian multipliers over the required time window
447 : }
448 5302 : if(time>=tend && tend >=0) { //setting tend<0 will disable this feature
449 98 : if(!done_average[i]) {
450 91 : avglambda[i]=avglambda[i]/avg_counter;
451 : done_average[i]=true;
452 91 : lambda[i]=avglambda[i];
453 : }
454 : else
455 7 : lambda[i]=avglambda[i]; //keep Lagrangian multipliers fixed to the previously computed average.
456 : }
457 21208 : work[i]+=(convert_lambda(type,lambda[i])-oldlambda[i])*getArgument(i); //compute the work performed in updating lambda
458 5302 : totalWork_+=work[i];
459 5302 : totalWork=totalWork_;
460 10604 : oldlambda[i]=convert_lambda(type,lambda[i]);
461 : };
462 550 : if(time>=tstart && time <=tend)
463 84 : avg_counter++;
464 550 : }
465 :
466 5050 : void MaxEnt::calculate() {
467 : double totf2=0.0;
468 : double ene=0.0;
469 5050 : double KbT=simtemp;
470 102414 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
471 194728 : getPntrToComponent(getPntrToArgument(i)->getName()+"_error")->set(expected_eps[i]);
472 146046 : getPntrToComponent(getPntrToArgument(i)->getName()+"_work")->set(work[i]);
473 48682 : valueWork->set(totalWork);
474 146046 : getPntrToComponent(getPntrToArgument(i)->getName()+"_coupling")->set(lambda[i]);
475 146046 : const double f=-KbT*convert_lambda(type,lambda[i])*apply_weights[myrep];
476 48682 : totf2+=f*f;
477 146046 : ene+=KbT*convert_lambda(type,lambda[i])*getArgument(i)*apply_weights[myrep];
478 48682 : setOutputForce(i,f);
479 : }
480 : setBias(ene);
481 5050 : valueForce2->set(totf2);
482 5050 : }
483 :
484 5050 : void MaxEnt::update() {
485 :
486 5050 : if(getStep()%stride_ == 0)
487 550 : WriteLagrangians(lambda,lagmultOfile_);
488 5050 : if(getStep()%pace_ == 0) {
489 550 : update_lambda();
490 550 : if(!no_broadcast) {
491 550 : if(comm.Get_rank()==0) //Comunicate Lagrangian multipliers from reference replica to higher ones
492 462 : multi_sim_comm.Bcast(lambda,learn_replica);
493 : }
494 550 : comm.Bcast(lambda,0);
495 : }
496 5050 : isFirstStep=false;
497 5050 : }
498 :
499 : }
500 :
501 5517 : }
|
