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