Skip to content
Snippets Groups Projects
Commit 82aac7b8 authored by Gareth Tribello's avatar Gareth Tribello
Browse files

Now possible to use moments with PAMM

parent 9826086a
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
regtest/multicolvar/rt-pamm-periodic/colvar.reference
+ 7
− 7
View file @ 82aac7b8
#! FIELDS time p.mean-1 p.between-1 p.between-2 p.highest
0.000000 0.9911 10.2441 0.5412 0.0002
1.000000 0.9980 10.0760 0.4862 0.0002
2.000000 0.9981 10.0747 0.4856 0.0003
3.000000 0.9959 10.1309 0.4965 0.0003
4.000000 0.9981 10.0696 0.4849 0.0019
5.000000 0.9954 10.1524 0.5005 0.0001
#! FIELDS time p.mean-1 p.between-1 p.between-2 p.highest p.moment-1-2 p.moment-1-3 p.moment-5-2 p.moment-5-3
0.000000 0.9911 10.2441 0.5412 0.0002 0.0009 -0.0001 0.0000 0.0000
1.000000 0.9980 10.0760 0.4862 0.0002 0.0000 -0.0000 0.0000 0.0000
2.000000 0.9981 10.0747 0.4856 0.0003 0.0000 -0.0000 0.0000 0.0000
3.000000 0.9959 10.1309 0.4965 0.0003 0.0001 -0.0000 0.0000 0.0000
4.000000 0.9981 10.0696 0.4849 0.0019 0.0000 -0.0000 0.0000 0.0000
5.000000 0.9954 10.1524 0.5005 0.0001 0.0001 -0.0000 0.0000 0.0000
regtest/multicolvar/rt-pamm-periodic/deriv.reference
+ 3079
− 3079
View file @ 82aac7b8
This diff is collapsed.
regtest/multicolvar/rt-pamm-periodic/plumed.dat
+ 2
− 2
View file @ 82aac7b8
......@@ -53,8 +53,8 @@ TORSIONS ...
DUMPMULTICOLVAR DATA=psi FILE=psi.xyz
DUMPMULTICOLVAR DATA=phi FILE=phi.xyz
PAMM DATA=phi,psi CLUSTERS=2D-testc-0.75.pammp MEAN1={COMPONENT=1} HISTOGRAM={GAUSSIAN NBINS=2 LOWER=0.0 UPPER=1.0 COMPONENT=2} HIGHEST={COMPONENT=4} LABEL=p
# PAMM DATA=phi,psi CLUSTERS=2D-testc-0.75.pammp MEAN1={COMPONENT=1} HISTOGRAM={GAUSSIAN NBINS=2 LOWER=0.0 UPPER=1.0 COMPONENT=2} HIGHEST={COMPONENT=4} NUMERICAL_DERIVATIVES LABEL=pn
PAMM DATA=phi,psi CLUSTERS=2D-testc-0.75.pammp MEAN1={COMPONENT=1} HISTOGRAM={GAUSSIAN NBINS=2 LOWER=0.0 UPPER=1.0 COMPONENT=2} HIGHEST={COMPONENT=4} MOMENTS1={MOMENTS=2-3 COMPONENT=1} MOMENTS2={MOMENTS=2-3 COMPONENT=5} LABEL=p
# PAMM DATA=phi,psi CLUSTERS=2D-testc-0.75.pammp MEAN1={COMPONENT=1} HISTOGRAM={GAUSSIAN NBINS=2 LOWER=0.0 UPPER=1.0 COMPONENT=2} HIGHEST={COMPONENT=4} MOMENTS1={MOMENTS=2-3 COMPONENT=1} MOMENTS2={MOMENTS=2-3 COMPONENT=5} NUMERICAL_DERIVATIVES LABEL=pn
PRINT ARG=p.* FILE=colvar FMT=%8.4f
DUMPDERIVATIVES ARG=p.* FILE=deriv FMT=%8.4f
......
src/multicolvar/PAMM.cpp
+ 1
− 2
View file @ 82aac7b8
......@@ -132,7 +132,7 @@ void PAMM::registerKeywords( Keywords& keys ){
keys.add("compulsory","CLUSTERS","the name of the file that contains the definitions of all the clusters");
keys.add("compulsory","REGULARISE","0.001","don't allow the denominator to be smaller then this value");
keys.use("MEAN"); keys.use("MORE_THAN"); keys.use("SUM"); keys.use("LESS_THAN"); keys.use("HISTOGRAM"); keys.use("HISTOGRAM");
keys.use("MIN"); keys.use("MAX"); keys.use("LOWEST"); keys.use("HIGHEST"); keys.use("ALT_MIN"); keys.use("BETWEEN");
keys.use("MIN"); keys.use("MAX"); keys.use("LOWEST"); keys.use("HIGHEST"); keys.use("ALT_MIN"); keys.use("BETWEEN"); keys.use("MOMENTS");
keys.setComponentsIntroduction("When the label of this action is used as the input for a second you are not referring to a scalar quantity as you are in "
"regular collective variables. The label is used to reference the full set of quantities calculated by "
"the action. This is usual when using \\ref multicolvarfunction. Generally when doing this the set of PAMM variables "
......@@ -202,7 +202,6 @@ MultiColvarFunction(ao)
// This builds the lists
buildSets();
// I think we can put this back for anything with not usespecies
// Now need to pass cutoff information to underlying distance multicolvars
// for(unsigned i=0;i<getNumberOfBaseMultiColvars();++i){
......
src/vesselbase/Moments.cpp
+ 58
− 23
View file @ 82aac7b8
......@@ -30,8 +30,10 @@ namespace vesselbase{
// The calculation of all the colvars is parallelized
// but the loops for calculating moments are not
// Feel free to reimplement this if you know how
class Moments : public StoreDataVessel {
class Moments : public Vessel {
private:
unsigned mycomponent;
StoreDataVessel* mystash;
std::vector<unsigned> powers;
std::vector<Value*> value_out;
public:
......@@ -40,6 +42,7 @@ public:
explicit Moments( const vesselbase::VesselOptions& da );
std::string description();
void resize();
bool calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const { return true; }
void finish( const std::vector<double>& buffer );
bool applyForce( std::vector<double>& forces );
};
......@@ -47,14 +50,19 @@ public:
PLUMED_REGISTER_VESSEL(Moments,"MOMENTS")
void Moments::registerKeywords( Keywords& keys ){
StoreDataVessel::registerKeywords( keys );
Vessel::registerKeywords( keys );
keys.remove("LABEL");
keys.add("compulsory","COMPONENT","1","the component of the vector for which to calculate this quantity");
keys.add("compulsory","MOMENTS","the list of moments that you would like to calculate");
}
void Moments::reserveKeyword( Keywords& keys ){
keys.reserve("optional","MOMENTS","calculate the moments of the distribution of collective variables. "
"The \\f$m\\f$th moment of a distribution is calculated using \\f$\\frac{1}{N} \\sum_{i=1}^N ( s_i - \\overline{s} )^m \\f$, where \\f$\\overline{s}\\f$ is "
"the average for the distribution. The moments keyword takes a lists of integers as input or a range. Each integer is a value of \\f$m\\f$. The final "
"calculated values can be referenced using moment-\\f$m\\f$.");
"calculated values can be referenced using moment-\\f$m\\f$. You can use the COMPONENT keyword in this action but the syntax is slightly different. "
"If you would like the 2nd and third moments of the 3rd component you would use MOMENTS={COMPONENT=3 MOMENTS=2-3}. The moments would then be refered to "
"using the labels moment-3-2 and moment-3-3. This syntax is also required if you are using numbered MOMENT keywords i.e. MOMENTS1, MOMENTS2...");
keys.reset_style("MOMENTS","vessel");
keys.addOutputComponent("moment","MOMENTS","the central moments of the distribution of values. The second moment "
"would be referenced elsewhere in the input file using "
......@@ -62,16 +70,28 @@ void Moments::reserveKeyword( Keywords& keys ){
}
Moments::Moments( const vesselbase::VesselOptions& da) :
StoreDataVessel(da)
Vessel(da)
{
// Build the data stashes in this vessel
mystash=getAction()->buildDataStashes( false, 0.0 );
ActionWithValue* a=dynamic_cast<ActionWithValue*>( getAction() );
plumed_massert(a,"cannot create passable values as base action does not inherit from ActionWithValue");
std::vector<std::string> moments=Tools::getWords(getAllInput(),"\t\n ,");
Tools::interpretRanges(moments); unsigned nn;
std::vector<std::string> moments; std::string valstr;
if( getNumericalLabel()==0 ){
valstr = "moment-";
moments=Tools::getWords(getAllInput(),"\t\n ,");
Tools::interpretRanges(moments); mycomponent=1;
} else {
std::string numstr; parse("COMPONENT",mycomponent);
Tools::convert( mycomponent, numstr); valstr = "moment-" + numstr + "-";
parseVector("MOMENTS",moments); Tools::interpretRanges(moments);
}
unsigned nn;
for(unsigned i=0;i<moments.size();++i){
a->addComponentWithDerivatives( "moment-" + moments[i] );
a->componentIsNotPeriodic( "moment-" + moments[i] );
a->addComponentWithDerivatives( valstr + moments[i] );
a->componentIsNotPeriodic( valstr + moments[i] );
value_out.push_back( a->copyOutput( a->getNumberOfComponents()-1 ) );
Tools::convert( moments[i], nn );
if( nn<2 ) error("moments are only possible for m>=2" );
......@@ -80,58 +100,73 @@ StoreDataVessel(da)
}
void Moments::resize(){
StoreDataVessel::resize();
resizeBuffer(0);
if( getAction()->derivativesAreRequired() ){
for(unsigned i=0;i<value_out.size();++i) value_out[i]->resizeDerivatives( getAction()->getNumberOfDerivatives() );
}
}
std::string Moments::description(){
std::string descri, num;
Tools::convert(powers[0],num);
descri = "value " + getAction()->getLabel() + "." + "moment-" + num + " contains the " + num + "th moment of the distribution";
for(unsigned i=1;i<powers.size();++i){
Tools::convert(powers[i],num);
descri = descri + "\n value " + getAction()->getLabel() + "." + "moment-" + num + " contains the " + num + "th moment of the distribution";
if( getNumericalLabel()==0 ){
descri = "value " + getAction()->getLabel() + "." + "moment-" + num + " contains the " + num + "th moment of the distribution";
for(unsigned i=1;i<powers.size();++i){
Tools::convert(powers[i],num);
descri = descri + "\n value " + getAction()->getLabel() + "." + "moment-" + num + " contains the " + num + "th moment of the distribution";
}
} else {
std::string numlab; Tools::convert( mycomponent, numlab );
descri = "value " + getAction()->getLabel() + "." + "moment-" + numlab + "-" + num + " contains the " + num + "th moment for the distribution of component " + numlab;
for(unsigned i=1;i<powers.size();++i){
Tools::convert(powers[i],num);
descri = descri + "\n value " + getAction()->getLabel() + "." + "moment-" + numlab + "-" + num + " contains the " + num + "th moment for the distribution of component " + numlab;
}
}
return descri;
}
void Moments::finish( const std::vector<double>& buffer ){
StoreDataVessel::finish( buffer );
const double pi=3.141592653589793238462643383279502884197169399375105820974944592307;
unsigned nvals=getAction()->getFullNumberOfTasks();
std::vector<double> myvalues( getAction()->getNumberOfQuantities() );
double mean=0; Value myvalue;
if( getAction()->isPeriodic() ){
std::string str_min, str_max; getAction()->retrieveDomain( str_min, str_max );
double pfactor, min, max; Tools::convert(str_min,min); Tools::convert(str_max,max);
pfactor = 2*pi / ( max-min ); myvalue.setDomain( str_min, str_max );
double sinsum=0, cossum=0, val;
for(unsigned i=0;i<nvals;++i){ val=pfactor*( buffer[i*nspace*vecsize+nspace] - min ); sinsum+=sin(val); cossum+=cos(val); }
double sinsum=0, cossum=0, val;
for(unsigned i=0;i<nvals;++i){ mystash->retrieveValue( i, false, myvalues ); val=pfactor*( myvalues[mycomponent] - min ); sinsum+=sin(val); cossum+=cos(val); }
mean = 0.5 + atan2( sinsum / static_cast<double>( nvals ) , cossum / static_cast<double>( nvals ) ) / (2*pi);
mean = min + (max-min)*mean;
} else {
for(unsigned i=0;i<nvals;++i) mean+=buffer[i*nspace*vecsize+nspace];
for(unsigned i=0;i<nvals;++i){ mystash->retrieveValue( i, false, myvalues ); mean+=myvalues[mycomponent]; }
mean/=static_cast<double>( nvals ); myvalue.setNotPeriodic();
}
for(unsigned npow=0;npow<powers.size();++npow){
double dev1=0;
if( value_out[0]->getNumberOfDerivatives()>0 ){
for(unsigned i=0;i<nvals;++i) dev1+=pow( myvalue.difference( mean, buffer[i*nspace*vecsize+nspace] ), powers[npow] - 1 );
for(unsigned i=0;i<nvals;++i){
mystash->retrieveValue( i, false, myvalues );
dev1+=pow( myvalue.difference( mean, myvalues[mycomponent] ), powers[npow] - 1 );
}
dev1/=static_cast<double>( nvals );
}
double moment=0;
MultiValue myvals( getNumberOfComponents(), getAction()->getNumberOfDerivatives() ); myvals.clearAll();
MultiValue myvals( getAction()->getNumberOfQuantities(), getAction()->getNumberOfDerivatives() ); myvals.clearAll();
for(unsigned i=0;i<nvals;++i){
double tmp=myvalue.difference( mean, buffer[i*nspace*vecsize+nspace] );
mystash->retrieveValue( i, false, myvalues );
double tmp=myvalue.difference( mean, myvalues[mycomponent] );
moment+=pow( tmp, powers[npow] );
if( value_out[npow]->getNumberOfDerivatives() ){
double pref=pow( tmp, powers[npow] - 1 ) - dev1;
retrieveDerivatives( i, false, myvals );
mystash->retrieveDerivatives( i, false, myvals );
for(unsigned j=0;j<myvals.getNumberActive();++j){
unsigned jatom=myvals.getActiveIndex(j);
value_out[npow]->addDerivative(jatom, pref*myvals.getDerivative( 1, jatom ) );
value_out[npow]->addDerivative(jatom, pref*myvals.getDerivative( mycomponent, jatom ) );
}
myvals.clearAll();
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment