diff --git a/regtest/basic/rt-make-extracv/COLVARX.reference b/regtest/basic/rt-make-extracv/COLVARX.reference new file mode 100644 index 0000000000000000000000000000000000000000..4e16cc2c16273656472754872c9d1b83fbe66a13 --- /dev/null +++ b/regtest/basic/rt-make-extracv/COLVARX.reference @@ -0,0 +1,11 @@ +#! FIELDS time e + 0.000000 0.000000 + 0.000000 1.000000 + 0.000000 2.000000 + 0.000000 3.000000 + 0.000000 4.000000 + 0.000000 5.000000 + 0.000000 6.000000 + 0.000000 7.000000 + 0.000000 8.000000 + 0.000000 9.000000 diff --git a/regtest/basic/rt-make-extracv/Makefile b/regtest/basic/rt-make-extracv/Makefile new file mode 100644 index 0000000000000000000000000000000000000000..3703b27cea227aa053fb6d1d73f861e4384dbcee --- /dev/null +++ b/regtest/basic/rt-make-extracv/Makefile @@ -0,0 +1 @@ +include ../../scripts/test.make diff --git a/regtest/basic/rt-make-extracv/config b/regtest/basic/rt-make-extracv/config new file mode 100644 index 0000000000000000000000000000000000000000..df1f95bf3ee289aa8367431334c00cf144754ddf --- /dev/null +++ b/regtest/basic/rt-make-extracv/config @@ -0,0 +1 @@ +type=make diff --git a/regtest/basic/rt-make-extracv/main.cpp b/regtest/basic/rt-make-extracv/main.cpp new file mode 100644 index 0000000000000000000000000000000000000000..33ce4647e4695b718e3e413bf439be219333f274 --- /dev/null +++ b/regtest/basic/rt-make-extracv/main.cpp @@ -0,0 +1,83 @@ +#include "plumed/wrapper/Plumed.h" +#include <vector> +#include <fstream> + +using namespace PLMD; + +int main(){ + Plumed* plumed=new Plumed; + + int natoms=10; + + std::vector<double> positions(3*natoms,0.0); + for(unsigned i=0;i<natoms;i++) positions[i]=i; + std::vector<double> masses(natoms,1.0); + std::vector<double> forces(3*natoms,0.0); + std::vector<double> box(9,0.0); + std::vector<double> virial(9,0.0); + + plumed->cmd("setNatoms",&natoms); + plumed->cmd("setLogFile","test.log"); + plumed->cmd("init"); + plumed->cmd("readInputLine","d: DISTANCE ATOMS=1,2"); + plumed->cmd("readInputLine","e: EXTRACV NAME=extra"); + plumed->cmd("readInputLine","e2: EXTRACV NAME=extra2"); + plumed->cmd("readInputLine","PRINT ARG=e FILE=COLVARX"); + plumed->cmd("readInputLine","RESTRAINT ARG=e AT=0 KAPPA=1"); + plumed->cmd("readInputLine","RESTRAINT ARG=d AT=0 KAPPA=1"); + plumed->cmd("readInputLine","RESTRAINT ARG=e2 AT=0 KAPPA=-1"); + + std::ofstream ofs("output"); + + for(int step=0;step<10;step++){ + double extracv=step; + double extracvf=0.0; + double extracv2=step; + double extracvf2=-step; + plumed->cmd("setStep",&step); + plumed->cmd("setPositions",&positions[0]); + plumed->cmd("setBox",&box[0]); + plumed->cmd("setForces",&forces[0]); + plumed->cmd("setVirial",&virial[0]); + plumed->cmd("setExtraCV extra",&extracv); + plumed->cmd("setExtraCVForce extra",&extracvf); + plumed->cmd("setExtraCV extra2",&extracv2); + plumed->cmd("setExtraCVForce extra2",&extracvf2); + plumed->cmd("setMasses",&masses[0]); +// first compute using modified positions: + positions[0]=0.5; + extracv2=100; + for(auto & f:forces) f=0.0; + extracvf=0.0; + extracvf2=0.0; + plumed->cmd("prepareCalc"); + plumed->cmd("performCalcNoUpdate"); + double bias; + plumed->cmd("getBias",&bias); + ofs<<"bias_pre: "<<bias<<"\n"; + ofs<<"extracvf_pre: "<<extracvf<<"\n"; + ofs<<"extracvf2_pre: "<<extracvf2<<"\n"; + ofs<<"f_pre:"; + for(auto & f:forces) ofs<<" "<<f; + ofs<<"\n"; +// then compute using regular positions: + positions[0]=0; + extracv2=0; + for(auto & f:forces) f=0.0; + extracvf=0.0; + extracvf2=0.0; + plumed->cmd("prepareCalc"); + plumed->cmd("performCalcNoUpdate"); + plumed->cmd("getBias",&bias); + ofs<<"bias: "<<bias<<"\n"; + ofs<<"extracvf: "<<extracvf<<"\n"; + ofs<<"extracvf2: "<<extracvf2<<"\n"; + ofs<<"f:"; + for(auto & f:forces) ofs<<" "<<f; + ofs<<"\n"; + plumed->cmd("update"); + } + + delete plumed; + return 0; +} diff --git a/regtest/basic/rt-make-extracv/output.reference b/regtest/basic/rt-make-extracv/output.reference new file mode 100644 index 0000000000000000000000000000000000000000..635ba78fc5e684ed12ee45e0aabbea24d6ea838d --- /dev/null +++ b/regtest/basic/rt-make-extracv/output.reference @@ -0,0 +1,80 @@ +bias_pre: -4987.88 +extracvf_pre: 0 +extracvf2_pre: 100 +f_pre: 2.5 3 3 -2.5 -3 -3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +bias: 13.5 +extracvf: 0 +extracvf2: 0 +f: 3 3 3 -3 -3 -3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +bias_pre: -4987.38 +extracvf_pre: -1 +extracvf2_pre: 100 +f_pre: 2.5 3 3 -2.5 -3 -3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +bias: 14 +extracvf: -1 +extracvf2: 0 +f: 3 3 3 -3 -3 -3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +bias_pre: -4985.88 +extracvf_pre: -2 +extracvf2_pre: 100 +f_pre: 2.5 3 3 -2.5 -3 -3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +bias: 15.5 +extracvf: -2 +extracvf2: 0 +f: 3 3 3 -3 -3 -3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +bias_pre: -4983.38 +extracvf_pre: -3 +extracvf2_pre: 100 +f_pre: 2.5 3 3 -2.5 -3 -3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +bias: 18 +extracvf: -3 +extracvf2: 0 +f: 3 3 3 -3 -3 -3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +bias_pre: -4979.88 +extracvf_pre: -4 +extracvf2_pre: 100 +f_pre: 2.5 3 3 -2.5 -3 -3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +bias: 21.5 +extracvf: -4 +extracvf2: 0 +f: 3 3 3 -3 -3 -3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +bias_pre: -4975.38 +extracvf_pre: -5 +extracvf2_pre: 100 +f_pre: 2.5 3 3 -2.5 -3 -3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +bias: 26 +extracvf: -5 +extracvf2: 0 +f: 3 3 3 -3 -3 -3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +bias_pre: -4969.88 +extracvf_pre: -6 +extracvf2_pre: 100 +f_pre: 2.5 3 3 -2.5 -3 -3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +bias: 31.5 +extracvf: -6 +extracvf2: 0 +f: 3 3 3 -3 -3 -3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +bias_pre: -4963.38 +extracvf_pre: -7 +extracvf2_pre: 100 +f_pre: 2.5 3 3 -2.5 -3 -3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +bias: 38 +extracvf: -7 +extracvf2: 0 +f: 3 3 3 -3 -3 -3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +bias_pre: -4955.88 +extracvf_pre: -8 +extracvf2_pre: 100 +f_pre: 2.5 3 3 -2.5 -3 -3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +bias: 45.5 +extracvf: -8 +extracvf2: 0 +f: 3 3 3 -3 -3 -3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +bias_pre: -4947.38 +extracvf_pre: -9 +extracvf2_pre: 100 +f_pre: 2.5 3 3 -2.5 -3 -3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +bias: 54 +extracvf: -9 +extracvf2: 0 +f: 3 3 3 -3 -3 -3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 diff --git a/regtest/basic/rt-make-extracv/plumed.dat b/regtest/basic/rt-make-extracv/plumed.dat new file mode 100644 index 0000000000000000000000000000000000000000..28977ccf003f3805102b35395552024feb7caa9e --- /dev/null +++ b/regtest/basic/rt-make-extracv/plumed.dat @@ -0,0 +1,3 @@ +THIS IS A CRAPPY INPUT + +USED HERE TO BE SURE THAT PLUMED DOES NOT READ IT BY DEFAULT diff --git a/src/colvar/ExtraCV.cpp b/src/colvar/ExtraCV.cpp new file mode 100644 index 0000000000000000000000000000000000000000..1a3374b65c3a9a345644a10d56faa2c8f61b62aa --- /dev/null +++ b/src/colvar/ExtraCV.cpp @@ -0,0 +1,100 @@ +/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + Copyright (c) 2017 The plumed team + (see the PEOPLE file at the root of the distribution for a list of names) + + See http://www.plumed.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 <http://www.gnu.org/licenses/>. ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ +#include "Colvar.h" +#include "ActionRegister.h" +#include "core/PlumedMain.h" +#include "core/Atoms.h" + +#include <string> +#include <cmath> + +namespace PLMD { +namespace colvar { + +//+PLUMEDOC COLVAR EXTRACV +/* +XX + +\par Examples + +XX + +*/ +//+ENDPLUMEDOC + + +class ExtraCV : public Colvar { + std::string name; +public: + explicit ExtraCV(const ActionOptions&); +// active methods: + void prepare(); + virtual void calculate(); + unsigned getNumberOfDerivatives(); + static void registerKeywords( Keywords& keys ); +}; + + +using namespace std; + + +PLUMED_REGISTER_ACTION(ExtraCV,"EXTRACV") + +ExtraCV::ExtraCV(const ActionOptions&ao): + PLUMED_COLVAR_INIT(ao) +{ + addValueWithDerivatives(); setNotPeriodic(); + getPntrToValue()->resizeDerivatives(1); + parse("NAME",name); + log<<" name: "<<name<<"\n"; + isExtraCV=true; + setExtraCV(name); +} + +void ExtraCV::registerKeywords( Keywords& keys ) { + Action::registerKeywords( keys ); + ActionAtomistic::registerKeywords( keys ); + ActionWithValue::registerKeywords( keys ); + keys.remove("NUMERICAL_DERIVATIVES"); + keys.add("compulsory","NAME","name of the CV as computed by the MD engine"); +} + +unsigned ExtraCV::getNumberOfDerivatives() { + return 1; +} + +void ExtraCV::prepare() { +/// \todo: notify Atoms that this is requested +} + +// calculator +void ExtraCV::calculate() { + double value=plumed.getAtoms().getExtraCV(name); + setValue( value ); + getPntrToComponent(0)->addDerivative(0,1.0); +} + +} +} + + + diff --git a/src/core/ActionAtomistic.cpp b/src/core/ActionAtomistic.cpp index baba14a59ff9e54a273549c83712ed6e85cb7ef6..1b145f0e0d84675b77487ac489453a298c02cf29 100644 --- a/src/core/ActionAtomistic.cpp +++ b/src/core/ActionAtomistic.cpp @@ -262,6 +262,7 @@ void ActionAtomistic::applyForces() { for(unsigned j=0; j<indexes.size(); j++) f[indexes[j].index()]+=forces[j]; v+=virial; atoms.forceOnEnergy+=forceOnEnergy; + if(extraCV.length()>0) atoms.updateExtraCVForce(extraCV,forceOnExtraCV); } void ActionAtomistic::clearOutputForces() { @@ -269,6 +270,7 @@ void ActionAtomistic::clearOutputForces() { if(donotforce) return; for(unsigned i=0; i<forces.size(); ++i)forces[i].zero(); forceOnEnergy=0.0; + forceOnExtraCV=0.0; } diff --git a/src/core/ActionAtomistic.h b/src/core/ActionAtomistic.h index 12acea4db122825f5736adecafc0b3535616160a..793b47402f7bfb2f8370e84390697d21b045cf56 100644 --- a/src/core/ActionAtomistic.h +++ b/src/core/ActionAtomistic.h @@ -29,6 +29,7 @@ #include "tools/ForwardDecl.h" #include <vector> #include <set> +#include <map> namespace PLMD { @@ -58,6 +59,10 @@ class ActionAtomistic : std::vector<Vector> forces; // forces on the needed atoms double forceOnEnergy; + double forceOnExtraCV; + + std::string extraCV; + bool lockRequestAtoms; // forbid changes to request atoms bool donotretrieve; @@ -66,6 +71,8 @@ class ActionAtomistic : protected: Atoms& atoms; + void setExtraCV(const std::string &name); + public: /// Request an array of atoms. /// This method is used to ask for a list of atoms. Atoms @@ -118,6 +125,8 @@ public: Tensor & modifyVirial(); /// Get a reference to force on energy double & modifyForceOnEnergy(); +/// Get a reference to force on extraCV + double & modifyForceOnExtraCV(); /// Get number of available atoms unsigned getNumberOfAtoms()const {return indexes.size();} /// Compute the pbc distance between two positions @@ -261,6 +270,11 @@ double & ActionAtomistic::modifyForceOnEnergy() { return forceOnEnergy; } +inline +double & ActionAtomistic::modifyForceOnExtraCV() { + return forceOnExtraCV; +} + inline const Pbc & ActionAtomistic::getPbc() const { return pbc; @@ -296,6 +310,11 @@ Pbc & ActionAtomistic::modifyGlobalPbc() { return atoms.pbc; } +inline +void ActionAtomistic::setExtraCV(const std::string &name) { + extraCV=name; +} + } diff --git a/src/core/Atoms.cpp b/src/core/Atoms.cpp index ae844b20e108036a729485906eda6b241ea9e969..3680912bf892791ce38c9afb7e4296f35f1feafc 100644 --- a/src/core/Atoms.cpp +++ b/src/core/Atoms.cpp @@ -581,4 +581,20 @@ void Atoms::getLocalMDForces(std::vector<Vector>& localForces) { } } +void Atoms::setExtraCV(const std::string &name,void*p) { + mdatoms->setExtraCV(name,p); +} + +void Atoms::setExtraCVForce(const std::string &name,void*p) { + mdatoms->setExtraCVForce(name,p); +} + +double Atoms::getExtraCV(const std::string &name) { + return mdatoms->getExtraCV(name); +} + +void Atoms::updateExtraCVForce(const std::string &name,double f) { + mdatoms->updateExtraCVForce(name,f); +} + } diff --git a/src/core/Atoms.h b/src/core/Atoms.h index 75ee16f5f6641d1e71fcbcd4a26403c5b31ef4c9..39fc0821abf43f3d8ae4fc98833be23bc64bf25e 100644 --- a/src/core/Atoms.h +++ b/src/core/Atoms.h @@ -232,6 +232,11 @@ public: bool usingNaturalUnits()const; void setNaturalUnits(bool n) {naturalUnits=n;} void setMDNaturalUnits(bool n) {MDnaturalUnits=n;} + + void setExtraCV(const std::string &name,void*p); + void setExtraCVForce(const std::string &name,void*p); + double getExtraCV(const std::string &name); + void updateExtraCVForce(const std::string &name,double f); }; inline diff --git a/src/core/Colvar.cpp b/src/core/Colvar.cpp index d04592e67a69e33abf94af1fbce50a8328ef79f5..20464ee3bcecf137555729c4462f3006e38e282f 100644 --- a/src/core/Colvar.cpp +++ b/src/core/Colvar.cpp @@ -31,7 +31,8 @@ Colvar::Colvar(const ActionOptions&ao): Action(ao), ActionAtomistic(ao), ActionWithValue(ao), - isEnergy(false) + isEnergy(false), + isExtraCV(false) { } @@ -67,7 +68,7 @@ void Colvar::apply() { unsigned nt=OpenMP::getNumThreads(); if(nt>ncp/(4*stride)) nt=1; - if(!isEnergy) { + if(!isEnergy && !isExtraCV) { #pragma omp parallel num_threads(nt) { vector<Vector> omp_f(fsz); @@ -107,6 +108,9 @@ void Colvar::apply() { } else if( isEnergy ) { vector<double> forces(1); if(getPntrToComponent(0)->applyForce(forces)) modifyForceOnEnergy()+=forces[0]; + } else if( isExtraCV ) { + vector<double> forces(1); + if(getPntrToComponent(0)->applyForce(forces)) modifyForceOnExtraCV()+=forces[0]; } } diff --git a/src/core/Colvar.h b/src/core/Colvar.h index 914f1bbfceaa022cf48fa10c4674e91bdd2bcefb..e5f0804f1eb4d359be7cb7a2ff5f891965052dcc 100644 --- a/src/core/Colvar.h +++ b/src/core/Colvar.h @@ -43,6 +43,7 @@ class Colvar : private: protected: bool isEnergy; + bool isExtraCV; void requestAtoms(const std::vector<AtomNumber> & a); // Set the derivatives for a particular atom equal to the input Vector // This routine is called setAtomsDerivatives because not because you diff --git a/src/core/MDAtoms.cpp b/src/core/MDAtoms.cpp index 7dccc583c98c54d42088a967c761e9d13869236e..4e8d1fc82193240f90543e8ee8e167796512e227 100644 --- a/src/core/MDAtoms.cpp +++ b/src/core/MDAtoms.cpp @@ -26,6 +26,7 @@ #include "tools/Units.h" #include <algorithm> #include <string> +#include <map> using namespace std; @@ -48,6 +49,8 @@ class MDAtomsTyped: T *fx; T *fy; T *fz; T *box; T *virial; + std::map<std::string,T*> extraCV; + std::map<std::string,T*> extraCVForce; public: MDAtomsTyped(); void setm(void*m); @@ -59,6 +62,18 @@ public: void setp(void*p,int i); void setf(void*f,int i); void setUnits(const Units&,const Units&); + void setExtraCV(const std::string &name,void*p) { + extraCV[name]=static_cast<T*>(p); + } + void setExtraCVForce(const std::string &name,void*p) { + extraCVForce[name]=static_cast<T*>(p); + } + double getExtraCV(const std::string &name) { + return static_cast<double>(*extraCV[name]); + } + void updateExtraCVForce(const std::string &name,double f) { + *extraCVForce[name]+=static_cast<T>(f); + } void MD2double(const void*m,double&d)const { d=double(*(static_cast<const T*>(m))); } diff --git a/src/core/MDAtoms.h b/src/core/MDAtoms.h index eaa012c826925330bd56cbf9d23ca8ab8b45f6e9..67fd6d9989fefbd4f3e5c413ee417ced25cec6ad 100644 --- a/src/core/MDAtoms.h +++ b/src/core/MDAtoms.h @@ -108,6 +108,16 @@ public: /// Rescale all the forces, including the virial. /// It is applied to all atoms with local index going from 0 to index.size()-1 virtual void rescaleForces(const std::vector<int>&index,double factor)=0; + +/// Set a pointer to an extra CV. + virtual void setExtraCV(const std::string &name,void*p)=0; +/// Set a pointer to an extra CV force. + virtual void setExtraCVForce(const std::string &name,void*p)=0; +/// Retrieve the value of an extra CV. + virtual double getExtraCV(const std::string &name)=0; +/// Update the value of an extra CV force. +/// \todo check if this should also be scaled when acting on total energy + virtual void updateExtraCVForce(const std::string &name,double f)=0; }; } diff --git a/src/core/PlumedMain.cpp b/src/core/PlumedMain.cpp index 46d5bc09623d3db907a4c069355c2330c61dcdf5..69983877f82e48936ce922353d9a02aae4c4edd1 100644 --- a/src/core/PlumedMain.cpp +++ b/src/core/PlumedMain.cpp @@ -437,6 +437,16 @@ void PlumedMain::cmd(const std::string & word,void*val) { plumed_assert(nw==2); *(static_cast<int*>(val))=(actionRegister().check(words[1]) ? 1:0); break; + case cmd_setExtraCV: + CHECK_NOTNULL(val,word); + plumed_assert(nw==2); + atoms.setExtraCV(words[1],val); + break; + case cmd_setExtraCVForce: + CHECK_NOTNULL(val,word); + plumed_assert(nw==2); + atoms.setExtraCVForce(words[1],val); + break; case cmd_GREX: if(!grex) grex.reset(new GREX(*this)); plumed_massert(grex,"error allocating grex"); diff --git a/user-doc/Miscelaneous.md b/user-doc/Miscelaneous.md index 4708dbd6c35e1cc1ecfc103b7b1e9910c3ce612c..3742fca0a22ab86585ea034c1eb685b14356e15c 100644 --- a/user-doc/Miscelaneous.md +++ b/user-doc/Miscelaneous.md @@ -6,6 +6,7 @@ - \subpage BashAutocompletion - \subpage includes - \subpage load +- \subpage embed - \subpage degub - \subpage exchange-patterns - \subpage mymodules