diff --git a/src/Action.cpp b/src/Action.cpp index 5021afe7e4f6cd303d70496a3bcbc9c631695814..930d974e0f9310fc886e6cd8d1175133dd4fd7db 100644 --- a/src/Action.cpp +++ b/src/Action.cpp @@ -26,6 +26,7 @@ Action::Action(const ActionOptions&ao): std::string s; Tools::convert(plumed.getActionSet().size(),s); label="@"+s; } + assert(!plumed.getActionSet().selectWithLabel<Action*>(label)); log.printf(" with label %s\n",label.c_str()); } diff --git a/src/ActionWithArguments.h b/src/ActionWithArguments.h index 9f4e9925410c6b5aa9fcb690e747d0644b15b2c7..839c91889ced349bcd5be6f0fbd78cd345dc1659 100644 --- a/src/ActionWithArguments.h +++ b/src/ActionWithArguments.h @@ -27,6 +27,8 @@ protected: double getArgument(int)const; /// Returns the number of arguments unsigned getNumberOfArguments()const; +/// Takes the difference taking into account pbc for arg i + double difference(int,double,double)const; public: }; @@ -46,6 +48,11 @@ unsigned ActionWithArguments::getNumberOfArguments()const{ return arguments.size(); } +inline +double ActionWithArguments::difference(int i,double d1,double d2)const{ + return arguments[i]->difference(d1,d2); +} + } diff --git a/src/BiasMovingRestraint.cpp b/src/BiasMovingRestraint.cpp index bea3ed798ff5773aab3ed5c78832412a9e5a4565..29ca5d69da4d7f910e902484e1b750f135b6a067 100644 --- a/src/BiasMovingRestraint.cpp +++ b/src/BiasMovingRestraint.cpp @@ -99,7 +99,7 @@ void BiasMovingRestraint::calculate(){ for(unsigned j=0;j<narg;j++) aa[j]=(c1*at[i-1][j]+c2*at[i][j]); } for(unsigned i=0;i<narg;++i){ - const double cv=(getArgument(i)-aa[i]); + const double cv=difference(i,aa[i],getArgument(i)); const double k=kk[i]; const double f=-k*cv; ene+=0.5*k*cv*cv; diff --git a/src/BiasRestraint.cpp b/src/BiasRestraint.cpp index 0483b06c80f1a62eab56986796cd29fd754a3540..952b7502b18a2e5196a72e8d5e7f1190b263b8dd 100644 --- a/src/BiasRestraint.cpp +++ b/src/BiasRestraint.cpp @@ -67,7 +67,7 @@ void BiasRestraint::calculate(){ double ene=0.0; double totf2=0.0; for(unsigned i=0;i<getNumberOfArguments();++i){ - const double cv=(getArgument(i)-at[i]); + const double cv=difference(i,at[i],getArgument(i)); const double k=kappa[i]; const double f=-k*cv; ene+=0.5*k*cv*cv; diff --git a/src/ColvarDistance.cpp b/src/ColvarDistance.cpp index 87057408151f4209f239e98a8bc2cdab5cd091fd..1db02c69dc7af37d009f083ab26c6c976a4ff29f 100644 --- a/src/ColvarDistance.cpp +++ b/src/ColvarDistance.cpp @@ -62,12 +62,17 @@ pbc(true) if(pbc) log.printf(" using periodic boundary conditions\n"); else log.printf(" without periodic boundary conditions\n"); + addValueWithDerivatives(""); + getValue("")->setPeriodicity(false); if(components){ addValueWithDerivatives("x"); + getValue("x")->setPeriodicity(false); addValueWithDerivatives("y"); + getValue("y")->setPeriodicity(false); addValueWithDerivatives("z"); + getValue("z")->setPeriodicity(false); } requestAtoms(atoms); diff --git a/src/ColvarEnergy.cpp b/src/ColvarEnergy.cpp index 0c91c0d0c949d68d058a9dfdd52bab3c31a11b35..feeb6a059fa68878fb40821f53ab8980dc95110d 100644 --- a/src/ColvarEnergy.cpp +++ b/src/ColvarEnergy.cpp @@ -44,6 +44,7 @@ components(false) requestAtoms(atoms); isEnergy=true; addValueWithDerivatives(""); + getValue("")->setPeriodicity(false); } diff --git a/src/ColvarVolume.cpp b/src/ColvarVolume.cpp index d702eed0c3730084c8af36980bfb06ede6767162..9d8b150cae244cc11d353e9b43145fb73acd1a8c 100644 --- a/src/ColvarVolume.cpp +++ b/src/ColvarVolume.cpp @@ -42,6 +42,8 @@ components(false) // todo } addValueWithDerivatives(""); + getValue("")->setPeriodicity(false); + requestAtoms(atoms); } diff --git a/src/FunctionCombine.cpp b/src/FunctionCombine.cpp index 0d2c562cd831a6ea028781e621fe0170a6469fea..243d146d5dc7b69d03d04de78c395b7b708481ff 100644 --- a/src/FunctionCombine.cpp +++ b/src/FunctionCombine.cpp @@ -42,9 +42,22 @@ powers(getNumberOfArguments(),1.0) assert(coefficients.size()==static_cast<unsigned>(getNumberOfArguments())); parseVector("POWERS",powers); assert(powers.size()==static_cast<unsigned>(getNumberOfArguments())); - checkRead(); addValueWithDerivatives(""); + vector<string> period; + + double min(0),max(0); + parseVector("PERIODIC",period); + if(period.size()==0){ + }else if(period.size()==1 && period[0]=="NO"){ + getValue("")->setPeriodicity(false); + } else if(period.size()==2 && Tools::convert(period[0],min) && Tools::convert(period[1],min)){ + getValue("")->setPeriodicity(true); + getValue("")->setDomain(min,max); + } + + checkRead(); + log.printf(" with coefficients:"); for(unsigned i=0;i<coefficients.size();i++) log.printf(" %f",coefficients[i]); log.printf("\n"); diff --git a/src/FunctionMatheval.cpp b/src/FunctionMatheval.cpp index 803827775c1431b450faf2f3b21159f70af97d9f..482b982e08ae989437eea2a4dd913ddea6d30621 100644 --- a/src/FunctionMatheval.cpp +++ b/src/FunctionMatheval.cpp @@ -63,6 +63,17 @@ names(getNumberOfArguments()) } assert(var.size()==getNumberOfArguments()); parse("FUNC",func); + addValueWithDerivatives(""); + vector<string> period; + double min(0),max(0); + parseVector("PERIODIC",period); + if(period.size()==0){ + }else if(period.size()==1 && period[0]=="NO"){ + getValue("")->setPeriodicity(false); + } else if(period.size()==2 && Tools::convert(period[0],min) && Tools::convert(period[1],min)){ + getValue("")->setPeriodicity(true); + getValue("")->setDomain(min,max); + } checkRead(); evaluator=evaluator_create(const_cast<char*>(func.c_str())); @@ -82,7 +93,6 @@ names(getNumberOfArguments()) for(unsigned i=0;i<getNumberOfArguments();i++) evaluator_deriv[i]=evaluator_derivative(evaluator,const_cast<char*>(var[i].c_str())); - addValueWithDerivatives(""); log.printf(" with function : %s\n",func.c_str()); log.printf(" with variables :"); diff --git a/src/Value.cpp b/src/Value.cpp index f95b594e0bc027a9cb60a9501602ee7ef826af82..618a1bb30bde323d4091e7bde4561b1d844ed895 100644 --- a/src/Value.cpp +++ b/src/Value.cpp @@ -3,6 +3,27 @@ using namespace PLMD; +Value::Value(ActionWithValue&action,const std::string& name): + action(action), + value(0.0), + name(name), + deriv(false), + periodicity(unset), + min(0.0), + max(0.0) +{} + +void Value::setPeriodicity(bool p){ + if(p) periodicity=periodic; + else periodicity=notperiodic; +} + +void Value::setDomain(double min,double max){ + this->min=min; + this->max=max; +} + + const std::string Value::getFullName()const{ return action.getLabel()+"."+name; } @@ -12,4 +33,18 @@ void Value::enableDerivatives() deriv=true;derivatives.resize(action.getNumberOfParameters()); } +double Value::difference(double d1,double d2)const{ + assert(periodicity!=unset); + if(periodicity==periodic){ + double s=(d2-d1)/(max-min); + s=Tools::pbc(s); + return s*(max-min); + }else{ + return d2-d1; + } +} + + + + diff --git a/src/Value.h b/src/Value.h index b252bdf1993b030202e17420dea1d7797a900612..88b04aa3ca65a24e6123cb349cad255ec361bafa 100644 --- a/src/Value.h +++ b/src/Value.h @@ -22,24 +22,98 @@ class Value{ std::vector<double> derivatives; std::string name; bool deriv; + enum {unset,periodic,notperiodic} periodicity; + double min,max; public: - Value(ActionWithValue&action,const std::string& name):action(action),value(0.0),name(name),deriv(false){}; - void set(double v){value=v;} - double get()const{return value;} - const std::string& getName()const{return name;} + Value(ActionWithValue&action,const std::string& name); + void set(double); + double get()const; + void setPeriodicity(bool); + void setDomain(double,double); + const std::string& getName()const; const std::string getFullName()const; void enableDerivatives(); - bool hasDerivatives(){return deriv;}; - void setNumberOfParameters(int n){if(deriv)derivatives.resize(n);} - void setDerivatives(int i,double d){derivatives[i]=d;} - void clearInputForce(){inputForce=0.0;} - void clearDerivatives(){derivatives.assign(derivatives.size(),0.0);} - double getForce()const{return inputForce;} - void addForce(double f){assert(hasDerivatives());inputForce+=f;} - const std::vector<double> & getDerivatives()const{return derivatives;} - ActionWithValue& getAction(){return action;}; + bool hasDerivatives()const; + void setNumberOfParameters(int n); + void setDerivatives(int i,double d); + void clearInputForce(); + void clearDerivatives(); + double getForce()const; + void addForce(double f); + const std::vector<double> & getDerivatives()const; + ActionWithValue& getAction(); + + double difference(double)const; + double difference(double,double)const; }; +inline +void Value::set(double v){ + value=v; +} + +inline +double Value::get()const{ + return value; +} + +inline +const std::string& Value::getName()const{ + return name; +} + +inline +ActionWithValue& Value::getAction(){ + return action; +} + +inline +double Value::getForce()const{ + return inputForce; +} + +inline +void Value::addForce(double f){ + assert(hasDerivatives()); + inputForce+=f; +} + +inline +const std::vector<double> & Value::getDerivatives()const{ + return derivatives; +} + +inline +bool Value::hasDerivatives()const{ + return deriv; +} + +inline +void Value::setNumberOfParameters(int n){ + if(deriv)derivatives.resize(n); +} + +inline +void Value::setDerivatives(int i,double d){ + derivatives[i]=d; +} + +inline +void Value::clearInputForce(){ + inputForce=0.0; +} + +inline +void Value::clearDerivatives(){ + derivatives.assign(derivatives.size(),0.0); +} + +inline +double Value::difference(double d)const{ + return difference(get(),d); +} + + } #endif diff --git a/test/simplemd/in b/test/simplemd/in index 75eb5f934559ea4a0f3be2f1f3b29f6bc5a63dc7..c0b90d649cdb9d0e742f9f66fd4a8e6c8b808d7f 100644 --- a/test/simplemd/in +++ b/test/simplemd/in @@ -5,6 +5,6 @@ tstep 0.005 friction 1 forcecutoff 2.5 listcutoff 3.0 -nstep 2000 +nstep 4000 nconfig 10 trajectory.xyz nstat 10 energies.dat diff --git a/test/simplemd/plumed.dat b/test/simplemd/plumed.dat index 59fe653f3ea6193dc843367982639c6452b70e08..f3d09701adb18659c3cf09b2755400244828c086 100644 --- a/test/simplemd/plumed.dat +++ b/test/simplemd/plumed.dat @@ -1,16 +1,21 @@ DISTANCE LABEL=C1 ATOMS=0,1 COMPONENTS +COMBINE LABEL=X ARG=C1.x,C1.y,C1.z COEFFICIENTS=1.0,1.0,1.0 POWERS=2.0,2.0,2.0 PERIODIC=0,1 +MATHEVAL LABEL=Y ARG=C1.x,C1.y FUNC=atan(y/x) PERIODIC=-1.5707963267949,+1.5707963267949 -RESTRAINT ... - KAPPA=10.0 - ARG=C1 - STRIDE=5 - AT=10.0 +MOVINGRESTRAINT ... + ARG=Y + STRIDE=1 + KAPPA0=100.0 + AT0=0.0 + AT1=6.0 + STEP0=0 + STEP1=4000 LABEL=rest -... RESTRAINT +... MOVINGRESTRAINT PRINT ... - STRIDE=25 - ARG=C1,C1.x,C1.y,C1.z,rest.* + STRIDE=5 + ARG=C1.x,C1.y,C1.z,Y,rest.Energy FILE=COLVAR ... PRINT