diff --git a/developer-doc/mdTemplate.txt b/developer-doc/mdTemplate.txt index 045da17dfe731c0da6b9608dd42d340fb66b40d9..98b7cf765c9bff6fe8604ca5f9e2e3478b6e1db1 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 eecaa01867e921f3f5964ad5ca4c61d5de248b12..b8b534d84b06ce7791f986b0a9f9dc28f32f43c5 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 3318fe1f4df0cea7fc1a7b995ca8e103cada66f6..576a3596d3200bccd18cbab4a4490382ae583ef8 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 1ca30cb3bfffdfba2e0b36ea8be0e1897f8b23b6..2d44a37582935c6bea8d39f165e050cb1c92e421 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 2c08c2678dd2b6df9029c0149058493c7a4df07f..f65aa9fe76630eb495850de1b48e368db2998984 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 0f2510e29663b1913965cc6931c827f54ee3e52f..839a08c11145ffb93666b39953476527072668ac 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 36d7600fd10aeef96bc67b3350dfed4ee95fa2ee..461a0665e1acc43cda1af9082c50868cb173a003 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 b9cded09331fda9c4fe63ee2657c478c9ea20076..0f41717d73ec1d26ef80e2057e710a95b6ef451f 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; +} + + }