From 1c87ddaeb68ea2c521da31ae54b9aaa257ce05fa Mon Sep 17 00:00:00 2001
From: Giovanni Bussi <giovanni.bussi@gmail.com>
Date: Fri, 17 Jun 2011 10:50:55 +0200
Subject: [PATCH] Tentative implementation of periodic variables
Implementation of periodicity as we discussed at the meeting.
In the constructor of a collective variable or function is it possible
to declare, for each value, if it is periodic or not. E.g., torsions
will be periodic, distances not. The difference between two values
(hypothetically, a stored value and the present one) can than
be computed automatically using a Value::difference() method.
If difference() is called for undeclared values, an error is raised.
Thus, functions of unpredictable periodicity (e.g. MATHEVAL) should
allow the user to specify the periodicity (or the non-periodicity)
in the input line.
See how I implemented the setting of periodicity in the available
Colvars. The syntax is a bit baroque (I will probably simplify it),
but it seems to work.
I also added a check on the fact that two labels can not be equal.
---
src/Action.cpp | 1 +
src/ActionWithArguments.h | 7 +++
src/BiasMovingRestraint.cpp | 2 +-
src/BiasRestraint.cpp | 2 +-
src/ColvarDistance.cpp | 5 ++
src/ColvarEnergy.cpp | 1 +
src/ColvarVolume.cpp | 2 +
src/FunctionCombine.cpp | 15 +++++-
src/FunctionMatheval.cpp | 12 ++++-
src/Value.cpp | 35 +++++++++++++
src/Value.h | 100 +++++++++++++++++++++++++++++++-----
test/simplemd/in | 2 +-
test/simplemd/plumed.dat | 21 +++++---
13 files changed, 179 insertions(+), 26 deletions(-)
diff --git a/src/Action.cpp b/src/Action.cpp
index 5021afe7e..930d974e0 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 9f4e99254..839c91889 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 bea3ed798..29ca5d69d 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 0483b06c8..952b7502b 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 870574081..1db02c69d 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 0c91c0d0c..feeb6a059 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 d702eed0c..9d8b150ca 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 0d2c562cd..243d146d5 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 803827775..482b982e0 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 f95b594e0..618a1bb30 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 b252bdf19..88b04aa3c 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 75eb5f934..c0b90d649 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 59fe653f3..f3d09701a 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
--
GitLab