Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2017-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 : /*
23 :
24 : */
25 : #include "bias/Bias.h"
26 : #include "bias/ActionRegister.h"
27 : #include "core/PlumedMain.h"
28 : #include "core/Atoms.h"
29 : #include "core/Value.h"
30 : #include "tools/File.h"
31 : #include "tools/Random.h"
32 : #include <ctime>
33 :
34 : namespace PLMD {
35 : namespace isdb {
36 :
37 : //+PLUMEDOC ISDB_BIAS RESCALE
38 : /*
39 : Scales the value of an another action, being a Collective Variable or a Bias.
40 :
41 : The rescaling factor is determined by a parameter defined on a logarithmic grid of dimension NBIN in the range
42 : from 1 to MAX_RESCALE. The current value of the rescaling parameter is stored and shared across
43 : other actions using a \ref SELECTOR. A Monte Carlo procedure is used to update the value
44 : of the rescaling factor every MC_STRIDE steps of molecular dynamics. Well-tempered metadynamics, defined by the
45 : parameters W0 and BIASFACTOR, is used to enhance the sampling in the space of the rescaling factor.
46 : The well-tempered metadynamics bias potential is written to the file BFILE every BSTRIDE steps and read
47 : when restarting the simulation using the directive \ref RESTART.
48 :
49 : \note
50 : Additional arguments not to be scaled, one for each bin in the rescaling parameter ladder, can be
51 : provided at the end of the ARG list. The number of such arguments is specified by the option NOT_RESCALED.
52 : These arguments will be not be scaled, but they will be
53 : considered as bias potentials and used in the computation of the Metropolis
54 : acceptance probability when proposing a move in the rescaling parameter. See example below.
55 :
56 : \note
57 : If PLUMED is running in a multiple-replica framework (for example using the -multi option in GROMACS),
58 : the arguments will be summed across replicas, unless the NOT_SHARED option is used. Also, the value of the
59 : \ref SELECTOR will be shared and thus will be the same in all replicas.
60 :
61 : \par Examples
62 :
63 : In this example we use \ref RESCALE to implement a simulated-tempering like approach.
64 : The total potential energy of the system is scaled by a parameter defined on a logarithmic grid
65 : of 5 bins in the range from 1 to 1.5.
66 : A well-tempered metadynamics bias potential is used to ensure diffusion in the space of the rescaling
67 : parameter.
68 :
69 : \plumedfile
70 : ene: ENERGY
71 :
72 : SELECTOR NAME=GAMMA VALUE=0
73 :
74 : RESCALE ...
75 : LABEL=res ARG=ene TEMP=300
76 : SELECTOR=GAMMA MAX_RESCALE=1.5 NBIN=5
77 : W0=1000 BIASFACTOR=100.0 BSTRIDE=2000 BFILE=bias.dat
78 : ...
79 :
80 : PRINT FILE=COLVAR ARG=* STRIDE=100
81 : \endplumedfile
82 :
83 : In this second example, we add to the simulated-tempering approach introduced above
84 : one Parallel Bias metadynamics simulation (see \ref PBMETAD) for each value of the rescaling parameter.
85 : At each moment of the simulation, only one of the \ref PBMETAD
86 : actions is activated, based on the current value of the associated \ref SELECTOR.
87 : The \ref PBMETAD bias potentials are not scaled, but just used in the calculation of
88 : the Metropolis acceptance probability when proposing a move in the rescaling parameter.
89 :
90 : \plumedfile
91 : ene: ENERGY
92 : d: DISTANCE ATOMS=1,2
93 :
94 : SELECTOR NAME=GAMMA VALUE=0
95 :
96 : pbmetad0: PBMETAD ARG=d SELECTOR=GAMMA SELECTOR_ID=0 SIGMA=0.1 PACE=500 HEIGHT=1 BIASFACTOR=8 FILE=HILLS.0
97 : pbmetad1: PBMETAD ARG=d SELECTOR=GAMMA SELECTOR_ID=1 SIGMA=0.1 PACE=500 HEIGHT=1 BIASFACTOR=8 FILE=HILLS.1
98 : pbmetad2: PBMETAD ARG=d SELECTOR=GAMMA SELECTOR_ID=2 SIGMA=0.1 PACE=500 HEIGHT=1 BIASFACTOR=8 FILE=HILLS.2
99 : pbmetad3: PBMETAD ARG=d SELECTOR=GAMMA SELECTOR_ID=3 SIGMA=0.1 PACE=500 HEIGHT=1 BIASFACTOR=8 FILE=HILLS.3
100 : pbmetad4: PBMETAD ARG=d SELECTOR=GAMMA SELECTOR_ID=4 SIGMA=0.1 PACE=500 HEIGHT=1 BIASFACTOR=8 FILE=HILLS.4
101 :
102 : RESCALE ...
103 : LABEL=res TEMP=300
104 : ARG=ene,pbmetad0.bias,pbmetad1.bias,pbmetad2.bias,pbmetad3.bias,pbmetad4.bias
105 : SELECTOR=GAMMA MAX_RESCALE=1.5 NOT_RESCALED=5 NBIN=5
106 : W0=1000 BIASFACTOR=100.0 BSTRIDE=2000 BFILE=bias.dat
107 : ...
108 :
109 : PRINT FILE=COLVAR ARG=* STRIDE=100
110 : \endplumedfile
111 :
112 :
113 :
114 : */
115 : //+ENDPLUMEDOC
116 :
117 : class Rescale : public bias::Bias {
118 : // gamma parameter
119 : std::vector<double> gamma_;
120 : double w0_;
121 : double biasf_;
122 : std::vector<double> bias_;
123 : std::vector<double> expo_;
124 : std::vector<unsigned> shared_;
125 : unsigned nores_;
126 : // bias
127 : unsigned int Biasstride_;
128 : unsigned int Biaspace_;
129 : std::string Biasfilename_;
130 : bool first_bias_;
131 : OFile Biasfile_;
132 : // temperature in kbt
133 : double kbt_;
134 : // Monte Carlo stuff
135 : unsigned MCsteps_;
136 : unsigned MCstride_;
137 : long long int MCfirst_;
138 : long long unsigned MCaccgamma_;
139 : // replica stuff
140 : unsigned nrep_;
141 : unsigned replica_;
142 : // selector
143 : std::string selector_;
144 :
145 : // Monte Carlo
146 : void doMonteCarlo(unsigned igamma, double oldE, std::vector<double> args, std::vector<double> bargs);
147 : unsigned proposeMove(unsigned x, unsigned xmin, unsigned xmax);
148 : bool doAccept(double oldE, double newE);
149 : // read and print bias
150 : void read_bias();
151 : void print_bias(long long int step);
152 :
153 : public:
154 : explicit Rescale(const ActionOptions&);
155 : ~Rescale();
156 : void calculate();
157 : static void registerKeywords(Keywords& keys);
158 : };
159 :
160 :
161 13785 : PLUMED_REGISTER_ACTION(Rescale,"RESCALE")
162 :
163 4 : void Rescale::registerKeywords(Keywords& keys) {
164 4 : Bias::registerKeywords(keys);
165 4 : keys.use("ARG");
166 8 : keys.add("compulsory","TEMP","temperature");
167 8 : keys.add("compulsory","SELECTOR", "name of the SELECTOR used for rescaling");
168 8 : keys.add("compulsory","MAX_RESCALE","maximum values for rescaling");
169 8 : keys.add("compulsory","NBIN","number of bins for gamma grid");
170 8 : keys.add("compulsory","W0", "initial bias height");
171 8 : keys.add("compulsory","BIASFACTOR", "bias factor");
172 8 : keys.add("compulsory","BSTRIDE", "stride for writing bias");
173 8 : keys.add("compulsory","BFILE", "file name for bias");
174 8 : keys.add("optional","NOT_SHARED", "list of arguments (from 1 to N) not summed across replicas");
175 8 : keys.add("optional","NOT_RESCALED", "these last N arguments will not be scaled");
176 8 : keys.add("optional","MC_STEPS","number of MC steps");
177 8 : keys.add("optional","MC_STRIDE","MC stride");
178 8 : keys.add("optional","PACE", "Pace for adding bias, in MC stride unit");
179 4 : componentsAreNotOptional(keys);
180 8 : keys.addOutputComponent("igamma", "default","gamma parameter");
181 8 : keys.addOutputComponent("accgamma","default","MC acceptance for gamma");
182 8 : keys.addOutputComponent("wtbias", "default","well-tempered bias");
183 4 : }
184 :
185 0 : Rescale::Rescale(const ActionOptions&ao):
186 : PLUMED_BIAS_INIT(ao),
187 0 : nores_(0), Biaspace_(1), first_bias_(true),
188 0 : MCsteps_(1), MCstride_(1), MCfirst_(-1), MCaccgamma_(0) {
189 : // set up replica stuff
190 0 : if(comm.Get_rank()==0) {
191 0 : nrep_ = multi_sim_comm.Get_size();
192 0 : replica_ = multi_sim_comm.Get_rank();
193 : } else {
194 0 : nrep_ = 0;
195 0 : replica_ = 0;
196 : }
197 0 : comm.Sum(&nrep_,1);
198 0 : comm.Sum(&replica_,1);
199 :
200 : // wt-parameters
201 0 : parse("W0", w0_);
202 0 : parse("BIASFACTOR", biasf_);
203 :
204 : // selector name
205 0 : parse("SELECTOR", selector_);
206 :
207 : // number of bins for gamma ladder
208 : unsigned nbin;
209 0 : parse("NBIN", nbin);
210 :
211 : // number of bias
212 0 : parse("NOT_RESCALED", nores_);
213 0 : if(nores_>0 && nores_!=nbin) {
214 0 : error("The number of non scaled arguments must be equal to either 0 or the number of bins");
215 : }
216 :
217 : // maximum value of rescale
218 : std::vector<double> max_rescale;
219 0 : parseVector("MAX_RESCALE", max_rescale);
220 : // check dimension of max_rescale
221 0 : if(max_rescale.size()!=(getNumberOfArguments()-nores_)) {
222 0 : error("Size of MAX_RESCALE array must be equal to the number of arguments that will to be scaled");
223 : }
224 :
225 : // calculate exponents
226 0 : double igamma_max = static_cast<double>(nbin);
227 0 : for(unsigned i=0; i<max_rescale.size(); ++i) {
228 0 : expo_.push_back(std::log(max_rescale[i])/std::log(igamma_max));
229 : }
230 :
231 : // allocate gamma grid and set bias to zero
232 0 : for(unsigned i=0; i<nbin; ++i) {
233 : // bias grid
234 0 : bias_.push_back(0.0);
235 : // gamma ladder
236 0 : double gamma = std::exp( static_cast<double>(i) / static_cast<double>(nbin-1) * std::log(igamma_max) );
237 0 : gamma_.push_back(gamma);
238 : }
239 : // print bias to file
240 0 : parse("BSTRIDE", Biasstride_);
241 0 : parse("BFILE", Biasfilename_);
242 :
243 : // create vectors of shared arguments
244 : // by default they are all shared
245 0 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
246 0 : shared_.push_back(1);
247 : }
248 : // share across replicas or not
249 : std::vector<unsigned> not_shared;
250 0 : parseVector("NOT_SHARED", not_shared);
251 : // and change the non-shared
252 0 : for(unsigned i=0; i<not_shared.size(); ++i) {
253 0 : if((not_shared[i]-1)>=(getNumberOfArguments()-nores_) && nrep_>1) {
254 0 : error("NOT_RESCALED args must always be shared when using multiple replicas");
255 : }
256 0 : if((not_shared[i]-1)>=getNumberOfArguments()) {
257 0 : error("NOT_SHARED args should be lower than total number of arguments");
258 : }
259 0 : shared_[not_shared[i]-1] = 0;
260 : }
261 :
262 : // monte carlo stuff
263 0 : parse("MC_STEPS",MCsteps_);
264 0 : parse("MC_STRIDE",MCstride_);
265 : // adjust for multiple-time steps
266 0 : MCstride_ *= getStride();
267 : // read bias deposition pace
268 0 : parse("PACE", Biaspace_);
269 : // multiply by MCstride
270 0 : Biaspace_ *= MCstride_;
271 :
272 : // get temperature
273 0 : double temp=0.0;
274 0 : parse("TEMP",temp);
275 0 : if(temp>0.0) {
276 0 : kbt_=plumed.getAtoms().getKBoltzmann()*temp;
277 : } else {
278 0 : kbt_=plumed.getAtoms().getKbT();
279 : }
280 :
281 0 : checkRead();
282 :
283 0 : log.printf(" temperature of the system in energy unit %f\n",kbt_);
284 0 : log.printf(" name of the SELECTOR use for this action %s\n",selector_.c_str());
285 0 : log.printf(" number of bins in grid %u\n",nbin);
286 0 : log.printf(" number of arguments that will not be scaled %u\n",nores_);
287 0 : if(nrep_>1) {
288 0 : log<<" number of arguments that will not be summed across replicas "<<not_shared.size()<<"\n";
289 : }
290 0 : log.printf(" biasfactor %f\n",biasf_);
291 0 : log.printf(" initial hills height %f\n",w0_);
292 0 : log.printf(" stride to write bias to file %u\n",Biasstride_);
293 0 : log.printf(" write bias to file : %s\n",Biasfilename_.c_str());
294 0 : log.printf(" number of replicas %u\n",nrep_);
295 0 : log.printf(" number of MC steps %d\n",MCsteps_);
296 0 : log.printf(" do MC every %d steps\n", MCstride_);
297 0 : log.printf("\n");
298 :
299 0 : log << " Bibliography" << plumed.cite("Bonomi, Camilloni, Bioinformatics, 33, 3999 (2017)") << "\n";
300 :
301 :
302 : // add components
303 0 : addComponent("igamma");
304 0 : componentIsNotPeriodic("igamma");
305 0 : addComponent("accgamma");
306 0 : componentIsNotPeriodic("accgamma");
307 0 : addComponent("wtbias");
308 0 : componentIsNotPeriodic("wtbias");
309 :
310 : // initialize random seed
311 0 : srand (time(NULL));
312 :
313 : // read bias if restarting
314 0 : if(getRestart()) {
315 0 : read_bias();
316 : }
317 0 : }
318 :
319 0 : Rescale::~Rescale() {
320 0 : Biasfile_.close();
321 0 : }
322 :
323 0 : void Rescale::read_bias() {
324 : // open file
325 0 : auto ifile=Tools::make_unique<IFile>();
326 0 : ifile->link(*this);
327 0 : if(ifile->FileExist(Biasfilename_)) {
328 0 : ifile->open(Biasfilename_);
329 : // read all the lines, store last value of bias
330 : double MDtime;
331 0 : while(ifile->scanField("MD_time",MDtime)) {
332 0 : for(unsigned i=0; i<bias_.size(); ++i) {
333 : // convert i to string
334 0 : std::stringstream ss;
335 : ss << i;
336 : // label
337 0 : std::string label = "b" + ss.str();
338 : // read entry
339 0 : ifile->scanField(label, bias_[i]);
340 0 : }
341 : // new line
342 0 : ifile->scanField();
343 : }
344 0 : ifile->close();
345 : } else {
346 0 : error("Cannot find bias file "+Biasfilename_+"\n");
347 : }
348 0 : }
349 :
350 0 : unsigned Rescale::proposeMove(unsigned x, unsigned xmin, unsigned xmax) {
351 0 : int xmin_i = static_cast<int>(xmin);
352 0 : int xmax_i = static_cast<int>(xmax);
353 : int dx;
354 0 : int r = rand() % 2;
355 0 : if( r % 2 == 0 ) {
356 : dx = +1;
357 : } else {
358 : dx = -1;
359 : }
360 : // new index, integer
361 0 : int x_new = static_cast<int>(x) + dx;
362 : // check boundaries
363 0 : if(x_new >= xmax_i) {
364 0 : x_new = xmax_i-1;
365 : }
366 : if(x_new < xmin_i) {
367 : x_new = xmin_i;
368 : }
369 0 : return static_cast<unsigned>(x_new);
370 : }
371 :
372 0 : bool Rescale::doAccept(double oldE, double newE) {
373 : bool accept = false;
374 : // calculate delta energy
375 0 : double delta = ( newE - oldE ) / kbt_;
376 : // if delta is negative always accept move
377 0 : if( delta < 0.0 ) {
378 : accept = true;
379 : } else {
380 : // otherwise extract random number
381 0 : double s = static_cast<double>(rand()) / RAND_MAX;
382 0 : if( s < std::exp(-delta) ) {
383 : accept = true;
384 : }
385 : }
386 0 : return accept;
387 : }
388 :
389 0 : void Rescale::doMonteCarlo(unsigned igamma, double oldE,
390 : std::vector<double> args, std::vector<double> bargs) {
391 : double oldB, newB;
392 :
393 : // cycle on MC steps
394 0 : for(unsigned i=0; i<MCsteps_; ++i) {
395 : // propose move in igamma
396 0 : unsigned new_igamma = proposeMove(igamma, 0, gamma_.size());
397 : // calculate new energy
398 : double newE = 0.0;
399 0 : for(unsigned j=0; j<args.size(); ++j) {
400 : // calculate energy term
401 0 : double fact = 1.0/pow(gamma_[new_igamma], expo_[j]) - 1.0;
402 0 : newE += args[j] * fact;
403 : }
404 : // calculate contributions from non-rescaled terms
405 0 : if(bargs.size()>0) {
406 0 : oldB = bias_[igamma]+bargs[igamma];
407 0 : newB = bias_[new_igamma]+bargs[new_igamma];
408 : } else {
409 0 : oldB = bias_[igamma];
410 0 : newB = bias_[new_igamma];
411 : }
412 : // accept or reject
413 0 : bool accept = doAccept(oldE+oldB, newE+newB);
414 0 : if(accept) {
415 0 : igamma = new_igamma;
416 : oldE = newE;
417 0 : MCaccgamma_++;
418 : }
419 : }
420 : // send values of gamma to all replicas
421 0 : if(comm.Get_rank()==0) {
422 0 : if(multi_sim_comm.Get_rank()!=0) {
423 0 : igamma = 0;
424 : }
425 0 : multi_sim_comm.Sum(&igamma, 1);
426 : } else {
427 0 : igamma = 0;
428 : }
429 : // local communication
430 0 : comm.Sum(&igamma, 1);
431 :
432 : // set the value of gamma into passMap
433 0 : plumed.passMap[selector_]=static_cast<double>(igamma);
434 0 : }
435 :
436 0 : void Rescale::print_bias(long long int step) {
437 : // if first time open the file
438 0 : if(first_bias_) {
439 0 : first_bias_ = false;
440 0 : Biasfile_.link(*this);
441 0 : Biasfile_.open(Biasfilename_);
442 : Biasfile_.setHeavyFlush();
443 0 : Biasfile_.fmtField("%30.5f");
444 : }
445 :
446 : // write fields
447 0 : double MDtime = static_cast<double>(step)*getTimeStep();
448 0 : Biasfile_.printField("MD_time", MDtime);
449 0 : for(unsigned i=0; i<bias_.size(); ++i) {
450 : // convert i to string
451 0 : std::stringstream ss;
452 : ss << i;
453 : // label
454 0 : std::string label = "b" + ss.str();
455 : // print entry
456 0 : Biasfile_.printField(label, bias_[i]);
457 0 : }
458 0 : Biasfile_.printField();
459 0 : }
460 :
461 0 : void Rescale::calculate() {
462 : // get the current value of the selector
463 0 : unsigned igamma = static_cast<unsigned>(plumed.passMap[selector_]);
464 :
465 : // collect data from other replicas
466 0 : std::vector<double> all_args(getNumberOfArguments(), 0.0);
467 : // first calculate arguments
468 0 : for(unsigned i=0; i<all_args.size(); ++i) {
469 0 : double arg = getArgument(i);
470 : // sum shared arguments across replicas
471 0 : if(shared_[i]==1) {
472 0 : if(comm.Get_rank()==0) {
473 0 : multi_sim_comm.Sum(arg);
474 : } else {
475 0 : arg = 0.0;
476 : }
477 0 : if(comm.Get_size()>1) {
478 0 : comm.Sum(arg);
479 : }
480 : }
481 : // put into all_args
482 0 : all_args[i] = arg;
483 : }
484 :
485 : // now separate terms that should be rescaled
486 : std::vector<double> args;
487 0 : if(getNumberOfArguments()-nores_>0) {
488 0 : args.resize(getNumberOfArguments()-nores_);
489 : }
490 0 : for(unsigned i=0; i<args.size(); ++i) {
491 0 : args[i] = all_args[i];
492 : }
493 : // and terms that should not
494 : std::vector<double> bargs;
495 0 : if(nores_>0) {
496 0 : bargs.resize(nores_);
497 : }
498 0 : for(unsigned i=0; i<bargs.size(); ++i) {
499 0 : bargs[i] = all_args[i+args.size()];
500 : }
501 :
502 : // calculate energy and forces, only on rescaled terms
503 : double ene = 0.0;
504 0 : for(unsigned i=0; i<args.size(); ++i) {
505 : // calculate energy term
506 0 : double fact = 1.0/pow(gamma_[igamma], expo_[i]) - 1.0;
507 0 : ene += args[i] * fact;
508 : // add force
509 0 : setOutputForce(i, -fact);
510 : }
511 :
512 : // set to zero on the others
513 0 : for(unsigned i=0; i<bargs.size(); ++i) {
514 0 : setOutputForce(i+args.size(), 0.0);
515 : }
516 :
517 : // set value of the bias
518 : setBias(ene);
519 : // set value of the wt-bias
520 0 : getPntrToComponent("wtbias")->set(bias_[igamma]);
521 : // set values of gamma
522 0 : getPntrToComponent("igamma")->set(igamma);
523 : // get time step
524 0 : long long int step = getStep();
525 0 : if(MCfirst_==-1) {
526 0 : MCfirst_=step;
527 : }
528 : // calculate gamma acceptance
529 0 : double MCtrials = std::floor(static_cast<double>(step-MCfirst_) / static_cast<double>(MCstride_))+1.0;
530 0 : double accgamma = static_cast<double>(MCaccgamma_) / static_cast<double>(MCsteps_) / MCtrials;
531 0 : getPntrToComponent("accgamma")->set(accgamma);
532 :
533 : // do MC at the right time step
534 0 : if(step%MCstride_==0&&!getExchangeStep()) {
535 0 : doMonteCarlo(igamma, ene, args, bargs);
536 : }
537 :
538 : // add well-tempered like bias
539 0 : if(step%Biaspace_==0) {
540 : // get updated igamma
541 0 : unsigned igamma = static_cast<unsigned>(plumed.passMap[selector_]);
542 : // add "Gaussian"
543 0 : double kbDT = kbt_ * ( biasf_ - 1.0 );
544 0 : bias_[igamma] += w0_ * std::exp(-bias_[igamma] / kbDT);
545 : }
546 :
547 : // print bias
548 0 : if(step%Biasstride_==0) {
549 0 : print_bias(step);
550 : }
551 :
552 0 : }
553 :
554 :
555 : }
556 : }
557 :
|
