/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Copyright (c) 2011-2015 The plumed team
(see the PEOPLE file at the root of the distribution for a list of names)
See http://www.plumed-code.org for more information.
This file is part of plumed, version 2.
plumed is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
plumed is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with plumed. If not, see .
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
#include "Bias.h"
#include "ActionRegister.h"
#include "core/ActionSet.h"
#include "tools/Grid.h"
#include "core/PlumedMain.h"
#include "core/Atoms.h"
#include "tools/Exception.h"
#include "core/FlexibleBin.h"
#include "tools/Matrix.h"
#include "tools/Random.h"
#include
#include
#include "tools/File.h"
#include
#include
#define DP2CUTOFF 6.25
using namespace std;
namespace PLMD{
namespace bias{
//+PLUMEDOC BIAS METAD
/*
Used to performed MetaDynamics on one or more collective variables.
In a metadynamics simulations a history dependent bias composed of
intermittently added Gaussian functions is added to the potential \cite metad.
\f[
V(\vec{s},t) = \sum_{ k \tau < t} W(k \tau)
\exp\left(
-\sum_{i=1}^{d} \frac{(s_i-s_i^{(0)}(k \tau))^2}{2\sigma_i^2}
\right).
\f]
This potential forces the system away from the kinetic traps in the potential energy surface
and out into the unexplored parts of the energy landscape. Information on the Gaussian
functions from which this potential is composed is output to a file called HILLS, which
is used both the restart the calculation and to reconstruct the free energy as a function of the CVs.
The free energy can be reconstructed from a metadynamics calculation because the final bias is given
by:
\f[
V(\vec{s}) = -F(\vec(s))
\f]
During post processing the free energy can be calculated in this way using the \ref sum_hills
utility.
In the simplest possible implementation of a metadynamics calculation the expense of a metadynamics
calculation increases with the length of the simulation as one has to, at every step, evaluate
the values of a larger and larger number of Gaussians. To avoid this issue you can
store the bias on a grid. This approach is similar to that proposed in \cite babi08jcp but has the
advantage that the grid spacing is independent on the Gaussian width.
Notice that you should
provide either the number of bins for every collective variable (GRID_BIN) or
the desired grid spacing (GRID_SPACING). In case you provide both PLUMED will use
the most conservative choice (highest number of bins) for each dimension.
In case you do not provide any information about bin size (neither GRID_BIN nor GRID_SPACING)
and if Gaussian width is fixed PLUMED will use 1/5 of the Gaussian width as grid spacing.
This default choice should be reasonable for most applications.
Another option that is available in plumed is well-tempered metadynamics \cite Barducci:2008. In this
varient of metadynamics the heights of the Gaussian hills are rescaled at each step so the bias is now
given by:
\f[
V({s},t)= \sum_{t'=0,\tau_G,2\tau_G,\dots}^{t' sw, the history dependent potential is still updated according to the above
equations but the metadynamics force is set to zero for s < sw. Notice that Gaussians are added also
if s < sw, as the tails of these Gaussians influence VG in the relevant region s > sw. In this way, the
force on the system in the region s > sw comes from both metadynamics and the force field, in the region
s < sw only from the latter. This approach allows obtaining a history-dependent bias potential VG that
fluctuates around a stable estimator, equal to the negative of the free energy far enough from the
boundaries. Note that:
- It works only for one-dimensional biases;
- It works both with and without GRID;
- The interval limit sw in a region where the free energy derivative is not large;
- If in the region outside the limit sw the system has a free energy minimum, the INTERVAL keyword should
be used together with a \ref UPPER_WALLS or \ref LOWER_WALLS at sw.
As a final note, since version 2.0.2 when the system is outside of the selected interval the force
is set to zero and the bias value to the value at the corresponding boundary. This allows acceptances
for replica exchange methods to be computed correctly.
Multiple walkers \cite multiplewalkers can also be used. See below the examples.
The c(t) reweighting factor can also be calculated on the fly using the equations
presented in \cite Tiwary_jp504920s.
The expression used to calculate c(t) follows directly from using Eq. 12 in
Eq. 3 in \cite Tiwary_jp504920s and gives smoother results than equivalent Eqs. 13
and Eqs. 14 in that paper. The c(t) is given by the rct component while the bias
normalized by c(t) is given by the rbias component (rbias=bias-ct) which can be used
to obtain a reweighted histogram.
The calculation of c(t) is enabled by using the keyword REWEIGHTING_NGRID where the grid used for the
calculation is specified. This grid should have a size that is equal or larger than the grid given in GRID_BIN./
By default c(t) is updated every 50 Gaussian hills but this
can be changed by using the REWEIGHTING_NHILLS keyword.
This option can only be employed together with Well-Tempered Metadynamics and requires that
a grid is used.
Additional material and examples can be also found in the tutorials:
- \ref belfast-6
- \ref belfast-7
- \ref belfast-8
Notice that at variance with PLUMED 1.3 it is now straightforward to apply concurrent metadynamics
as done e.g. in Ref. \cite gil2015enhanced . This indeed can be obtained by using the METAD
action multiple times in the same input file.
\par Examples
The following input is for a standard metadynamics calculation using as
collective variables the distance between atoms 3 and 5
and the distance between atoms 2 and 4. The value of the CVs and
the metadynamics bias potential are written to the COLVAR file every 100 steps.
\verbatim
DISTANCE ATOMS=3,5 LABEL=d1
DISTANCE ATOMS=2,4 LABEL=d2
METAD ARG=d1,d2 SIGMA=0.2,0.2 HEIGHT=0.3 PACE=500 LABEL=restraint
PRINT ARG=d1,d2,restraint.bias STRIDE=100 FILE=COLVAR
\endverbatim
(See also \ref DISTANCE \ref PRINT).
\par
If you use adaptive Gaussians, with diffusion scheme where you use
a Gaussian that should cover the space of 20 timesteps in collective variables.
Note that in this case the histogram correction is needed when summing up hills.
\verbatim
DISTANCE ATOMS=3,5 LABEL=d1
DISTANCE ATOMS=2,4 LABEL=d2
METAD ARG=d1,d2 SIGMA=20 HEIGHT=0.3 PACE=500 LABEL=restraint ADAPTIVE=DIFF
PRINT ARG=d1,d2,restraint.bias STRIDE=100 FILE=COLVAR
\endverbatim
\par
If you use adaptive Gaussians, with geometrical scheme where you use
a Gaussian that should cover the space of 0.05 nm in Cartesian space.
Note that in this case the histogram correction is needed when summing up hills.
\verbatim
DISTANCE ATOMS=3,5 LABEL=d1
DISTANCE ATOMS=2,4 LABEL=d2
METAD ARG=d1,d2 SIGMA=0.05 HEIGHT=0.3 PACE=500 LABEL=restraint ADAPTIVE=GEOM
PRINT ARG=d1,d2,restraint.bias STRIDE=100 FILE=COLVAR
\endverbatim
\par
When using adaptive Gaussians you might want to limit how the hills width can change.
You can use SIGMA_MIN and SIGMA_MAX keywords.
The sigmas should specified in terms of CV so you should use the CV units.
Note that if you use a negative number, this means that the limit is not set.
Note also that in this case the histogram correction is needed when summing up hills.
\verbatim
DISTANCE ATOMS=3,5 LABEL=d1
DISTANCE ATOMS=2,4 LABEL=d2
METAD ...
ARG=d1,d2 SIGMA=0.05 HEIGHT=0.3 PACE=500 LABEL=restraint ADAPTIVE=GEOM
SIGMA_MIN=0.2,0.1 SIGMA_MAX=0.5,1.0
... METAD
PRINT ARG=d1,d2,restraint.bias STRIDE=100 FILE=COLVAR
\endverbatim
\par
Multiple walkers can be also use as in \cite multiplewalkers
These are enabled by setting the number of walker used, the id of the
current walker which interprets the input file, the directory where the
hills containing files resides, and the frequency to read the other walkers.
Here is an example
\verbatim
DISTANCE ATOMS=3,5 LABEL=d1
METAD ...
ARG=d1 SIGMA=0.05 HEIGHT=0.3 PACE=500 LABEL=restraint
WALKERS_N=10
WALKERS_ID=3
WALKERS_DIR=../
WALKERS_RSTRIDE=100
... METAD
\endverbatim
where WALKERS_N is the total number of walkers, WALKERS_ID is the
id of the present walker (starting from 0 ) and the WALKERS_DIR is the directory
where all the walkers are located. WALKERS_RSTRIDE is the number of step between
one update and the other.
\par
The c(t) reweighting factor can be calculated on the fly using the equations
presented in \cite Tiwary_jp504920s as described above.
This is enabled by using the keyword REWEIGHTING_NGRID where the grid used for
the calculation is set. The number of grid points given in REWEIGHTING_NGRID
should be equal or larger than the number of grid points given in GRID_BIN.
\verbatim
METAD ...
LABEL=metad
ARG=phi,psi SIGMA=0.20,0.20 HEIGHT=1.20 BIASFACTOR=5 TEMP=300.0 PACE=500
GRID_MIN=-pi,-pi GRID_MAX=pi,pi GRID_BIN=150,150
REWEIGHTING_NGRID=150,150
REWEIGHTING_NHILLS=20
... METAD
\endverbatim
Here we have asked that the calculation is performed every 20 hills by using
REWEIGHTING_NHILLS keyword. If this keyword is not given the calculation will
by default be performed every 50 hills. The c(t) reweighting factor will be given
in the rct component while the instantaneous value of the bias potential
normalized using the c(t) reweighting factor is given in the rbias component
[rbias=bias-c(t)] which can be used to obtain a reweighted histogram or
free energy surface using the \ref HISTOGRAM analysis.
\par
The kinetics of the transitions between basins can also be analysed on the fly as
in \cite PRL230602. The flag ACCELERATION turn on accumulation of the acceleration
factor that can then be used to determine the rate. This method can be used together
with \ref COMMITTOR analysis to stop the simulation when the system get to the target basin.
It must be used together with Well-Tempered Metadynamics.
*/
//+ENDPLUMEDOC
class MetaD : public Bias{
private:
struct Gaussian {
vector center;
vector sigma;
double height;
bool multivariate; // this is required to discriminate the one dimensional case
vector invsigma;
Gaussian(const vector & center,const vector & sigma,double height, bool multivariate ):
center(center),sigma(sigma),height(height),multivariate(multivariate),invsigma(sigma){
// to avoid troubles from zero element in flexible hills
for(unsigned i=0;i1.e-20?invsigma[i]=1.0/invsigma[i]:0.;
}
};
vector sigma0_;
vector sigma0min_;
vector sigma0max_;
vector hills_;
OFile hillsOfile_;
OFile gridfile_;
Grid* BiasGrid_;
Grid* ExtGrid_;
bool storeOldGrids_;
std::string gridfilename_,gridreadfilename_;
int wgridstride_;
bool grid_,hasextgrid_;
double height0_;
double biasf_;
double kbt_;
int stride_;
bool welltemp_;
double* dp_;
int adaptive_;
FlexibleBin *flexbin;
int mw_n_;
string mw_dir_;
int mw_id_;
int mw_rstride_;
bool walkers_mpi;
bool acceleration;
double acc;
vector ifiles;
vector ifilesnames;
double uppI_;
double lowI_;
bool doInt_;
bool isFirstStep;
double reweight_factor;
std::vector rewf_grid_;
unsigned int rewf_ustride_;
/// accumulator for work
double work_;
long int last_step_warn_grid;
int altruistic_; // altruistic
double altfactor_; // altruistic
double alt_minor_koef_; // altruistic_mpi
void readGaussians(IFile*);
bool readChunkOfGaussians(IFile *ifile, unsigned n);
void writeGaussian(const Gaussian&,OFile&);
void addGaussian(const Gaussian&);
double getHeight(const vector&);
double getBiasAndDerivatives(const vector&,double* der=NULL);
double evaluateGaussian(const vector&, const Gaussian&,double* der=NULL);
void finiteDifferenceGaussian(const vector&, const Gaussian&);
double getGaussianNormalization( const Gaussian& );
vector getGaussianSupport(const Gaussian&);
bool scanOneHill(IFile *ifile, vector &v, vector ¢er, vector &sigma, double &height, bool &multivariate );
void computeReweightingFactor();
std::string fmt;
public:
explicit MetaD(const ActionOptions&);
~MetaD();
void calculate();
void update();
static void registerKeywords(Keywords& keys);
bool checkNeedsGradients()const{if(adaptive_==FlexibleBin::geometry){return true;}else{return false;}}
};
PLUMED_REGISTER_ACTION(MetaD,"METAD")
void MetaD::registerKeywords(Keywords& keys){
Bias::registerKeywords(keys);
componentsAreNotOptional(keys);
keys.addOutputComponent("bias","default","the instantaneous value of the bias potential");
keys.addOutputComponent("rbias","REWEIGHTING_NGRID","the instantaneous value of the bias normalized using the \\f$c(t)\\f$ reweighting factor [rbias=bias-c(t)]. This component can be used to obtain a reweighted histogram.");
keys.addOutputComponent("rct","REWEIGHTING_NGRID","the reweighting factor \\f$c(t)\\f$.");
keys.addOutputComponent("work","default","accumulator for work");
keys.addOutputComponent("acc","ACCELERATION","the metadynamics acceleration factor");
keys.use("ARG");
keys.add("compulsory","SIGMA","the widths of the Gaussian hills");
keys.add("compulsory","PACE","the frequency for hill addition");
keys.add("compulsory","FILE","HILLS","a file in which the list of added hills is stored");
keys.add("optional","HEIGHT","the heights of the Gaussian hills. Compulsory unless TAU, TEMP and BIASFACTOR are given");
keys.add("optional","FMT","specify format for HILLS files (useful for decrease the number of digits in regtests)");
keys.add("optional","BIASFACTOR","use well tempered metadynamics and use this biasfactor. Please note you must also specify temp");
keys.add("optional","TEMP","the system temperature - this is only needed if you are doing well-tempered metadynamics");
keys.add("optional","TAU","in well tempered metadynamics, sets height to (kb*DeltaT*pace*timestep)/tau");
keys.add("optional","GRID_MIN","the lower bounds for the grid");
keys.add("optional","GRID_MAX","the upper bounds for the grid");
keys.add("optional","GRID_BIN","the number of bins for the grid");
keys.add("optional","GRID_SPACING","the approximate grid spacing (to be used as an alternative or together with GRID_BIN)");
keys.add("optional","REWEIGHTING_NGRID","calculate the c(t) reweighting factor and use that to obtain the normalized bias [rbias=bias-c(t)]. Here you should specify the number of grid points required in each dimension. The number of grid points should be equal or larger to the number of grid points given in GRID_BIN." "This method is not compatible with metadynamics not on a grid.");
keys.add("optional","REWEIGHTING_NHILLS","how many Gaussian hills should be deposited between calculating the c(t) reweighting factor. The default is to do this every 50 hills.");
keys.addFlag("GRID_SPARSE",false,"use a sparse grid to store hills");
keys.addFlag("GRID_NOSPLINE",false,"don't use spline interpolation with grids");
keys.add("optional","GRID_WSTRIDE","write the grid to a file every N steps");
keys.add("optional","GRID_WFILE","the file on which to write the grid");
keys.addFlag("STORE_GRIDS",false,"store all the grid files the calculation generates. They will be deleted if this keyword is not present");
keys.add("optional","ADAPTIVE","use a geometric (=GEOM) or diffusion (=DIFF) based hills width scheme. Sigma is one number that has distance units or timestep dimensions");
keys.add("optional","WALKERS_ID", "walker id");
keys.add("optional","WALKERS_N", "number of walkers");
keys.add("optional","WALKERS_DIR", "shared directory with the hills files from all the walkers");
keys.add("optional","WALKERS_RSTRIDE","stride for reading hills files");
keys.add("optional","INTERVAL","monodimensional lower and upper limits, outside the limits the system will not feel the biasing force.");
keys.add("optional","GRID_RFILE","a grid file from which the bias should be read at the initial step of the simulation");
keys.add("optional","SIGMA_MAX","the upper bounds for the sigmas (in CV units) when using adaptive hills. Negative number means no bounds ");
keys.add("optional","SIGMA_MIN","the lower bounds for the sigmas (in CV units) when using adaptive hills. Negative number means no bounds ");
keys.add("optional","ALTFACTOR","altfactor for altruistic MTD"); // altruistic
keys.add("optional","ALTTYPE","version of height modificators for altruistic MTD"); // altruistic
keys.addFlag("WALKERS_MPI",false,"Switch on MPI version of multiple walkers - not compatible with other WALKERS_* options");
keys.addFlag("ACCELERATION",false,"Set to TRUE if you want to compute the metadynamics acceleration factor.");
keys.use("RESTART");
keys.use("UPDATE_FROM");
keys.use("UPDATE_UNTIL");
}
MetaD::~MetaD(){
if(flexbin) delete flexbin;
if(BiasGrid_) delete BiasGrid_;
hillsOfile_.close();
if(wgridstride_>0) gridfile_.close();
delete [] dp_;
// close files
for(int i=0;iisOpen()) ifiles[i]->close();
delete ifiles[i];
}
}
MetaD::MetaD(const ActionOptions& ao):
PLUMED_BIAS_INIT(ao),
// Grid stuff initialization
BiasGrid_(NULL),ExtGrid_(NULL), wgridstride_(0), grid_(false), hasextgrid_(false),
// Metadynamics basic parameters
height0_(std::numeric_limits::max()), biasf_(1.0), kbt_(0.0),
stride_(0), welltemp_(false),
altfactor_(-1.0), altruistic_(0), alt_minor_koef_(-1.0), // altruistic
// Other stuff
dp_(NULL), adaptive_(FlexibleBin::none),
flexbin(NULL),
// Multiple walkers initialization
mw_n_(1), mw_dir_("./"), mw_id_(0), mw_rstride_(1),
walkers_mpi(false),
acceleration(false), acc(0.0),
// Interval initialization
uppI_(-1), lowI_(-1), doInt_(false),
isFirstStep(true),
reweight_factor(0.0),
rewf_ustride_(1),
work_(0),
last_step_warn_grid(0)
{
// parse the flexible hills
string adaptiveoption;
adaptiveoption="NONE";
parse("ADAPTIVE",adaptiveoption);
if (adaptiveoption=="GEOM"){
log.printf(" Uses Geometry-based hills width: sigma must be in distance units and only one sigma is needed\n");
adaptive_=FlexibleBin::geometry;
}else if (adaptiveoption=="DIFF"){
log.printf(" Uses Diffusion-based hills width: sigma must be in timesteps and only one sigma is needed\n");
adaptive_=FlexibleBin::diffusion;
}else if (adaptiveoption=="NONE"){
adaptive_=FlexibleBin::none;
}else{
error("I do not know this type of adaptive scheme");
}
// parse the sigma
parseVector("SIGMA",sigma0_);
parse("FMT",fmt);
if (adaptive_==FlexibleBin::none){
// if you use normal sigma you need one sigma per argument
if( sigma0_.size()!=getNumberOfArguments() ) error("number of arguments does not match number of SIGMA parameters");
}else{
// if you use flexible hills you need one sigma
if(sigma0_.size()!=1){
error("If you choose ADAPTIVE you need only one sigma according to your choice of type (GEOM/DIFF)");
}
// if adaptive then the number must be an integer
if(adaptive_==FlexibleBin::diffusion){
if(int(sigma0_[0])-sigma0_[0]>1.e-9 || int(sigma0_[0])-sigma0_[0] <-1.e-9 || int(sigma0_[0])<1 ){
error("In case of adaptive hills with diffusion, the sigma must be an integer which is the number of timesteps\n");
}
}
// here evtl parse the sigma min and max values
parseVector("SIGMA_MIN",sigma0min_);
if(sigma0min_.size()>0 && sigma0min_.size()0 && sigma0max_.size()0.0) kbt_=plumed.getAtoms().getKBoltzmann()*temp;
else kbt_=plumed.getAtoms().getKbT();
if(biasf_>1.0){
if(kbt_==0.0) error("Unless the MD engine passes the temperature to plumed, with well-tempered metad you must specify it using TEMP");
welltemp_=true;
}
double tau=0.0;
parse("TAU",tau);
if(tau==0.0){
if(height0_==std::numeric_limits::max()) error("At least one between HEIGHT and TAU should be specified");
// if tau is not set, we compute it here from the other input parameters
if(welltemp_) tau=(kbt_*(biasf_-1.0))/height0_*getTimeStep()*stride_;
} else {
if(!welltemp_)error("TAU only makes sense in well-tempered metadynamics");
if(height0_!=std::numeric_limits::max()) error("At most one between HEIGHT and TAU should be specified");
height0_=(kbt_*(biasf_-1.0))/tau*getTimeStep()*stride_;
}
// Grid Stuff
vector gmin(getNumberOfArguments());
parseVector("GRID_MIN",gmin);
if(gmin.size()!=getNumberOfArguments() && gmin.size()!=0) error("not enough values for GRID_MIN");
vector gmax(getNumberOfArguments());
parseVector("GRID_MAX",gmax);
if(gmax.size()!=getNumberOfArguments() && gmax.size()!=0) error("not enough values for GRID_MAX");
vector gbin(getNumberOfArguments());
vector gspacing;
parseVector("GRID_BIN",gbin);
if(gbin.size()!=getNumberOfArguments() && gbin.size()!=0) error("not enough values for GRID_BIN");
parseVector("GRID_SPACING",gspacing);
if(gspacing.size()!=getNumberOfArguments() && gspacing.size()!=0) error("not enough values for GRID_SPACING");
if(gmin.size()!=gmax.size()) error("GRID_MAX and GRID_MIN should be either present or absent");
if(gspacing.size()!=0 && gmin.size()==0) error("If GRID_SPACING is present also GRID_MIN should be present");
if(gbin.size()!=0 && gmin.size()==0) error("If GRID_SPACING is present also GRID_MIN should be present");
if(gmin.size()!=0){
if(gbin.size()==0 && gspacing.size()==0){
if(adaptive_==FlexibleBin::none){
log<<" Binsize not specified, 1/5 of sigma will be be used\n";
plumed_assert(sigma0_.size()==getNumberOfArguments());
gspacing.resize(getNumberOfArguments());
for(unsigned i=0;i0){grid_=true;}
parse("GRID_WSTRIDE",wgridstride_);
parse("GRID_WFILE",gridfilename_);
parseFlag("STORE_GRIDS",storeOldGrids_);
if(grid_ && gridfilename_.length()>0){
if(wgridstride_==0 ) error("frequency with which to output grid not specified use GRID_WSTRIDE");
}
if(grid_ && wgridstride_>0){
if(gridfilename_.length()==0) error("grid filename not specified use GRID_WFILE");
}
parse("GRID_RFILE",gridreadfilename_);
if(grid_){
parseVector("REWEIGHTING_NGRID",rewf_grid_);
if( rewf_grid_.size()>0 && rewf_grid_.size()!=getNumberOfArguments() ){
error("size mismatch for REWEIGHTING_NGRID keyword");
} else if( rewf_grid_.size()==getNumberOfArguments() ){
for(unsigned j=0;jisPeriodic() ) rewf_grid_[j] += 1;
}
}
if( adaptive_==FlexibleBin::diffusion || adaptive_==FlexibleBin::geometry) warning("reweighting has not been proven to work with adaptive Gaussians");
rewf_ustride_=50; parse("REWEIGHTING_NHILLS",rewf_ustride_);
}
// Multiple walkers
parse("WALKERS_N",mw_n_);
parse("WALKERS_ID",mw_id_);
if(mw_n_<=mw_id_) error("walker ID should be a numerical value less than the total number of walkers");
parse("WALKERS_DIR",mw_dir_);
parse("WALKERS_RSTRIDE",mw_rstride_);
// MPI version
parseFlag("WALKERS_MPI",walkers_mpi);
parse("ALTFACTOR",altfactor_); // altruistic
parse("ALTTYPE",altruistic_); // altruistic
if(altfactor_>=0.0) { // altruistic
if(!altruistic_) altruistic_ = 3; // altruistic
if(altfactor_>1.0) error("altfactor must be 0-1"); // altruistic
int nw = walkers_mpi ? multi_sim_comm.Get_size() : mw_n_;
if(nw < 2) error("altruistic MTD needs walkers"); // altruistic
switch(altruistic_){ // altruistic
case 1: height0_ /= 2; // altruistic
case 4: height0_ *= (2-altfactor_); // altruistic
case 3: alt_minor_koef_ = altfactor_/(2.0-altfactor_)/(nw-1.0);break;
case 2: alt_minor_koef_ = altfactor_/(nw-1.0);break; // altruistic
default: error("Unknown ALTTYPE - specify 1-4"); // altruistic
} // altruistic
} // altruistic
if(altruistic_
&& welltemp_
&& !walkers_mpi) error("file-based altruistic MTD with Welltempered is not implemented"); // altruistic
// Inteval keyword
vector tmpI(2);
parseVector("INTERVAL",tmpI);
if(tmpI.size()!=2&&tmpI.size()!=0) error("both a lower and an upper limits must be provided with INTERVAL");
else if(tmpI.size()==2) {
lowI_=tmpI.at(0);
uppI_=tmpI.at(1);
if(getNumberOfArguments()!=1) error("INTERVAL limits correction works only for monodimensional metadynamics!");
if(uppI_0){log.printf(" Grid is written on file %s with stride %d\n",gridfilename_.c_str(),wgridstride_);}
}
if(gridreadfilename_.length()>0){
log.printf(" Reading an additional bias from grid in file %s \n",gridreadfilename_.c_str());
}
if(mw_n_>1){
if(walkers_mpi) error("MPI version of multiple walkers is not compatible with filesystem version of multiple walkers");
log.printf(" %d multiple walkers active\n",mw_n_);
log.printf(" walker id %d\n",mw_id_);
log.printf(" reading stride %d\n",mw_rstride_);
log.printf(" directory with hills files %s\n",mw_dir_.c_str());
} else {
if(walkers_mpi) log.printf(" Multiple walkers active using MPI communnication\n");
}
if(altruistic_){ //altruistic
log.printf(" Running altruistic MTD with %d walkers and altfactor %f (major=%f minor=%f)\n",
walkers_mpi?multi_sim_comm.Get_size():mw_n_, altfactor_, height0_, alt_minor_koef_*height0_); //altruistic
}else{ log.printf(" ALTRUISTIC not engaged");} //altruistic
addComponent("bias"); componentIsNotPeriodic("bias");
if( rewf_grid_.size()>0 ){
addComponent("rbias"); componentIsNotPeriodic("rbias");
addComponent("rct"); componentIsNotPeriodic("rct");
}
addComponent("work"); componentIsNotPeriodic("work");
if(acceleration) {
if(!welltemp_) error("The calculation of the acceleration works only if Well-Tempered Metadynamics is on");
log.printf(" calculation on the fly of the acceleration factor");
addComponent("acc"); componentIsNotPeriodic("acc");
}
// for performance
dp_ = new double[getNumberOfArguments()];
// initializing and checking grid
if(grid_){
// check for adaptive and sigma_min
if(sigma0min_.size()==0&&adaptive_!=FlexibleBin::none) error("When using Adaptive Gaussians on a grid SIGMA_MIN must be specified");
// check for mesh and sigma size
for(unsigned i=0;i0.5*sigma0_[i]) log<<" WARNING: Using a METAD with a Grid Spacing larger than half of the Gaussians width can produce artifacts\n";
}
std::string funcl=getLabel() + ".bias";
if(!sparsegrid){BiasGrid_=new Grid(funcl,getArguments(),gmin,gmax,gbin,spline,true);}
else{BiasGrid_=new SparseGrid(funcl,getArguments(),gmin,gmax,gbin,spline,true);}
std::vector actualmin=BiasGrid_->getMin();
std::vector actualmax=BiasGrid_->getMax();
for(unsigned i=0;i0){
gridfile_.link(*this);
gridfile_.open(gridfilename_);
}
// initializing external grid
if(gridreadfilename_.length()>0){
hasextgrid_=true;
// read the grid in input, find the keys
IFile gridfile; gridfile.open(gridreadfilename_);
std::string funcl=getLabel() + ".bias";
ExtGrid_=Grid::create(funcl,getArguments(),gridfile,false,false,true);
gridfile.close();
if(ExtGrid_->getDimension()!=getNumberOfArguments()) error("mismatch between dimensionality of input grid and number of arguments");
for(unsigned i=0;iisPeriodic()!=ExtGrid_->getIsPeriodic()[i] ) error("periodicity mismatch between arguments and input bias");
}
}
// Tiwary-Parrinello reweighting factor
if(rewf_grid_.size()>0){
log.printf(" the c(t) reweighting factor will be calculated every %u hills\n",rewf_ustride_);
getPntrToComponent("rct")->set(reweight_factor);
}
// creating vector of ifile* for hills reading
// open all files at the beginning and read Gaussians if restarting
for(int i=0;i1) {
stringstream out; out << i;
fname = mw_dir_+"/"+hillsfname+"."+out.str();
} else {
fname = hillsfname;
}
IFile *ifile = new IFile();
ifile->link(*this);
ifiles.push_back(ifile);
ifilesnames.push_back(fname);
if(ifile->FileExist(fname)){
ifile->open(fname);
if(getRestart()){
log.printf(" Restarting from %s:",ifilesnames[i].c_str());
readGaussians(ifiles[i]);
}
ifiles[i]->reset(false);
// close only the walker own hills file for later writing
if(i==mw_id_) ifiles[i]->close();
}
}
comm.Barrier();
// this barrier is needed when using walkers_mpi
// to be sure that all files have been read before
// backing them up
// it should not be used when walkers_mpi is false otherwise
// it would introduce troubles when using replicas without METAD
// (e.g. in bias exchange with a neutral replica)
// see issue #168 on github
if(comm.Get_rank()==0 && walkers_mpi) multi_sim_comm.Barrier();
// Calculate the Tiwary-Parrinello reweighting factor if we are restarting from previous hills
if(plumed.getRestart() && rewf_grid_.size()>0 ){computeReweightingFactor();}
// open hills file for writing
hillsOfile_.link(*this);
//if(walkers_mpi){ // altruistic_mpi
//int r=0; // altruistic_mpi
//if(comm.Get_rank()==0) r=multi_sim_comm.Get_rank(); // altruistic_mpi
//comm.Bcast(r,0); // altruistic_mpi
//if(r>0) ifilesnames[mw_id_]="/dev/null"; // altruistic_mpi
//hillsOfile_.enforceSuffix(""); // altruistic_mpi
//}
hillsOfile_.open(ifilesnames[mw_id_]);
if(fmt.length()>0) hillsOfile_.fmtField(fmt);
hillsOfile_.addConstantField("multivariate");
if(doInt_) {
hillsOfile_.addConstantField("lower_int").printField("lower_int",lowI_);
hillsOfile_.addConstantField("upper_int").printField("upper_int",uppI_);
}
hillsOfile_.setHeavyFlush();
// output periodicities of variables
for(unsigned i=0;i(*p)){ concurrent=true; break; }
if(concurrent) log<<" You are using concurrent metadynamics\n";
log<<" Bibliography "<1||walkers_mpi) log<0) log< center(ncv);
vector sigma(ncv);
double height;
int nhills=0;
bool multivariate=false;
std::vector tmpvalues;
for(unsigned j=0;jgetName(), false ) );
while(scanOneHill(ifile,tmpvalues,center,sigma,height,multivariate)){;
nhills++;
if(welltemp_){height*=(biasf_-1.0)/biasf_;}
addGaussian(Gaussian(center,sigma,height,multivariate));
}
log.printf(" %d Gaussians read\n",nhills);
}
bool MetaD::readChunkOfGaussians(IFile *ifile, unsigned n)
{
unsigned ncv=getNumberOfArguments();
vector center(ncv);
vector sigma(ncv);
double height;
unsigned nhills=0;
bool multivariate=false;
std::vector tmpvalues;
for(unsigned j=0;jgetName(), false ) );
while(scanOneHill(ifile,tmpvalues,center,sigma,height,multivariate)){;
if(welltemp_){height*=(biasf_-1.0)/biasf_;}
addGaussian(Gaussian(center,sigma,height,multivariate));
if(nhills==n){
log.printf(" %u Gaussians read\n",nhills);
return true;
}
nhills++;
}
log.printf(" %u Gaussians read\n",nhills);
return false;
}
void MetaD::writeGaussian(const Gaussian& hill, OFile&file){
unsigned ncv=getNumberOfArguments();
file.printField("time",getTimeStep()*getStep());
for(unsigned i=0;i mymatrix(ncv,ncv);
unsigned k=0;
for(unsigned i=0;i invmatrix(ncv,ncv);
Invert(mymatrix,invmatrix);
// enforce symmetry
for(unsigned i=0;i lower(ncv,ncv);
cholesky(invmatrix,lower); // now this , in band form , is similar to the sigmas
// loop in band form
for (unsigned i=0;igetName()+"_"+getPntrToArgument(j)->getName(),lower(j+i,j));
}
}
} else {
hillsOfile_.printField("multivariate","false");
for(unsigned i=0;igetName(),hill.sigma[i]);
}
double height=hill.height;
if(altruistic_ && !walkers_mpi){ // altruistic_mpi
height = alt_minor_koef_*hill.height; // altruistic_mpi
} // altruistic_mpi
if(welltemp_){height*=biasf_/(biasf_-1.0);}
file.printField("height",height).printField("biasf",biasf_);
if(mw_n_>1) file.printField("clock",int(time(0)));
file.printField();
}
void MetaD::addGaussian(const Gaussian& hill){
if(!grid_){hills_.push_back(hill);}
else{
unsigned ncv=getNumberOfArguments();
vector nneighb=getGaussianSupport(hill);
vector neighbors=BiasGrid_->getNeighbors(hill.center,nneighb);
vector der(ncv);
vector xx(ncv);
if(comm.Get_size()==1){
for(unsigned i=0;igetPoint(ineigh,xx);
double bias=evaluateGaussian(xx,hill,&der[0]);
BiasGrid_->addValueAndDerivatives(ineigh,bias,der);
}
} else {
unsigned stride=comm.Get_size();
unsigned rank=comm.Get_rank();
vector allder(ncv*neighbors.size(),0.0);
vector allbias(neighbors.size(),0.0);
for(unsigned i=rank;igetPoint(ineigh,xx);
allbias[i]=evaluateGaussian(xx,hill,&allder[ncv*i]);
}
comm.Sum(allbias);
comm.Sum(allder);
for(unsigned i=0;iaddValueAndDerivatives(ineigh,allbias[i],der);
}
}
}
}
vector MetaD::getGaussianSupport(const Gaussian& hill)
{
vector nneigh;
if(doInt_){
double cutoff=sqrt(2.0*DP2CUTOFF)*hill.sigma[0];
if(hill.center[0]+cutoff > uppI_ || hill.center[0]-cutoff < lowI_) {
// in this case, we updated the entire grid to avoid problems
return BiasGrid_->getNbin();
} else {
nneigh.push_back( static_cast(ceil(cutoff/BiasGrid_->getDx()[0])) );
return nneigh;
}
}
// traditional or flexible hill?
if(hill.multivariate){
unsigned ncv=getNumberOfArguments();
unsigned k=0;
//log<<"------- GET GAUSSIAN SUPPORT --------\n";
Matrix mymatrix(ncv,ncv);
for(unsigned i=0;i myinv(ncv,ncv);
Invert(mymatrix,myinv);
Matrix myautovec(ncv,ncv);
vector myautoval(ncv); //should I take this or their square root?
diagMat(myinv,myautoval,myautovec);
double maxautoval;maxautoval=0.;
unsigned ind_maxautoval;ind_maxautoval=ncv;
for (unsigned i=0;imaxautoval){maxautoval=myautoval[i];ind_maxautoval=i;}
}
for (unsigned i=0;i(ceil(cutoff/BiasGrid_->getDx()[i])) );
}
}else{
for(unsigned i=0;i(ceil(cutoff/BiasGrid_->getDx()[i])) );
}
}
//log<<"------- END GET GAUSSIAN SUPPORT --------\n";
return nneigh;
}
double MetaD::getBiasAndDerivatives(const vector& cv, double* der)
{
double bias=0.0;
if(!grid_){
if(hills_.size()>10000 && (getStep()-last_step_warn_grid)>10000){
std::string msg;
Tools::convert(hills_.size(),msg);
msg="You have accumulated "+msg+" hills, you should enable GRIDs to avoid serious performance hits";
warning(msg);
last_step_warn_grid=getStep();
}
unsigned stride=comm.Get_size();
unsigned rank=comm.Get_rank();
for(unsigned i=rank;i vder(getNumberOfArguments());
bias=BiasGrid_->getValueAndDerivatives(cv,vder);
for(unsigned i=0;igetValue(cv);
}
}
if(hasextgrid_){
if(der){
vector vder(getNumberOfArguments());
bias+=ExtGrid_->getValueAndDerivatives(cv,vder);
for(unsigned i=0;igetValue(cv);
}
}
return bias;
}
double MetaD::getGaussianNormalization( const Gaussian& hill ){
double norm=1; unsigned ncv=hill.center.size();
if(hill.multivariate){
// recompose the full sigma from the upper diag cholesky
unsigned k=0;
Matrix mymatrix(ncv,ncv);
for(unsigned i=0;i(ncv)/2.0);
}
double MetaD::evaluateGaussian
(const vector& cv, const Gaussian& hill, double* der)
{
double dp2=0.0;
double bias=0.0;
// I use a pointer here because cv is const (and should be const)
// but when using doInt it is easier to locally replace cv[0] with
// the upper/lower limit in case it is out of range
const double *pcv=NULL; // pointer to cv
double tmpcv[1]; // tmp array with cv (to be used with doInt_)
if(cv.size()>0) pcv=&cv[0];
if(doInt_){
plumed_dbg_assert(cv.size()==1);
pcv=&(tmpcv[0]);
tmpcv[0]=cv[0];
if(cv[0]uppI_) tmpcv[0]=uppI_;
}
if(hill.multivariate){
unsigned k=0;
unsigned ncv=cv.size();
// recompose the full sigma from the upper diag cholesky
Matrix mymatrix(ncv,ncv);
for(unsigned i=0;iuppI_) && der ) for(unsigned i=0;i& cv)
{
double height=height0_;
if(welltemp_){
double vbias=getBiasAndDerivatives(cv);
if(altruistic_ && walkers_mpi){ // altruistic mpiWT
if(comm.Get_rank()==0){ // altruistic mpiWT
multi_sim_comm.Sum(vbias); // altruistic mpiWT
} // altruistic mpiWT
//vbias/=multi_sim_comm.Get_size(); // altruistic mpiWT
comm.Bcast(vbias,0); // altruistic mpiWT
} // altruistic mpiWT
height=height0_*exp(-vbias/(kbt_*(biasf_-1.0)));
}
return height;
}
void MetaD::calculate()
{
// this is because presently there is no way to properly pass information
// on adaptive hills (diff) after exchanges:
if(adaptive_==FlexibleBin::diffusion && getExchangeStep()) error("ADAPTIVE=DIFF is not compatible with replica exchange");
unsigned ncv=getNumberOfArguments();
vector cv(ncv);
for(unsigned i=0;iset(ene);
if( rewf_grid_.size()>0 ) getPntrToComponent("rbias")->set(ene - reweight_factor);
// calculate the acceleration factor
if(acceleration&&!isFirstStep) {
acc += exp(ene/(kbt_));
double mean_acc = acc/((double) getStep());
getPntrToComponent("acc")->set(mean_acc);
}
getPntrToComponent("work")->set(work_);
// set Forces
for(unsigned i=0;i cv(getNumberOfArguments());
vector thissigma;
bool multivariate;
// adding hills criteria (could be more complex though)
bool nowAddAHill;
if(getStep()%stride_==0 && !isFirstStep ){nowAddAHill=true;}else{nowAddAHill=false;isFirstStep=false;}
for(unsigned i=0;iupdate(nowAddAHill);
multivariate=true;
}else{
multivariate=false;
};
if(nowAddAHill){ // probably this can be substituted with a signal
// add a Gaussian
double height=getHeight(cv);
// use normal sigma or matrix form?
if (adaptive_!=FlexibleBin::none){
thissigma=flexbin->getInverseMatrix(); // returns upper diagonal inverse
//cerr<<"ADDING HILLS "< all_cv(nw*cv.size(),0.0);
std::vector all_sigma(nw*thissigma.size(),0.0);
std::vector all_height(nw,0.0);
std::vector all_multivariate(nw,0);
if(comm.Get_rank()==0){
// Communicate (only root)
multi_sim_comm.Allgather(cv,all_cv);
multi_sim_comm.Allgather(thissigma,all_sigma);
if(altruistic_){ // altruistic_mpi
multi_sim_comm.Allgather(alt_minor_koef_*height,all_height); // altruistic_mpi
}else{ // altruistic_mpi
multi_sim_comm.Allgather(height,all_height);
} // altruistic_mpi
multi_sim_comm.Allgather(int(multivariate),all_multivariate);
}
// Share info with group members
comm.Bcast(all_cv,0);
comm.Bcast(all_sigma,0);
comm.Bcast(all_height,0);
comm.Bcast(all_multivariate,0);
for(int i=0;i cv_now(cv.size());
std::vector sigma_now(thissigma.size());
for(unsigned j=0;j0&&getStep()%wgridstride_==0){
// in case old grids are stored, a sequence of grids should appear
// this call results in a repetition of the header:
if(storeOldGrids_) gridfile_.clearFields();
// in case only latest grid is stored, file should be rewound
// this will overwrite previously written grids
else gridfile_.rewind();
BiasGrid_->writeToFile(gridfile_);
// if a single grid is stored, it is necessary to flush it, otherwise
// the file might stay empty forever (when a single grid is not large enough to
// trigger flushing from the operating system).
// on the other hand, if grids are stored one after the other this is
// no necessary, and we leave the flushing control to the user as usual
// (with FLUSH keyword)
if(!storeOldGrids_) gridfile_.flush();
}
// if multiple walkers and time to read Gaussians
if(mw_n_>1 && getStep()%mw_rstride_==0){
for(int i=0;iisOpen())){
// check if it exists now and open it!
if(ifiles[i]->FileExist(ifilesnames[i])) {
ifiles[i]->open(ifilesnames[i]);
ifiles[i]->reset(false);
}
// otherwise read the new Gaussians
} else {
log.printf(" Reading hills from %s:",ifilesnames[i].c_str());
readGaussians(ifiles[i]);
ifiles[i]->reset(false);
}
}
}
if(getStep()%(stride_*rewf_ustride_)==0 && nowAddAHill && rewf_grid_.size()>0 ){computeReweightingFactor();}
}
void MetaD::finiteDifferenceGaussian
(const vector& cv, const Gaussian& hill)
{
log<<"--------- finiteDifferenceGaussian: size "< oldder(cv.size());
vector der(cv.size());
vector mycv(cv.size());
mycv=cv;
double step=1.e-6;
Random random;
// just displace a tiny bit
for(unsigned i=0;i &tmpvalues, vector ¢er, vector &sigma, double &height , bool &multivariate ){
double dummy;
multivariate=false;
if(ifile->scanField("time",dummy)){
unsigned ncv; ncv=tmpvalues.size();
for(unsigned i=0;iscanField( &tmpvalues[i] );
if( tmpvalues[i].isPeriodic() && ! getPntrToArgument(i)->isPeriodic() ){
error("in hills file periodicity for variable " + tmpvalues[i].getName() + " does not match periodicity in input");
} else if( tmpvalues[i].isPeriodic() ){
std::string imin, imax; tmpvalues[i].getDomain( imin, imax );
std::string rmin, rmax; getPntrToArgument(i)->getDomain( rmin, rmax );
if( imin!=rmin || imax!=rmax ){
error("in hills file periodicity for variable " + tmpvalues[i].getName() + " does not match periodicity in input");
}
}
center[i]=tmpvalues[i].get();
}
// scan for multivariate label: record the actual file position so to eventually rewind
std::string sss;
ifile->scanField("multivariate",sss);
if(sss=="true") multivariate=true;
else if(sss=="false") multivariate=false;
else plumed_merror("cannot parse multivariate = "+ sss);
if(multivariate){
sigma.resize(ncv*(ncv+1)/2);
Matrix upper(ncv,ncv);
Matrix lower(ncv,ncv);
for (unsigned i=0;iscanField("sigma_"+getPntrToArgument(j+i)->getName()+"_"+getPntrToArgument(j)->getName(),lower(j+i,j));
upper(j,j+i)=lower(j+i,j);
}
}
Matrix mymult(ncv,ncv);
Matrix invmatrix(ncv,ncv);
//log<<"Lower \n";
//matrixOut(log,lower);
//log<<"Upper \n";
//matrixOut(log,upper);
mult(lower,upper,mymult);
//log<<"Mult \n";
//matrixOut(log,mymult);
// now invert and get the sigmas
Invert(mymult,invmatrix);
//log<<"Invert \n";
//matrixOut(log,invmatrix);
// put the sigmas in the usual order: upper diagonal (this time in normal form and not in band form)
unsigned k=0;
for (unsigned i=0;iscanField("sigma_"+getPntrToArgument(i)->getName(),sigma[i]);
}
}
ifile->scanField("height",height);
ifile->scanField("biasf",dummy);
if(ifile->FieldExist("clock")) ifile->scanField("clock",dummy);
if(ifile->FieldExist("lower_int")) ifile->scanField("lower_int",dummy);
if(ifile->FieldExist("upper_int")) ifile->scanField("upper_int",dummy);
ifile->scanField();
return true;
}else{
return false;
};
}
void MetaD::computeReweightingFactor(){
if( !welltemp_ ) error("cannot compute the c(t) reweighting factors for non well-tempered metadynamics");
// Recover the minimum values for the grid
unsigned ncv=getNumberOfArguments();
double grid_size=1.0; unsigned ntotgrid=1;
std::vector dmin( ncv ),dmax( ncv ), grid_spacing( ncv ), vals( ncv );
for(unsigned j=0;jgetMin()[j], dmin[j] );
Tools::convert( BiasGrid_->getMax()[j], dmax[j] );
grid_spacing[j] = ( dmax[j] - dmin[j] ) / static_cast( rewf_grid_[j] );
grid_size *= grid_spacing[j];
if( !getPntrToArgument(j)->isPeriodic() ) dmax[j] += grid_spacing[j];
ntotgrid *= rewf_grid_[j];
}
// Now sum over whole grid
reweight_factor=0.0; double* der=new double[ncv]; std::vector t_index( ncv );
double sum1=0.0; double sum2=0.0;
double afactor = biasf_ / (kbt_*(biasf_-1.0)); double afactor2 = 1.0 / (kbt_*(biasf_-1.0));
unsigned rank=comm.Get_rank(), stride=comm.Get_size();
for(unsigned i=rank;i=2 ) t_index[ncv-1]=((kk-t_index[ncv-1])/rewf_grid_[ncv-2]);
for(unsigned j=0;jset(reweight_factor);
}
}
}