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