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 35f6086d23e944bb56b862d1fd99807ab2301f17..cfb16c8fd165f0038f91715e5fed800ec1b3de3e 100644 --- a/src/core/ActionAtomistic.cpp +++ b/src/core/ActionAtomistic.cpp @@ -251,6 +251,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() { @@ -258,6 +259,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 46ce9ebc2bcef2f7edc7dc487bc1c0fb3c56dc62..a4285bebcddf819660ed1682a12b6451543e8fdb 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 @@ -115,6 +122,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 @@ -258,6 +267,11 @@ double & ActionAtomistic::modifyForceOnEnergy() { return forceOnEnergy; } +inline +double & ActionAtomistic::modifyForceOnExtraCV() { + return forceOnExtraCV; +} + inline const Pbc & ActionAtomistic::getPbc() const { return pbc; @@ -293,6 +307,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 67b8e7ad61b3cc6f165e386d9b48d64ec21f0e92..b63c8e5d15729503975cc98d2aa62c69aca2f49a 100644 --- a/src/core/Atoms.cpp +++ b/src/core/Atoms.cpp @@ -569,4 +569,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 4aa74757a307907ae5b7cbc301a84d04ffaa3ac3..8d91c582f9196bd46095e811d9047c4579ce9abc 100644 --- a/src/core/Atoms.h +++ b/src/core/Atoms.h @@ -226,6 +226,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 963d236338d49321593097e439be25c72ed0f6df..1dcf8ee2a1f1077023dce3d6e8779da8b93f0e08 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/(2.*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 31bd861f22611f916b8abcaab60eaa2ac19380b3..7fc484f31749c7a52364f167523ef64e5a936bbd 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 9f40013cb972acb1e00b551b71bc93878bccdb45..590fd0091b80a60f98bd058b8c5cccbadeb5beaf 100644 --- a/src/core/MDAtoms.cpp +++ b/src/core/MDAtoms.cpp @@ -25,6 +25,7 @@ #include "tools/Exception.h" #include <algorithm> #include <string> +#include <map> using namespace std; @@ -47,6 +48,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); @@ -58,6 +61,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 1ec82c1fc0be38ff8a72118ec0fc63bf1d116d7b..2ed264c8f10895d2a0c4715d590032d95d9cc078 100644 --- a/src/core/MDAtoms.h +++ b/src/core/MDAtoms.h @@ -107,6 +107,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 afb83393bf4b1042274dd7ae8fb095d2f0044a82..452f13f3d29f11ebe82161b55ae39ebf664258db 100644 --- a/src/core/PlumedMain.cpp +++ b/src/core/PlumedMain.cpp @@ -443,6 +443,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");