All Pamm stuff now works using one implementation only

src/pamm/HBPammObject.cpp

@@ -25,54 +25,17 @@
 namespace PLMD {
 namespace pamm {
 
-HBPammObject::HBPammObject():
-mymulti(NULL),
-regulariser(0.001)
-{
-}
-
-HBPammObject::HBPammObject( const HBPammObject& in ):
-mymulti(in.mymulti),
-regulariser(in.regulariser)
-{
-  for(unsigned i=0;i<in.kernels.size();++i) kernels.push_back( new KernelFunctions( in.kernels[i] ) ); 
-}
-
-HBPammObject::~HBPammObject(){
-  for(unsigned i=0;i<kernels.size();++i) delete kernels[i];
-}
-
 void HBPammObject::setup( const std::string& filename, const double& reg, multicolvar::MultiColvarBase* mybase, std::string& errorstr ){
-  IFile ifile; regulariser=reg; mymulti=mybase;
-  if( !ifile.FileExist(filename) ){
-     errorstr = "could not find file named " + filename;
-     return;
-  }
-
-  std::vector<std::string> valnames(3); 
-  valnames[0]="ptc"; valnames[1]="ssc"; valnames[2]="adc";
-
-  std::vector<Value*> pos; 
-  for(unsigned i=0;i<3;++i){
-    pos.push_back( new Value() ); pos[i]->setNotPeriodic();
-  }
-  
-  ifile.open(filename); ifile.allowIgnoredFields(); kernels.resize(0);
-  for(unsigned k=0;;++k){
-      KernelFunctions* kk = KernelFunctions::read( &ifile, false, valnames );
-      if( !kk ) break ;
-      kk->normalize( pos );
-      kernels.push_back( kk );
-      ifile.scanField();
-  }
-  ifile.close();
-  for(unsigned i=0;i<3;++i) delete pos[i];
+  mymulti=mybase; std::vector<std::string> valnames(3);
+  valnames[0]="ptc"; valnames[1]="ssc"; valnames[2]="adc"; 
+  std::vector<std::string> min(3), max(3); std::vector<bool> pbc(3, false);
+  mypamm.setup( filename, reg, valnames, pbc, min, max, errorstr );
 }
 
 double HBPammObject::get_cutoff() const {
   double sfmax=0;
-  for(unsigned k=0;k<kernels.size();++k){
-      double rcut = kernels[k]->getCenter()[2] + kernels[k]->getContinuousSupport()[2];
+  for(unsigned k=0;k<mypamm.getNumberOfKernels();++k){
+      double rcut = mypamm.getKernelCenter(k)[2] + mypamm.getKernelSupport(k)[2]; 
       if( rcut>sfmax ){ sfmax=rcut; }
   }
   return sfmax;
@@ -83,47 +46,29 @@ double HBPammObject::evaluate( const unsigned& dno, const unsigned& ano, const u
   Vector d_dh = mymulti->getSeparation( myatoms.getPosition(dno), myatoms.getPosition(hno) ); double md_dh = d_dh.modulo(); // hydrogen - donor
   Vector d_ah = mymulti->getSeparation( myatoms.getPosition(ano), myatoms.getPosition(hno) ); double md_ah = d_ah.modulo(); // hydrogen - acceptor
 
-  std::vector<Value*> pos;
-  for(unsigned i=0;i<3;++i){
-    pos.push_back( new Value() ); pos[i]->setNotPeriodic();
-  }
-  // Proton transfer coordinate
-  pos[0]->set( md_dh - md_ah );
-  // Symmetric stretch coordinate
-  pos[1]->set( md_dh + md_ah );
-  // Acceptor donor distance 
-  pos[2]->set( md_da );
-
-  std::vector<double> tderiv(3), tmderiv(3), dderiv(3,0);
-
-  // Now evaluate all kernels  
-  double denom=regulariser, numer=0;
-  for(unsigned i=0;i<kernels.size();++i){
-      double val=kernels[i]->evaluate( pos, tmderiv ); denom+=val;
-      if(i==0){ numer=val; for(unsigned j=0;j<3;++j) tderiv[j]=tmderiv[j]; }
-      for(unsigned j=0;j<3;++j) dderiv[j] += tmderiv[j];
-  }
+  // Create some vectors locally for pamm evaluation
+  std::vector<double> invals( 3 ), outvals( mypamm.getNumberOfKernels() ); 
+  std::vector<std:: vector<double> > der( mypamm.getNumberOfKernels() );
+  for(unsigned i=0;i<der.size();++i) der[i].resize(3);
+  
+  // Evaluate the pamm object
+  invals[0]=md_dh - md_ah; invals[1]=md_dh+md_ah; invals[2]=md_da;
+  mypamm.evaluate( invals, outvals, der );
 
   if( !mymulti->doNotCalculateDerivatives() ){
-      double denom2 = denom*denom, pref;
-      pref = tderiv[0] / denom - numer*dderiv[0]/denom2;
-      mymulti->addAtomDerivatives( 1, dno, ((-pref)/md_dh)*d_dh, myatoms );
-      mymulti->addAtomDerivatives( 1, ano, ((+pref)/md_ah)*d_ah, myatoms  );
-      mymulti->addAtomDerivatives( 1, hno, ((+pref)/md_dh)*d_dh - ((+pref)/md_ah)*d_ah, myatoms );
-      myatoms.addBoxDerivatives( 1, ((-pref)/md_dh)*Tensor(d_dh,d_dh) - ((-pref)/md_ah)*Tensor(d_ah,d_ah) );
-      pref = tderiv[1] / denom - numer*dderiv[1]/denom2;
-      mymulti->addAtomDerivatives( 1, dno, ((-pref)/md_dh)*d_dh, myatoms );
-      mymulti->addAtomDerivatives( 1, ano, ((-pref)/md_ah)*d_ah, myatoms );
-      mymulti->addAtomDerivatives( 1, hno, ((+pref)/md_dh)*d_dh + ((+pref)/md_ah)*d_ah, myatoms );
-      myatoms.addBoxDerivatives( 1, ((-pref)/md_dh)*Tensor(d_dh,d_dh) + ((-pref)/md_ah)*Tensor(d_ah,d_ah) );
-      pref = tderiv[2] / denom - numer*dderiv[2]/denom2;
-      mymulti->addAtomDerivatives( 1, dno, ((-pref)/md_da)*d_da, myatoms );
-      mymulti->addAtomDerivatives( 1, ano, ((+pref)/md_da)*d_da, myatoms );
-      myatoms.addBoxDerivatives( 1, ((-pref)/md_da)*Tensor(d_da,d_da) );
+      mymulti->addAtomDerivatives( 1, dno, ((-der[0][0])/md_dh)*d_dh, myatoms );
+      mymulti->addAtomDerivatives( 1, ano, ((+der[0][0])/md_ah)*d_ah, myatoms  );
+      mymulti->addAtomDerivatives( 1, hno, ((+der[0][0])/md_dh)*d_dh - ((+der[0][0])/md_ah)*d_ah, myatoms );
+      myatoms.addBoxDerivatives( 1, ((-der[0][0])/md_dh)*Tensor(d_dh,d_dh) - ((-der[0][0])/md_ah)*Tensor(d_ah,d_ah) );
+      mymulti->addAtomDerivatives( 1, dno, ((-der[0][1])/md_dh)*d_dh, myatoms );
+      mymulti->addAtomDerivatives( 1, ano, ((-der[0][1])/md_ah)*d_ah, myatoms );
+      mymulti->addAtomDerivatives( 1, hno, ((+der[0][1])/md_dh)*d_dh + ((+der[0][1])/md_ah)*d_ah, myatoms );
+      myatoms.addBoxDerivatives( 1, ((-der[0][1])/md_dh)*Tensor(d_dh,d_dh) + ((-der[0][1])/md_ah)*Tensor(d_ah,d_ah) );
+      mymulti->addAtomDerivatives( 1, dno, ((-der[0][2])/md_da)*d_da, myatoms );
+      mymulti->addAtomDerivatives( 1, ano, ((+der[0][2])/md_da)*d_da, myatoms );
+      myatoms.addBoxDerivatives( 1, ((-der[0][2])/md_da)*Tensor(d_da,d_da) );
   }
-  for(unsigned i=0;i<3;++i) delete pos[i];
-  pos.resize(0);
-  return numer/denom;
+  return outvals[0];
 
 }
 

src/pamm/HBPammObject.h

@@ -22,28 +22,20 @@
 #ifndef __PLUMED_pamm_HBPammObject_h
 #define __PLUMED_pamm_HBPammObject_h
 
-#include <vector>
 #include "tools/Vector.h"
-#include "core/Value.h"
 #include "multicolvar/AtomValuePack.h"
-#include "tools/KernelFunctions.h"
+#include "PammObject.h"
 
 namespace PLMD {
 namespace pamm {
 
 class HBPammObject {
 private:
+/// The Pamm object underlying this HBPamm calculation
+  PammObject mypamm; 
 /// Pointer to base class in multicolvar
   multicolvar::MultiColvarBase* mymulti;
-/// Regularisation parameter to use
-  double regulariser;
-/// List of kernel functions involved
-  std::vector<KernelFunctions*> kernels;
 public:
-// Explicit definitions for constructor, copy constructor and destructor
-  HBPammObject();
-  HBPammObject( const HBPammObject& );
-  ~HBPammObject();
 /// Setup the HBPamm object
   void setup( const std::string& filename, const double& reg, multicolvar::MultiColvarBase* mybase, std::string& errorstr );
 /// Get the cutoff to use throughout

src/pamm/PAMM.cpp

@@ -23,6 +23,7 @@
 #include "tools/KernelFunctions.h"
 #include "tools/IFile.h"
 #include "multicolvar/MultiColvarFunction.h"
+#include "PammObject.h"
 
 //+PLUMEDOC MCOLVARF PAMM
 /*
@@ -108,13 +109,10 @@ namespace pamm {
 
 class PAMM : public multicolvar::MultiColvarFunction {
 private:
-  double regulariser;
-  std::vector<KernelFunctions*> kernels;
-  std::vector<Value*> pos;
+  PammObject mypamm;
 public:
   static void registerKeywords( Keywords& keys );
   PAMM(const ActionOptions&);
-  ~PAMM();
 /// We have to overwrite this here
   unsigned getNumberOfQuantities();
 /// Calculate the weight of this object ( average of input weights )
@@ -162,48 +160,27 @@ MultiColvarFunction(ao)
    }
 
    bool mixed=getBaseMultiColvar(0)->isPeriodic();
+   std::vector<bool> pbc( getNumberOfBaseMultiColvars() ); 
+   std::vector<std::string> valnames( getNumberOfBaseMultiColvars() );
+   std::vector<std::string> min( getNumberOfBaseMultiColvars() ), max( getNumberOfBaseMultiColvars() );
    for(unsigned i=0;i<getNumberOfBaseMultiColvars();++i){
-       pos.push_back( new Value() );
        if( getBaseMultiColvar(i)->isPeriodic()!=mixed ) warning("mixing of periodic and aperiodic base variables in pamm is untested");
-       if( !getBaseMultiColvar(i)->isPeriodic() ){
-           pos[i]->setNotPeriodic();
-       } else {
-           std::string min, max;
-           getBaseMultiColvar(i)->retrieveDomain( min, max );
-           pos[i]->setDomain( min, max );
-       }
-   }
-  
-   // Get labels for base multicolvars
-   std::vector<std::string> valnames( getNumberOfBaseMultiColvars() );
-   for(unsigned i=0;i<valnames.size();++i) valnames[i]=getBaseMultiColvar(i)->getLabel();
-
-   std::string filename; parse("CLUSTERS",filename); IFile ifile; 
-   if( !ifile.FileExist(filename) ) error("could not find file named " + filename);
-   ifile.open(filename); ifile.allowIgnoredFields(); 
- 
-   for(unsigned i=0;;++i){
-       KernelFunctions* kk = KernelFunctions::read( &ifile, false, valnames );
-       if( !kk ) break ;
-       kk->normalize( pos );
-       kernels.push_back( kk );
-       ifile.scanField(); 
+       pbc[i]=getBaseMultiColvar(i)->isPeriodic();
+       if( pbc[i] ) getBaseMultiColvar(i)->retrieveDomain( min[i], max[i] );
+       valnames[i]=getBaseMultiColvar(i)->getLabel();
    }
-   ifile.close();  
+
+   double regulariser; parse("REGULARISE",regulariser);
+   std::string errorstr, filename; parse("CLUSTERS",filename); 
+   mypamm.setup( filename, regulariser, valnames, pbc, min, max, errorstr );
+   if( errorstr.length()>0 ) error( errorstr );
 
    // This builds the lists
    buildSets( false );
-   // Read the regularisation parameter
-   parse("REGULARISE",regulariser);
-}
-
-PAMM::~PAMM(){
-   for(unsigned i=0;i<pos.size();++i) delete pos[i];
-   for(unsigned i=0;i<kernels.size();++i) delete kernels[i];
 }
 
 unsigned PAMM::getNumberOfQuantities(){
-   return 1 + kernels.size();    
+   return 1 + mypamm.getNumberOfKernels();    
 }
 
 void PAMM::calculateWeight( multicolvar::AtomValuePack& myatoms ){
@@ -230,40 +207,34 @@ void PAMM::calculateWeight( multicolvar::AtomValuePack& myatoms ){
 
 double PAMM::compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const {
    unsigned nvars = getNumberOfBaseMultiColvars();    
-   std::vector<std::vector<double> > tderiv( kernels.size() );
-   for(unsigned i=0;i<kernels.size();++i) tderiv[i].resize( nvars );
-   std::vector<double> tval(2), vals( kernels.size() ), dderiv( kernels.size(), 0 );
+   std::vector<std::vector<double> > tderiv( mypamm.getNumberOfKernels() );
+   for(unsigned i=0;i<tderiv.size();++i) tderiv[i].resize( nvars );
+   std::vector<double> tval(2), invals( nvars ), vals( mypamm.getNumberOfKernels() );
 
    for(unsigned i=0;i<nvars;++i){
-       getVectorForTask( myatoms.getIndex(i), false, tval ); pos[i]->set( tval[1] );
+       getVectorForTask( myatoms.getIndex(i), false, tval ); invals[i]=tval[1];
    }
+   mypamm.evaluate( invals, vals, tderiv );
 
-   // Evaluate the set of kernels 
-   double denom=regulariser;
-   for(unsigned i=0;i<kernels.size();++i){ 
-       vals[i]=kernels[i]->evaluate( pos, tderiv[i] );
-       denom+=vals[i]; 
-       for(unsigned j=0;j<nvars;++j) dderiv[j] += tderiv[i][j];
-   }
    // Now set all values other than the first one 
    // This is because of some peverse choices in multicolvar
-   for(unsigned i=1;i<kernels.size();++i) myatoms.setValue( 1+i, vals[i]/denom );
+   for(unsigned i=1;i<vals.size();++i) myatoms.setValue( 1+i, vals[i] );
 
    if( !doNotCalculateDerivatives() ){
-       double denom2=denom*denom; std::vector<double> mypref( 1 + kernels.size() );
+       std::vector<double> mypref( 1 + vals.size() );
        for(unsigned ivar=0;ivar<nvars;++ivar){
            // Get the values of the derivatives
            MultiValue& myder = getVectorDerivatives( myatoms.getIndex(ivar), false );
            // And calculate the derivatives
-           for(unsigned i=0;i<kernels.size();++i) mypref[1+i] = tderiv[i][ivar]/denom - vals[i]*dderiv[ivar]/denom2;
+           for(unsigned i=0;i<vals.size();++i) mypref[1+i] = tderiv[i][ivar];    
            // This is basically doing the chain rule to get the final derivatives
-           superChainRule( 1, 1, 1+kernels.size(), myatoms.getIndex(ivar), mypref, myder, myatoms );
+           superChainRule( 1, 1, 1+vals.size(), myatoms.getIndex(ivar), mypref, myder, myatoms );
            // And clear the derivatives
            myder.clearAll();
        }
    }
 
-   return vals[0] / denom;
+   return vals[0];
 }
 
 Vector PAMM::getCentralAtom(){

src/pamm/PammObject.cpp

+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+   Copyright (c) 2012-2015 The plumed team
+   (see the PEOPLE file at the root of the distribution for a list of names)
+
+   See http://www.plumed-code.org for more information.
+
+   This file is part of plumed, version 2.
+
+   plumed is free software: you can redistribute it and/or modify
+   it under the terms of the GNU Lesser General Public License as published by
+   the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
+
+   plumed is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with plumed.  If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "PammObject.h"
+#include "tools/IFile.h"
+
+namespace PLMD {
+namespace pamm {
+
+PammObject::PammObject():
+regulariser(0.001)
+{
+}
+
+PammObject::PammObject( const PammObject& in ):
+regulariser(in.regulariser),
+pbc(in.pbc),
+min(in.min),
+max(in.max)
+{
+  for(unsigned i=0;i<in.kernels.size();++i) kernels.push_back( new KernelFunctions( in.kernels[i] ) );
+}
+
+PammObject::~PammObject(){
+  for(unsigned i=0;i<kernels.size();++i) delete kernels[i];
+}
+
+void PammObject::setup( const std::string& filename, const double& reg, const std::vector<std::string>& valnames, 
+                        const std::vector<bool>& pbcin, const std::vector<std::string>& imin, const std::vector<std::string>& imax, 
+                        std::string& errorstr ){
+  IFile ifile; regulariser=reg; 
+  if( !ifile.FileExist(filename) ){
+     errorstr = "could not find file named " + filename;
+     return;
+  }
+
+  std::vector<Value*> pos; 
+  pbc.resize( valnames.size() ); 
+  min.resize( valnames.size() );
+  max.resize( valnames.size() );
+  for(unsigned i=0;i<valnames.size();++i){
+    pbc[i]=pbcin[i]; min[i]=imin[i]; max[i]=imax[i];
+    pos.push_back( new Value() ); 
+    if( !pbc[i] ) pos[i]->setNotPeriodic();
+    else pos[i]->setDomain( min[i], max[i] );
+  }
+  
+  ifile.open(filename); ifile.allowIgnoredFields(); kernels.resize(0);
+  for(unsigned k=0;;++k){
+      KernelFunctions* kk = KernelFunctions::read( &ifile, false, valnames );
+      if( !kk ) break ;
+      kk->normalize( pos );
+      kernels.push_back( kk );
+      ifile.scanField();
+  }
+  ifile.close();
+  for(unsigned i=0;i<valnames.size();++i) delete pos[i];
+}
+
+void PammObject::evaluate( const std::vector<double>& invar, std::vector<double>& outvals, std::vector<std::vector<double> >& der ) const {
+   std::vector<Value*> pos;
+   for(unsigned i=0;i<pbc.size();++i){
+      pos.push_back( new Value() ); 
+      if( !pbc[i] ) pos[i]->setNotPeriodic();
+      else pos[i]->setDomain( min[i], max[i] );
+      // And set the value 
+      pos[i]->set( invar[i] );
+   }
+
+   // Evaluate the set of kernels 
+   double denom=regulariser;
+   std::vector<double> dderiv( der[0].size(), 0 );
+   for(unsigned i=0;i<kernels.size();++i){
+       outvals[i]=kernels[i]->evaluate( pos, der[i] );
+       denom+=outvals[i];
+       for(unsigned j=0;j<der[i].size();++j) dderiv[j] += der[i][j];
+   }
+   // Evaluate the set of derivatives
+   for(unsigned i=0;i<kernels.size();++i){
+       outvals[i]/=denom; 
+       for(unsigned j=0;j<der[i].size();++j) der[i][j]=der[i][j]/denom - outvals[i]*dderiv[j]/denom;
+   }
+
+   for(unsigned i=0;i<pbc.size();++i) delete pos[i];
+}
+
+
+}
+}

src/pamm/PammObject.h

+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+   Copyright (c) 2012-2015 The plumed team
+   (see the PEOPLE file at the root of the distribution for a list of names)
+
+   See http://www.plumed-code.org for more information.
+
+   This file is part of plumed, version 2.
+
+   plumed is free software: you can redistribute it and/or modify
+   it under the terms of the GNU Lesser General Public License as published by
+   the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
+
+   plumed is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with plumed.  If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#ifndef __PLUMED_pamm_PammObject_h
+#define __PLUMED_pamm_PammObject_h
+
+#include <vector>
+#include "core/Value.h"
+#include "tools/KernelFunctions.h"
+
+namespace PLMD {
+namespace pamm {
+
+class PammObject {
+private:
+/// Regularisation parameter to use
+  double regulariser;
+/// Is the domain periodic
+  std::vector<bool> pbc;
+/// The domain of the function
+  std::vector<std::string> min, max;
+/// List of kernel functions involved
+  std::vector<KernelFunctions*> kernels;
+public:
+// Explicit definitions for constructor, copy constructor and destructor
+  PammObject();
+  PammObject( const PammObject& );
+  ~PammObject();  
+/// Setup the Pamm object
+  void setup( const std::string& filename, const double& reg, const std::vector<std::string>& valnames, 
+              const std::vector<bool>& pbcin, const std::vector<std::string>& imin, const std::vector<std::string>& imax,
+              std::string& errorstr );
+/// 
+  void evaluate( const std::vector<double>& invar, std::vector<double>& outvals, std::vector<std::vector<double> >& der ) const ;
+///
+  unsigned getNumberOfKernels() const ;
+///
+  std::vector<double> getKernelCenter( const unsigned& kno ) const ;
+///
+  std::vector<double> getKernelSupport( const unsigned& kno ) const ;
+};
+
+inline
+unsigned PammObject::getNumberOfKernels() const {
+  return kernels.size();
+}
+
+inline
+std::vector<double> PammObject::getKernelCenter( const unsigned& kno ) const {
+  return kernels[kno]->getCenter();
+}
+
+inline
+std::vector<double> PammObject::getKernelSupport( const unsigned& kno ) const {
+  return kernels[kno]->getContinuousSupport();
+}
+
+}
+}
+
+#endif
+