diff --git a/patches/gromacs-4.6.5.diff/src/kernel/md.c b/patches/gromacs-4.6.5.diff/src/kernel/md.c
index c896ffadd94a9bfac98caa32822d7f8482161473..a96f3bcd40e5f62081515b9a7980411be94a632c 100644
--- a/patches/gromacs-4.6.5.diff/src/kernel/md.c
+++ b/patches/gromacs-4.6.5.diff/src/kernel/md.c
@@ -726,6 +726,13 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
/* PLUMED */
if(plumedswitch){
+/* detect plumed API version */
+ int pversion=0;
+ plumed_cmd(plumedmain,"getApiVersion",&pversion);
+/* setting kbT is only implemented with api>1) */
+ real kbT=ir->opts.ref_t[0]*BOLTZ;
+ if(pversion>1) plumed_cmd(plumedmain,"setKbT",&kbT);
+
if(cr->ms && cr->ms->nsim>1) {
if(MASTER(cr)) plumed_cmd(plumedmain,"GREX setMPIIntercomm",&cr->ms->mpi_comm_masters);
if(PAR(cr)){
diff --git a/src/analysis/Analysis.cpp b/src/analysis/Analysis.cpp
index 9e3a6cc9f2b31d59ea77b33e78eb4c865f629334..425f51a8553e28e420190e890ec5fa1e1d013c20 100644
--- a/src/analysis/Analysis.cpp
+++ b/src/analysis/Analysis.cpp
@@ -150,14 +150,19 @@ argument_names(getNumberOfArguments())
}
}
- simtemp=0.; parse("TEMP",simtemp);
- if( simtemp==0 && !biases.empty() ) error("to reweight you must specify a temperature use TEMP");
- rtemp=0;
+ rtemp=0;
if( keywords.exists("REWEIGHT_TEMP") ) parse("REWEIGHT_TEMP",rtemp);
if( rtemp!=0 ){
- if( simtemp==0 ) error("to reweight you must specify a temperature use TEMP");
- needeng=true;
- log.printf(" reweighting simulation at %f to probabilities at temperature %f\n",simtemp,rtemp);
+ needeng=true;
+ log.printf(" reweighting simulation to probabilities at temperature %f\n",rtemp);
+ rtemp*=plumed.getAtoms().getKBoltzmann();
+ }
+ simtemp=0.; parse("TEMP",simtemp);
+ if( rtemp>0 || !biases.empty() ){
+ if(simtemp>0) simtemp*=plumed.getAtoms().getKBoltzmann();
+ else simtemp*=plumed.getAtoms().getKbT();
+
+ if(simtemp==0) error("The MD engine does not pass the temperature to plumed so you have to specify it using TEMP");
}
parseFlag("USE_ALL_DATA",single_run);
@@ -260,10 +265,10 @@ void Analysis::calculate(){
if(needeng){
double energy=plumed.getAtoms().getEnergy()+bias;
// Reweighting because of temperature difference
- ww=-( (1.0/rtemp) - (1.0/simtemp) )*(energy+bias) / plumed.getAtoms().getKBoltzmann();
+ ww=-( (1.0/rtemp) - (1.0/simtemp) )*(energy+bias);
}
// Reweighting because of biases
- if( !biases.empty() ) ww += bias/( plumed.getAtoms().getKBoltzmann()*simtemp );
+ if( !biases.empty() ) ww += bias/ simtemp;
// Get the arguments ready to transfer to reference configuration
for(unsigned i=0;i<getNumberOfArguments();++i) current_args[i]=getArgument(i);
diff --git a/src/analysis/Analysis.h b/src/analysis/Analysis.h
index 22c30f8e433912308e5a43b462c561a0bb851d34..5b4343a6785f0dee31d5fd235b0eac7d119ec29d 100644
--- a/src/analysis/Analysis.h
+++ b/src/analysis/Analysis.h
@@ -121,7 +121,7 @@ protected:
bool getPeriodicityInformation(const unsigned& i, std::string& dmin, std::string& dmax);
/// Return the normalization constant
double getNormalization() const;
-/// Return the set temperature
+/// Return the set temperature (N.B. k_B T is what is returned by this function)
double getTemp () const;
/// Are we analyzing each data block separately (if we are not this also returns the old normalization )
bool usingMemory() const;
diff --git a/src/analysis/Histogram.cpp b/src/analysis/Histogram.cpp
index 3fa0c092f993977428f402cc8a4f54bca8af8890..b4f984169f33fed41e4401406538e903e9bbe29f 100644
--- a/src/analysis/Histogram.cpp
+++ b/src/analysis/Histogram.cpp
@@ -258,7 +258,7 @@ void Histogram::performAnalysis(){
// Normalize the histogram
gg->scaleAllValuesAndDerivatives( 1.0 / getNormalization() );
if(fenergy) {
- gg->logAllValuesAndDerivatives( -getTemp() * plumed.getAtoms().getKBoltzmann() );
+ gg->logAllValuesAndDerivatives( -getTemp() );
gg->setMinToZero();
}
diff --git a/src/bias/MetaD.cpp b/src/bias/MetaD.cpp
index 3266464283f9b8be796400fd9d903ed9b0b0e3fe..f9f2619d0d3c1d3bb41a46de1d73c0fcd79f9c5f 100644
--- a/src/bias/MetaD.cpp
+++ b/src/bias/MetaD.cpp
@@ -245,7 +245,7 @@ private:
bool grid_,hasextgrid_;
double height0_;
double biasf_;
- double temp_;
+ double kbt_;
int stride_;
bool welltemp_;
double* dp_;
@@ -342,7 +342,7 @@ PLUMED_BIAS_INIT(ao),
// Grid stuff initialization
BiasGrid_(NULL),ExtGrid_(NULL), wgridstride_(0), grid_(false), hasextgrid_(false),
// Metadynamics basic parameters
-height0_(std::numeric_limits<double>::max()), biasf_(1.0), temp_(0.0),
+height0_(std::numeric_limits<double>::max()), biasf_(1.0), kbt_(0.0),
stride_(0), welltemp_(false),
// Other stuff
dp_(NULL), adaptive_(FlexibleBin::none),
@@ -417,9 +417,12 @@ isFirstStep(true)
parse("FILE",hillsfname);
parse("BIASFACTOR",biasf_);
if( biasf_<1.0 ) error("well tempered bias factor is nonsensical");
- parse("TEMP",temp_);
+ double temp;
+ parse("TEMP",temp);
+ if(temp>0.0) kbt_=plumed.getAtoms().getKBoltzmann()*temp;
+ else kbt_=plumed.getAtoms().getKbT();
if(biasf_>1.0){
- if(temp_==0.0) error("if you are doing well tempered metadynamics you must specify the temperature using TEMP");
+ 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;
@@ -427,11 +430,11 @@ isFirstStep(true)
if(tau==0.0){
if(height0_==std::numeric_limits<double>::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=(plumed.getAtoms().getKBoltzmann()*temp_*(biasf_-1.0))/height0_*getTimeStep()*stride_;
+ 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<double>::max()) error("At most one between HEIGHT and TAU should be specified");
- height0_=(plumed.getAtoms().getKBoltzmann()*temp_*(biasf_-1.0))/tau*getTimeStep()*stride_;
+ height0_=(kbt_*(biasf_-1.0))/tau*getTimeStep()*stride_;
}
// Grid Stuff
@@ -537,6 +540,7 @@ isFirstStep(true)
if(welltemp_){
log.printf(" Well-Tempered Bias Factor %f\n",biasf_);
log.printf(" Hills relaxation time (tau) %f\n",tau);
+ log.printf(" KbT %f\n",kbt_);
}
if(doInt_) log.printf(" Upper and Lower limits boundaries for the bias are activated at %f - %f\n", lowI_, uppI_);
if(grid_){
@@ -971,7 +975,7 @@ double MetaD::getHeight(const vector<double>& cv)
double height=height0_;
if(welltemp_){
double vbias=getBiasAndDerivatives(cv);
- height=height0_*exp(-vbias/(plumed.getAtoms().getKBoltzmann()*temp_*(biasf_-1.0)));
+ height=height0_*exp(-vbias/(kbt_*(biasf_-1.0)));
}
return height;
}
@@ -993,7 +997,7 @@ void MetaD::calculate()
getPntrToComponent("bias")->set(ene);
// calculate the acceleration factor
if(acceleration&&!isFirstStep) {
- acc += exp(ene/(plumed.getAtoms().getKBoltzmann()*temp_));
+ acc += exp(ene/(kbt_));
double mean_acc = acc/((double) getStep());
getPntrToComponent("acc")->set(mean_acc);
}
diff --git a/src/core/Atoms.cpp b/src/core/Atoms.cpp
index 2ce582a1e1b119f887eb14b57a9a7d278671aaaa..1fa1b04da19456ac633b89a68c49bdc38aa9823d 100644
--- a/src/core/Atoms.cpp
+++ b/src/core/Atoms.cpp
@@ -55,7 +55,8 @@ Atoms::Atoms(PlumedMain&plumed):
plumed(plumed),
naturalUnits(false),
timestep(0.0),
- forceOnEnergy(0.0)
+ forceOnEnergy(0.0),
+ kbT(0.0)
{
mdatoms=MDAtomsBase::create(sizeof(double));
}
@@ -356,6 +357,15 @@ double Atoms::getTimeStep()const{
return timestep/units.getTime()*MDUnits.getTime();
}
+void Atoms::setKbT(void*p){
+ MD2double(p,kbT);
+}
+
+double Atoms::getKbT()const{
+ return kbT/units.getEnergy()*MDUnits.getEnergy();
+}
+
+
void Atoms::createFullList(int*n){
vector<AtomNumber> fullListTmp;
for(unsigned i=0;i<actions.size();i++) if(actions[i]->isActive())
diff --git a/src/core/Atoms.h b/src/core/Atoms.h
index ee3daeae2ec2d3a28bc7638fb18686c51847c05c..3d5ca4030c8d3203e4d53a13839307d24510e7d6 100644
--- a/src/core/Atoms.h
+++ b/src/core/Atoms.h
@@ -90,6 +90,8 @@ class Atoms
double timestep;
double forceOnEnergy;
+ double kbT;
+
std::vector<const ActionAtomistic*> actions;
std::vector<int> gatindex;
@@ -137,6 +139,9 @@ public:
void setTimeStep(void*);
double getTimeStep()const;
+ void setKbT(void*);
+ double getKbT()const;
+
void setNatoms(int);
const int & getNatoms()const;
diff --git a/src/core/PlumedMain.cpp b/src/core/PlumedMain.cpp
index 5319d29fd5ce670ee2114f18081c854bbd93a2d5..e3a6c43d86f791a805b1cb0ac42fe300090eff65 100644
--- a/src/core/PlumedMain.cpp
+++ b/src/core/PlumedMain.cpp
@@ -46,7 +46,7 @@
using namespace std;
-enum { SETBOX, SETPOSITIONS, SETMASSES, SETCHARGES, SETPOSITIONSX, SETPOSITIONSY, SETPOSITIONSZ, SETVIRIAL, SETENERGY, SETFORCES, SETFORCESX, SETFORCESY, SETFORCESZ, CALC, PREPAREDEPENDENCIES, SHAREDATA, PREPARECALC, PERFORMCALC, SETSTEP, SETSTEPLONG, SETATOMSNLOCAL, SETATOMSGATINDEX, SETATOMSCONTIGUOUS, CREATEFULLLIST, GETFULLLIST, CLEARFULLLIST, READ, CLEAR, GETAPIVERSION, INIT, SETREALPRECISION, SETMDLENGTHUNITS, SETMDENERGYUNITS, SETMDTIMEUNITS, SETNATURALUNITS, SETNOVIRIAL, SETPLUMEDDAT, SETMPICOMM, SETMPIFCOMM, SETMPIMULTISIMCOMM, SETNATOMS, SETTIMESTEP, SETMDENGINE, SETLOG, SETLOGFILE, SETSTOPFLAG, GETEXCHANGESFLAG, SETEXCHANGESSEED, SETNUMBEROFREPLICAS, GETEXCHANGESLIST, RUNFINALJOBS, ISENERGYNEEDED, GETBIAS };
+enum { SETBOX, SETPOSITIONS, SETMASSES, SETCHARGES, SETPOSITIONSX, SETPOSITIONSY, SETPOSITIONSZ, SETVIRIAL, SETENERGY, SETFORCES, SETFORCESX, SETFORCESY, SETFORCESZ, CALC, PREPAREDEPENDENCIES, SHAREDATA, PREPARECALC, PERFORMCALC, SETSTEP, SETSTEPLONG, SETATOMSNLOCAL, SETATOMSGATINDEX, SETATOMSCONTIGUOUS, CREATEFULLLIST, GETFULLLIST, CLEARFULLLIST, READ, CLEAR, GETAPIVERSION, INIT, SETREALPRECISION, SETMDLENGTHUNITS, SETMDENERGYUNITS, SETMDTIMEUNITS, SETNATURALUNITS, SETNOVIRIAL, SETPLUMEDDAT, SETMPICOMM, SETMPIFCOMM, SETMPIMULTISIMCOMM, SETNATOMS, SETTIMESTEP, SETMDENGINE, SETLOG, SETLOGFILE, SETSTOPFLAG, GETEXCHANGESFLAG, SETEXCHANGESSEED, SETNUMBEROFREPLICAS, GETEXCHANGESLIST, RUNFINALJOBS, ISENERGYNEEDED, GETBIAS, SETKBT };
namespace PLMD{
@@ -130,6 +130,7 @@ PlumedMain::PlumedMain():
word_map["runFinalJobs"]=RUNFINALJOBS;
word_map["isEnergyNeeded"]=ISENERGYNEEDED;
word_map["getBias"]=GETBIAS;
+ word_map["setKbT"]=SETKBT;
}
PlumedMain::~PlumedMain(){
@@ -292,7 +293,7 @@ void PlumedMain::cmd(const std::string & word,void*val){
break;
case GETAPIVERSION:
CHECK_NULL(val,word);
- *(static_cast<int*>(val))=1;
+ *(static_cast<int*>(val))=2;
break;
// commands which can be used only before initialization:
case INIT:
@@ -362,6 +363,11 @@ void PlumedMain::cmd(const std::string & word,void*val){
CHECK_NULL(val,word);
atoms.setTimeStep(val);
break;
+ case SETKBT: /* ADDED WITH API==2 */
+ CHECK_NOTINIT(initialized,word);
+ CHECK_NULL(val,word);
+ atoms.setKbT(val);
+ break;
case SETMDENGINE:
CHECK_NOTINIT(initialized,word);
CHECK_NULL(val,word);
@@ -469,6 +475,12 @@ void PlumedMain::init(){
}
atoms.updateUnits();
log.printf("Timestep: %f\n",atoms.getTimeStep());
+ if(atoms.getKbT()>0.0)
+ log.printf("KbT: %f\n",atoms.getKbT());
+ else {
+ log.printf("KbT has not been set by the MD engine\n");
+ log.printf("It should be set by hand where needed\n");
+ }
log<<"Relevant bibliography:\n";
log<<citations;
log<<"Please read and cite where appropriate!\n";