diff --git a/regtest/rt9/plumed.dat b/regtest/rt9/plumed.dat index ca79e46ab2bf60e6a21e29e7e89c7c16492eeb0b..5e6ee5231a48323b34ffe7bde9ac7cef7fe802d1 100644 --- a/regtest/rt9/plumed.dat +++ b/regtest/rt9/plumed.dat @@ -6,7 +6,7 @@ LOWER_WALLS ARG=ang2 AT=0.2 KAPPA=40.0 EXP=4 EPS=0.3 OFFSET=0.1 LABEL=lwall PRINT ... STRIDE=10 - ARG=uwall.bias,lwall.bias,uwall.force2,lwall.force2 + ARG=*.bias,*.force2 FILE=COLVAR FMT=%6.3f ... PRINT diff --git a/src/ActionAtomistic.cpp b/src/ActionAtomistic.cpp index 60f2830649757e52d881876b94d6d1142f18573c..0a757e1e0519084b3644ab73dd6777ee344d4ee7 100644 --- a/src/ActionAtomistic.cpp +++ b/src/ActionAtomistic.cpp @@ -58,7 +58,7 @@ Vector ActionAtomistic::pbcDistance(const Vector &v1,const Vector &v2)const{ void ActionAtomistic::calculateNumericalDerivatives(){ ActionWithValue*a=dynamic_cast<ActionWithValue*>(this); plumed_massert(a,"only Actions with a value can be differentiated"); - const int nval=a->getNumberOfValues(); + const int nval=a->getNumberOfComponents(); const int natoms=getNumberOfAtoms(); std::vector<Vector> value(nval*natoms); std::vector<Tensor> valuebox(nval); @@ -70,8 +70,8 @@ void ActionAtomistic::calculateNumericalDerivatives(){ positions[i][k]=positions[i][k]+delta; calculate(); positions[i][k]=savedPositions[i][k]; - for(int j=0;j<nval;j++){ - value[j*natoms+i][k]=a->getValue(j)->get(); + for(unsigned j=0;j<nval;j++){ + value[j*natoms+i][k]=a->getOutputQuantity(j); } } for(int i=0;i<3;i++) for(int k=0;k<3;k++){ @@ -84,24 +84,24 @@ void ActionAtomistic::calculateNumericalDerivatives(){ box(i,k)=arg0; pbc.setBox(box); for(int j=0;j<natoms;j++) positions[j]=savedPositions[j]; - for(int j=0;j<nval;j++) valuebox[j](i,k)=a->getValue(j)->get(); + for(unsigned j=0;j<nval;j++) valuebox[j](i,k)=a->getOutputQuantity(j); } calculate(); - - for(int j=0;j<nval;j++){ - Value* v=a->getValue(j); + a->clearDerivatives(); + for(unsigned j=0;j<nval;j++){ + Value* v=a->copyOutput(j); double ref=v->get(); if(v->hasDerivatives()){ for(int i=0;i<natoms;i++) for(int k=0;k<3;k++) { double d=(value[j*natoms+i][k]-ref)/delta; - v->setDerivatives(3*i+k,d); + v->addDerivative(3*i+k,d); } Tensor virial; for(int i=0;i<3;i++) for(int k=0;k<3;k++)virial(i,k)= (valuebox[j](i,k)-ref)/delta; // BE CAREFUL WITH NON ORTHOROMBIC CELL virial=-1.0*matmul(box.transpose(),virial.transpose()); - for(int i=0;i<3;i++) for(int k=0;k<3;k++) v->setDerivatives(3*natoms+3*k+i,virial(i,k)); + for(int i=0;i<3;i++) for(int k=0;k<3;k++) v->addDerivative(3*natoms+3*k+i,virial(i,k)); } } } diff --git a/src/ActionWithArguments.cpp b/src/ActionWithArguments.cpp index 815a6968858d0d4ed20f1d686f228e27862912d1..020d179bd57fe7cb14f3772bb4d643a3a62351af 100644 --- a/src/ActionWithArguments.cpp +++ b/src/ActionWithArguments.cpp @@ -11,47 +11,57 @@ void ActionWithArguments::registerKeywords(Keywords& keys){ } void ActionWithArguments::parseArgumentList(const std::string&key,std::vector<Value*>&arg){ - vector<string> c; - arg.clear(); - parseVector(key,c); + vector<string> c; arg.clear(); parseVector(key,c); + for(unsigned i=0;i<c.size();i++){ - std::size_t dot=c[i].find_first_of('.'); - if(dot!=string::npos){ -// if it contains a dot: + std::size_t dot=c[i].find_first_of('.'); string a=c[i].substr(0,dot); string name=c[i].substr(dot+1); - if(a=="*"){ - plumed_massert(name=="*","arguments in the form *.something are not allowed, but for *.*"); -// this is *.*: all the values - std::vector<ActionWithValue*> all=plumed.getActionSet().select<ActionWithValue*>(); - for(unsigned j=0;j<all.size();j++){ - for(int k=0;k<all[j]->getNumberOfValues();++k) arg.push_back(all[j]->getValue(k)); - }; - } else { - ActionWithValue* action=plumed.getActionSet().selectWithLabel<ActionWithValue*>(a); - plumed_massert(action,"cannot find action named " + a); - if(name=="*"){ -// this is something.*: all the values in "something" - for(int k=0;k<action->getNumberOfValues();++k) arg.push_back(action->getValue(k)); + if(c[i].find(".")!=string::npos){ // if it contains a dot: + if(a=="*" && name=="*"){ + // Take all values from all actions + std::vector<ActionWithValue*> all=plumed.getActionSet().select<ActionWithValue*>(); + if( all.size()==0 ) error("your input file is not telling plumed to calculate anything"); + for(unsigned j=0;j<all.size();j++){ + for(int k=0;k<all[j]->getNumberOfComponents();++k) arg.push_back(all[j]->copyOutput(k)); + } + } else if ( name=="*"){ + // Take all the values from an action with a specific name + ActionWithValue* action=plumed.getActionSet().selectWithLabel<ActionWithValue*>(a); + if(!action) error("cannot find action named " + a); + for(int k=0;k<action->getNumberOfComponents();++k) arg.push_back(action->copyOutput(k)); + } else if ( a=="*" ){ + // Take components from all actions with a specific name + std::vector<ActionWithValue*> all=plumed.getActionSet().select<ActionWithValue*>(); + if( all.size()==0 ) error("your input file is not telling plumed to calculate anything"); + std::string lab; unsigned nval=0; + for(unsigned j=0;j<all.size();j++){ + std::string flab; flab=all[j]->getLabel() + "." + name; + if( all[j]->exists(flab) ){ arg.push_back(all[j]->copyOutput(flab)); nval++; } + } + if(nval==0) error("found no actions with a component called " + name ); + } else { + // Take values with a specific name + ActionWithValue* action=plumed.getActionSet().selectWithLabel<ActionWithValue*>(a); + if(!action) error("cannot find action named " + a); + if( !(action->exists(c[i])) ) error("action " + a + " has no component named " + name ); + arg.push_back(action->copyOutput(c[i])); + } + } else { // if it doesn't contain a dot + if(c[i]=="*"){ + // Take all values from all actions + std::vector<ActionWithValue*> all=plumed.getActionSet().select<ActionWithValue*>(); + if( all.size()==0 ) error("your input file is not telling plumed to calculate anything"); + for(unsigned j=0;j<all.size();j++){ + for(int k=0;k<all[j]->getNumberOfComponents();++k) arg.push_back(all[j]->copyOutput(k)); + } } else { -// this is something.x: take that component - plumed_massert(action->hasNamedValue(name),"action " + a + " has no component with name "+ name); - arg.push_back(action->getValue(name)); + ActionWithValue* action=plumed.getActionSet().selectWithLabel<ActionWithValue*>(c[i]); + if(!action) error("cannot find action named " + c[i]); + if( !(action->exists(c[i])) ) error("action " + c[i] + " has no component named " + c[i] ); + arg.push_back(action->copyOutput(c[i])); } } - } else if(c[i]=="*"){ -// this is *: all the values (same as *.*) - std::vector<ActionWithValue*> all=plumed.getActionSet().select<ActionWithValue*>(); - for(unsigned j=0;j<all.size();j++){ - for(int k=0;k<all[j]->getNumberOfValues();++k) arg.push_back(all[j]->getValue(k)); - }; - } else{ -// this is something: take the unnamed component - ActionWithValue* action=plumed.getActionSet().selectWithLabel<ActionWithValue*>(c[i]); - plumed_massert(action,"cannot find action named " +c[i]); - plumed_massert(action->hasNamedValue(""),"action "+c[i]+" has no default component (unnamed one)"); - arg.push_back(action->getValue("")); - } } } @@ -59,7 +69,19 @@ void ActionWithArguments::requestArguments(const vector<Value*> &arg){ plumed_massert(!lockRequestArguments,"requested argument list can only be changed in the prepare() method"); arguments=arg; clearDependencies(); - for(unsigned i=0;i<arguments.size();i++) addDependency(&arguments[i]->getAction()); + std::string fullname,name; + for(unsigned i=0;i<arguments.size();i++){ + fullname=arguments[i]->getName(); + if(fullname.find(".")!=string::npos){ + std::size_t dot=fullname.find_first_of('.'); + name=fullname.substr(0,dot); + } else { + name=fullname; + } + ActionWithValue* action=plumed.getActionSet().selectWithLabel<ActionWithValue*>(name); + if(!action) error("cannot find action named (in requestArguments - this is weird)" + name); + addDependency(action); + } } ActionWithArguments::ActionWithArguments(const ActionOptions&ao): @@ -71,7 +93,7 @@ ActionWithArguments::ActionWithArguments(const ActionOptions&ao): if(arg.size()>0){ log.printf(" with arguments"); - for(unsigned i=0;i<arg.size();i++) log.printf(" %s",arg[i]->getFullName().c_str()); + for(unsigned i=0;i<arg.size();i++) log.printf(" %s",arg[i]->getName().c_str()); log.printf("\n"); } @@ -82,23 +104,24 @@ ActionWithArguments::ActionWithArguments(const ActionOptions&ao): void ActionWithArguments::calculateNumericalDerivatives(){ ActionWithValue*a=dynamic_cast<ActionWithValue*>(this); plumed_massert(a,"cannot compute numerical derivatives for an action without values"); - const int nval=a->getNumberOfValues(); - const int npar=a->getNumberOfParameters(); + const int nval=a->getNumberOfComponents(); + const int npar=arguments.size(); std::vector<double> value (nval*npar); for(int i=0;i<npar;i++){ double arg0=arguments[i]->get(); arguments[i]->set(arg0+sqrt(epsilon)); calculate(); arguments[i]->set(arg0); - for(int j=0;j<nval;j++){ - value[i*nval+j]=a->getValue(j)->get(); + for(unsigned j=0;j<nval;j++){ + value[i*nval+j]=a->getOutputQuantity(j); } } calculate(); + a->clearDerivatives(); std::vector<double> value0(nval); - for(int j=0;j<nval;j++){ - Value* v=a->getValue(j); - if(v->hasDerivatives())for(int i=0;i<npar;i++) v->setDerivatives(i,(value[i*nval+j]-a->getValue(j)->get())/sqrt(epsilon)); + for(unsigned j=0;j<nval;j++){ + Value* v=a->copyOutput(j); + if( v->hasDerivatives() ) for(int i=0;i<npar;i++) v->addDerivative(i,(value[i*nval+j]-a->getOutputQuantity(j))/sqrt(epsilon)); } } diff --git a/src/ActionWithArguments.h b/src/ActionWithArguments.h index d97491daa1b305cbc8d47e84e9b5def7eed30792..c4da9fd3438ba9d4b7e5fc5fb884236b85ec33c9 100644 --- a/src/ActionWithArguments.h +++ b/src/ActionWithArguments.h @@ -24,40 +24,43 @@ class ActionWithArguments: bool lockRequestArguments; protected: - ActionWithArguments(const ActionOptions&); - virtual ~ActionWithArguments(){}; /// double getProjection(unsigned i,unsigned j)const; public: /// Registers the list of keywords static void registerKeywords( Keywords& keys ); - -/// Returns an array of pointers to the arguments - std::vector<Value*> & getArguments(); /// Returns the value of an argument - double getArgument(int)const; + double getArgument( const unsigned n ) const; +/// Return a pointer to specific argument + Value* getPntrToArgument( const unsigned n ); /// Returns the number of arguments - unsigned getNumberOfArguments()const; -/// - void calculateNumericalDerivatives(); + unsigned getNumberOfArguments() const ; /// Takes the difference taking into account pbc for arg i - double difference(int,double,double)const; + double difference(int, double, double) const; /// Parse a list of arguments - void parseArgumentList(const std::string&key,std::vector<Value*>&args); - void requestArguments(const std::vector<Value*> &arg); + void parseArgumentList(const std::string&key,std::vector<Value*>&args); +/// Setup the dependencies + void requestArguments(const std::vector<Value*> &arg); +public: + ActionWithArguments(const ActionOptions&); + virtual ~ActionWithArguments(){}; +/// Registers the list of keywords + static void registerKeywords( Keywords& keys ); +/// Calculate the numerical derivatives + void calculateNumericalDerivatives(); void lockRequests(); void unlockRequests(); }; inline -std::vector<Value*> & ActionWithArguments::getArguments(){ - return arguments; +Value* ActionWithArguments::getPntrToArgument( const unsigned n ){ + return arguments[n]; } inline -double ActionWithArguments::getArgument(int i)const{ - return arguments[i]->get(); +double ActionWithArguments::getArgument(const unsigned n) const { + return arguments[n]->get(); } inline diff --git a/src/ActionWithValue.cpp b/src/ActionWithValue.cpp index 4b22c1f1c353e8d6177ec040d343c19a9a481040..7444bc1b3f5db7842467e6236af30460f7eeb3ba 100644 --- a/src/ActionWithValue.cpp +++ b/src/ActionWithValue.cpp @@ -16,10 +16,7 @@ void ActionWithValue::noAnalyticalDerivatives(Keywords& keys){ ActionWithValue::ActionWithValue(const ActionOptions&ao): Action(ao), - numberOfParameters(0), - numericalDerivatives(false), - hasMultipleValues(false), - hasUnnamedValue(false) + numericalDerivatives(false) { if( keywords.exists("NUMERICAL_DERIVATIVES") ) parseFlag("NUMERICAL_DERIVATIVES",numericalDerivatives); if(numericalDerivatives) log.printf(" using numerical derivatives\n"); @@ -29,66 +26,101 @@ ActionWithValue::~ActionWithValue(){ for(unsigned i=0;i<values.size();++i)delete values[i]; } -void ActionWithValue::check(const std::string&name){ - assertUnique(name); - if(name==""){ - hasUnnamedValue=true; - plumed_massert(!hasMultipleValues,"cannot use unnamed components for a multicomponent Action"); - }else{ - hasMultipleValues=true; - plumed_massert(!hasUnnamedValue,"cannot use unnamed components for a multicomponent Action"); +void ActionWithValue::clearInputForces(){ + for(unsigned i=0;i<values.size();i++) values[i]->clearInputForce(); +} +void ActionWithValue::clearDerivatives(){ + for(unsigned i=0;i<values.size();i++) values[i]->clearDerivatives(); +} + +// -- These are the routine for copying the value pointers to other classes -- // + +bool ActionWithValue::exists( const std::string& name ) const { + for(unsigned i=0;i<values.size();++i){ + if (values[i]->name==name) return true; } + return false; } -void ActionWithValue::addValue(const std::string&name){ - check(name); - values.push_back(new Value(*this,name)); +Value* ActionWithValue::copyOutput( const std::string& name ) const { + for(unsigned i=0;i<values.size();++i){ + if (values[i]->name==name) return values[i]; + } + plumed_massert(0,"there is no pointer with name " + name); } -void ActionWithValue::addValueWithDerivatives(const std::string&name){ - check(name); - Value* v=new Value(*this,name); - v->enableDerivatives(); - values.push_back(v); +Value* ActionWithValue::copyOutput( const unsigned& n ) const { + plumed_massert(n<values.size(),"you have requested a pointer that is out of bounds"); + return values[n]; } -bool ActionWithValue::hasNamedValue(const std::string&name)const{ - for(unsigned i=0;i<values.size();++i){ - if(name==values[i]->getName()) return true; - } - return false; +// -- HERE WE HAVE THE STUFF FOR THE DEFAULT VALUE -- // + +void ActionWithValue::addValue(){ + plumed_massert(values.size()==0,"You have already added the default value for this action"); + values.push_back(new Value(getLabel(), false ) ); } -int ActionWithValue::getValueIndex(const std::string&name)const{ - for(unsigned i=0;i<values.size();++i) if(name==values[i]->getName()) return i; - plumed_merror("value not found" + name); - return -1; // otherwise the compiler complains +void ActionWithValue::addValueWithDerivatives(){ + plumed_massert(values.size()==0,"You have already added the default value for this action"); + values.push_back(new Value(getLabel(), true ) ); } -Value* ActionWithValue::getValue(const std::string&name)const{ - return values[getValueIndex(name)]; +void ActionWithValue::setNotPeriodic(){ + plumed_massert(values.size()==1,"The number of components is not equal to one"); + plumed_massert(values[0]->name==getLabel(), "The value you are trying to set is not the default"); + values[0]->min=0; values[0]->max=0; + values[0]->setupPeriodicity(); } -Value* ActionWithValue::getValue(int i)const{ - return values[i]; +void ActionWithValue::setPeriodic( const double min, const double max ){ + plumed_massert(values.size()==1,"The number of components is not equal to one"); + plumed_massert(values[0]->name==getLabel(), "The value you are trying to set is not the default"); + values[0]->min=min; values[0]->max=max; + values[0]->setupPeriodicity(); } -std::vector<std::string> ActionWithValue::getValueNames()const{ - std::vector<std::string> ret; - for(unsigned i=0;i<values.size();i++) ret.push_back(values[i]->getName()); - return ret; +Value* ActionWithValue::getPntrToValue(){ + plumed_massert(values.size()==1,"The number of components is not equal to one"); + plumed_massert(values[0]->name==getLabel(), "The value you are trying to retrieve is not the default"); + return values[0]; } -void ActionWithValue::setNumberOfParameters(int n){ - numberOfParameters=n; - for(unsigned i=0;i<values.size();i++) values[i]->setNumberOfParameters(n); +// -- HERE WE HAVE THE STUFF FOR NAMED VALUES / COMPONENTS -- // + +void ActionWithValue::addComponent( const std::string& name ){ + std::string thename; thename=getLabel() + "." + name; + for(unsigned i=0;i<values.size();++i){ + plumed_massert(values[i]->name!=getLabel(),"Cannot mix single values with components"); + plumed_massert(values[i]->name!=thename,"there is already a value with this name"); + } + values.push_back(new Value(thename, false ) ); } -void ActionWithValue::clearInputForces(){ - for(unsigned i=0;i<values.size();i++) values[i]->clearInputForce(); +void ActionWithValue::addComponentWithDerivatives( const std::string& name ){ + std::string thename; thename=getLabel() + "." + name; + for(unsigned i=0;i<values.size();++i){ + plumed_massert(values[i]->name!=getLabel(),"Cannot mix single values with components"); + plumed_massert(values[i]->name!=thename,"there is already a value with this name"); + } + values.push_back(new Value(thename, true ) ); } -void ActionWithValue::clearDerivatives(){ - for(unsigned i=0;i<values.size();i++) values[i]->clearDerivatives(); + +int ActionWithValue::getComponent( const std::string& name ) const { + plumed_massert( !exists( getLabel() ), "You should not be calling this routine if you are using a value"); + plumed_massert(name!=getLabel(),"You should never be calling this routine to retrieve the value"); + std::string thename; thename=getLabel() + "." + name; + for(unsigned i=0;i<values.size();++i){ + if (values[i]->name==thename) return i; + } + plumed_massert(0,"there is no component with name" + name); + return -1; +} + +void ActionWithValue::componentIsNotPeriodic( const std::string& name ){ + int kk=getComponent(name); + values[kk]->min=0; values[kk]->max=0; + values[kk]->setupPeriodicity(); } void ActionWithValue::setGradientsIfNeeded(){ @@ -97,3 +129,18 @@ void ActionWithValue::setGradientsIfNeeded(){ } } +void ActionWithValue::componentIsPeriodic( const std::string& name, const double min, const double max ){ + int kk=getComponent(name); + values[kk]->min=min; values[kk]->max=max; + values[kk]->setupPeriodicity(); +} + +Value* ActionWithValue::getPntrToComponent( const std::string& name ){ + int kk=getComponent(name); + return values[kk]; +} + +Value* ActionWithValue::getPntrToComponent( int n ){ + plumed_massert(n<values.size(),"you have requested a pointer that is out of bounds"); + return values[n]; +} diff --git a/src/ActionWithValue.h b/src/ActionWithValue.h index 5c8f721ef52e46f619ce5be10cca680798cbf064..0be87217f12b7da784429eb8c666ece8fb7447e7 100644 --- a/src/ActionWithValue.h +++ b/src/ActionWithValue.h @@ -15,22 +15,74 @@ This is used for PLMD::Bias, PLMD::Colvar and PLMD::Function */ //+ENDDEVELDOC -/// Action which can take one or more values. -/// This object contains an array of PLMD::Value, one for each component. -/// It also stores all the derivatives of these values wrt the parameters -/// Parameters are other values (from other Action s) or atomic positions. -class ActionWithValue: +/// The vast majority of the PLMD::Action objects that are implemented in +/// plumed calculate some quantity or a set of quantities. This could be +/// the value of a CV, the value of a function or the potential due to a bias. +/// PLMD::ActionWithValue provides the functionality for storing these quantities +/// and (in tandem with PLMD::ActionWithArguments) the functionality for passing +/// quantities between PLMD::Actions. When you are deciding what quantities +/// your new PLMD::Action will need to store using PLMD::ActionWithValue you must +/// ask yourself the following two questions: + +/// - Do I need to differentiate my output quantities +/// - Is my PLMD::Action calculating a single thing or does the output have multiple components + +/// If the answer to the first of these questions is yes then you must setup your values +/// you using either PLMD::ActionWithValue::addValueWithDerivatives() or +/// PLMD::ActionWithValue::addComponentWithDerivatives. If the answer is no you +/// can set up values using PLMD::ActionWithValue::addValue() or PLMD::ActionWithValue::addComponent(). +/// The precise routine you use to setup your values will depend on your answer to the +/// second question. As you are probably aware if the output of your PLMD::Action is a +/// single quantity you can reference that quantity in the input file using the label of the +/// PLMD::Action it was calculated in. If your action <b> outputs only one quantity </b> +/// we call that quantity the <b> value </b> of the Action. To set the <b> value </b> and get pointers to it +/// you should <b> use the set of routines that have the word value in the name </b>. If, by contrast, +/// your PLMD::Action calculates multiple quantities then these quantities are referenced in input using the +/// label.component syntax. We refer to these <b> multiple quantities </b> the <b> components </b> +/// of the PLMD::Action. Perhaps unsurprisingly, when you manipulate the <b> components </b> of an +/// PLMD::Action you should use <b> the routines with the word component in the name. </b> + +class ActionWithValue : public virtual Action { - int numberOfParameters; +private: +/// An array containing the values for this action std::vector<Value*> values; - void assertUnique(const std::string&name); - void check(const std::string&name); - int getValueIndex(const std::string&name)const; +/// Are we using numerical derivatives to differentiate bool numericalDerivatives; - bool hasMultipleValues; - bool hasUnnamedValue; +/// Return the index for the component named name + int getComponent( const std::string& name ) const; protected: + +// -------- The action has one value only ---------------- // + +/// Add a value with the name <label> + void addValue(); +/// Add a value with the name <label> that has derivatives + void addValueWithDerivatives(); +/// Set your default value to have no periodicity + void setNotPeriodic(); +/// Set the value to be periodic with a particular domain + void setPeriodic( const double min, const double max ); +/// Set the value of quantity with name <label> + void setValue( const double& f); +/// Get a pointer to the default value + Value* getPntrToValue(); + +// -------- The action has multiple components ---------- // + +/// Add a value with a name like <label>.name + void addComponent( const std::string& name ); +/// Add a value with a name like <label>.name that has derivatives + void addComponentWithDerivatives( const std::string& name ); +/// Set your value component to have no periodicity + void componentIsNotPeriodic( const std::string& name ); +/// Set the value to be periodic with a particular domain + void componentIsPeriodic( const std::string& name, const double min, const double max ); +/// Return a pointer to the component by index + Value* getPntrToComponent(int i); +/// Return a pointer to the value by name + Value* getPntrToComponent(const std::string& name); public: ActionWithValue(const ActionOptions&ao); ~ActionWithValue(); @@ -39,36 +91,29 @@ public: static void registerKeywords( Keywords& keys ); /// Insist that numerical derivatives should always be used for an action and make this fact appear in the manual static void noAnalyticalDerivatives(Keywords& keys); +/// Get the value of one of the components of the PLMD::Action + double getOutputQuantity( const unsigned j ) const ; +/// Get the value with a specific name (N.B. if there is no such value this returns zero) + double getOutputQuantity( const std::string& name ) const ; + +// --- Routines for passing stuff to ActionWithArguments -- // + +/// Check if a value with a particular name is present. This is only used in PLMD::ActionWithArguments. +/// You should not use it when manipulating components. + bool exists( const std::string& name ) const; +/// Return a pointer to the value with name (this is used to retrieve values in other PLMD::Actions) +/// You should NEVER use this routine to refer to the components of your PLMD::Action. Use +/// getPntrToComponent instead. + Value* copyOutput( const std::string&name ) const; +/// Return a pointer to the value with this number (this is used to retrieve values in other PLMD::Actions) +/// You should NEVER use this routine to refer to the components of your PLMD::Action. Use +/// getPntrToComponent instead. + Value* copyOutput( const unsigned& n ) const; + +// -- Routines for everything else -- // -/// Add a new value without derivatives. -/// This should be used for values which are only evaluated (e.g. for printing) -/// but for which we do not make derivatives available so that forces cannot -/// be applied - void addValue(const std::string&name); -/// Add a new value with derivatives. -/// This should be used for values for which we make derivatives available -/// so that forces can be applied - void addValueWithDerivatives(const std::string&name); -/// Check if a value with a given name is already used - bool hasNamedValue(const std::string&name)const; -/// Return a pointer to the value by name - Value* getValue(const std::string&name)const; -/// Return a pointer to the value by index -/// This should be an indexes growing for new inserted values. -/// E.g., the default value (no name) is number 0, ... - Value* getValue(int i)const; -/// Returns an array of strings with the names of the values - std::vector<std::string> getValueNames()const; /// Returns the number of values defined - int getNumberOfValues(); -/// Set the number of parameters on which this Action depends. -/// Example: for a Bias, this is the number of arguments, for a Colvar -/// is 3*Natoms+cell variables - void setNumberOfParameters(int n); -/// Returns the number of parameters on which this Action depends. - int getNumberOfParameters()const; -/// Returns the total force applied on i-th value - double getForce(int n); + int getNumberOfComponents() const ; /// Clear the forces on the values void clearInputForces(); /// Clear the derivatives of values wrt parameters @@ -80,41 +125,40 @@ public: /// Set the default value (the one without name) void setValue(double); /// Check if numerical derivatives should be used - bool checkNumericalDerivatives()const; + bool checkNumericalDerivatives() const ; }; inline -void ActionWithValue::setValue(Value*v,double d){ - v->set(d); +double ActionWithValue::getOutputQuantity(const unsigned j) const { + plumed_massert(j<values.size(),"index requested is out of bounds"); + return values[j]->get(); } inline -void ActionWithValue::setValue(double d){ - values[0]->set(d); +double ActionWithValue::getOutputQuantity( const std::string& name ) const { + std::string thename; thename=getLabel() + "." + name; + if( exists(thename) ){ + for(unsigned i=0;i<values.size();++i){ + if( values[i]->name==thename ) return values[i]->value; + } + } + return 0.0; } inline -double ActionWithValue::getForce(int n){ - return values[n]->getForce(); -} - -inline -void ActionWithValue::assertUnique(const std::string&name){ - plumed_massert(!hasNamedValue(name),"action " + getLabel() + " already has a value " + name); +void ActionWithValue::setValue(const double& d){ + plumed_massert(values.size()==1, "cannot use setValue in multi-component actions"); + plumed_massert(values[0]->name==getLabel(), "The value you are trying to set is not the default"); + values[0]->set(d); } inline -int ActionWithValue::getNumberOfValues(){ +int ActionWithValue::getNumberOfComponents() const { return values.size(); } inline -int ActionWithValue::getNumberOfParameters()const{ - return numberOfParameters; -} - -inline -bool ActionWithValue::checkNumericalDerivatives()const{ +bool ActionWithValue::checkNumericalDerivatives() const { return numericalDerivatives; } diff --git a/src/Bias.cpp b/src/Bias.cpp index e802b9bd16705c9edca3e8c29ac71b831f099929..f06e87abd2852acdec697108d0d12341d04baaa1 100644 --- a/src/Bias.cpp +++ b/src/Bias.cpp @@ -9,7 +9,7 @@ Action(ao), ActionPilot(ao), ActionWithValue(ao), ActionWithArguments(ao), -outputForces(getArguments().size(),0.0) +outputForces(getNumberOfArguments(),0.0) { } @@ -24,7 +24,7 @@ void Bias::registerKeywords( Keywords& keys ){ void Bias::apply(){ if(onStep()) for(unsigned i=0;i<getNumberOfArguments();++i){ - getArguments()[i]->addForce(getStride()*outputForces[i]); + getPntrToArgument(i)->addForce(getStride()*outputForces[i]); } } diff --git a/src/BiasExternal.cpp b/src/BiasExternal.cpp index ddc08c9c985a3fe570e4e0513a4d31d972bbdd44..b20aa37fdbc9d6e4b12ce4a0a9a0e699253f5af4 100644 --- a/src/BiasExternal.cpp +++ b/src/BiasExternal.cpp @@ -70,15 +70,15 @@ BiasGrid_(NULL) if(spline){log.printf(" External potential uses spline interpolation\n");} if(sparsegrid){log.printf(" External potential uses sparse grid\n");} - addValue("bias"); + addComponent("bias"); // read grid FILE* gridfile=fopen(filename.c_str(),"r"); BiasGrid_=Grid::create(gridfile,sparsegrid,spline,true); fclose(gridfile); - plumed_assert(BiasGrid_->getDimension()==getNumberOfArguments()); + if(BiasGrid_->getDimension()!=getNumberOfArguments()) error("mismatch between dimensionality of input grid and number of arguments"); for(unsigned i=0;i<getNumberOfArguments();++i){ - plumed_assert(getArguments()[i]->isPeriodic()==BiasGrid_->getIsPeriodic()[i]); + if( getPntrToArgument(i)->isPeriodic()!=BiasGrid_->getIsPeriodic()[i] ) error("periodicity mismatch between arguments and input bias"); } } @@ -91,8 +91,7 @@ void BiasExternal::calculate() double ene=BiasGrid_->getValueAndDerivatives(cv,der); - Value* value=getValue("bias"); - setValue(value,ene); + getPntrToComponent("bias")->set(ene); // set Forces for(unsigned i=0;i<ncv;++i){ diff --git a/src/BiasLWalls.cpp b/src/BiasLWalls.cpp index fc3d55ea914cd3190b06d4d21297589aeb8abdd5..dae30cc601be210c894f9cdb2b983a4e475d6185 100644 --- a/src/BiasLWalls.cpp +++ b/src/BiasLWalls.cpp @@ -97,8 +97,8 @@ offset(getNumberOfArguments(),0.0) for(unsigned i=0;i<eps.size();i++) log.printf(" %f",eps[i]); log.printf("\n"); - addValue("bias"); - addValue("force2"); + addComponent("bias"); + addComponent("force2"); } void BiasLWalls::calculate(){ @@ -118,9 +118,8 @@ void BiasLWalls::calculate(){ totf2+=f*f; } }; - Value* value; - value=getValue("bias"); setValue(value,ene); - value=getValue("force2"); setValue(value,totf2); + getPntrToComponent("bias")->set(ene); + getPntrToComponent("force2")->set(totf2); } } diff --git a/src/BiasMetaD.cpp b/src/BiasMetaD.cpp index a2314f882fef6d249663499b815b7f874cbe5e79..cd442db20311e13fdc72cb4d10f804049b4007f9 100644 --- a/src/BiasMetaD.cpp +++ b/src/BiasMetaD.cpp @@ -218,7 +218,7 @@ grid_(false) if(wgridstride_>0){log.printf(" Grid is written on file %s with stride %d\n",gridfname.c_str(),wgridstride_);} } - addValue("bias"); + addComponent("bias"); // for performance dp_ = new double[getNumberOfArguments()]; @@ -227,11 +227,11 @@ grid_(false) if(grid_){ vector<bool> pbc; for(unsigned i=0;i<getNumberOfArguments();++i){ - pbc.push_back(getArguments()[i]->isPeriodic()); + pbc.push_back(getPntrToArgument(i)->isPeriodic()); // if periodic, use CV domain for grid boundaries if(pbc[i]){ double dmin,dmax; - getArguments()[i]->getDomain(dmin,dmax); + getPntrToArgument(i)->getDomain(dmin,dmax); gmin[i]=dmin; gmax[i]=dmax; } @@ -309,7 +309,7 @@ vector<unsigned> BiasMetaD::getGaussianSupport(const Gaussian& hill) vector<unsigned> nneigh; for(unsigned i=0;i<getNumberOfArguments();++i){ double cutoff=sqrt(2.0*DP2CUTOFF)*hill.sigma[i]; - nneigh.push_back((unsigned)ceil(cutoff/BiasGrid_->getDx()[i])); + nneigh.push_back( static_cast<unsigned>(ceil(cutoff/BiasGrid_->getDx()[i])) ); } return nneigh; } @@ -372,9 +372,7 @@ void BiasMetaD::calculate() double* der=new double[ncv]; for(unsigned i=0;i<ncv;++i){der[i]=0.0;} double ene=getBiasAndDerivatives(cv,der); - - Value* value=getValue("bias"); - setValue(value,ene); + getPntrToComponent("bias")->set(ene); // set Forces for(unsigned i=0;i<ncv;++i){ diff --git a/src/BiasMovingRestraint.cpp b/src/BiasMovingRestraint.cpp index cf9b13f5a54e8fb3f58d409fb3d0fecec1669849..8b5fcff9f9caf529faf3d1c967264efd2d801745 100644 --- a/src/BiasMovingRestraint.cpp +++ b/src/BiasMovingRestraint.cpp @@ -131,9 +131,9 @@ verse(getNumberOfArguments()) log.printf("\n"); }; - addValue("bias"); - addValue("force2"); - addValue("work"); + addComponent("bias"); + addComponent("force2"); + addComponent("work"); } @@ -168,11 +168,8 @@ void BiasMovingRestraint::calculate(){ setOutputForces(i,f); totf2+=f*f; }; - Value* value; - value=getValue("bias"); - setValue(value,ene); - value=getValue("force2"); - setValue(value,totf2); + getPntrToComponent("bias")->set(ene); + getPntrToComponent("force2")->set(totf2); } } diff --git a/src/BiasRatchet.cpp b/src/BiasRatchet.cpp index cdc57175763aea9e49de0b6af3ef9aa969d197e5..8480557cb3a77fea1768244d7d49289206f07300 100644 --- a/src/BiasRatchet.cpp +++ b/src/BiasRatchet.cpp @@ -87,9 +87,9 @@ kappa(getNumberOfArguments(),0.0) for(unsigned i=0;i<kappa.size();i++) log.printf(" %f",kappa[i]); log.printf("\n"); - for(unsigned i=0;i<getNumberOfArguments();i++) {char str_min[6]; sprintf(str_min,"min_%u",i+1); addValue(str_min);} - addValue("bias"); - addValue("force2"); + for(unsigned i=0;i<getNumberOfArguments();i++) {char str_min[6]; sprintf(str_min,"min_%u",i+1); addComponent(str_min);} + addComponent("bias"); + addComponent("force2"); } @@ -111,11 +111,10 @@ void BiasRatchet::calculate(){ } char str_min[6]; sprintf(str_min,"min_%u",i+1); - value=getValue(str_min); - setValue(value,min[i]); + getPntrToComponent(str_min)->set(min[i]); }; - value=getValue("bias"); setValue(value,ene); - value=getValue("force2"); setValue(value,totf2); + getPntrToComponent("bias")->set(ene); + getPntrToComponent("force2")->set(totf2); } } diff --git a/src/BiasRestraint.cpp b/src/BiasRestraint.cpp index 074dbaeb7c8dbab3e8fb213f33d48e175001b8e1..48f025607dd8bef762a66cdbd562585500201047 100644 --- a/src/BiasRestraint.cpp +++ b/src/BiasRestraint.cpp @@ -80,8 +80,8 @@ slope(getNumberOfArguments(),0.0) for(unsigned i=0;i<slope.size();i++) log.printf(" %f",slope[i]); log.printf("\n"); - addValue("bias"); - addValue("force2"); + addComponent("bias"); + addComponent("force2"); } @@ -97,9 +97,8 @@ void BiasRestraint::calculate(){ setOutputForces(i,f); totf2+=f*f; }; - Value* value; - value=getValue("bias"); setValue(value,ene); - value=getValue("force2"); setValue(value,totf2); + getPntrToComponent("bias")->set(ene); + getPntrToComponent("force2")->set(totf2); } } diff --git a/src/BiasUWalls.cpp b/src/BiasUWalls.cpp index 287764121891d93e32f1cf244c18f00010413f22..41d6a45d5cbd86ef5527d9e7973cd8558ad0a7de 100644 --- a/src/BiasUWalls.cpp +++ b/src/BiasUWalls.cpp @@ -97,8 +97,8 @@ offset(getNumberOfArguments(),0.0) for(unsigned i=0;i<eps.size();i++) log.printf(" %f",eps[i]); log.printf("\n"); - addValue("bias"); - addValue("force2"); + addComponent("bias"); + addComponent("force2"); } void BiasUWalls::calculate(){ @@ -118,9 +118,8 @@ void BiasUWalls::calculate(){ totf2+=f*f; } }; - Value* value; - value=getValue("bias"); setValue(value,ene); - value=getValue("force2"); setValue(value,totf2); + getPntrToComponent("bias")->set(ene); + getPntrToComponent("force2")->set(totf2); } } diff --git a/src/Colvar.cpp b/src/Colvar.cpp index 28985fc22d724e2e5d8ae125421fda8222d8b5e0..276af58d7cbaaf6335c68d96cd209f0066f40105 100644 --- a/src/Colvar.cpp +++ b/src/Colvar.cpp @@ -23,13 +23,19 @@ void Colvar::registerKeywords( Keywords& keys ){ } void Colvar::requestAtoms(const vector<AtomNumber> & a){ + plumed_massert(!isEnergy,"request atoms should not be called if this is energy"); +// Tell actionAtomistic what atoms we are getting ActionAtomistic::requestAtoms(a); - setNumberOfParameters(3*a.size()+9); +// Resize the derivatives of all atoms + for(int i=0;i<getNumberOfComponents();++i) getPntrToComponent(i)->resizeDerivatives(3*a.size()+9); +// Set the size of the forces array + forces.resize(3*getNumberOfAtoms()+9); } void Colvar::apply(){ vector<Vector>& f(modifyForces()); Tensor& v(modifyVirial()); + unsigned nat=getNumberOfAtoms(); for(unsigned i=0;i<f.size();i++){ f[i][0]=0.0; @@ -38,31 +44,27 @@ void Colvar::apply(){ } v.clear(); - if(!isEnergy) - for(int i=0;i<getNumberOfValues();++i){ - if(!getValue(i)->checkForced())continue; - const vector<double> & derivatives(getValue(i)->getDerivatives()); - const unsigned nat=f.size(); - const double force=getValue(i)->getForce(); - for(unsigned j=0;j<nat;++j){ - f[j][0]+=force*derivatives[3*j+0]; - f[j][1]+=force*derivatives[3*j+1]; - f[j][2]+=force*derivatives[3*j+2]; + if(!isEnergy){ + for(int i=0;i<getNumberOfComponents();++i){ + if( getPntrToComponent(i)->applyForce( forces ) ){ + for(unsigned j=0;j<nat;++j){ + f[j][0]+=forces[3*j+0]; + f[j][1]+=forces[3*j+1]; + f[j][2]+=forces[3*j+2]; + } + v(0,0)+=forces[3*nat+0]; + v(0,1)+=forces[3*nat+1]; + v(0,2)+=forces[3*nat+2]; + v(1,0)+=forces[3*nat+3]; + v(1,1)+=forces[3*nat+4]; + v(1,2)+=forces[3*nat+5]; + v(2,0)+=forces[3*nat+6]; + v(2,1)+=forces[3*nat+7]; + v(2,2)+=forces[3*nat+8]; } - v(0,0)+=force*derivatives[3*nat+0]; - v(0,1)+=force*derivatives[3*nat+1]; - v(0,2)+=force*derivatives[3*nat+2]; - v(1,0)+=force*derivatives[3*nat+3]; - v(1,1)+=force*derivatives[3*nat+4]; - v(1,2)+=force*derivatives[3*nat+5]; - v(2,0)+=force*derivatives[3*nat+6]; - v(2,1)+=force*derivatives[3*nat+7]; - v(2,2)+=force*derivatives[3*nat+8]; + } + } else if( isEnergy ){ + forces.resize(1); + if( getPntrToComponent(0)->applyForce( forces ) ) modifyForceOnEnergy()+=forces[0]; } - if(isEnergy) modifyForceOnEnergy()+=getValue(0)->getForce(); } - - - - - diff --git a/src/Colvar.h b/src/Colvar.h index 60079ac79a79c8dce26b95e2928a39a896d1c9fe..799336bc2689610cab776b06a2cd2d50f4d5a3f2 100644 --- a/src/Colvar.h +++ b/src/Colvar.h @@ -20,6 +20,9 @@ class Colvar : public ActionAtomistic, public ActionWithValue { +private: +/// This is used by apply to retrive the forces on the atoms + std::vector<double> forces; protected: bool isEnergy; void requestAtoms(const std::vector<AtomNumber> & a); @@ -44,34 +47,34 @@ public: inline void Colvar::setAtomsDerivatives(Value*v,int i,const Vector&d){ - v->setDerivatives(3*i+0,d[0]); - v->setDerivatives(3*i+1,d[1]); - v->setDerivatives(3*i+2,d[2]); + v->addDerivative(3*i+0,d[0]); + v->addDerivative(3*i+1,d[1]); + v->addDerivative(3*i+2,d[2]); } inline void Colvar::setBoxDerivatives(Value* v,const Tensor&d){ unsigned nat=getNumberOfAtoms(); - v->setDerivatives(3*nat+0,d(0,0)); - v->setDerivatives(3*nat+1,d(0,1)); - v->setDerivatives(3*nat+2,d(0,2)); - v->setDerivatives(3*nat+3,d(1,0)); - v->setDerivatives(3*nat+4,d(1,1)); - v->setDerivatives(3*nat+5,d(1,2)); - v->setDerivatives(3*nat+6,d(2,0)); - v->setDerivatives(3*nat+7,d(2,1)); - v->setDerivatives(3*nat+8,d(2,2)); + v->addDerivative(3*nat+0,d(0,0)); + v->addDerivative(3*nat+1,d(0,1)); + v->addDerivative(3*nat+2,d(0,2)); + v->addDerivative(3*nat+3,d(1,0)); + v->addDerivative(3*nat+4,d(1,1)); + v->addDerivative(3*nat+5,d(1,2)); + v->addDerivative(3*nat+6,d(2,0)); + v->addDerivative(3*nat+7,d(2,1)); + v->addDerivative(3*nat+8,d(2,2)); } inline void Colvar::setAtomsDerivatives(int i,const Vector&d){ - setAtomsDerivatives(getValue(0),i,d); + setAtomsDerivatives(getPntrToValue(),i,d); } inline void Colvar::setBoxDerivatives(const Tensor&d){ - setBoxDerivatives(getValue(0),d); + setBoxDerivatives(getPntrToValue(),d); } } diff --git a/src/ColvarAngle.cpp b/src/ColvarAngle.cpp index 0c2d358b25d10e92b33ec81a5ffc6e8c209b1d51..1ebc5f12fa2f5bac47df4dc86347a9d7e3baeb0d 100644 --- a/src/ColvarAngle.cpp +++ b/src/ColvarAngle.cpp @@ -74,9 +74,7 @@ pbc(true) if(pbc) log.printf(" using periodic boundary conditions\n"); else log.printf(" without periodic boundary conditions\n"); - addValueWithDerivatives(""); - getValue("")->setPeriodicity(false); - + addValueWithDerivatives(); setNotPeriodic(); requestAtoms(atoms); } diff --git a/src/ColvarCoordination.cpp b/src/ColvarCoordination.cpp index e7239645980f1181a0a8e3c0ef787d8b11bb3721..a9e782b3d4213c39ce25c6f10d6a43dd06ed218f 100644 --- a/src/ColvarCoordination.cpp +++ b/src/ColvarCoordination.cpp @@ -114,9 +114,7 @@ reduceListAtNextStep(false) checkRead(); - addValueWithDerivatives(""); - getValue("")->setPeriodicity(false); - + addValueWithDerivatives(); setNotPeriodic(); if(doneigh) nl= new NeighborList(ga_lista,gb_lista,dopair,pbc,getPbc(),nl_cut,nl_st); else nl= new NeighborList(ga_lista,gb_lista,dopair,pbc,getPbc()); diff --git a/src/ColvarDipole.cpp b/src/ColvarDipole.cpp index 43a9122868675b275ff5cb94de82e661183fcb5b..94b6182a6f4d4ed59a000097637f85e800de9454 100644 --- a/src/ColvarDipole.cpp +++ b/src/ColvarDipole.cpp @@ -43,8 +43,7 @@ PLUMED_COLVAR_INIT(ao) { parseAtomList("GROUP",ga_lista); checkRead(); - addValueWithDerivatives(""); - getValue("")->setPeriodicity(false); + addValueWithDerivatives(); setNotPeriodic(); log.printf(" of %d atoms\n",ga_lista.size()); for(unsigned int i=0;i<ga_lista.size();++i){ diff --git a/src/ColvarDistance.cpp b/src/ColvarDistance.cpp index d954834cb859d52df5fb94a90d815fc7bf7baf48..28d429d777a112376c107f57bf2c134367aa110a 100644 --- a/src/ColvarDistance.cpp +++ b/src/ColvarDistance.cpp @@ -68,17 +68,12 @@ pbc(true) if(!components){ - addValueWithDerivatives(""); - getValue("")->setPeriodicity(false); + addValueWithDerivatives(); setNotPeriodic(); - }else{ - - addValueWithDerivatives("x"); - getValue("x")->setPeriodicity(false); - addValueWithDerivatives("y"); - getValue("y")->setPeriodicity(false); - addValueWithDerivatives("z"); - getValue("z")->setPeriodicity(false); + } else{ + addComponentWithDerivatives("x"); componentIsNotPeriodic("x"); + addComponentWithDerivatives("y"); componentIsNotPeriodic("y"); + addComponentWithDerivatives("z"); componentIsNotPeriodic("z"); } requestAtoms(atoms); @@ -106,24 +101,24 @@ void ColvarDistance::calculate(){ }else{ - Value* valuex=getValue("x"); - Value* valuey=getValue("y"); - Value* valuez=getValue("z"); + Value* valuex=getPntrToComponent("x"); + Value* valuey=getPntrToComponent("y"); + Value* valuez=getPntrToComponent("z"); setAtomsDerivatives (valuex,0,Vector(-1,0,0)); setAtomsDerivatives (valuex,1,Vector(+1,0,0)); setBoxDerivatives (valuex,Tensor(distance,Vector(-1,0,0))); - setValue (valuex,distance[0]); + valuex->set(distance[0]); setAtomsDerivatives (valuey,0,Vector(0,-1,0)); setAtomsDerivatives (valuey,1,Vector(0,+1,0)); setBoxDerivatives (valuey,Tensor(distance,Vector(0,-1,0))); - setValue (valuey,distance[1]); + valuey->set(distance[1]); setAtomsDerivatives (valuez,0,Vector(0,0,-1)); setAtomsDerivatives (valuez,1,Vector(0,0,+1)); setBoxDerivatives (valuez,Tensor(distance,Vector(0,0,-1))); - setValue (valuez,distance[2]); + valuez->set(distance[2]); }; } diff --git a/src/ColvarEnergy.cpp b/src/ColvarEnergy.cpp index acc530aa2e7fa827706a6cb469b1b8d695b5ae45..aa5a305fcf9ba0e1c25974b3a2fbbee7b9699481 100644 --- a/src/ColvarEnergy.cpp +++ b/src/ColvarEnergy.cpp @@ -47,11 +47,9 @@ PLUMED_COLVAR_INIT(ao), components(false) { assert(!checkNumericalDerivatives()); - std::vector<AtomNumber> atoms; - requestAtoms(atoms); isEnergy=true; - addValueWithDerivatives(""); - getValue("")->setPeriodicity(false); + addValueWithDerivatives(); setNotPeriodic(); + getPntrToValue()->resizeDerivatives(1); } void ColvarEnergy::registerKeywords( Keywords& keys ){ @@ -64,7 +62,8 @@ void ColvarEnergy::registerKeywords( Keywords& keys ){ // calculator void ColvarEnergy::calculate(){ - setValue(getEnergy()); + setValue( getEnergy() ); + getPntrToComponent(0)->addDerivative(0,1.0); } } diff --git a/src/ColvarRMSD.cpp b/src/ColvarRMSD.cpp index e62e63d1221cd8394c0136f5d894b225dae90cc5..31615d7ec43b32b6f4e9688a4667caa9066f44ce 100644 --- a/src/ColvarRMSD.cpp +++ b/src/ColvarRMSD.cpp @@ -69,9 +69,7 @@ PLUMED_COLVAR_INIT(ao),rmsd(log) checkRead(); - addValueWithDerivatives(""); - getValue("")->setPeriodicity(false); - + addValueWithDerivatives(); setNotPeriodic(); PDB pdb; // read everything in ang and transform to nm diff --git a/src/ColvarTemplate.cpp b/src/ColvarTemplate.cpp index 32555c46d04458cc319f0a93a99b82a9d915cd82..9b9bac5bcb79aa4bc60f4e8d0cc5a823cc72f21c 100644 --- a/src/ColvarTemplate.cpp +++ b/src/ColvarTemplate.cpp @@ -60,7 +60,7 @@ pbc(true) if(pbc) log.printf(" using periodic boundary conditions\n"); else log.printf(" without periodic boundary conditions\n"); - addValueWithDerivatives(""); + addValueWithDerivatives(); setNotPeriodic(); requestAtoms(atoms); } diff --git a/src/ColvarTorsion.cpp b/src/ColvarTorsion.cpp index e12c461a5edbdded0c9244c3219dec94ba44780c..1c748b0da97fd9bb2b7a73cf2f7f28b38481d570 100644 --- a/src/ColvarTorsion.cpp +++ b/src/ColvarTorsion.cpp @@ -93,10 +93,7 @@ pbc(true) if(pbc) log.printf(" using periodic boundary conditions\n"); else log.printf(" without periodic boundary conditions\n"); - addValueWithDerivatives(""); - getValue("")->setPeriodicity(true); - getValue("")->setDomain(-pi,pi); - + addValueWithDerivatives(); setPeriodic(-pi,pi); requestAtoms(atoms); } diff --git a/src/ColvarVolume.cpp b/src/ColvarVolume.cpp index a70943fb3098040214e550e1abf264aab58fa9b8..c0ddfc693bdb79b5fca0961c762a5504a09ec921 100644 --- a/src/ColvarVolume.cpp +++ b/src/ColvarVolume.cpp @@ -49,9 +49,7 @@ components(false) if(components){ // todo } - addValueWithDerivatives(""); - getValue("")->setPeriodicity(false); - + addValueWithDerivatives(); setNotPeriodic(); requestAtoms(atoms); } diff --git a/src/Function.cpp b/src/Function.cpp index 4bc04b4e0722cb4f0047e9e1b0376d6fa79983a4..b77092de7a323527175048fbdde92424fda9221a 100644 --- a/src/Function.cpp +++ b/src/Function.cpp @@ -9,7 +9,18 @@ Action(ao), ActionWithValue(ao), ActionWithArguments(ao) { - setNumberOfParameters(getNumberOfArguments()); +} + +void Function::addValueWithDerivatives(){ + plumed_massert( getNumberOfArguments()!=0, "for functions you must requestArguments before adding values"); + ActionWithValue::addValueWithDerivatives(); + getPntrToValue()->resizeDerivatives(getNumberOfArguments()); +} + +void Function::addComponentWithDerivatives( const std::string& name ){ + plumed_massert( getNumberOfArguments()!=0, "for functions you must requestArguments before adding values"); + ActionWithValue::addComponentWithDerivatives(name); + getPntrToComponent(name)->resizeDerivatives(getNumberOfArguments()); } void Function::apply(){ @@ -17,16 +28,16 @@ void Function::apply(){ vector<double> f(getNumberOfArguments(),0.0); bool at_least_one_forced=false; - for(int i=0;i<getNumberOfValues();++i){ - if(!getValue(i)->checkForced()) continue; - at_least_one_forced=true; - const vector<double> & derivatives(getValue(i)->getDerivatives()); - for(unsigned j=0;j<derivatives.size();j++){ - f[j]+=getForce(i)*derivatives[j]; + std::vector<double> forces( getNumberOfArguments() ); + for(int i=0;i<getNumberOfComponents();++i){ + if( getPntrToComponent(i)->applyForce( forces ) ){ + at_least_one_forced=true; + for(unsigned j=0;j<forces.size();j++){ f[j]+=forces[j]; } } } + if(at_least_one_forced) for(unsigned i=0;i<getNumberOfArguments();++i){ - getArguments()[i]->addForce(f[i]); + getPntrToArgument(i)->addForce(f[i]); } } @@ -34,7 +45,7 @@ void Function::registerKeywords(Keywords& keys){ Action::registerKeywords(keys); ActionWithValue::registerKeywords(keys); ActionWithArguments::registerKeywords(keys); - keys.add("compulsory","PERIODIC","if the output of your function is periodic then you should specify the periodicity of the function. If the output is not periodic you must state this using PERIODIC=NO"); + keys.reserve("compulsory","PERIODIC","if the output of your function is periodic then you should specify the periodicity of the function. If the output is not periodic you must state this using PERIODIC=NO"); } diff --git a/src/Function.h b/src/Function.h index b982f3d6a02a3e21125d0842fb13cbab14aba35c..5bf74446ddbd08210f4572388a66364d3e82df35 100644 --- a/src/Function.h +++ b/src/Function.h @@ -19,23 +19,26 @@ class Function: public ActionWithValue, public ActionWithArguments { +protected: + void setDerivatives(int,double); + void setDerivatives(Value*,int,double); + void addValueWithDerivatives(); + void addComponentWithDerivatives( const std::string& name ); public: Function(const ActionOptions&); virtual ~Function(){}; void apply(); - void setDerivatives(int,double); - void setDerivatives(Value*,int,double); static void registerKeywords(Keywords&); }; inline void Function::setDerivatives(Value*v,int i,double d){ - v->setDerivatives(i,d); + v->addDerivative(i,d); } inline void Function::setDerivatives(int i,double d){ - setDerivatives(getValue(0),i,d); + setDerivatives(getPntrToValue(),i,d); } } diff --git a/src/FunctionCombine.cpp b/src/FunctionCombine.cpp index d9b2e286b743df6582a14a0c0e760c310c8fc75d..e9ea9e02180fee7a7262d15860aa7028956e5dcb 100644 --- a/src/FunctionCombine.cpp +++ b/src/FunctionCombine.cpp @@ -46,6 +46,7 @@ PLUMED_REGISTER_ACTION(FunctionCombine,"COMBINE") void FunctionCombine::registerKeywords(Keywords& keys){ Function::registerKeywords(keys); + keys.use("ARG"); keys.use("PERIODIC"); keys.add("compulsory","COEFFICIENTS","1.0","the coefficients of the arguments in your function"); keys.add("compulsory","POWERS","1.0","the powers to which you are raising each of the arguments in your function"); keys.addFlag("NORMALIZE",false,"normalize all the coefficents so that in total they are equal to one"); @@ -72,14 +73,17 @@ powers(getNumberOfArguments(),1.0) for(unsigned i=0;i<coefficients.size();i++) coefficients[i]*=(1.0/n); } - addValueWithDerivatives(""); +<<<<<<< HEAD + addValueWithDerivatives(); getPntrToComponent(0)->resizeDerivatives(getNumberOfArguments()); +======= + addValueWithDerivatives(); +>>>>>>> 6166262... Made it so that we don't have to resize the values in the individual double min(0),max(0); std::vector<std::string> period; parseVector("PERIODIC",period); if(period.size()==1 && period[0]=="NO"){ - getValue("")->setPeriodicity(false); + setNotPeriodic(); } else if(period.size()==2 && Tools::convert(period[0],min) && Tools::convert(period[1],max)){ - getValue("")->setPeriodicity(true); - getValue("")->setDomain(min,max); + setPeriodic(min,max); } else error("missing PERIODIC keyword"); checkRead(); diff --git a/src/FunctionMatheval.cpp b/src/FunctionMatheval.cpp index ae77c2fa690cfdf9630b789d042fd17ca10e2b0d..970c51abc04c6f36c1f96045aec935efca986479 100644 --- a/src/FunctionMatheval.cpp +++ b/src/FunctionMatheval.cpp @@ -62,6 +62,7 @@ PLUMED_REGISTER_ACTION(FunctionMatheval,"MATHEVAL") void FunctionMatheval::registerKeywords(Keywords& keys){ Function::registerKeywords(keys); + keys.use("ARG"); keys.use("PERIODIC"); keys.add("compulsory","FUNC","the function you wish to evaluate"); keys.add("optional","VAR","the names to give each of the arguments in the function. If you have up to three arguments in your function you can use x, y and z to refer to them. Otherwise you must use this flag to give your variables names."); @@ -84,14 +85,17 @@ names(getNumberOfArguments()) } assert(var.size()==getNumberOfArguments()); parse("FUNC",func); - addValueWithDerivatives(""); +<<<<<<< HEAD + addValueWithDerivatives(); getPntrToComponent(0)->resizeDerivatives(getNumberOfArguments()); +======= + addValueWithDerivatives(); +>>>>>>> 6166262... Made it so that we don't have to resize the values in the individual double min(0),max(0); std::vector<std::string> period; parseVector("PERIODIC",period); if(period.size()==1 && period[0]=="NO"){ - getValue("")->setPeriodicity(false); + setNotPeriodic(); } else if(period.size()==2 && Tools::convert(period[0],min) && Tools::convert(period[1],max)){ - getValue("")->setPeriodicity(true); - getValue("")->setDomain(min,max); + setPeriodic(min,max); } else error("missing PERIODIC keyword"); checkRead(); diff --git a/src/GenericDumpDerivatives.cpp b/src/GenericDumpDerivatives.cpp index 8678b555d9a3e6c551fe43a10ed410a6916fac90..127a687b9bc39c82429a238df4ef3c5d4507334f 100644 --- a/src/GenericDumpDerivatives.cpp +++ b/src/GenericDumpDerivatives.cpp @@ -73,16 +73,16 @@ fp(NULL) fp=fopen(file.c_str(),"wa"); log.printf(" on file %s\n",file.c_str()); log.printf(" with format %s\n",fmt.c_str()); - fprintf(fp,"%s","#! FIELDS time parameter"); - const std::vector<Value*>& arguments(getArguments()); - assert(arguments.size()>0); - unsigned npar=arguments[0]->getDerivatives().size(); - assert(npar>0); - for(unsigned i=1;i<arguments.size();i++){ - assert(npar==arguments[i]->getDerivatives().size()); + unsigned nargs=getNumberOfArguments(); + if( nargs==0 ) error("no arguments specified"); + unsigned npar=getPntrToArgument(0)->getNumberOfDerivatives(); + if( npar==0 ) error("one or more arguments has no derivatives"); + for(unsigned i=1;i<nargs;i++){ + if( npar!=getPntrToArgument(i)->getNumberOfDerivatives() ) error("the number of derivatives must be the same in all values being dumped"); } - for(unsigned i=0;i<arguments.size();i++){ - fprintf(fp," %s",arguments[i]->getFullName().c_str()); + fprintf(fp,"%s","#! FIELDS time parameter"); + for(unsigned i=0;i<nargs;i++){ + fprintf(fp," %s",getPntrToArgument(i)->getName().c_str()); }; fprintf(fp,"%s","\n"); } @@ -92,13 +92,12 @@ fp(NULL) void GenericDumpDerivatives::update(){ if(comm.Get_rank()!=0)return; - const std::vector<Value*>& arguments(getArguments()); - unsigned npar=arguments[0]->getDerivatives().size(); + unsigned npar=getPntrToArgument(0)->getNumberOfDerivatives(); for(unsigned ipar=0;ipar<npar;ipar++){ fprintf(fp," %f",getTime()); fprintf(fp," %u",ipar); for(unsigned i=0;i<getNumberOfArguments();i++){ - fprintf(fp,fmt.c_str(),arguments[i]->getDerivatives()[ipar]); + fprintf( fp, fmt.c_str(),getPntrToArgument(i)->getDerivative(ipar) ); }; fprintf(fp,"\n"); } diff --git a/src/GenericDumpForces.cpp b/src/GenericDumpForces.cpp index e5c6e57d63c1c883fd738eee964679d7406cf73b..2cde9f2458710c0fd24780527f48f74a27b350a6 100644 --- a/src/GenericDumpForces.cpp +++ b/src/GenericDumpForces.cpp @@ -68,11 +68,10 @@ fp(NULL) if(comm.Get_rank()==0){ fp=fopen(file.c_str(),"wa"); log.printf(" on file %s\n",file.c_str()); + if( getNumberOfArguments()==0 ) error("no arguments have been specified"); fprintf(fp,"%s","#! FIELDS time parameter"); - const std::vector<Value*>& arguments(getArguments()); - assert(arguments.size()>0); - for(unsigned i=0;i<arguments.size();i++){ - fprintf(fp," %s",arguments[i]->getFullName().c_str()); + for(unsigned i=0;i<getNumberOfArguments();i++){ + fprintf(fp," %s",getPntrToArgument(i)->getName().c_str()); }; fprintf(fp,"%s","\n"); } @@ -82,10 +81,9 @@ fp(NULL) void GenericDumpForces::update(){ if(comm.Get_rank()!=0)return; - const std::vector<Value*>& arguments(getArguments()); fprintf(fp," %f",getTime()); for(unsigned i=0;i<getNumberOfArguments();i++){ - fprintf(fp," %15.10f",arguments[i]->getForce()); + fprintf(fp," %15.10f",getPntrToArgument(i)->getForce()); }; fprintf(fp,"\n"); } diff --git a/src/GenericPrint.cpp b/src/GenericPrint.cpp index 7d1607e9ec137a9da097904a0ff925547e16d6cd..5dbab1f0e9fc8825ac295bb901daa5d2ca7a6b42 100644 --- a/src/GenericPrint.cpp +++ b/src/GenericPrint.cpp @@ -80,9 +80,8 @@ rotate(0) fp=fopen(file.c_str(),"a"); log.printf(" on file %s\n",file.c_str()); fprintf(fp,"#! FIELDS time"); - const std::vector<Value*>& arguments(getArguments()); - for(unsigned i=0;i<arguments.size();i++){ - fprintf(fp," %s",arguments[i]->getFullName().c_str()); + for(unsigned i=0;i<getNumberOfArguments();i++){ + fprintf(fp," %s",getPntrToArgument(i)->getName().c_str()); }; fprintf(fp,"\n"); } @@ -99,7 +98,7 @@ rotate(0) parse("_ROTATE",rotate); if(rotate>0){ rotateCountdown=rotate; - rotateArguments=getArguments(); + for(unsigned i=0;i<getNumberOfArguments();++i) rotateArguments.push_back( getPntrToArgument(i) ); vector<Value*> a(1,rotateArguments[0]); requestArguments(vector<Value*>(1,rotateArguments[0])); rotateLast=0; diff --git a/src/PlumedMain.cpp b/src/PlumedMain.cpp index afda4e33700ca53f0fcef429cf7c5ae2c71e32f7..113e0437792849f34dbc84190539bd51be54c714 100644 --- a/src/PlumedMain.cpp +++ b/src/PlumedMain.cpp @@ -454,9 +454,8 @@ void PlumedMain::justCalculate(){ if((*p)->isActive()){ if((*p)->checkNumericalDerivatives()) (*p)->calculateNumericalDerivatives(); else (*p)->calculate(); - if(av)for(int i=0;i<av->getNumberOfValues();++i){ - if(av->getValue(i)->getName()=="bias") bias+=av->getValue(i)->get(); - } + // This retrieves components called bias + if(av) bias+=av->getOutputQuantity("bias"); if(av)av->setGradientsIfNeeded(); ActionWithVirtualAtom*avv=dynamic_cast<ActionWithVirtualAtom*>(*p); if(avv)avv->setGradientsIfNeeded(); diff --git a/src/Value.cpp b/src/Value.cpp index f3ea75151ae3b2bac1c8cbf9bb5bc7da2a1cf8fd..b704ea6c1c8e1def3b7afde7b45121d920351085 100644 --- a/src/Value.cpp +++ b/src/Value.cpp @@ -9,50 +9,47 @@ using namespace PLMD; -Value::Value(ActionWithValue&action,const std::string& name): - action(action), +Value::Value(const std::string& name, const bool withderiv): value(0.0), name(name), - deriv(false), + hasDeriv(withderiv), periodicity(unset), min(0.0), max(0.0), max_minus_min(0.0), inv_max_minus_min(0.0) -{} - -bool Value::isPeriodic()const{ - plumed_massert(periodicity!=unset,"periodicity should be set"); - return periodicity==periodic; -} - -void Value::getDomain(double&min,double&max)const{ - min=this->min; - max=this->max; +{ } -void Value::setPeriodicity(bool p){ - if(p) periodicity=periodic; - else periodicity=notperiodic; +void Value::setupPeriodicity(){ + if( min==0 && max==0 ){ + periodicity=notperiodic; + } else { + periodicity=periodic; + max_minus_min=max-min; + plumed_massert(max_minus_min>0, "your function has a very strange domain?"); + inv_max_minus_min=1.0/max_minus_min; + } } -void Value::setDomain(double min,double max){ - this->min=min; - this->max=max; - max_minus_min=max-min; - if(max_minus_min!=0.0) inv_max_minus_min=1.0/max_minus_min; +bool Value::isPeriodic()const{ + plumed_massert(periodicity!=unset,"periodicity should be set"); + return periodicity==periodic; } - -const std::string Value::getFullName()const{ - if(name.length()==0) return action.getLabel(); - else return action.getLabel()+"."+name; +bool Value::applyForce(std::vector<double>& forces ) const { + plumed_massert( derivatives.size()==forces.size()," forces array has wrong size" ); + if( !hasForce ) return false; + for(unsigned i=0;i<derivatives.size();++i) forces[i]=inputForce*derivatives[i]; + return true; } -void Value::enableDerivatives() -{ - deriv=true;derivatives.resize(action.getNumberOfParameters()); +void Value::getDomain(double&min,double&max) const { + plumed_massert(periodicity==periodic,"function should be periodic"); + min=this->min; + max=this->max; } +<<<<<<< HEAD void Value::setGradients(){ gradients.clear(); diff --git a/src/Value.h b/src/Value.h index 5963b6f3b7986dd25af93951c56dc906d75f7156..09e1d62fc78ac4628ce053b53e63a959ce422014 100644 --- a/src/Value.h +++ b/src/Value.h @@ -13,53 +13,83 @@ namespace PLMD{ class ActionWithValue; -/// Class containing a value which can be addressed by PLMD::ActionWithArguments. -/// It also contains the derivative of this value with respect to -/// an arbitrary number of parameters. -/// Typically, an object of type PLMD::ActionWithValue will contain one or -/// more objects of type PLUMD::Value, one per component. +//+DEVELDOC TOOLBOX Value +/** +A class for holding the value of a function together with its derivatives. +*/ +//+ENDDEVELDOC + +/// Typically, an object of type PLMD::ActionWithValue will contain one +/// object of type PLUMD::Value that will be named after the label. If the +/// PLMD::ActionWithValue is part of a class that calculates multiple components +/// then the class will contain multiple that will be called label.component-name +/// This class is used to pass information between different PLMD::Action +/// objects. However, if you find a use for a tempory PLMD::Value in some method +/// you are implementing please feel free to use it. class Value{ - ActionWithValue&action; +friend class ActionWithValue; +private: +/// The value of the quantity double value; +/// The force acting on this quantity double inputForce; - bool forced; +/// A flag telling us we have a force acting on this quantity + bool hasForce; +/// The derivatives of the quantity stored in value std::vector<double> derivatives; std::map<AtomNumber,Vector> gradients; +/// The name of this quantiy std::string name; - bool deriv; +/// Does this quanity have derivatives + bool hasDeriv; +/// Is this quantity periodic enum {unset,periodic,notperiodic} periodicity; +/// Various quantities that describe the domain of this value double min,max; double max_minus_min; double inv_max_minus_min; +/// Complete the setup of the periodicity + void setupPeriodicity(); public: void setGradients(); - Value(ActionWithValue&action,const std::string& name); + Value(const std::string& name, const bool withderiv); +/// Set the value of the function void set(double); - double get()const; +/// Get the value of the function + double get() const; void setPeriodicity(bool); void setDomain(double,double); - bool isPeriodic()const; - void getDomain(double&,double&)const; - const std::string& getName()const; - const std::string getFullName()const; - void enableDerivatives(); +/// Check if the value is periodic + bool isPeriodic() const; +/// Get the domain of the quantity + void getDomain(double&,double&) const; +/// Get the name of the quantity + const std::string& getName() const; +/// Check whether or not this particular quantity has derivatives bool hasDerivatives()const; - void setNumberOfParameters(int n); - void setDerivatives(int i,double d); - void clearInputForce(); +/// Get the number of derivatives that this particular value has + unsigned getNumberOfDerivatives() const; +/// Set the number of derivatives + void resizeDerivatives(int n); +/// Set all the derivatives to zero void clearDerivatives(); - double getForce()const; +/// Add some derivative to the ith component of the derivatives array + void addDerivative(int i,double d); +/// Get the derivative with respect to component n + double getDerivative(const unsigned n) const; +/// Clear the input force on the variable + void clearInputForce(); +/// Add some force on this value void addForce(double f); - const std::vector<double> & getDerivatives()const; - ActionWithValue& getAction(); +/// Get the value of the force on this colvar + double getForce() const ; +/// Apply the forces to the derivatives using the chain rule (if there are no forces this routine returns false + bool applyForce( std::vector<double>& forces ) const ; +/// Calculate the difference between the instantaneous value of the function and some other point double difference(double)const; +/// Calculate the difference between two values of this function double difference(double,double)const; - -/// check if a force has been added at this step - bool checkForced()const; - static double projection(const Value&,const Value&); - }; inline @@ -78,38 +108,37 @@ const std::string& Value::getName()const{ } inline -ActionWithValue& Value::getAction(){ - return action; -} - -inline -double Value::getForce()const{ - return inputForce; +unsigned Value::getNumberOfDerivatives() const { + plumed_massert(derivatives.size()!=0,"the derivatives array for this value has zero size"); + return derivatives.size(); } inline -const std::vector<double> & Value::getDerivatives()const{ - return derivatives; +double Value::getDerivative(const unsigned n) const { + plumed_massert(n<derivatives.size(),"you are asking for a derivative that is out of bounds"); + return derivatives[n]; } inline -bool Value::hasDerivatives()const{ - return deriv; +bool Value::hasDerivatives() const { + return (derivatives.size()!=0); } inline -void Value::setNumberOfParameters(int n){ - if(deriv)derivatives.resize(n); +void Value::resizeDerivatives(int n){ + plumed_massert(hasDeriv,"cannot resize derivatives in values that have not got derivatives"); + derivatives.resize(n); } inline -void Value::setDerivatives(int i,double d){ +void Value::addDerivative(int i,double d){ + plumed_massert(i<derivatives.size(),"derivative is out of bounds"); derivatives[i]=d; } inline void Value::clearInputForce(){ - forced=false; + hasForce=false; inputForce=0.0; } @@ -118,18 +147,18 @@ void Value::clearDerivatives(){ derivatives.assign(derivatives.size(),0.0); } -inline -bool Value::checkForced()const{ - return forced; -} - inline void Value::addForce(double f){ plumed_massert(hasDerivatives(),"forces can only be added to values with derivatives"); - forced=true; + hasForce=true; inputForce+=f; } +inline +double Value::getForce() const { + return inputForce; +} + inline double Value::difference(double d1,double d2)const{ if(periodicity==notperiodic){