From 7d48a57c272e2284bd35d3749fd6c21739d641d5 Mon Sep 17 00:00:00 2001
From: Giovanni Bussi <giovanni.bussi@gmail.com>
Date: Tue, 5 Jan 2016 13:33:56 +0100
Subject: [PATCH] Added units for masses and charges
I added the possibility to personalize charge and mass units
in the same way as length, energy, and time units were treated so far.
Notice that since two new cmd() strings have been added to set
units from an MD code, it is necessary to increase the API version
to 4.
I also added the appropriate developer documentation.
Fixes #179
---
developer-doc/mdTemplate.txt | 7 ++++++
src/cltools/Driver.cpp | 8 +++++++
src/core/Atoms.h | 2 ++
src/core/MDAtoms.cpp | 11 +++++++--
src/core/PlumedMain.cpp | 18 +++++++++++++--
src/setup/Units.cpp | 14 +++++++++++
src/tools/Units.cpp | 42 ++++++++++++++++++++++++++++++++-
src/tools/Units.h | 45 ++++++++++++++++++++++++++++++++++++
8 files changed, 142 insertions(+), 5 deletions(-)
diff --git a/developer-doc/mdTemplate.txt b/developer-doc/mdTemplate.txt
index 045da17df..98b7cf765 100644
--- a/developer-doc/mdTemplate.txt
+++ b/developer-doc/mdTemplate.txt
@@ -190,6 +190,13 @@ plumed_cmd(plumedmain,"setRealPrecision",&real_precision); // Pass a pointer
plumed_cmd(plumedmain,"setMDEnergyUnits",&energyUnits); // Pass a pointer to the conversion factor between the energy unit used in your code and kJ mol-1
plumed_cmd(plumedmain,"setMDLengthUnits",&lengthUnits); // Pass a pointer to the conversion factor between the length unit used in your code and nm
plumed_cmd(plumedmain,"setMDTimeUnits",&timeUnits); // Pass a pointer to the conversion factor between the time unit used in your code and ps
+
+// This is valid only if API VERSION > 3
+plumed_cmd(plumedmain,"setMDChargeUnits",&chargeUnits); // Pass a pointer to the conversion factor between the charge unit used in your code and e
+
+// This is valid only if API VERSION > 3
+plumed_cmd(plumedmain,"setMDMassUnits",&massUnits); // Pass a pointer to the conversion factor between the mass unit used in your code and amu
+
plumed_cmd(plumedmain,"setPlumedDat",&plumedInput); // Pass the name of the plumed input file from the md code to plumed
plumed_cmd(plumedmain,"setMPIComm",&MPI_COMM_WORLD); // Pass a pointer to the MPI communicator to plumed
// notice that from fortran the command "setMPIFComm" should be used instead
diff --git a/src/cltools/Driver.cpp b/src/cltools/Driver.cpp
index eecaa0186..b8b534d84 100644
--- a/src/cltools/Driver.cpp
+++ b/src/cltools/Driver.cpp
@@ -187,6 +187,8 @@ void Driver<real>::registerKeywords( Keywords& keys ){
keys.add("atoms","--itrr","the trajectory in trr format (xdrfile implementation)");
#endif
keys.add("optional","--length-units","units for length, either as a string or a number");
+ keys.add("optional","--mass-units","units for mass in pdb and mc file, either as a string or a number");
+ keys.add("optional","--charge-units","units for charge in pdb and mc file, either as a string or a number");
keys.add("optional","--dump-forces","dump the forces on a file");
keys.add("optional","--dump-forces-fmt","( default=%%f ) the format to use to dump the forces");
keys.addFlag("--dump-full-virial",false,"with --dump-forces, it dumps the 9 components of the virial");
@@ -375,6 +377,10 @@ int Driver<real>::main(FILE* in,FILE*out,Communicator& pc){
}
string lengthUnits(""); parse("--length-units",lengthUnits);
if(lengthUnits.length()>0) units.setLength(lengthUnits);
+ string chargeUnits(""); parse("--charge-units",chargeUnits);
+ if(chargeUnits.length()>0) units.setCharge(chargeUnits);
+ string massUnits(""); parse("--mass-units",massUnits);
+ if(massUnits.length()>0) units.setMass(massUnits);
parse("--pdb",pdbfile);
if(pdbfile.length()>0){
@@ -423,6 +429,8 @@ int Driver<real>::main(FILE* in,FILE*out,Communicator& pc){
p.cmd("setMPIComm",&intracomm.Get_comm());
}
p.cmd("setMDLengthUnits",&units.getLength());
+ p.cmd("setMDChargeUnits",&units.getCharge());
+ p.cmd("setMDMassUnits",&units.getMass());
p.cmd("setMDEngine","driver");
p.cmd("setTimestep",×tep);
p.cmd("setPlumedDat",plumedFile.c_str());
diff --git a/src/core/Atoms.h b/src/core/Atoms.h
index 3318fe1f4..576a3596d 100644
--- a/src/core/Atoms.h
+++ b/src/core/Atoms.h
@@ -194,6 +194,8 @@ public:
void setMDEnergyUnits(double d){MDUnits.setEnergy(d);}
void setMDLengthUnits(double d){MDUnits.setLength(d);}
void setMDTimeUnits(double d){MDUnits.setTime(d);}
+ void setMDChargeUnits(double d){MDUnits.setCharge(d);}
+ void setMDMassUnits(double d){MDUnits.setMass(d);}
const Units& getMDUnits(){return MDUnits;}
void setUnits(const Units&u){units=u;}
const Units& getUnits(){return units;}
diff --git a/src/core/MDAtoms.cpp b/src/core/MDAtoms.cpp
index 1ca30cb3b..2d44a3758 100644
--- a/src/core/MDAtoms.cpp
+++ b/src/core/MDAtoms.cpp
@@ -39,6 +39,7 @@ public MDAtomsBase
{
T scalep,scalef;
T scaleb,scalev;
+ T scalec,scalem; // factor to scale charges and masses
int stride;
T *m;
T *c;
@@ -79,12 +80,16 @@ template <class T>
void MDAtomsTyped<T>::setUnits(const Units& units,const Units& MDUnits){
double lscale=units.getLength()/MDUnits.getLength();
double escale=units.getEnergy()/MDUnits.getEnergy();
+ double cscale=units.getCharge()/MDUnits.getCharge();
+ double mscale=units.getMass()/MDUnits.getMass();
// scalep and scaleb are used to convert MD to plumed
scalep=1.0/lscale;
scaleb=1.0/lscale;
// scalef and scalev are used to convert plumed to MD
scalef=escale/lscale;
scalev=escale;
+ scalec=1.0/cscale;
+ scalem=1.0/mscale;
}
template <class T>
@@ -127,13 +132,13 @@ void MDAtomsTyped<T>::getLocalPositions(vector<Vector>&positions)const{
template <class T>
void MDAtomsTyped<T>::getMasses(const vector<int>&index,vector<double>&masses)const{
- if(m) for(unsigned i=0;i<index.size();++i) masses[index[i]]=m[i];
+ if(m) for(unsigned i=0;i<index.size();++i) masses[index[i]]=scalem*m[i];
else for(unsigned i=0;i<index.size();++i) masses[index[i]]=0.0;
}
template <class T>
void MDAtomsTyped<T>::getCharges(const vector<int>&index,vector<double>&charges)const{
- if(c) for(unsigned i=0;i<index.size();++i) charges[index[i]]=c[i];
+ if(c) for(unsigned i=0;i<index.size();++i) charges[index[i]]=scalec*c[i];
else for(unsigned i=0;i<index.size();++i) charges[index[i]]=0.0;
}
@@ -238,6 +243,8 @@ MDAtomsTyped<T>::MDAtomsTyped():
scalef(1.0),
scaleb(1.0),
scalev(1.0),
+ scalec(1.0),
+ scalem(1.0),
stride(0),
m(NULL),
c(NULL),
diff --git a/src/core/PlumedMain.cpp b/src/core/PlumedMain.cpp
index 2c08c2678..f65aa9fe7 100644
--- a/src/core/PlumedMain.cpp
+++ b/src/core/PlumedMain.cpp
@@ -47,7 +47,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, SETATOMSFGATINDEX, 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, SETRESTART };
+enum { SETBOX, SETPOSITIONS, SETMASSES, SETCHARGES, SETPOSITIONSX, SETPOSITIONSY, SETPOSITIONSZ, SETVIRIAL, SETENERGY, SETFORCES, SETFORCESX, SETFORCESY, SETFORCESZ, CALC, PREPAREDEPENDENCIES, SHAREDATA, PREPARECALC, PERFORMCALC, SETSTEP, SETSTEPLONG, SETATOMSNLOCAL, SETATOMSGATINDEX, SETATOMSFGATINDEX, SETATOMSCONTIGUOUS, CREATEFULLLIST, GETFULLLIST, CLEARFULLLIST, READ, CLEAR, GETAPIVERSION, INIT, SETREALPRECISION, SETMDLENGTHUNITS, SETMDENERGYUNITS, SETMDTIMEUNITS, SETMDCHARGEUNITS, SETMDMASSUNITS, SETNATURALUNITS, SETNOVIRIAL, SETPLUMEDDAT, SETMPICOMM, SETMPIFCOMM, SETMPIMULTISIMCOMM, SETNATOMS, SETTIMESTEP, SETMDENGINE, SETLOG, SETLOGFILE, SETSTOPFLAG, GETEXCHANGESFLAG, SETEXCHANGESSEED, SETNUMBEROFREPLICAS, GETEXCHANGESLIST, RUNFINALJOBS, ISENERGYNEEDED, GETBIAS, SETKBT, SETRESTART };
namespace PLMD{
@@ -114,6 +114,8 @@ PlumedMain::PlumedMain():
word_map["setMDLengthUnits"]=SETMDLENGTHUNITS;
word_map["setMDEnergyUnits"]=SETMDENERGYUNITS;
word_map["setMDTimeUnits"]=SETMDTIMEUNITS;
+ word_map["setMDMassUnits"]=SETMDMASSUNITS;
+ word_map["setMDChargeUnits"]=SETMDCHARGEUNITS;
word_map["setNaturalUnits"]=SETNATURALUNITS;
word_map["setNoVirial"]=SETNOVIRIAL;
word_map["setPlumedDat"]=SETPLUMEDDAT;
@@ -301,7 +303,7 @@ void PlumedMain::cmd(const std::string & word,void*val){
break;
case GETAPIVERSION:
CHECK_NULL(val,word);
- *(static_cast<int*>(val))=3;
+ *(static_cast<int*>(val))=4;
break;
// commands which can be used only before initialization:
case INIT:
@@ -319,6 +321,18 @@ void PlumedMain::cmd(const std::string & word,void*val){
atoms.MD2double(val,d);
atoms.setMDLengthUnits(d);
break;
+ case SETMDCHARGEUNITS:
+ CHECK_NOTINIT(initialized,word);
+ CHECK_NULL(val,word);
+ atoms.MD2double(val,d);
+ atoms.setMDChargeUnits(d);
+ break;
+ case SETMDMASSUNITS:
+ CHECK_NOTINIT(initialized,word);
+ CHECK_NULL(val,word);
+ atoms.MD2double(val,d);
+ atoms.setMDMassUnits(d);
+ break;
case SETMDENERGYUNITS:
CHECK_NOTINIT(initialized,word);
CHECK_NULL(val,word);
diff --git a/src/setup/Units.cpp b/src/setup/Units.cpp
index 0f2510e29..839a08c11 100644
--- a/src/setup/Units.cpp
+++ b/src/setup/Units.cpp
@@ -73,6 +73,8 @@ void Units::registerKeywords( Keywords& keys ){
keys.add("optional","LENGTH","the units of lengths. Either specify a conversion factor from the default, nm, or A (for angstroms) or um");
keys.add("optional","ENERGY","the units of energy. Either specify a conversion factor from the default, kj/mol, or use j/mol or kcal/mol");
keys.add("optional","TIME","the units of time. Either specify a conversion factor from the default, ps, or use ns or fs");
+ keys.add("optional","MASS","the units of masses. Specify a conversion factor from the default, amu");
+ keys.add("optional","CHARGE","the units of charges. Specify a conversion factor from the default, e");
keys.addFlag("NATURAL",false,"use natural units");
}
@@ -102,6 +104,18 @@ ActionSetup(ao)
if(u.getTimeString().length()>0) log.printf(" time: %s\n",u.getTimeString().c_str());
else log.printf(" time: %f ps\n",u.getTime());
+ s="";
+ parse("CHARGE",s);
+ if(s.length()>0) u.setCharge(s);
+ if(u.getChargeString().length()>0) log.printf(" time: %s\n",u.getChargeString().c_str());
+ else log.printf(" time: %f e\n",u.getCharge());
+
+ s="";
+ parse("MASS",s);
+ if(s.length()>0) u.setMass(s);
+ if(u.getMassString().length()>0) log.printf(" time: %s\n",u.getMassString().c_str());
+ else log.printf(" time: %f amu\n",u.getMass());
+
bool natural=false;
parseFlag("NATURAL",natural);
plumed.getAtoms().setNaturalUnits(natural);
diff --git a/src/tools/Units.cpp b/src/tools/Units.cpp
index 36d7600fd..461a0665e 100644
--- a/src/tools/Units.cpp
+++ b/src/tools/Units.cpp
@@ -32,7 +32,11 @@ Units::Units():
length(1.0),
lengthString("nm"),
time(1.0),
- timeString("ps")
+ timeString("ps"),
+ charge(1.0),
+ chargeString("e"),
+ mass(1.0),
+ massString("amu")
{
}
@@ -86,6 +90,30 @@ void Units::setTime(const std::string &s){
}
}
+void Units::setCharge(const std::string &s){
+ chargeString=s;
+ if(s=="e"){
+ charge=1.0;
+ } else {
+ charge=-1.0;
+ chargeString="";
+ Tools::convert(s,charge);
+ plumed_massert(charge>0.0,"charge units should be positive");
+ }
+}
+
+void Units::setMass(const std::string &s){
+ massString=s;
+ if(s=="amu"){
+ mass=1.0;
+ } else {
+ mass=-1.0;
+ massString="";
+ Tools::convert(s,mass);
+ plumed_massert(mass>0.0,"mass units should be positive");
+ }
+}
+
void Units::setEnergy(const double s){
energyString="";
energy=s;
@@ -101,5 +129,17 @@ void Units::setTime(const double s){
time=s;
}
+void Units::setCharge(const double s){
+ chargeString="";
+ charge=s;
+}
+
+void Units::setMass(const double s){
+ massString="";
+ mass=s;
+}
+
+
+
}
diff --git a/src/tools/Units.h b/src/tools/Units.h
index b9cded093..0f41717d7 100644
--- a/src/tools/Units.h
+++ b/src/tools/Units.h
@@ -48,6 +48,12 @@ class Units{
/// Units for time, expressed in ps (e.g. 0.001 means fs)
double time;
std::string timeString;
+/// Units for charges, expressed in proton charge (e.g. 1./18.2223 are sqrt(kcal/mol*A), as used in Amber)
+ double charge;
+ std::string chargeString;
+/// Units for masses, expressed in amu
+ double mass;
+ std::string massString;
public:
/// Constructor, setting default values (1.0)
Units();
@@ -63,6 +69,10 @@ public:
/// Also understands the following strings:
/// nm, A, um.
void setLength(const std::string &);
+/// Set charge units from string.
+ void setCharge(const std::string &);
+/// Set mass units from string.
+ void setMass(const std::string &);
/// Set energy units from double.
/// Should be specified in units of kj/mol (e.g. 4.184 means kcal/mol)
void setEnergy(double);
@@ -72,18 +82,32 @@ public:
/// Set lenght units from double.
/// Should be specified in units of nm (e.g. 0.1 means A)
void setLength(double);
+/// Set charge units from double.
+/// Should be specified in units of proton charge.
+ void setCharge(double);
+/// Set mass units from double.
+/// Should be specified in units of amu.
+ void setMass(double);
/// Get energy units as double.
const double & getEnergy()const;
/// Get length units as double.
const double & getLength()const;
/// Get time units as double.
const double & getTime()const;
+/// Get charge units as double.
+ const double & getCharge()const;
+/// Get mass units as double.
+ const double & getMass()const;
/// Get energy units as string.
const std::string & getEnergyString()const;
/// Get length units as string.
const std::string & getLengthString()const;
/// Get time units as string.
const std::string & getTimeString()const;
+/// Get charge units as string.
+ const std::string & getChargeString()const;
+/// Get mass units as string.
+ const std::string & getMassString()const;
};
inline
@@ -101,6 +125,16 @@ const double & Units::getTime()const{
return time;
}
+inline
+const double & Units::getCharge()const{
+ return charge;
+}
+
+inline
+const double & Units::getMass()const{
+ return mass;
+}
+
inline
const std::string & Units::getEnergyString()const{
return energyString;
@@ -116,6 +150,17 @@ const std::string & Units::getTimeString()const{
return timeString;
}
+inline
+const std::string & Units::getChargeString()const{
+ return chargeString;
+}
+
+inline
+const std::string & Units::getMassString()const{
+ return massString;
+}
+
+
}
--
GitLab