diff --git a/src/analysis/AnalysisWithLandmarks.cpp b/src/analysis/AnalysisWithLandmarks.cpp
index 32df489ef81ee3acbd4677d0629e7c01bb6f91be..4a4548b884a11ec9317ba71c904d336be0426fbc 100644
--- a/src/analysis/AnalysisWithLandmarks.cpp
+++ b/src/analysis/AnalysisWithLandmarks.cpp
@@ -71,7 +71,7 @@ void AnalysisWithLandmarks::performAnalysis(){
analyzeLandmarks();
}
-void AnalysisWithLandmarks::performTask( const unsigned& taskIndex, const unsigned& current, vesselbase::MultiValue& myvals ){
+void AnalysisWithLandmarks::performTask( const unsigned& taskIndex, const unsigned& current, MultiValue& myvals ) const {
plumed_merror("Should not be here");
}
diff --git a/src/analysis/AnalysisWithLandmarks.h b/src/analysis/AnalysisWithLandmarks.h
index bd9842ddd8e5ad488ee684671b545284440a3acc..86426fc7961c612b86af1e9c0bf6ad25cd96991b 100644
--- a/src/analysis/AnalysisWithLandmarks.h
+++ b/src/analysis/AnalysisWithLandmarks.h
@@ -53,7 +53,7 @@ public:
void performAnalysis();
virtual void analyzeLandmarks()=0;
/// This does nothing
- void performTask( const unsigned& , const unsigned& , vesselbase::MultiValue& );
+ void performTask( const unsigned& , const unsigned& , MultiValue& ) const ;
};
}
diff --git a/src/analysis/Histogram.cpp b/src/analysis/Histogram.cpp
index a6fb5ff3ccaaa9bcc5ef9dc095276b947bbe9135..c3a3899947aa00bec6b03ac13d01dc202ec8ab55 100644
--- a/src/analysis/Histogram.cpp
+++ b/src/analysis/Histogram.cpp
@@ -124,7 +124,7 @@ public:
static void registerKeywords( Keywords& keys );
Histogram(const ActionOptions&ao);
void performAnalysis();
- void performTask( const unsigned& , const unsigned& , vesselbase::MultiValue& );
+ void performTask( const unsigned& , const unsigned& , MultiValue& ) const ;
};
PLUMED_REGISTER_ACTION(Histogram,"HISTOGRAM")
@@ -212,7 +212,7 @@ unnormalized(false)
log.printf("\n");
}
-void Histogram::performTask( const unsigned& , const unsigned& , vesselbase::MultiValue& ){ plumed_error(); }
+void Histogram::performTask( const unsigned& , const unsigned& , MultiValue& ) const { plumed_error(); }
void Histogram::performAnalysis(){
// Back up old histogram files
diff --git a/src/colvar/DRMSD.cpp b/src/colvar/DRMSD.cpp
index 4d45a32977717841b99bc4c604a4f3663473decb..b5693e74e04b6ba42570508396cdfdb22c70ef20 100644
--- a/src/colvar/DRMSD.cpp
+++ b/src/colvar/DRMSD.cpp
@@ -75,9 +75,11 @@ DRMSD REFERENCE=file.pdb LOWER_CUTOFF=0.1 UPPER_CUTOFF=0.8
class DRMSD : public Colvar {
-
+
+ bool pbc_;
+ MultiValue myvals;
+ ReferenceValuePack mypack;
PLMD::DRMSD* drmsd_;
- bool pbc_;
public:
DRMSD(const ActionOptions&);
@@ -96,7 +98,7 @@ void DRMSD::registerKeywords(Keywords& keys){
}
DRMSD::DRMSD(const ActionOptions&ao):
-PLUMED_COLVAR_INIT(ao), pbc_(true)
+PLUMED_COLVAR_INIT(ao), pbc_(true), myvals(1,0), mypack(0,0,myvals)
{
string reference;
parse("REFERENCE",reference);
@@ -125,9 +127,14 @@ PLUMED_COLVAR_INIT(ao), pbc_(true)
std::vector<AtomNumber> atoms;
drmsd_->getAtomRequests( atoms );
- drmsd_->setNumberOfAtoms( atoms.size() );
+// drmsd_->setNumberOfAtoms( atoms.size() );
requestAtoms( atoms );
+ // Setup the derivative pack
+ myvals.resize( 1, 3*atoms.size()+9 );
+ mypack.resize( 0, atoms.size() ); mypack.setValIndex(0);
+ for(unsigned i=0;i<atoms.size();++i) mypack.setAtomIndex( i, i );
+
log.printf(" reference from file %s\n",reference.c_str());
log.printf(" which contains %d atoms\n",getNumberOfAtoms());
@@ -140,14 +147,15 @@ DRMSD::~DRMSD(){
// calculator
void DRMSD::calculate(){
- double drmsd; Tensor virial;
+ double drmsd; Tensor virial; mypack.clear();
- drmsd=drmsd_->calculate(getPositions(), getPbc(), false);
+ drmsd=drmsd_->calculate(getPositions(), getPbc(), mypack, false);
setValue(drmsd);
- for(unsigned i=0;i<getNumberOfAtoms();++i) { setAtomsDerivatives( i, drmsd_->getAtomDerivative(i) ); }
-
- drmsd_->getVirial( virial ); setBoxDerivatives(virial);
+ for(unsigned i=0;i<getNumberOfAtoms();++i) { setAtomsDerivatives( i, mypack.getAtomDerivative(i) ); }
+ setBoxDerivatives( mypack.getBoxDerivatives() );
+
+ // drmsd_->getVirial( virial ); setBoxDerivatives(virial);
}
diff --git a/src/colvar/MultiRMSD.cpp b/src/colvar/MultiRMSD.cpp
index 977151dc3c254d79b3b6ecf72f6344354d171861..a52aa688ee29b263d8dc717e56243abeb0b000aa 100644
--- a/src/colvar/MultiRMSD.cpp
+++ b/src/colvar/MultiRMSD.cpp
@@ -37,6 +37,8 @@ class MultiRMSD : public Colvar {
PLMD::MultiDomainRMSD* rmsd;
bool squared;
+ MultiValue myvals;
+ ReferenceValuePack mypack;
public:
MultiRMSD(const ActionOptions&);
@@ -124,7 +126,7 @@ void MultiRMSD::registerKeywords(Keywords& keys){
}
MultiRMSD::MultiRMSD(const ActionOptions&ao):
-PLUMED_COLVAR_INIT(ao),squared(false)
+PLUMED_COLVAR_INIT(ao),squared(false),myvals(1,0), mypack(0,0,myvals)
{
string reference;
parse("REFERENCE",reference);
@@ -147,9 +149,13 @@ PLUMED_COLVAR_INIT(ao),squared(false)
std::vector<AtomNumber> atoms;
rmsd->getAtomRequests( atoms );
- rmsd->setNumberOfAtoms( atoms.size() );
+// rmsd->setNumberOfAtoms( atoms.size() );
requestAtoms( atoms );
+ myvals.resize( 1, 3*atoms.size()+9 );
+ mypack.resize( 0, atoms.size() ); mypack.setValIndex(0);
+ for(unsigned i=0;i<atoms.size();++i) mypack.setAtomIndex( i, i );
+
log.printf(" reference from file %s\n",reference.c_str());
log.printf(" which contains %d atoms\n",getNumberOfAtoms());
log.printf(" method for alignment : %s \n",type.c_str() );
@@ -163,14 +169,13 @@ MultiRMSD::~MultiRMSD(){
// calculator
void MultiRMSD::calculate(){
- double r=rmsd->calculate( getPositions(), getPbc(), squared );
+ mypack.clear(); double r=rmsd->calculate( getPositions(), getPbc(), mypack, squared );
setValue(r);
- for(unsigned i=0;i<getNumberOfAtoms();i++) setAtomsDerivatives( i, rmsd->getAtomDerivative(i) );
+ for(unsigned i=0;i<getNumberOfAtoms();i++) setAtomsDerivatives( i, mypack.getAtomDerivative(i) );
- Tensor virial;
- if( !rmsd->getVirial( virial ) ) setBoxDerivativesNoPbc();
- else setBoxDerivatives( virial );
+ if( !mypack.virialWasSet() ) setBoxDerivativesNoPbc();
+ else setBoxDerivatives( mypack.getBoxDerivatives() );
}
}
diff --git a/src/colvar/RMSD.cpp b/src/colvar/RMSD.cpp
index 33d040994a6167ab9d9992c2c4a0f4b1bdb1b113..8e9333e42bcc1654109ef0e8d77f54e0b89ed6cb 100644
--- a/src/colvar/RMSD.cpp
+++ b/src/colvar/RMSD.cpp
@@ -35,6 +35,8 @@ namespace colvar{
class RMSD : public Colvar {
+ MultiValue myvals;
+ ReferenceValuePack mypack;
PLMD::RMSDBase* rmsd;
bool squared;
@@ -150,7 +152,7 @@ void RMSD::registerKeywords(Keywords& keys){
}
RMSD::RMSD(const ActionOptions&ao):
-PLUMED_COLVAR_INIT(ao),squared(false)
+PLUMED_COLVAR_INIT(ao),myvals(1,0), mypack(0,0,myvals),squared(false)
{
string reference;
parse("REFERENCE",reference);
@@ -173,9 +175,14 @@ PLUMED_COLVAR_INIT(ao),squared(false)
std::vector<AtomNumber> atoms;
rmsd->getAtomRequests( atoms );
- rmsd->setNumberOfAtoms( atoms.size() );
+// rmsd->setNumberOfAtoms( atoms.size() );
requestAtoms( atoms );
+ // Setup the derivative pack
+ myvals.resize( 1, 3*atoms.size()+9 );
+ mypack.resize( 0, atoms.size() ); mypack.setValIndex(0);
+ for(unsigned i=0;i<atoms.size();++i) mypack.setAtomIndex( i, i );
+
log.printf(" reference from file %s\n",reference.c_str());
log.printf(" which contains %d atoms\n",getNumberOfAtoms());
log.printf(" method for alignment : %s \n",type.c_str() );
@@ -189,12 +196,12 @@ RMSD::~RMSD(){
// calculator
void RMSD::calculate(){
- double r=rmsd->calculate( getPositions(), squared );
+ mypack.clear(); double r=rmsd->calculate( getPositions(), mypack, squared );
setValue(r);
- for(unsigned i=0;i<getNumberOfAtoms();i++) setAtomsDerivatives( i, rmsd->getAtomDerivative(i) );
+ for(unsigned i=0;i<getNumberOfAtoms();i++) setAtomsDerivatives( i, mypack.getAtomDerivative(i) );
- Tensor virial; plumed_dbg_assert( !rmsd->getVirial(virial) );
+ Tensor virial; plumed_dbg_assert( !mypack.virialWasSet() );
setBoxDerivativesNoPbc();
}
diff --git a/src/core/ActionWithArguments.h b/src/core/ActionWithArguments.h
index 3db42230bb8016533b6dceb660dc67a072980964..feef08a6da5f054e6c3d043e1ec1ebffb81d07cd 100644
--- a/src/core/ActionWithArguments.h
+++ b/src/core/ActionWithArguments.h
@@ -78,7 +78,7 @@ public:
void lockRequests();
void unlockRequests();
/// Returns an array of pointers to the arguments
- std::vector<Value*> & getArguments();
+ const std::vector<Value*> & getArguments() const ;
/// Convert a list of argument names into a list of pointers to the values
void interpretArgumentList(const std::vector<std::string>& c, std::vector<Value*>&arg);
};
@@ -120,7 +120,7 @@ void ActionWithArguments::unlockRequests(){
}
inline
-std::vector<Value*> & ActionWithArguments::getArguments(){
+const std::vector<Value*> & ActionWithArguments::getArguments() const {
return arguments;
}
diff --git a/src/crystallization/DFSClustering.cpp b/src/crystallization/DFSClustering.cpp
index 174ebb4433fd01575ff1dc4b60ab987d4afe81b7..3ba1a0d4a3d203adfb9327e0f7dc827787f88ff4 100644
--- a/src/crystallization/DFSClustering.cpp
+++ b/src/crystallization/DFSClustering.cpp
@@ -66,7 +66,7 @@ public:
void completeCalculation();
/// Derivatives of elements of adjacency matrix are unimportant. We thus
/// overwrite this routine as this makes the code faster
- void updateActiveAtoms( multicolvar::AtomValuePack& myatoms ){}
+ void updateActiveAtoms( multicolvar::AtomValuePack& myatoms ) const {}
};
PLUMED_REGISTER_ACTION(DFSClustering,"DFSCLUSTERING")
@@ -151,8 +151,8 @@ void DFSClustering::completeCalculation(){
// Get size for buffer
unsigned bsize=0; std::vector<double> buffer( getSizeOfBuffer( bsize ), 0.0 );
std::vector<double> vals( getNumberOfQuantities() ); std::vector<unsigned> der_index;
- vesselbase::MultiValue myvals( getNumberOfQuantities(), getNumberOfDerivatives() );
- vesselbase::MultiValue bvals( getNumberOfQuantities(), getNumberOfDerivatives() );
+ MultiValue myvals( getNumberOfQuantities(), getNumberOfDerivatives() );
+ MultiValue bvals( getNumberOfQuantities(), getNumberOfDerivatives() );
// Get rid of bogus derivatives
clearDerivatives(); getAdjacencyVessel()->setFinishedTrue();
for(unsigned j=rank;j<myatoms.size();j+=size){
diff --git a/src/crystallization/Fccubic.cpp b/src/crystallization/Fccubic.cpp
index a7d40a05e8743a93c805826e00e97f8b6f59c116..3511b6c37fa6e0ded8eb980c867ae575df38873e 100644
--- a/src/crystallization/Fccubic.cpp
+++ b/src/crystallization/Fccubic.cpp
@@ -47,7 +47,7 @@ public:
static void registerKeywords( Keywords& keys );
Fccubic(const ActionOptions&);
// active methods:
- virtual double compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms );
+ virtual double compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const ;
/// Returns the number of coordinates of the field
bool isPeriodic(){ return false; }
};
@@ -95,7 +95,7 @@ PLUMED_MULTICOLVAR_INIT(ao)
checkRead();
}
-double Fccubic::compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ){
+double Fccubic::compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const {
double value=0, norm=0, dfunc; Vector distance;
// Calculate the coordination number
diff --git a/src/crystallization/Gradient.cpp b/src/crystallization/Gradient.cpp
index 498664c3ab45b829d6e16eef27a02bbbde4db79f..ea4f1a389d79a6c04cce80518399e1468a5aa835 100644
--- a/src/crystallization/Gradient.cpp
+++ b/src/crystallization/Gradient.cpp
@@ -89,9 +89,7 @@ nbins(3)
std::transform( functype.begin(), functype.end(), functype.begin(), tolower );
log.printf(" calculating gradient of %s in %s direction \n",functype.c_str(), direction.c_str() );
- parse("SIGMA",sigma);
- bead.isNotPeriodic(); std::string kerneltype;
- parse("KERNEL",kerneltype); bead.setKernelType( kerneltype );
+ parse("SIGMA",sigma); parse("KERNEL",kerneltype);
checkRead(); requestAtoms(atom);
// And setup the vessel
@@ -104,7 +102,10 @@ void Gradient::setupRegions(){
// if( !getPbc().isOrthorombic() ) error("cell must be orthorhombic when using gradient");
}
-void Gradient::calculateAllVolumes( const unsigned& curr, vesselbase::MultiValue& outvals ){
+void Gradient::calculateAllVolumes( const unsigned& curr, MultiValue& outvals ) const {
+ // Setup the bead
+ HistogramBead bead; bead.isNotPeriodic(); bead.setKernelType( kerneltype );
+
Vector cpos = pbcDistance( getPosition(0), getPntrToMultiColvar()->getCentralAtomPos( curr ) );
// Note we use the pbc from base multicolvar so that we get numerical derivatives correct
Vector oderiv, fpos = (getPntrToMultiColvar()->getPbc()).realToScaled( cpos );
diff --git a/src/crystallization/Gradient.h b/src/crystallization/Gradient.h
index b16d0376d9d5dceb29675a01e0892136aee328b8..7a8b6a8e1754e029c96f05b2e51547231c990822 100644
--- a/src/crystallization/Gradient.h
+++ b/src/crystallization/Gradient.h
@@ -36,8 +36,8 @@ private:
unsigned vend, nquantities;
/// Number of bins in each direction
std::vector<unsigned> nbins;
-/// The bead for the histogram
- HistogramBead bead;
+/// The type of kernel for the histogram
+ std::string kerneltype;
public:
static void registerKeywords( Keywords& keys );
Gradient(const ActionOptions&);
@@ -46,7 +46,7 @@ public:
/// Check on pbc - is it orthorhombic
void setupRegions();
/// Calculate whats in the volume
- void calculateAllVolumes( const unsigned& curr, vesselbase::MultiValue& outvals );
+ void calculateAllVolumes( const unsigned& curr, MultiValue& outvals ) const ;
};
inline
diff --git a/src/crystallization/GradientVessel.cpp b/src/crystallization/GradientVessel.cpp
index fbbb2adbced0d282228130a6fb31d712d52d6c2d..31d9c34e163efa44f54f5f3e448d188533c5ec5a 100644
--- a/src/crystallization/GradientVessel.cpp
+++ b/src/crystallization/GradientVessel.cpp
@@ -40,7 +40,7 @@ public:
GradientVessel( const vesselbase::VesselOptions& da );
std::string function_description();
void resize();
- bool calculate( const unsigned& current, vesselbase::MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const ;
+ bool calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const ;
void finish( const std::vector<double>& buffer );
};
@@ -99,7 +99,7 @@ void GradientVessel::resize(){
}
}
-bool GradientVessel::calculate( const unsigned& current, vesselbase::MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const {
+bool GradientVessel::calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const {
unsigned nder=getAction()->getNumberOfDerivatives();
unsigned wstart, cstart; if( ncomponents==1 ){ cstart=1; wstart=2; } else { cstart=2; wstart=2+ncomponents; }
diff --git a/src/crystallization/MoleculeOrientation.cpp b/src/crystallization/MoleculeOrientation.cpp
index 9517bd0ca6b8a50ba2c584fc0439af62bd4c002e..c47a9143593bdba079ef587cf8980d5429bcd76a 100644
--- a/src/crystallization/MoleculeOrientation.cpp
+++ b/src/crystallization/MoleculeOrientation.cpp
@@ -54,7 +54,7 @@ private:
public:
static void registerKeywords( Keywords& keys );
MoleculeOrientation( const ActionOptions& ao );
- void calculateVector( multicolvar::AtomValuePack& myatoms );
+ void calculateVector( multicolvar::AtomValuePack& myatoms ) const;
};
PLUMED_REGISTER_ACTION(MoleculeOrientation,"MOLECULES")
@@ -87,7 +87,7 @@ VectorMultiColvar(ao)
}
}
-void MoleculeOrientation::calculateVector( multicolvar::AtomValuePack& myatoms ){
+void MoleculeOrientation::calculateVector( multicolvar::AtomValuePack& myatoms ) const {
Vector distance; distance=getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
myatoms.addAtomsDerivatives( 2, 0, Vector(-1.0,0,0) );
diff --git a/src/crystallization/MoleculePlane.cpp b/src/crystallization/MoleculePlane.cpp
index b38129161bcc46b792fff08698136c1da69a512e..e9c0c74d1dd4a282d75fc3dfdc23a6d46c6ab718 100644
--- a/src/crystallization/MoleculePlane.cpp
+++ b/src/crystallization/MoleculePlane.cpp
@@ -41,7 +41,7 @@ private:
public:
static void registerKeywords( Keywords& keys );
MoleculePlane( const ActionOptions& ao );
- void calculateVector( multicolvar::AtomValuePack& myatoms );
+ void calculateVector( multicolvar::AtomValuePack& myatoms ) const ;
};
PLUMED_REGISTER_ACTION(MoleculePlane,"PLANES")
@@ -71,7 +71,7 @@ VectorMultiColvar(ao)
setVectorDimensionality( 3, natoms );
}
-void MoleculePlane::calculateVector( multicolvar::AtomValuePack& myatoms ){
+void MoleculePlane::calculateVector( multicolvar::AtomValuePack& myatoms ) const {
Vector d1, d2, cp;
if( myatoms.getNumberOfAtoms()==3 ){
d1=getSeparation( myatoms.getPosition(1), myatoms.getPosition(0) );
diff --git a/src/crystallization/OrientationSphere.cpp b/src/crystallization/OrientationSphere.cpp
index a6d8042b881c795b5fe48d6032eb0e38f231e564..4bdf3807d47182d220e14cd50d2c9f519b5892c3 100644
--- a/src/crystallization/OrientationSphere.cpp
+++ b/src/crystallization/OrientationSphere.cpp
@@ -65,20 +65,20 @@ MultiColvarFunction(ao)
buildSymmetryFunctionLists();
}
-double OrientationSphere::compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ){
+double OrientationSphere::compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const {
// Make sure derivatives for central atom are only calculated once
VectorMultiColvar* vv = dynamic_cast<VectorMultiColvar*>( getBaseMultiColvar(0) );
vv->firstcall=true;
double d2, sw, value=0, denom=0, dot, f_dot, dot_df, dfunc; Vector distance;
unsigned ncomponents=getBaseMultiColvar(0)->getNumberOfQuantities();
- unsigned nder=getNumberOfDerivatives();
+ unsigned nder=myatoms.getNumberOfDerivatives();
std::vector<double> catom_orient( ncomponents ), this_orient( ncomponents ), catom_der( ncomponents );
Vector catom_pos = myatoms.getPosition(0);
getVectorForTask( myatoms.getIndex(0), true, catom_orient );
multicolvar::CatomPack atom0;
- vesselbase::MultiValue myder0(ncomponents,nder), myder1(ncomponents,nder);
+ MultiValue myder0(ncomponents,nder), myder1(ncomponents,nder);
if( !doNotCalculateDerivatives() ){
atom0=getCentralAtomPackFromInput( myatoms.getIndex(0) );
getVectorDerivatives( myatoms.getIndex(0), true, myder0 );
@@ -121,7 +121,7 @@ double OrientationSphere::compute( const unsigned& tindex, multicolvar::AtomValu
// Now divide everything
double denom2=denom*denom;
- updateActiveAtoms( myatoms ); vesselbase::MultiValue& myvals=myatoms.getUnderlyingMultiValue();
+ updateActiveAtoms( myatoms ); MultiValue& myvals=myatoms.getUnderlyingMultiValue();
for(unsigned i=0;i<myvals.getNumberActive();++i){
unsigned ider=myvals.getActiveIndex(i);
myvals.setDerivative( 1, ider, myvals.getDerivative(1,ider)/denom - (value*myvals.getDerivative(0,ider))/denom2 );
diff --git a/src/crystallization/OrientationSphere.h b/src/crystallization/OrientationSphere.h
index 1b137910ae3e1c8b14efc11b6baf38d3d52902e1..53473b5ef5341796f73547ee273f3c6ecd3f195e 100644
--- a/src/crystallization/OrientationSphere.h
+++ b/src/crystallization/OrientationSphere.h
@@ -38,13 +38,13 @@ private:
public:
static void registerKeywords( Keywords& keys );
OrientationSphere(const ActionOptions&);
- double compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms );
- virtual double transformDotProduct( const double& dot, double& df );
+ double compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const ;
+ virtual double transformDotProduct( const double& dot, double& df ) const ;
bool isPeriodic(){ return false; }
};
inline
-double OrientationSphere::transformDotProduct( const double& dot, double& df ){
+double OrientationSphere::transformDotProduct( const double& dot, double& df ) const {
df=1.0; return dot;
}
diff --git a/src/crystallization/SMAC.cpp b/src/crystallization/SMAC.cpp
index be67b9f581464cf30499059f915c8a8525540b04..9437d0288cdd0aaf2dec1ac2ca6689678033f446 100644
--- a/src/crystallization/SMAC.cpp
+++ b/src/crystallization/SMAC.cpp
@@ -29,13 +29,10 @@ namespace crystallization {
class SMAC : public OrientationSphere {
private:
std::vector<KernelFunctions> kernels;
- std::vector<double> deriv;
- std::vector<Value*> pos;
public:
static void registerKeywords( Keywords& keys );
SMAC(const ActionOptions& ao);
- ~SMAC();
- double transformDotProduct( const double& dot, double& df );
+ double transformDotProduct( const double& dot, double& df ) const ;
};
PLUMED_REGISTER_ACTION(SMAC,"SMAC")
@@ -48,8 +45,7 @@ void SMAC::registerKeywords( Keywords& keys ){
SMAC::SMAC(const ActionOptions& ao):
Action(ao),
-OrientationSphere(ao),
-deriv(1)
+OrientationSphere(ao)
{
std::string kernelinpt;
for(int i=1;;i++){
@@ -58,21 +54,17 @@ deriv(1)
kernels.push_back( mykernel );
}
if( kernels.size()==0 ) error("no kernels defined");
-
- pos.push_back( new Value() );
- pos[0]->setNotPeriodic();
-}
-
-SMAC::~SMAC(){
- delete pos[0];
}
-double SMAC::transformDotProduct( const double& dot, double& df ){
- double ans=0; df=0; pos[0]->set( acos( dot ) ); double dcos=-1./sqrt( 1. - dot*dot );
+double SMAC::transformDotProduct( const double& dot, double& df ) const {
+ std::vector<Value*> pos; pos.push_back( new Value() ); std::vector<double> deriv(1);
+ pos[0]->setNotPeriodic(); pos[0]->set( acos( dot ) );
+ double ans=0; df=0; double dcos=-1./sqrt( 1. - dot*dot );
for(unsigned i=0;i<kernels.size();++i){
ans += kernels[i].evaluate( pos, deriv );
df += deriv[0]*dcos;
}
+ delete pos[0];
return ans;
}
diff --git a/src/crystallization/SimpleCubic.cpp b/src/crystallization/SimpleCubic.cpp
index 09ab4b73fc83cc8eab7c7f434c499afaf1b01115..f53a1d4e87e6c10086fb32e47be10eb173e20048 100644
--- a/src/crystallization/SimpleCubic.cpp
+++ b/src/crystallization/SimpleCubic.cpp
@@ -64,7 +64,7 @@ public:
static void registerKeywords( Keywords& keys );
SimpleCubic(const ActionOptions&);
// active methods:
- virtual double compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms );
+ virtual double compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const ;
/// Returns the number of coordinates of the field
bool isPeriodic(){ return false; }
};
@@ -112,7 +112,7 @@ PLUMED_MULTICOLVAR_INIT(ao)
checkRead();
}
-double SimpleCubic::compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ){
+double SimpleCubic::compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const {
double d2, value=0, norm=0, dfunc; Vector distance;
// Calculate the coordination number
diff --git a/src/crystallization/Steinhardt.cpp b/src/crystallization/Steinhardt.cpp
index 61382e07519f5c290397e9a969f970a8fbe968ae..30c372a413a0ebac40152bd25e6d71080b5aa370 100644
--- a/src/crystallization/Steinhardt.cpp
+++ b/src/crystallization/Steinhardt.cpp
@@ -64,7 +64,7 @@ void Steinhardt::setAngularMomentum( const unsigned& ang ){
tmom=ang; setVectorDimensionality( 2*(2*ang + 1), 2 );
}
-void Steinhardt::calculateVector( multicolvar::AtomValuePack& myatoms ){
+void Steinhardt::calculateVector( multicolvar::AtomValuePack& myatoms ) const {
double dfunc, dpoly_ass, md, tq6, itq6, real_z, imag_z;
Vector distance, dz, myrealvec, myimagvec, real_dz, imag_dz;
// The square root of -1
@@ -163,7 +163,7 @@ void Steinhardt::calculateVector( multicolvar::AtomValuePack& myatoms ){
myatoms.getUnderlyingMultiValue().clear(1);
}
-double Steinhardt::deriv_poly( const unsigned& m, const double& val, double& df ){
+double Steinhardt::deriv_poly( const unsigned& m, const double& val, double& df ) const {
double fact=1.0;
for(unsigned j=1;j<=m;++j) fact=fact*j;
double res=coeff_poly[m]*fact;
diff --git a/src/crystallization/Steinhardt.h b/src/crystallization/Steinhardt.h
index 542ccc3582b076691fe2f5bec3090dd1600a8145..06d06618c57b1a3e57ed9b752edd9ac24ff9c9fd 100644
--- a/src/crystallization/Steinhardt.h
+++ b/src/crystallization/Steinhardt.h
@@ -41,8 +41,8 @@ protected:
public:
static void registerKeywords( Keywords& keys );
Steinhardt( const ActionOptions& ao );
- void calculateVector( multicolvar::AtomValuePack& myatoms );
- double deriv_poly( const unsigned&, const double&, double& );
+ void calculateVector( multicolvar::AtomValuePack& myatoms ) const ;
+ double deriv_poly( const unsigned&, const double&, double& ) const ;
};
}
diff --git a/src/crystallization/Tetrahedral.cpp b/src/crystallization/Tetrahedral.cpp
index 886075a556fb74ca2b26bb1d0e26b33a7cea5489..ba0ecda48ca2203ace1b084c61f0fb441245f53f 100644
--- a/src/crystallization/Tetrahedral.cpp
+++ b/src/crystallization/Tetrahedral.cpp
@@ -50,7 +50,7 @@ public:
static void registerKeywords( Keywords& keys );
Tetrahedral(const ActionOptions&);
// active methods:
- virtual double compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms );
+ virtual double compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const ;
/// Returns the number of coordinates of the field
bool isPeriodic(){ return false; }
};
@@ -98,7 +98,7 @@ PLUMED_MULTICOLVAR_INIT(ao)
checkRead();
}
-double Tetrahedral::compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ){
+double Tetrahedral::compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const {
double value=0, norm=0, dfunc; Vector distance;
// Calculate the coordination number
diff --git a/src/crystallization/VectorMean.cpp b/src/crystallization/VectorMean.cpp
index 965ab55eb9195b73608060f2974225c016946e44..c4c99db2b8a0db622e353f64b4dec0619497bfb4 100644
--- a/src/crystallization/VectorMean.cpp
+++ b/src/crystallization/VectorMean.cpp
@@ -35,7 +35,7 @@ public:
VectorMean( const vesselbase::VesselOptions& da );
std::string function_description();
void resize();
- bool calculate( const unsigned& current, vesselbase::MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const ;
+ bool calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const ;
void finish( const std::vector<double>& buffer );
};
@@ -73,7 +73,7 @@ void VectorMean::resize(){
}
}
-bool VectorMean::calculate( const unsigned& current, vesselbase::MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const {
+bool VectorMean::calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const {
unsigned ncomp=getAction()->getNumberOfQuantities()-2, nder=getAction()->getNumberOfDerivatives();
double weight=myvals.get(0); plumed_dbg_assert( weight>=getTolerance() );
diff --git a/src/crystallization/VectorMultiColvar.cpp b/src/crystallization/VectorMultiColvar.cpp
index 414cf8962cb9a162cdac30b540d191e8bd37744e..39927875ecbfe2509d1e8ddc1422031701a62e01 100644
--- a/src/crystallization/VectorMultiColvar.cpp
+++ b/src/crystallization/VectorMultiColvar.cpp
@@ -43,23 +43,13 @@ void VectorMultiColvar::setVectorDimensionality( const unsigned& ncomp, const in
ncomponents = ncomp;
// Read in the atoms if we are using multicolvar reading
int natoms=nat; readAtoms( natoms );
- // Create the store vector object
-// std::string param; vesselbase::VesselOptions da("","",0,param,this);
-// Keywords keys; StoreVectorsVessel::registerKeywords( keys );
-// vesselbase::VesselOptions da2(da,keys);
-// vecs = new StoreVectorsVessel(da2);
-// // Add the vessel to the base
-// addVessel(vecs);
-// // Read in any vessels
-// readVesselKeywords();
- // Resize a holder for the derivatives of the norm of the vector
}
void VectorMultiColvar::doNotCalculateDirector(){
store_director=false; // Need a sanity check in here so that you don't use the same instance of Q4 to calcualte vectors and directors
}
-double VectorMultiColvar::doCalculation( const unsigned& taskIndex, multicolvar::AtomValuePack& myatoms ){
+double VectorMultiColvar::doCalculation( const unsigned& taskIndex, multicolvar::AtomValuePack& myatoms ) const {
// Now calculate the vector
calculateVector( myatoms );
// Sort out the active derivatives
@@ -74,139 +64,16 @@ double VectorMultiColvar::doCalculation( const unsigned& taskIndex, multicolvar:
double inorm = 1.0 / norm; std::vector<double> dervec( ncomponents );
for(unsigned i=0;i<ncomponents;++i) dervec[i] = inorm*myatoms.getValue(2+i);
- vesselbase::MultiValue& myvals=myatoms.getUnderlyingMultiValue();
+ MultiValue& myvals=myatoms.getUnderlyingMultiValue();
for(unsigned j=0;j<myvals.getNumberActive();++j){
unsigned jder=myvals.getActiveIndex(j);
for(unsigned i=0;i<ncomponents;++i) myvals.addDerivative( 1, jder, dervec[i]*myvals.getDerivative( 2+i, jder ) );
}
}
- // if(complexvec){
- // for(unsigned i=0;i<ncomponents;++i) norm += getComponent(i)*getComponent(i) + getImaginaryComponent(i)*getImaginaryComponent(i);
- // norm=sqrt(norm); inorm = 1.0 / norm;
- // for(unsigned i=0;i<ncomponents;++i){ dervec[i] = inorm*getComponent(i); dervec[ncomponents+i] = inorm*getImaginaryComponent(i); }
- // } else {
- // for(unsigned i=0;i<ncomponents;++i) norm += getComponent(i)*getComponent(i);
- // norm=sqrt(norm); inorm = 1.0 / norm;
- // for(unsigned i=0;i<ncomponents;++i) dervec[i] = inorm*getComponent(i);
- // }
-
- // if( !doNotCalculateDerivatives() ){
- // if( usingLowMem() ){
- // vecs->storeDerivativesLowMem( 0 );
- // vecs->chainRule( 0, dervec );
- // } else {
- // //vecs->storeDerivativesHighMem( getCurrentPositionInTaskList() );
- // vecs->chainRule( getCurrentPositionInTaskList(), dervec );
- // }
-
- // // Add derivatives to base multicolvars
- // Vector tmpd;
- // for(unsigned i=0;i<atoms_with_derivatives.getNumberActive();++i){
- // unsigned k=atoms_with_derivatives[i];
- // tmpd[0]=vecs->getFinalDerivative(3*i+0);
- // tmpd[1]=vecs->getFinalDerivative(3*i+1);
- // tmpd[2]=vecs->getFinalDerivative(3*i+2);
- // MultiColvarBase::addAtomsDerivatives( 0, k, tmpd );
- // }
- // unsigned vvbase=3*atoms_with_derivatives.getNumberActive(); Tensor tmpv;
- // for(unsigned i=0;i<3;++i){
- // for(unsigned j=0;j<3;++j){
- // tmpv(i,j) = vecs->getFinalDerivative( vvbase+3*i+j );
- // }
- // }
- // MultiColvarBase::addBoxDerivatives( 0, tmpv );
- // }
return norm;
}
-// vesselbase::StoreDataVessel* VectorMultiColvar::buildDataStashes( const bool& allow_wcutoff, const double& wtol ){
-// // Build everyting for the multicolvar
-// vesselbase::StoreDataVessel* vsv=MultiColvarBase::buildDataStashes( allow_wcutoff, wtol );
-// if( allow_wcutoff ) vsv->setHardCutoffOnWeight( wtol );
-// // Resize the variable
-// vecs->resize();
-// // And make sure we set up the vector storage correctly
-// vv1.resize( 1 ); vv2.resize( getNumberOfQuantities() - 5 );
-// // And return
-// return vsv;
-// }
-//
-// void VectorMultiColvar::getValueForTask( const unsigned& iatom, std::vector<double>& vals ){
-// plumed_dbg_assert( vecs && vals.size()==(getNumberOfQuantities()-4) );
-// MultiColvarBase::getValueForTask( iatom, vv1 ); vecs->getVector( iatom, vv2 );
-// vals[0]=vv1[0]; for(unsigned i=0;i<vv2.size();++i) vals[i+1]=vv2[i];
-// }
-
-// void VectorMultiColvar::addWeightedValueDerivatives( const unsigned& iatom, const unsigned& base_cv_no, const double& weight, multicolvar::MultiColvarFunction* func ){
-// if( usingLowMem() ){
-// vecs->recompute( iatom, 1 );
-// for(unsigned j=0;j<getNumberOfQuantities()-5;++j) vecs->chainRuleForComponent( 1, j, 5+j, base_cv_no, weight, func );
-// } else {
-// for(unsigned j=0;j<getNumberOfQuantities()-5;++j) vecs->chainRuleForComponent( iatom, j, 5+j, base_cv_no, weight, func );
-// }
-// }
-
-// void VectorMultiColvar::finishWeightedAverageCalculation( multicolvar::MultiColvarFunction* func ){
-// // And calculate the norm of the vector
-// double norm=0, inorm; std::vector<unsigned> tmpindices( 1 + func->getNumberOfDerivatives() );
-// if(complexvec){
-// for(unsigned i=0;i<ncomponents;++i){
-// // Calculate average vector
-// func->quotientRule(5+i, 1, 5+i); func->quotientRule(5+ncomponents+i, 1, 5+ncomponents+i);
-// // Calculate length of vector
-// norm += func->getElementValue(5+i)*func->getElementValue(5+i) + func->getElementValue(5+ncomponents+i)*func->getElementValue(5+ncomponents+i);
-// }
-// norm=sqrt(norm); inorm = 1.0 / norm;
-// for(unsigned i=0;i<ncomponents;++i){
-// dervec[i] = inorm*func->getElementValue(5+i); dervec[ncomponents+i] = inorm*func->getElementValue(5+ncomponents+i);
-// }
-// func->getIndexList( 1, 0, func->getNumberOfDerivatives(), tmpindices );
-// unsigned nder = func->getNumberOfDerivatives();
-// for(unsigned i=0;i<tmpindices[0];++i){
-// unsigned ind = tmpindices[1+i];
-// for(unsigned j=0;j<ncomponents;++j){
-// func->addElementDerivative( ind, dervec[j]*func->getElementDerivative(nder*(5+j) + ind) );
-// func->addElementDerivative( ind, dervec[ncomponents+j]*func->getElementDerivative(nder*(5+ncomponents+j) + ind) );
-// }
-// }
-// } else {
-// for(unsigned i=0;i<ncomponents;++i){
-// // Calculate average vector
-// func->quotientRule(5+i, 1, 5+i);
-// // Calculate length of vector
-// norm += func->getElementValue(5+i)*func->getElementValue(5+i);
-// }
-// norm=sqrt(norm); inorm = 1.0 / norm;
-// for(unsigned i=0;i<ncomponents;++i) dervec[i] = inorm*func->getElementValue(5+i);
-// func->getIndexList( 1, 0, func->getNumberOfDerivatives(), tmpindices );
-// // And set derivatives given magnitude of the vector
-// unsigned nder = func->getNumberOfDerivatives();
-// for(unsigned i=0;i<tmpindices[0];++i){
-// unsigned ind = tmpindices[1+i];
-// for(unsigned j=0;j<ncomponents;++j){
-// func->addElementDerivative( ind, dervec[j]*func->getElementDerivative(nder*(5+j) + ind) );
-// }
-// }
-// }
-// func->setElementValue( 0, norm );
-// }
-
-// void VectorMultiColvar::addOrientationDerivativesToBase( const unsigned& iatom, const unsigned& jstore, const unsigned& base_cv_no,
-// const std::vector<double>& der, multicolvar::MultiColvarFunction* func ){
-// if( usingLowMem() ){
-// if(jstore==1){
-// if(firstcall){ vecs->recompute( iatom, jstore ); firstcall=false; }
-// vecs->chainRuleForVector( jstore, 0, base_cv_no, der, func );
-// } else {
-// vecs->recompute( iatom, jstore );
-// vecs->chainRuleForVector( jstore, 0, base_cv_no, der, func );
-// }
-// } else {
-// vecs->chainRuleForVector( iatom, 0, base_cv_no, der, func );
-// }
-// }
-
void VectorMultiColvar::addForcesOnAtoms( const std::vector<double>& inforces ){
plumed_dbg_assert( inforces.size()==getNumberOfDerivatives() );
std::vector<double> oldforces( getNumberOfDerivatives() );
@@ -215,25 +82,5 @@ void VectorMultiColvar::addForcesOnAtoms( const std::vector<double>& inforces ){
setForcesOnAtoms( oldforces );
}
-// void VectorMultiColvar::copyElementsToBridgedColvar( multicolvar::BridgedMultiColvarFunction* func ){
-// MultiColvarBase::copyElementsToBridgedColvar( func );
-//
-// for(unsigned icomp=5;icomp<getNumberOfQuantities();++icomp){
-// func->setElementValue( icomp-4, getElementValue(icomp) );
-// unsigned nbase = icomp * getNumberOfDerivatives();
-// unsigned nbasev = (icomp-4) * func->getNumberOfDerivatives();
-// for(unsigned jatom=0;jatom<atoms_with_derivatives.getNumberActive();++jatom){
-// unsigned n=atoms_with_derivatives[jatom], nx=nbase + 3*n, ny=nbasev + 3*n;
-// func->addElementDerivative( ny+0, getElementDerivative(nx+0) );
-// func->addElementDerivative( ny+1, getElementDerivative(nx+1) );
-// func->addElementDerivative( ny+2, getElementDerivative(nx+2) );
-// }
-// unsigned nwvir=nbase + 3*getNumberOfAtoms(), nwvirv=nbasev + 3*getNumberOfAtoms();
-// for(unsigned k=0;k<9;++k){
-// func->addElementDerivative( nwvirv, getElementDerivative(nwvir) ); nwvir++; nwvirv++;
-// }
-// }
-// }
-
}
}
diff --git a/src/crystallization/VectorMultiColvar.h b/src/crystallization/VectorMultiColvar.h
index e0324a27383193aa7bfce83204d9d7c08e4b8288..37db089aa83b4632c155717a30c34b08ae867bd0 100644
--- a/src/crystallization/VectorMultiColvar.h
+++ b/src/crystallization/VectorMultiColvar.h
@@ -44,36 +44,6 @@ private:
protected:
/// Set the dimensionality of the vector
void setVectorDimensionality( const unsigned&, const int& );
-// /// Add some value to the ith component of the vector
-// void addComponent( const unsigned&, const double& );
-// /// Get the ith component
-// double getComponent( const unsigned& ) const ;
-// /// Set the ith component
-// void setComponent( const unsigned&, const double& );
-// /// Add derivatives of ith component of vector with repect to jth atom
-// void addAtomsDerivative( const unsigned&, const unsigned&, const Vector& );
-// /// Add atomic derivatives to all components of matrix (note iatom is treated literally here - cf above)
-// void addAtomDerivativeToAllRealComponents( const unsigned& iatom, const std::vector<double>& vec, const Vector& avec );
-// /// Add derivatives of ith component of vector with respect to the box
-// void addBoxDerivatives( const unsigned&, const Tensor& );
-// /// Add box derivatives to all components of matrix
-// void addBoxDerivativesToAllRealComponents( const std::vector<double>& vec, const Tensor& avec );
-// /// Add some value to the imaginary part of the ith component of the vector
-// void addImaginaryComponent( const unsigned&, const double& );
-// /// Get the ith component
-// double getImaginaryComponent( const unsigned& ) const ;
-// /// Set the ith component
-// void setImaginaryComponent( const unsigned&, const double& );
-// /// Add derivatives of the imaginary part of the ith component of vector with repect to jth atom
-// void addImaginaryAtomsDerivative( const unsigned&, const unsigned&, const Vector& );
-// /// Add atomic derivatives to all components of matrix (note iatom is treated literally here - cf above)
-// void addAtomDerivativeToAllImagComponents( const unsigned& iatom, const std::vector<double>& vec, const Vector& avec );
-// /// Add derivatives of the imaginary part of the ith component of vector with respect to the box
-// void addImaginaryBoxDerivatives( const unsigned&, const Tensor& );
-// /// Add box derivatives to all components of matrix
-// void addBoxDerivativesToAllImagComponents( const std::vector<double>& vec, const Tensor& avec );
-// /// This can be used to accumulate derivative from a store of vectors
-// void accumulateDerivativesFromVector( const unsigned& ivec, const unsigned& base_cv_no, const double& weight, StoreVectorsVessel* vectors );
/// Used in vector average to add forces from vector the the forces from here
void addForcesOnAtoms( const std::vector<double>& inforces );
public:
@@ -83,32 +53,19 @@ public:
/// The norm of a vector is not periodic
virtual bool isPeriodic(){ return false; }
/// Calculate the multicolvar
- double doCalculation( const unsigned& taskIndex, multicolvar::AtomValuePack& myatoms );
+ double doCalculation( const unsigned& taskIndex, multicolvar::AtomValuePack& myatoms ) const ;
/// This shouldn't do anything
- double compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ){ plumed_error(); }
+ double compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const { plumed_error(); }
/// Calculate the vector
- virtual void calculateVector( multicolvar::AtomValuePack& myatoms )=0;
+ virtual void calculateVector( multicolvar::AtomValuePack& myatoms ) const=0;
/// Get the number of components in the vector
unsigned getNumberOfComponentsInVector() const ;
/// Get the number of quantities we are calculating per step
unsigned getNumberOfQuantities();
-/// Create places to store the data
-// vesselbase::StoreDataVessel* buildDataStashes( const bool& allow_wcutoff, const double& wtol );
-/// Get the vector
-// void getValueForTask( const unsigned& iatom, std::vector<double>& vals );
-/// Used to accumulate values
-// void addWeightedValueDerivatives( const unsigned& iatom, const unsigned& base_cv_no, const double& weight, multicolvar::MultiColvarFunction* func );
-/// Used for calculating weighted averages
-// void finishWeightedAverageCalculation( multicolvar::MultiColvarFunction* func );
-/// Used in functions to add derivatives to the orientation vector
-// void addOrientationDerivativesToBase( const unsigned& iatom, const unsigned& jstore, const unsigned& base_cv_no,
-// const std::vector<double>& der, multicolvar::MultiColvarFunction* func );
/// Can we differentiate the orientation - yes we can the multicolvar is a vector
bool hasDifferentiableOrientation() const { return true; }
/// This makes sure we are not calculating the director when we do LocalAverage
virtual void doNotCalculateDirector();
-/// Used by ActionVolume
-// void copyElementsToBridgedColvar( multicolvar::BridgedMultiColvarFunction* func );
};
inline
@@ -116,67 +73,6 @@ unsigned VectorMultiColvar::getNumberOfComponentsInVector() const {
return ncomponents;
}
-// inline
-// void VectorMultiColvar::addComponent( const unsigned& icomp, const double& val ){
-// plumed_dbg_assert( icomp<ncomponents );
-// addElementValue( 5 + icomp, val );
-// }
-//
-// inline
-// void VectorMultiColvar::setComponent( const unsigned& icomp, const double& val ){
-// plumed_dbg_assert( icomp<ncomponents );
-// setElementValue( 5 + icomp, val );
-// }
-//
-// inline
-// double VectorMultiColvar::getComponent( const unsigned& icomp ) const {
-// plumed_dbg_assert( icomp<ncomponents );
-// return getElementValue( 5 + icomp );
-// }
-//
-//
-// inline
-// void VectorMultiColvar::addAtomsDerivative( const unsigned& icomp, const unsigned& jatom, const Vector& der ){
-// plumed_dbg_assert( icomp<ncomponents && jatom<getNAtoms() );
-// MultiColvarBase::addAtomsDerivatives( 5 + icomp, current_atoms[jatom], der );
-// }
-
-// inline
-// void VectorMultiColvar::addBoxDerivatives( const unsigned& icomp, const Tensor& vir ){
-// plumed_dbg_assert( icomp<ncomponents );
-// MultiColvarBase::addBoxDerivatives( 5 + icomp, vir );
-// }
-//
-// inline
-// void VectorMultiColvar::addImaginaryComponent( const unsigned& icomp, const double& val ){
-// plumed_dbg_assert( icomp<ncomponents && complexvec );
-// addElementValue( 5 + ncomponents + icomp, val );
-// }
-//
-// inline
-// void VectorMultiColvar::setImaginaryComponent( const unsigned& icomp, const double& val ){
-// plumed_dbg_assert( icomp<ncomponents && complexvec );
-// setElementValue( 5 + ncomponents + icomp, val );
-// }
-//
-// inline
-// double VectorMultiColvar::getImaginaryComponent( const unsigned& icomp ) const {
-// plumed_dbg_assert( icomp<ncomponents && complexvec );
-// return getElementValue( 5 + ncomponents + icomp );
-// }
-//
-// inline
-// void VectorMultiColvar::addImaginaryAtomsDerivative( const unsigned& icomp, const unsigned& jatom, const Vector& der){
-// plumed_dbg_assert( icomp<ncomponents && complexvec && jatom<getNAtoms() );
-// MultiColvarBase::addAtomsDerivatives( 5 + ncomponents + icomp, current_atoms[jatom], der );
-// }
-//
-// inline
-// void VectorMultiColvar::addImaginaryBoxDerivatives( const unsigned& icomp, const Tensor& vir ){
-// plumed_dbg_assert( icomp<ncomponents && complexvec );
-// MultiColvarBase::addBoxDerivatives( 5 + ncomponents + icomp, vir );
-// }
-
inline
unsigned VectorMultiColvar::getNumberOfQuantities(){
return 2 + ncomponents;
diff --git a/src/crystallization/VectorSum.cpp b/src/crystallization/VectorSum.cpp
index 1c6b3bf7d8c5727213994653dfb2329baa9c451b..52a6e3f4439226adbcf92cbf6747987aec92326f 100644
--- a/src/crystallization/VectorSum.cpp
+++ b/src/crystallization/VectorSum.cpp
@@ -35,7 +35,7 @@ public:
VectorSum( const vesselbase::VesselOptions& da );
std::string function_description();
void resize();
- bool calculate( const unsigned& current, vesselbase::MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const ;
+ bool calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const ;
void finish( const std::vector<double>& buffer );
};
@@ -73,7 +73,7 @@ void VectorSum::resize(){
}
}
-bool VectorSum::calculate( const unsigned& current, vesselbase::MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const {
+bool VectorSum::calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const {
unsigned ncomp=getAction()->getNumberOfQuantities()-2, nder=getAction()->getNumberOfDerivatives();
double weight=myvals.get(0);
diff --git a/src/function/Target.cpp b/src/function/Target.cpp
index 2f416fc4ad0199a71305454cabecaabc2e6c2036..b9eef4c2bc006ae9238b5a8903a0784f7ea49e6f 100644
--- a/src/function/Target.cpp
+++ b/src/function/Target.cpp
@@ -45,6 +45,8 @@ set of collective variables.
class Target : public Function {
private:
+ MultiValue myvals;
+ ReferenceValuePack mypack;
PLMD::ArgumentOnlyDistance* target;
public:
Target(const ActionOptions&);
@@ -67,7 +69,9 @@ void Target::registerKeywords(Keywords& keys){
Target::Target(const ActionOptions&ao):
Action(ao),
-Function(ao)
+Function(ao),
+myvals(1,0),
+mypack(0,0,myvals)
{
std::string type; parse("TYPE",type);
std::string reference; parse("REFERENCE",reference);
@@ -84,13 +88,17 @@ Function(ao)
// Get the argument names
std::vector<std::string> args_to_retrieve;
target->getArgumentRequests( args_to_retrieve, false );
- target->setNumberOfArguments( args_to_retrieve.size() );
// Get the arguments
std::vector<Value*> myargs;
interpretArgumentList( args_to_retrieve, myargs );
requestArguments( myargs );
+ // Now create packs
+ myvals.resize( 1, myargs.size() );
+ mypack.resize( myargs.size(), 0 );
+ mypack.setValIndex( 0 );
+
// Create the value
addValueWithDerivatives(); setNotPeriodic();
}
@@ -100,8 +108,8 @@ Target::~Target(){
}
void Target::calculate(){
- double r=target->calculate( getArguments(), false ); setValue(r);
- for(unsigned i=0;i<getNumberOfArguments();i++) setDerivative( i, target->getArgumentDerivative(i) );
+ mypack.clear(); double r=target->calculate( getArguments(), mypack, false ); setValue(r);
+ for(unsigned i=0;i<getNumberOfArguments();i++) setDerivative( i, mypack.getArgumentDerivative(i) );
}
}
diff --git a/src/manyrestraints/ManyRestraintsBase.cpp b/src/manyrestraints/ManyRestraintsBase.cpp
index f2bc0b9981e9042c25e2bc5d9afe0016ee2a1fc6..bc7a999f3a1019589a0983a87218a9ff2af4259c 100644
--- a/src/manyrestraints/ManyRestraintsBase.cpp
+++ b/src/manyrestraints/ManyRestraintsBase.cpp
@@ -65,7 +65,7 @@ void ManyRestraintsBase::doJobsRequiredBeforeTaskList(){
ActionWithValue::clearDerivatives();
}
-void ManyRestraintsBase::transformBridgedDerivatives( const unsigned& current, vesselbase::MultiValue& invals, vesselbase::MultiValue& outvals ){
+void ManyRestraintsBase::transformBridgedDerivatives( const unsigned& current, MultiValue& invals, MultiValue& outvals ) const {
outvals.setValue( 0, invals.get(0) );
// Get the potential
diff --git a/src/manyrestraints/ManyRestraintsBase.h b/src/manyrestraints/ManyRestraintsBase.h
index 468e01508690def8f0b60a05d0e74e2f65fea503..2cf2b38eefb2f62c2c17c6112a43a2fc0e9dc79d 100644
--- a/src/manyrestraints/ManyRestraintsBase.h
+++ b/src/manyrestraints/ManyRestraintsBase.h
@@ -53,13 +53,13 @@ public:
/// Do jobs required before tasks are undertaken
void doJobsRequiredBeforeTaskList();
/// This actually does the calculation
- void transformBridgedDerivatives( const unsigned& current, vesselbase::MultiValue& invals, vesselbase::MultiValue& outvals );
+ void transformBridgedDerivatives( const unsigned& current, MultiValue& invals, MultiValue& outvals ) const ;
/// Calculate the potential
- virtual double calcPotential( const double& val, double& df )=0;
+ virtual double calcPotential( const double& val, double& df ) const=0;
// Calculate does nothing
void calculate(){};
/// This should never be called
- void performTask( const unsigned& , const unsigned& , vesselbase::MultiValue& ){ plumed_error(); }
+ void performTask( const unsigned& , const unsigned& , MultiValue& ) const { plumed_error(); }
/// Deactivate task now does nothing
void apply();
void applyBridgeForces( const std::vector<double>& bb ){ plumed_assert( bb.size()==0 ); }
diff --git a/src/manyrestraints/UWalls.cpp b/src/manyrestraints/UWalls.cpp
index 7292a7872853f4776c57393a961ec21865e0dbee..0fda53ec788da377f11c9e87e0cd161a1544b6c3 100644
--- a/src/manyrestraints/UWalls.cpp
+++ b/src/manyrestraints/UWalls.cpp
@@ -69,7 +69,7 @@ private:
public:
static void registerKeywords( Keywords& keys );
UWalls( const ActionOptions& );
- double calcPotential( const double& val, double& df );
+ double calcPotential( const double& val, double& df ) const ;
};
PLUMED_REGISTER_ACTION(UWalls,"UWALLS")
@@ -95,7 +95,7 @@ ManyRestraintsBase(ao)
checkRead();
}
-double UWalls::calcPotential( const double& val, double& df ){
+double UWalls::calcPotential( const double& val, double& df ) const {
double uscale = (val - at + offset)/eps;
if( uscale > 0. ){
double power = pow( uscale, exp );
diff --git a/src/mapping/Mapping.cpp b/src/mapping/Mapping.cpp
index 25fc1a0dfd98f3849048d80f3d8c8d231bd0b312..be3bc0aadf27ce7d8a1acf4b19184e0d02467d4b 100644
--- a/src/mapping/Mapping.cpp
+++ b/src/mapping/Mapping.cpp
@@ -22,6 +22,7 @@
#include "Mapping.h"
#include "vesselbase/Vessel.h"
#include "reference/MetricRegister.h"
+#include "reference/ReferenceAtoms.h"
#include "tools/PDB.h"
#include "tools/Matrix.h"
#include "core/PlumedMain.h"
@@ -106,10 +107,10 @@ ActionWithVessel(ao)
interpretArgumentList( args, req_args ); requestArguments( req_args );
// Duplicate all frames (duplicates are used by sketch-map)
mymap->duplicateFrameList();
- fframes.resize( 2*nfram, 0.0 ); dfframes.resize( 2*nfram, 0.0 );
+ // fframes.resize( 2*nfram, 0.0 ); dfframes.resize( 2*nfram, 0.0 );
plumed_assert( !mymap->mappingNeedsSetup() );
// Resize all derivative arrays
- mymap->setNumberOfAtomsAndArguments( atoms.size(), args.size() );
+ // mymap->setNumberOfAtomsAndArguments( atoms.size(), args.size() );
// Resize forces array
if( getNumberOfAtoms()>0 ){
forcesToApply.resize( 3*getNumberOfAtoms() + 9 + getNumberOfArguments() );
@@ -138,10 +139,10 @@ void Mapping::prepare(){
mymap->duplicateFrameList();
// Get the number of frames in the path
unsigned nfram=getNumberOfReferencePoints();
- fframes.resize( 2*nfram, 0.0 ); dfframes.resize( 2*nfram, 0.0 );
+ // fframes.resize( 2*nfram, 0.0 ); dfframes.resize( 2*nfram, 0.0 );
plumed_assert( !mymap->mappingNeedsSetup() );
// Resize all derivative arrays
- mymap->setNumberOfAtomsAndArguments( atoms.size(), args.size() );
+ // mymap->setNumberOfAtomsAndArguments( atoms.size(), args.size() );
// Resize forces array
forcesToApply.resize( 3*getNumberOfAtoms() + 9 + getNumberOfArguments() );
}
@@ -165,12 +166,28 @@ std::string Mapping::getArgumentName( unsigned& iarg ){
return "pos" + atnum + "z";
}
-double Mapping::calculateDistanceFunction( const unsigned& ifunc, const bool& squared ){
+void Mapping::finishPackSetup( const unsigned& ifunc, ReferenceValuePack& mypack ) const {
+ ReferenceConfiguration* myref=mymap->getFrame(ifunc); mypack.setValIndex(0);
+ unsigned nargs2=myref->getNumberOfReferenceArguments(); unsigned nat2=myref->getNumberOfReferencePositions();
+ if( mypack.getNumberOfAtoms()!=nat2 || mypack.getNumberOfArguments()!=nargs2 ) mypack.resize( nargs2, nat2 );
+ if( nat2>0 ){
+ ReferenceAtoms* myat2=dynamic_cast<ReferenceAtoms*>( myref ); plumed_dbg_assert( myat2 );
+ for(unsigned i=0;i<nat2;++i) mypack.setAtomIndex( i, myat2->getAtomIndex(i) );
+ }
+}
+
+double Mapping::calculateDistanceFunction( const unsigned& ifunc, ReferenceValuePack& myder, const bool& squared ) const {
// Calculate the distance
- double dd = mymap->calcDistanceFromConfiguration( ifunc, getPositions(), getPbc(), getArguments(), squared );
+ double dd = mymap->calcDistanceFromConfiguration( ifunc, getPositions(), getPbc(), getArguments(), myder, squared );
// Transform distance by whatever
- fframes[ifunc]=transformHD( dd, dfframes[ifunc] );
- return fframes[ifunc];
+ double df, ff=transformHD( dd, df ); myder.scaleAllDerivatives( df );
+ // And the virial
+ if( !myder.virialWasSet() ){
+ Tensor tvir; tvir.zero();
+ for(unsigned i=0;i<myder.getNumberOfAtoms();++i) tvir +=-1.0*Tensor( getPosition( myder.getAtomIndex(i) ), myder.getAtomDerivative(i) );
+ myder.addBoxDerivatives( tvir );
+ }
+ return ff;
}
void Mapping::calculateNumericalDerivatives( ActionWithValue* a ){
@@ -189,36 +206,6 @@ void Mapping::calculateNumericalDerivatives( ActionWithValue* a ){
}
}
-void Mapping::transferDerivatives( const unsigned& ider, const unsigned& fno, const unsigned& cur, vesselbase::MultiValue& myvals ){
- if( !derivativesAreRequired() ) return;
-
- unsigned frameno=fno*getNumberOfReferencePoints() + cur;
- for(unsigned i=0;i<getNumberOfArguments();++i) myvals.addDerivative( ider, i, dfframes[frameno]*mymap->getArgumentDerivative(cur,i) );
-
- if( getNumberOfAtoms()>0 ){
- Vector ader; Tensor tmpvir; tmpvir.zero();
- unsigned n=getNumberOfArguments();
- for(unsigned i=0;i<getNumberOfAtoms();++i){
- ader=mymap->getAtomDerivatives( cur, i );
- myvals.addDerivative( ider, n, dfframes[frameno]*ader[0] ); n++;
- myvals.addDerivative( ider, n, dfframes[frameno]*ader[1] ); n++;
- myvals.addDerivative( ider, n, dfframes[frameno]*ader[2] ); n++;
- tmpvir += -1.0*Tensor( getPosition(i), ader );
- }
- Tensor vir;
- if( !mymap->getVirial( cur, vir ) ) vir=tmpvir;
- myvals.addDerivative( ider, n, dfframes[frameno]*vir(0,0) ); n++;
- myvals.addDerivative( ider, n, dfframes[frameno]*vir(0,1) ); n++;
- myvals.addDerivative( ider, n, dfframes[frameno]*vir(0,2) ); n++;
- myvals.addDerivative( ider, n, dfframes[frameno]*vir(1,0) ); n++;
- myvals.addDerivative( ider, n, dfframes[frameno]*vir(1,1) ); n++;
- myvals.addDerivative( ider, n, dfframes[frameno]*vir(1,2) ); n++;
- myvals.addDerivative( ider, n, dfframes[frameno]*vir(2,0) ); n++;
- myvals.addDerivative( ider, n, dfframes[frameno]*vir(2,1) ); n++;
- myvals.addDerivative( ider, n, dfframes[frameno]*vir(2,2) );
- }
-}
-
void Mapping::apply(){
if( getForcesFromVessels( forcesToApply ) ){
addForcesOnArguments( forcesToApply );
diff --git a/src/mapping/Mapping.h b/src/mapping/Mapping.h
index 25ab9dad96eb7d0dcb5d3eeb09dc8c8133338d1e..f02769fed90d301fbeb33f43c24f90b2cf25ff45 100644
--- a/src/mapping/Mapping.h
+++ b/src/mapping/Mapping.h
@@ -53,14 +53,14 @@ protected:
std::vector<double> fframes;
/// Get the number of frames in the path
unsigned getNumberOfReferencePoints() const ;
+/// Finish the setup of the referenceValuePack by transfering atoms and args
+ void finishPackSetup( const unsigned& ifunc, ReferenceValuePack& mypack ) const ;
/// Calculate the value of the distance from the ith frame
- double calculateDistanceFunction( const unsigned& ifunc, const bool& squared );
+ double calculateDistanceFunction( const unsigned& ifunc, ReferenceValuePack& myder, const bool& squared ) const ;
/// Store the distance function
void storeDistanceFunction( const unsigned& ifunc );
/// Get the value of the weight
double getWeight( const unsigned& weight ) const ;
-/// Transfer the derivatives to a MultiValue object
- void transferDerivatives( const unsigned& ider, const unsigned& fno, const unsigned& cur, vesselbase::MultiValue& myvals );
public:
static void registerKeywords( Keywords& keys );
Mapping(const ActionOptions&);
@@ -77,7 +77,7 @@ public:
/// Get the value of lambda for paths and property maps
virtual double getLambda();
/// This does the transformation of the distance by whatever function is required
- virtual double transformHD( const double& dist, double& df )=0;
+ virtual double transformHD( const double& dist, double& df ) const=0;
/// Get the number of properties we are projecting onto
unsigned getNumberOfProperties() const ;
/// Get the name of the ith property we are projecting
diff --git a/src/mapping/PathBase.cpp b/src/mapping/PathBase.cpp
index 789edae8e310fcb4b3577b8aae1cf93eefa3b7db..00ae362019c90604717dd9a65ce2ec443d066b67 100644
--- a/src/mapping/PathBase.cpp
+++ b/src/mapping/PathBase.cpp
@@ -55,16 +55,18 @@ void PathBase::calculate(){
runAllTasks();
}
-void PathBase::performTask( const unsigned& task_index, const unsigned& current, vesselbase::MultiValue& myvals ){
+void PathBase::performTask( const unsigned& task_index, const unsigned& current, MultiValue& myvals ) const {
+ // This builds a pack to hold the derivatives
+ ReferenceValuePack mypack( getNumberOfArguments(), getNumberOfAtoms(), myvals );
+ finishPackSetup( current, mypack );
// Calculate the distance from the frame
- double val=calculateDistanceFunction( current, true );
+ double val=calculateDistanceFunction( current, mypack, true );
// Put the element value in element zero
myvals.setValue( 0, val ); myvals.setValue( 1, 1.0 );
- transferDerivatives( 0, 0, current, myvals );
return;
}
-double PathBase::transformHD( const double& dist, double& df ){
+double PathBase::transformHD( const double& dist, double& df ) const {
double val = exp( -dist*lambda );
df = -lambda*val;
return val;
diff --git a/src/mapping/PathBase.h b/src/mapping/PathBase.h
index 238a80d07b074e4984d4a02fb13733b31851a296..da36cd41ab7845a45d701e0bb044661362ad815c 100644
--- a/src/mapping/PathBase.h
+++ b/src/mapping/PathBase.h
@@ -35,8 +35,8 @@ public:
PathBase(const ActionOptions&);
double getLambda();
void calculate();
- void performTask( const unsigned& , const unsigned& , vesselbase::MultiValue& );
- double transformHD( const double& dist, double& df );
+ void performTask( const unsigned& , const unsigned& , MultiValue& ) const ;
+ double transformHD( const double& dist, double& df ) const ;
};
}
diff --git a/src/mapping/SpathVessel.cpp b/src/mapping/SpathVessel.cpp
index 2ca3ad40c908c680c5e8dcf95a44f568d695f24a..0bf7045814701328b0d03e3716e69b4c8a897cb8 100644
--- a/src/mapping/SpathVessel.cpp
+++ b/src/mapping/SpathVessel.cpp
@@ -38,7 +38,7 @@ public:
SpathVessel( const vesselbase::VesselOptions& da );
std::string function_description();
void prepare();
- bool calculate( const unsigned& current, vesselbase::MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_index ) const ;
+ bool calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_index ) const ;
};
PLUMED_REGISTER_VESSEL(SpathVessel,"SPATH")
@@ -74,7 +74,7 @@ void SpathVessel::prepare(){
foundoneclose=false;
}
-bool SpathVessel::calculate( const unsigned& current, vesselbase::MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_index ) const {
+bool SpathVessel::calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_index ) const {
double pp=mymap->getPropertyValue( current, mycoordnumber ), weight=myvals.get(0);
if( weight<getTolerance() ) return false;
buffer[bufstart] += weight*pp; buffer[bufstart+1+nderiv] += weight;
diff --git a/src/multicolvar/ActionVolume.cpp b/src/multicolvar/ActionVolume.cpp
index 2ef73aec2bce36a17fca94e101c9071ea730cddf..cdfbde1152a27f37bdcf846fa99fdc24f4de0979 100644
--- a/src/multicolvar/ActionVolume.cpp
+++ b/src/multicolvar/ActionVolume.cpp
@@ -49,9 +49,7 @@ VolumeGradientBase(ao)
log.printf(" calculating %s inside region of insterest\n",functype.c_str() );
parseFlag("OUTSIDE",not_in); parse("SIGMA",sigma);
- bead.isNotPeriodic();
- std::string kerneltype; parse("KERNEL",kerneltype);
- bead.setKernelType( kerneltype );
+ parse("KERNEL",kerneltype);
if( getPntrToMultiColvar()->isDensity() ){
std::string input;
@@ -61,11 +59,11 @@ VolumeGradientBase(ao)
}
}
-void ActionVolume::calculateAllVolumes( const unsigned& curr, vesselbase::MultiValue& outvals ){
+void ActionVolume::calculateAllVolumes( const unsigned& curr, MultiValue& outvals ) const {
Vector catom_pos=getPntrToMultiColvar()->getCentralAtomPos( curr );
double weight; Vector wdf; Tensor vir; std::vector<Vector> refders( getNumberOfAtoms() );
- weight=calculateNumberInside( catom_pos, bead, wdf, vir, refders );
+ weight=calculateNumberInside( catom_pos, wdf, vir, refders );
if( not_in ){
weight = 1.0 - weight; wdf *= -1.; vir *=-1;
for(unsigned i=0;i<refders.size();++i) refders[i]*=-1;
diff --git a/src/multicolvar/ActionVolume.h b/src/multicolvar/ActionVolume.h
index 10a459b341b19d594fd36fc20c409b03a0d93f1a..e1868802e3c7800523b696bb486b40c2d678c7cc 100644
--- a/src/multicolvar/ActionVolume.h
+++ b/src/multicolvar/ActionVolume.h
@@ -43,18 +43,19 @@ private:
double sigma;
/// Are we interested in the area outside the colvar
bool not_in;
-/// The bead for the histogram
- HistogramBead bead;
+/// The kernel type for this histogram
+ std::string kerneltype;
protected:
double getSigma() const ;
+ std::string getKernelType() const ;
public:
static void registerKeywords( Keywords& keys );
ActionVolume(const ActionOptions&);
/// Get the number of quantities that are calculated each time
virtual unsigned getNumberOfQuantities();
/// Calculate whats in the volume
- void calculateAllVolumes( const unsigned& curr, vesselbase::MultiValue& outvals );
- virtual double calculateNumberInside( const Vector& cpos, HistogramBead& bead, Vector& derivatives, Tensor& vir, std::vector<Vector>& refders )=0;
+ void calculateAllVolumes( const unsigned& curr, MultiValue& outvals ) const ;
+ virtual double calculateNumberInside( const Vector& cpos, Vector& derivatives, Tensor& vir, std::vector<Vector>& refders ) const=0;
unsigned getCentralAtomElementIndex();
};
@@ -68,6 +69,11 @@ double ActionVolume::getSigma() const {
return sigma;
}
+inline
+std::string ActionVolume::getKernelType() const {
+ return kerneltype;
+}
+
inline
unsigned ActionVolume::getCentralAtomElementIndex(){
return 1;
diff --git a/src/multicolvar/AdjacencyMatrixAction.cpp b/src/multicolvar/AdjacencyMatrixAction.cpp
index 0b2664b6467f8a0b8f802dc66befe025a004554b..159cacfc7dac7163bd2a7d073bd5c33018bf95aa 100644
--- a/src/multicolvar/AdjacencyMatrixAction.cpp
+++ b/src/multicolvar/AdjacencyMatrixAction.cpp
@@ -83,7 +83,6 @@ tmpdf(1)
// Build active elements array
for(unsigned i=0;i<getFullNumberOfTasks();++i) active_elements.addIndexToList( i );
- active_elements.setupMPICommunication( comm );
// Find the largest sf cutoff
double sfmax=switchingFunction(0,0).get_dmax();
@@ -101,6 +100,8 @@ tmpdf(1)
Keywords keys; AdjacencyMatrixVessel::registerKeywords( keys );
vesselbase::VesselOptions da2(da,keys);
mat = new AdjacencyMatrixVessel(da2);
+ // Set a cutoff for clustering
+ mat->setHardCutoffOnWeight( getTolerance() );
// Add the vessel to the base
addVessel( mat );
// And resize everything
@@ -111,20 +112,20 @@ void AdjacencyMatrixAction::doJobsRequiredBeforeTaskList(){
// Do jobs required by ActionWithVessel
ActionWithVessel::doJobsRequiredBeforeTaskList();
// Make it possible only to go through loop once
- gathered=false; active_elements.deactivateAll();
+ gathered=false;
// Dont calculate derivatives on first loop
if( usingLowMem() ) dertime=false;
else dertime=true;
}
-void AdjacencyMatrixAction::calculateWeight( AtomValuePack& myatoms ){
+void AdjacencyMatrixAction::calculateWeight( AtomValuePack& myatoms ) const {
Vector distance = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
double dfunc, sw = switchingFunction( getBaseColvarNumber( myatoms.getIndex(0) ),getBaseColvarNumber( myatoms.getIndex(1) ) ).calculate( distance.modulo(), dfunc );
myatoms.setValue(0,sw);
}
-double AdjacencyMatrixAction::compute( const unsigned& tindex, AtomValuePack& myatoms ){
- active_elements.activate( tindex );
+double AdjacencyMatrixAction::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
+// active_elements.activate( tindex );
double f_dot, dot_df;
unsigned ncomp=getBaseMultiColvar(0)->getNumberOfQuantities();
@@ -160,9 +161,9 @@ double AdjacencyMatrixAction::compute( const unsigned& tindex, AtomValuePack& my
for(unsigned k=2;k<orient0.size();++k){
orient0[k]*=sw*dot_df; orient1[k]*=sw*dot_df;
}
- vesselbase::MultiValue myder0(0,0); getVectorDerivatives( myatoms.getIndex(0), true, myder0 );
+ MultiValue myder0(0,0); getVectorDerivatives( myatoms.getIndex(0), true, myder0 );
mergeVectorDerivatives( 1, 2, orient1.size(), myatoms.getIndex(0), orient1, myder0, myatoms );
- vesselbase::MultiValue myder1(0,0); getVectorDerivatives( myatoms.getIndex(1), true, myder1 );
+ MultiValue myder1(0,0); getVectorDerivatives( myatoms.getIndex(1), true, myder1 );
mergeVectorDerivatives( 1, 2, orient0.size(), myatoms.getIndex(1), orient0, myder1, myatoms );
}
}
@@ -171,8 +172,13 @@ double AdjacencyMatrixAction::compute( const unsigned& tindex, AtomValuePack& my
void AdjacencyMatrixAction::retrieveMatrix( Matrix<double>& mymatrix ){
// Gather active elements in matrix
- if(!gathered) active_elements.mpi_gatherActiveMembers( comm );
- gathered=true;
+ if(!gathered){
+ active_elements.deactivateAll();
+ for(unsigned i=0;i<getFullNumberOfTasks();++i){
+ if( mat->storedValueIsActive(i) ) active_elements.activate(i);
+ }
+ active_elements.updateActiveMembers(); gathered=true;
+ }
std::vector<unsigned> myatoms(2); std::vector<double> vals(2);
for(unsigned i=0;i<active_elements.getNumberActive();++i){
@@ -187,8 +193,13 @@ void AdjacencyMatrixAction::retrieveAdjacencyLists( std::vector<unsigned>& nneig
plumed_dbg_assert( nneigh.size()==getFullNumberOfBaseTasks() && adj_list.nrows()==getFullNumberOfBaseTasks() &&
adj_list.ncols()==getFullNumberOfBaseTasks() );
// Gather active elements in matrix
- if(!gathered) active_elements.mpi_gatherActiveMembers( comm );
- gathered=true;
+ if(!gathered){
+ active_elements.deactivateAll();
+ for(unsigned i=0;i<getFullNumberOfTasks();++i){
+ if( mat->storedValueIsActive(i) ) active_elements.activate(i);
+ }
+ active_elements.updateActiveMembers(); gathered=true;
+ }
// Currently everything has zero neighbors
for(unsigned i=0;i<nneigh.size();++i) nneigh[i]=0;
diff --git a/src/multicolvar/AdjacencyMatrixAction.h b/src/multicolvar/AdjacencyMatrixAction.h
index e397e89c0a62a1f684695a9c682f7f10c4d394ad..39a96a3a672eb180ec5bd02cb94e916d4fe78dec 100644
--- a/src/multicolvar/AdjacencyMatrixAction.h
+++ b/src/multicolvar/AdjacencyMatrixAction.h
@@ -64,8 +64,8 @@ protected:
public:
static void registerKeywords( Keywords& keys );
AdjacencyMatrixAction(const ActionOptions&);
- double compute( const unsigned& tindex, AtomValuePack& myatoms );
- void calculateWeight( AtomValuePack& myatoms );
+ double compute( const unsigned& tindex, AtomValuePack& myatoms ) const ;
+ void calculateWeight( AtomValuePack& myatoms ) const ;
void doJobsRequiredBeforeTaskList();
/// Finish the calculation
virtual void completeCalculation()=0;
diff --git a/src/multicolvar/AlphaBeta.cpp b/src/multicolvar/AlphaBeta.cpp
index 085fd49c123d609300a0976aeb85ccc9365a2c5f..b203e5b55d79f26b5799fcf7db668aa977869719 100644
--- a/src/multicolvar/AlphaBeta.cpp
+++ b/src/multicolvar/AlphaBeta.cpp
@@ -98,7 +98,7 @@ private:
public:
static void registerKeywords( Keywords& keys );
AlphaBeta(const ActionOptions&);
- virtual double compute( const unsigned& tindex, AtomValuePack& myatoms );
+ virtual double compute( const unsigned& tindex, AtomValuePack& myatoms ) const ;
bool isPeriodic(){ return false; }
};
@@ -148,7 +148,7 @@ PLUMED_MULTICOLVAR_INIT(ao)
checkRead();
}
-double AlphaBeta::compute( const unsigned& tindex, AtomValuePack& myatoms ){
+double AlphaBeta::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
Vector d0,d1,d2;
d0=getSeparation(myatoms.getPosition(1),myatoms.getPosition(0));
d1=getSeparation(myatoms.getPosition(2),myatoms.getPosition(1));
diff --git a/src/multicolvar/Angles.cpp b/src/multicolvar/Angles.cpp
index 5d412c9fe8329f5751e6fcfc26c5816d9ce22edd..cb640f9974b901809fb156386a8f71a22f519415 100644
--- a/src/multicolvar/Angles.cpp
+++ b/src/multicolvar/Angles.cpp
@@ -87,16 +87,15 @@ class Angles : public MultiColvar {
private:
bool use_sf;
double rcut2_1, rcut2_2;
- Vector dij, dik;
SwitchingFunction sf1;
SwitchingFunction sf2;
public:
static void registerKeywords( Keywords& keys );
Angles(const ActionOptions&);
/// Updates neighbor list
- virtual double compute( const unsigned& tindex, AtomValuePack& );
+ virtual double compute( const unsigned& tindex, AtomValuePack& ) const ;
/// Returns the number of coordinates of the field
- void calculateWeight( AtomValuePack& );
+ void calculateWeight( AtomValuePack& ) const ;
bool isPeriodic(){ return false; }
};
@@ -169,10 +168,10 @@ use_sf(false)
setAtomsForCentralAtom( catom_ind );
}
-void Angles::calculateWeight( AtomValuePack& myatoms ){
- dij=getSeparation( myatoms.getPosition(0), myatoms.getPosition(2) );
- dik=getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
+void Angles::calculateWeight( AtomValuePack& myatoms ) const {
if(!use_sf){ myatoms.setValue( 0, 1.0 ); return; }
+ Vector dij=getSeparation( myatoms.getPosition(0), myatoms.getPosition(2) );
+ Vector dik=getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
double w1, w2, dw1, dw2, wtot;
double ldij = dij.modulo2(), ldik = dik.modulo2();
@@ -192,7 +191,10 @@ void Angles::calculateWeight( AtomValuePack& myatoms ){
myatoms.addBoxDerivatives( 0, (-dw1)*Tensor(dij,dij) + (-dw2)*Tensor(dik,dik) );
}
-double Angles::compute( const unsigned& tindex, AtomValuePack& myatoms ){
+double Angles::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
+ Vector dij=getSeparation( myatoms.getPosition(0), myatoms.getPosition(2) );
+ Vector dik=getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
+
Vector ddij,ddik; PLMD::Angle a;
double angle=a.compute(dij,dik,ddij,ddik);
diff --git a/src/multicolvar/AtomValuePack.cpp b/src/multicolvar/AtomValuePack.cpp
index a42dab9a8946ec6a60b8ffd7faf057f0e43689dd..651ee3de12da33c0b0902f92d83b4367597b5c0b 100644
--- a/src/multicolvar/AtomValuePack.cpp
+++ b/src/multicolvar/AtomValuePack.cpp
@@ -25,11 +25,11 @@
namespace PLMD {
namespace multicolvar {
-AtomValuePack::AtomValuePack( vesselbase::MultiValue& vals, MultiColvarBase* mcolv ):
+AtomValuePack::AtomValuePack( MultiValue& vals, MultiColvarBase const * mcolv ):
myvals(vals),
mycolv(mcolv)
{
- indices.resize( mycolv->getNumberOfDerivatives() );
+ indices.resize( vals.getNumberOfDerivatives() ); // mycolv->getNumberOfDerivatives() );
}
void AtomValuePack::updateUsingIndices(){
diff --git a/src/multicolvar/AtomValuePack.h b/src/multicolvar/AtomValuePack.h
index ffe0bbe0f6be42e639b764ba71ca6cd2077b0b82..c7b7d9a8a8a00920f12c8ec658cec63c2e9f27e4 100644
--- a/src/multicolvar/AtomValuePack.h
+++ b/src/multicolvar/AtomValuePack.h
@@ -22,7 +22,7 @@
#ifndef __PLUMED_multicolvar_AtomValuePack_h
#define __PLUMED_multicolvar_AtomValuePack_h
-#include "vesselbase/MultiValue.h"
+#include "tools/MultiValue.h"
#include "MultiColvarBase.h"
namespace PLMD {
@@ -33,15 +33,15 @@ class CatomPack;
class AtomValuePack {
private:
/// Copy of the values that we are adding to
- vesselbase::MultiValue& myvals;
+ MultiValue& myvals;
/// Copy of the underlying multicolvar
- MultiColvarBase* mycolv;
+ MultiColvarBase const * mycolv;
/// Number of atoms at the moment
unsigned natoms;
/// Atom indices
std::vector<unsigned> indices;
public:
- AtomValuePack( vesselbase::MultiValue& vals, MultiColvarBase* mcolv );
+ AtomValuePack( MultiValue& vals, MultiColvarBase const * mcolv );
/// Set the number of atoms
void setNumberOfAtoms( const unsigned& );
/// Set the index for one of the atoms
@@ -50,6 +50,8 @@ public:
unsigned getIndex( const unsigned& j ) const ;
///
unsigned getNumberOfAtoms() const ;
+///
+ unsigned getNumberOfDerivatives() const ;
/// Get the position of the ith atom
Vector getPosition( const unsigned& );
///
@@ -71,7 +73,7 @@ public:
///
void addComDerivatives( const unsigned& , const Vector& , CatomPack& );
///
- vesselbase::MultiValue& getUnderlyingMultiValue();
+ MultiValue& getUnderlyingMultiValue();
};
inline
@@ -84,6 +86,11 @@ unsigned AtomValuePack::getNumberOfAtoms() const {
return natoms;
}
+inline
+unsigned AtomValuePack::getNumberOfDerivatives() const {
+ return myvals.getNumberOfDerivatives();
+}
+
inline
void AtomValuePack::setIndex( const unsigned& j, const unsigned& ind ){
plumed_dbg_assert( j<natoms ); indices[j]=ind;
@@ -141,7 +148,7 @@ void AtomValuePack::updateDynamicList(){
}
inline
-vesselbase::MultiValue& AtomValuePack::getUnderlyingMultiValue(){
+MultiValue& AtomValuePack::getUnderlyingMultiValue(){
return myvals;
}
diff --git a/src/multicolvar/Bridge.cpp b/src/multicolvar/Bridge.cpp
index 21ed24ecdb48d138c053e2dde56c5dbee1fb5c8c..9acb1484c669b0aee46ecd047929e66da042365e 100644
--- a/src/multicolvar/Bridge.cpp
+++ b/src/multicolvar/Bridge.cpp
@@ -68,8 +68,8 @@ public:
static void registerKeywords( Keywords& keys );
Bridge(const ActionOptions&);
// active methods:
- virtual double compute( const unsigned& tindex, AtomValuePack& myatoms );
- void calculateWeight( AtomValuePack& myatoms );
+ virtual double compute( const unsigned& tindex, AtomValuePack& myatoms ) const ;
+ void calculateWeight( AtomValuePack& myatoms ) const ;
bool isPeriodic(){ return false; }
};
@@ -138,7 +138,7 @@ PLUMED_MULTICOLVAR_INIT(ao)
checkRead();
}
-void Bridge::calculateWeight( AtomValuePack& myatoms ){
+void Bridge::calculateWeight( AtomValuePack& myatoms ) const {
Vector dij=getSeparation( myatoms.getPosition(0), myatoms.getPosition(2) );
double ldij = dij.modulo2();
if( ldij>rcut2 ) { myatoms.setValue(0,0); return; }
@@ -150,7 +150,7 @@ void Bridge::calculateWeight( AtomValuePack& myatoms ){
myatoms.addBoxDerivatives( 0, (-dw)*Tensor(dij,dij) );
}
-double Bridge::compute( const unsigned& tindex, AtomValuePack& myatoms ){
+double Bridge::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
Vector dik=getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
double dw, w=sf1.calculateSqr( dik.modulo2(), dw );
diff --git a/src/multicolvar/BridgedMultiColvarFunction.cpp b/src/multicolvar/BridgedMultiColvarFunction.cpp
index 6b82ce32298ac15bf99387eef3e65c8f9095b699..83e032e08531afcf848ebc7d317bd1f3fd3df28e 100644
--- a/src/multicolvar/BridgedMultiColvarFunction.cpp
+++ b/src/multicolvar/BridgedMultiColvarFunction.cpp
@@ -51,7 +51,7 @@ MultiColvarBase(ao)
for(unsigned i=0;i<mycolv->getFullNumberOfTasks();++i) addTaskToList( mycolv->getTaskCode(i) );
}
-void BridgedMultiColvarFunction::transformBridgedDerivatives( const unsigned& current, vesselbase::MultiValue& invals, vesselbase::MultiValue& outvals ){
+void BridgedMultiColvarFunction::transformBridgedDerivatives( const unsigned& current, MultiValue& invals, MultiValue& outvals ) const {
completeTask( current, invals, outvals );
// Now update the outvals derivatives lists
@@ -66,8 +66,8 @@ void BridgedMultiColvarFunction::transformBridgedDerivatives( const unsigned& cu
outvals.sortActiveList();
}
-void BridgedMultiColvarFunction::performTask( const unsigned& taskIndex, const unsigned& current, vesselbase::MultiValue& myvals ){
- vesselbase::MultiValue invals( mycolv->getNumberOfQuantities(), mycolv->getNumberOfDerivatives() );
+void BridgedMultiColvarFunction::performTask( const unsigned& taskIndex, const unsigned& current, MultiValue& myvals ) const {
+ MultiValue invals( mycolv->getNumberOfQuantities(), mycolv->getNumberOfDerivatives() );
mycolv->performTask( taskIndex, current, invals );
transformBridgedDerivatives( taskIndex, invals, myvals );
}
diff --git a/src/multicolvar/BridgedMultiColvarFunction.h b/src/multicolvar/BridgedMultiColvarFunction.h
index 192c16c2272425d7544c6556d04456f180975076..e59d7dcce6ab1c4768c1c4d2aa78e1f42e330aff 100644
--- a/src/multicolvar/BridgedMultiColvarFunction.h
+++ b/src/multicolvar/BridgedMultiColvarFunction.h
@@ -64,9 +64,9 @@ public:
void deactivate_task( const unsigned& taskno );
void calculate(){}
/// This does the task
- void transformBridgedDerivatives( const unsigned& current, vesselbase::MultiValue& invals, vesselbase::MultiValue& outvals );
- void performTask( const unsigned& , const unsigned& , vesselbase::MultiValue& );
- virtual void completeTask( const unsigned& curr, vesselbase::MultiValue& invals, vesselbase::MultiValue& outvals )=0;
+ void transformBridgedDerivatives( const unsigned& current, MultiValue& invals, MultiValue& outvals ) const ;
+ void performTask( const unsigned& , const unsigned& , MultiValue& ) const ;
+ virtual void completeTask( const unsigned& curr, MultiValue& invals, MultiValue& outvals ) const=0;
/// Get the central atom position
Vector retrieveCentralAtomPos();
/// We need our own calculate numerical derivatives here
@@ -76,9 +76,9 @@ public:
bool isCurrentlyActive( const unsigned& );
/// This should not be called
Vector calculateCentralAtomPosition(){ plumed_error(); }
- double compute( const unsigned& tindex, AtomValuePack& myvals ){ plumed_error(); }
- Vector getPositionOfAtomForLinkCells( const unsigned& iatom ){ plumed_error(); }
- void updateActiveAtoms(){ plumed_error(); }
+ double compute( const unsigned& tindex, AtomValuePack& myvals ) const { plumed_error(); }
+ Vector getPositionOfAtomForLinkCells( const unsigned& iatom ) const ;
+ void updateActiveAtoms( AtomValuePack& myatoms ) const { plumed_error(); }
void getIndexList( const unsigned& ntotal, const unsigned& jstore, const unsigned& maxder, std::vector<unsigned>& indices );
void applyBridgeForces( const std::vector<double>& bb );
};
@@ -104,7 +104,7 @@ unsigned BridgedMultiColvarFunction::getSizeOfAtomsWithDerivatives(){
}
inline
-Vector BridgedMultiColvarFunction::getPositionOfAtomForLinkCells( const unsigned& iatom ){
+Vector BridgedMultiColvarFunction::getPositionOfAtomForLinkCells( const unsigned& iatom ) const {
return mycolv->getPositionOfAtomForLinkCells(iatom);
}
diff --git a/src/multicolvar/CoordinationNumbers.cpp b/src/multicolvar/CoordinationNumbers.cpp
index 7545b773a21c6790e39173deaf3ed915b6a0f68f..1738046c487349af167300b3e7f3d21828494528 100644
--- a/src/multicolvar/CoordinationNumbers.cpp
+++ b/src/multicolvar/CoordinationNumbers.cpp
@@ -71,7 +71,7 @@ public:
static void registerKeywords( Keywords& keys );
CoordinationNumbers(const ActionOptions&);
// active methods:
- virtual double compute( const unsigned& tindex, AtomValuePack& myatoms );
+ virtual double compute( const unsigned& tindex, AtomValuePack& myatoms ) const ;
/// Returns the number of coordinates of the field
bool isPeriodic(){ return false; }
};
@@ -119,7 +119,7 @@ PLUMED_MULTICOLVAR_INIT(ao)
checkRead();
}
-double CoordinationNumbers::compute( const unsigned& tindex, AtomValuePack& myatoms ){
+double CoordinationNumbers::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
double value=0, dfunc; Vector distance;
// Calculate the coordination number
diff --git a/src/multicolvar/Density.cpp b/src/multicolvar/Density.cpp
index 149db001ce2133a266a32f20300bd8e97df39203..32794ad573b485dabf71cbe8119380793401e4de 100644
--- a/src/multicolvar/Density.cpp
+++ b/src/multicolvar/Density.cpp
@@ -54,10 +54,10 @@ public:
static void registerKeywords( Keywords& keys );
Density(const ActionOptions&);
// active methods:
- virtual double compute( const unsigned& tindex, AtomValuePack& myatoms );
+ virtual double compute( const unsigned& tindex, AtomValuePack& myatoms ) const ;
/// Returns the number of coordinates of the field
bool isPeriodic(){ return false; }
- bool isDensity(){ return true; }
+ bool isDensity() const { return true; }
bool hasDifferentiableOrientation() const { return true; }
// void addOrientationDerivativesToBase( const unsigned& iatom, const unsigned& jstore, const unsigned& base_cv_no,
// const std::vector<double>& weight, MultiColvarFunction* func ){}
@@ -81,7 +81,7 @@ PLUMED_MULTICOLVAR_INIT(ao)
checkRead();
}
-double Density::compute( const unsigned& tindex, AtomValuePack& myvals ){
+double Density::compute( const unsigned& tindex, AtomValuePack& myvals ) const {
return 1.0;
}
diff --git a/src/multicolvar/DihedralCorrelation.cpp b/src/multicolvar/DihedralCorrelation.cpp
index 579f737b51446ae510f2e84f47c6e7ed59ad91f2..14b72c4339c169c3b16acc0fcb2bb9764f313997 100644
--- a/src/multicolvar/DihedralCorrelation.cpp
+++ b/src/multicolvar/DihedralCorrelation.cpp
@@ -85,7 +85,7 @@ private:
public:
static void registerKeywords( Keywords& keys );
DihedralCorrelation(const ActionOptions&);
- virtual double compute( const unsigned& tindex, AtomValuePack& myatoms );
+ virtual double compute( const unsigned& tindex, AtomValuePack& myatoms ) const ;
bool isPeriodic(){ return false; }
};
@@ -117,7 +117,7 @@ PLUMED_MULTICOLVAR_INIT(ao)
checkRead();
}
-double DihedralCorrelation::compute( const unsigned& tindex, AtomValuePack& myatoms ){
+double DihedralCorrelation::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
Vector d10,d11,d12;
d10=getSeparation(myatoms.getPosition(1),myatoms.getPosition(0));
d11=getSeparation(myatoms.getPosition(2),myatoms.getPosition(1));
diff --git a/src/multicolvar/Distances.cpp b/src/multicolvar/Distances.cpp
index 207f74029b859ba83c5785e19e816c27a0393612..420fc3fb5d4ca6c636affab4da49e31b9f99cea5 100644
--- a/src/multicolvar/Distances.cpp
+++ b/src/multicolvar/Distances.cpp
@@ -82,7 +82,7 @@ public:
static void registerKeywords( Keywords& keys );
Distances(const ActionOptions&);
// active methods:
- virtual double compute( const unsigned& tindex, AtomValuePack& myatoms );
+ virtual double compute( const unsigned& tindex, AtomValuePack& myatoms ) const ;
/// Returns the number of coordinates of the field
bool isPeriodic(){ return false; }
};
@@ -138,7 +138,7 @@ PLUMED_MULTICOLVAR_INIT(ao)
}
}
-double Distances::compute( const unsigned& tindex, AtomValuePack& myatoms ){
+double Distances::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
Vector distance;
distance=getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
const double value=distance.modulo();
diff --git a/src/multicolvar/FilterBetween.cpp b/src/multicolvar/FilterBetween.cpp
index 2359c7b42b871ca85986c6bc2cf5a30f6922a678..e2e2822a7004346fa67694f5adafe81396b5c508 100644
--- a/src/multicolvar/FilterBetween.cpp
+++ b/src/multicolvar/FilterBetween.cpp
@@ -42,7 +42,7 @@ private:
public:
static void registerKeywords( Keywords& keys );
FilterBetween(const ActionOptions& ao);
- double applyFilter( const double& val, double& df );
+ double applyFilter( const double& val, double& df ) const ;
};
PLUMED_REGISTER_ACTION(FilterBetween,"MFILTER_BETWEEN")
@@ -86,7 +86,7 @@ MultiColvarFilter(ao)
checkRead();
}
-double FilterBetween::applyFilter( const double& val, double& df ){
+double FilterBetween::applyFilter( const double& val, double& df ) const {
double f = hb.calculate( val, df );
return f;
}
diff --git a/src/multicolvar/FilterLessThan.cpp b/src/multicolvar/FilterLessThan.cpp
index a9cb12ebb3f4edfa37336e76b6dce141a7f3a017..17acee9182a5435e064160d7296d93161f467539 100644
--- a/src/multicolvar/FilterLessThan.cpp
+++ b/src/multicolvar/FilterLessThan.cpp
@@ -42,7 +42,7 @@ private:
public:
static void registerKeywords( Keywords& keys );
FilterLess(const ActionOptions& ao);
- double applyFilter( const double& val, double& df );
+ double applyFilter( const double& val, double& df ) const ;
};
PLUMED_REGISTER_ACTION(FilterLess,"MFILTER_LESS")
@@ -79,7 +79,7 @@ MultiColvarFilter(ao)
checkRead();
}
-double FilterLess::applyFilter( const double& val, double& df ){
+double FilterLess::applyFilter( const double& val, double& df ) const {
double f = sf.calculate( val, df ); df*=val;
return f;
}
diff --git a/src/multicolvar/FilterMoreThan.cpp b/src/multicolvar/FilterMoreThan.cpp
index 7bf3c4b888950a107619b8c2cd3292d5f177cb67..d22fb6c968d58ab40a3953ee3b58b7d101af4ecc 100644
--- a/src/multicolvar/FilterMoreThan.cpp
+++ b/src/multicolvar/FilterMoreThan.cpp
@@ -42,7 +42,7 @@ private:
public:
static void registerKeywords( Keywords& keys );
FilterMore(const ActionOptions& ao);
- double applyFilter( const double& val, double& df );
+ double applyFilter( const double& val, double& df ) const ;
};
PLUMED_REGISTER_ACTION(FilterMore,"MFILTER_MORE")
@@ -79,7 +79,7 @@ MultiColvarFilter(ao)
checkRead();
}
-double FilterMore::applyFilter( const double& val, double& df ){
+double FilterMore::applyFilter( const double& val, double& df ) const {
double f = 1.0 - sf.calculate( val, df ); df*=-val;
return f;
}
diff --git a/src/multicolvar/LocalAverage.cpp b/src/multicolvar/LocalAverage.cpp
index ea0536c559f3fc9476ebda02ef215572f5a8721c..aa1359af53f0946f104348af4cbdc84c14c6bc41 100644
--- a/src/multicolvar/LocalAverage.cpp
+++ b/src/multicolvar/LocalAverage.cpp
@@ -92,7 +92,7 @@ public:
/// We have to overwrite this here
unsigned getNumberOfQuantities();
/// Actually do the calculation
- double compute( const unsigned& tindex, AtomValuePack& myatoms );
+ double compute( const unsigned& tindex, AtomValuePack& myatoms ) const ;
/// Is the variable periodic
bool isPeriodic(){ return false; }
};
@@ -140,10 +140,10 @@ unsigned LocalAverage::getNumberOfQuantities(){
return getBaseMultiColvar(0)->getNumberOfQuantities();
}
-double LocalAverage::compute( const unsigned& tindex, AtomValuePack& myatoms ){
+double LocalAverage::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
Vector distance; double d2, sw, dfunc, nbond=1; CatomPack atom0, atom1;
std::vector<double> values( getBaseMultiColvar(0)->getNumberOfQuantities() );
- vesselbase::MultiValue myder(values.size(), getNumberOfDerivatives());
+ MultiValue myder(values.size(), myatoms.getNumberOfDerivatives());
getVectorForTask( myatoms.getIndex(0), false, values );
if( values.size()>2 ){
@@ -219,7 +219,7 @@ double LocalAverage::compute( const unsigned& tindex, AtomValuePack& myatoms ){
myatoms.setValue( 0, nbond ); updateActiveAtoms( myatoms );
if( values.size()>2){
double norm=0;
- vesselbase::MultiValue& myvals=myatoms.getUnderlyingMultiValue();
+ MultiValue& myvals=myatoms.getUnderlyingMultiValue();
for(unsigned i=2;i<values.size();++i){
myvals.quotientRule( i, 0, i );
// Calculate length of vector
diff --git a/src/multicolvar/MultiColvar.cpp b/src/multicolvar/MultiColvar.cpp
index d44c61a1e593f83f6bd34bd1a9dcd3af5f9536f0..618778204ce175dfbfff8300c311e49919a77446 100644
--- a/src/multicolvar/MultiColvar.cpp
+++ b/src/multicolvar/MultiColvar.cpp
@@ -363,7 +363,7 @@ void MultiColvar::calculate(){
runAllTasks();
}
-void MultiColvar::updateActiveAtoms( AtomValuePack& myatoms ){
+void MultiColvar::updateActiveAtoms( AtomValuePack& myatoms ) const {
myatoms.updateUsingIndices();
}
diff --git a/src/multicolvar/MultiColvar.h b/src/multicolvar/MultiColvar.h
index 3b1c728e233588157f5f6c54667e2cb3b3ebda81..5bc5f2f63e93322cb40756313fbf1fe3d4c7597a 100644
--- a/src/multicolvar/MultiColvar.h
+++ b/src/multicolvar/MultiColvar.h
@@ -66,15 +66,15 @@ public:
/// Calculate the multicolvar
virtual void calculate();
/// Update the atoms that have derivatives
- void updateActiveAtoms( AtomValuePack& myatoms );
+ void updateActiveAtoms( AtomValuePack& myatoms ) const ;
/// This is used in MultiColvarBase only - it is used to setup the link cells
- Vector getPositionOfAtomForLinkCells( const unsigned& iatom );
+ Vector getPositionOfAtomForLinkCells( const unsigned& iatom ) const ;
/// Atoms are always active
bool isCurrentlyActive( const unsigned& code ){ return true; }
};
inline
-Vector MultiColvar::getPositionOfAtomForLinkCells( const unsigned& iatom ){
+Vector MultiColvar::getPositionOfAtomForLinkCells( const unsigned& iatom ) const {
return ActionAtomistic::getPosition( iatom );
}
diff --git a/src/multicolvar/MultiColvarBase.cpp b/src/multicolvar/MultiColvarBase.cpp
index 62e9ef8a559cd27d0c27eaee865dbf3a5b1e3ca1..38f4a27fc8c7d4911a3af40f60d85fff254c90b4 100644
--- a/src/multicolvar/MultiColvarBase.cpp
+++ b/src/multicolvar/MultiColvarBase.cpp
@@ -193,7 +193,7 @@ void MultiColvarBase::setupLinkCells(){
}
}
-void MultiColvarBase::decodeIndexToAtoms( const unsigned& taskCode, std::vector<unsigned>& atoms ){
+void MultiColvarBase::decodeIndexToAtoms( const unsigned& taskCode, std::vector<unsigned>& atoms ) const {
plumed_dbg_assert( atoms.size()==ablocks.size() && !usespecies && ablocks.size()<4 );
unsigned scode = taskCode;
for(unsigned i=0;i<ablocks.size();++i){
@@ -203,7 +203,7 @@ void MultiColvarBase::decodeIndexToAtoms( const unsigned& taskCode, std::vector<
}
}
-bool MultiColvarBase::setupCurrentAtomList( const unsigned& taskCode, AtomValuePack& myatoms ){
+bool MultiColvarBase::setupCurrentAtomList( const unsigned& taskCode, AtomValuePack& myatoms ) const {
if( usespecies ){
unsigned natomsper=1; std::vector<unsigned> current_atoms( getNumberOfAtoms() );
if( isDensity() ) return true;
@@ -223,7 +223,7 @@ bool MultiColvarBase::setupCurrentAtomList( const unsigned& taskCode, AtomValueP
return true;
}
-void MultiColvarBase::performTask( const unsigned& task_index, const unsigned& current, vesselbase::MultiValue& myvals ){
+void MultiColvarBase::performTask( const unsigned& task_index, const unsigned& current, MultiValue& myvals ) const {
AtomValuePack myatoms( myvals, this );
// Retrieve the atom list
@@ -242,11 +242,11 @@ void MultiColvarBase::performTask( const unsigned& task_index, const unsigned& c
return;
}
-void MultiColvarBase::calculateWeight( AtomValuePack& myatoms ){
+void MultiColvarBase::calculateWeight( AtomValuePack& myatoms ) const {
myatoms.setValue( 0, 1.0 );
}
-double MultiColvarBase::doCalculation( const unsigned& taskIndex, AtomValuePack& myatoms ){
+double MultiColvarBase::doCalculation( const unsigned& taskIndex, AtomValuePack& myatoms ) const {
double val=compute( taskIndex, myatoms ); updateActiveAtoms( myatoms );
return val;
}
diff --git a/src/multicolvar/MultiColvarBase.h b/src/multicolvar/MultiColvarBase.h
index 06b2cdf93620e878de06abc8df6cbb4d0188f850..ecd60f3f2875ac435b218a2fa74b9d48195858cd 100644
--- a/src/multicolvar/MultiColvarBase.h
+++ b/src/multicolvar/MultiColvarBase.h
@@ -27,8 +27,6 @@
#include "tools/DynamicList.h"
#include "tools/LinkCells.h"
#include "vesselbase/ActionWithVessel.h"
-// #include "StoreColvarVessel.h"
-// #include "StoreCentralAtomsVessel.h"
#include <vector>
namespace PLMD {
@@ -92,9 +90,9 @@ protected:
/// Get the number of atoms in this particular colvar
unsigned getNAtoms() const;
/// This sets up the list of atoms that are involved in this colvar
- bool setupCurrentAtomList( const unsigned& taskCode, AtomValuePack& myatoms );
+ bool setupCurrentAtomList( const unsigned& taskCode, AtomValuePack& myatoms ) const ;
/// Decode indices if there are 2 or 3 atoms involved
- void decodeIndexToAtoms( const unsigned& taskCode, std::vector<unsigned>& atoms );
+ void decodeIndexToAtoms( const unsigned& taskCode, std::vector<unsigned>& atoms ) const ;
public:
MultiColvarBase(const ActionOptions&);
~MultiColvarBase(){}
@@ -103,15 +101,15 @@ public:
virtual void turnOnDerivatives();
/// Prepare for the calculation
/// Perform one of the tasks
- virtual void performTask( const unsigned& , const unsigned& , vesselbase::MultiValue& );
+ virtual void performTask( const unsigned& , const unsigned& , MultiValue& ) const ;
/// This gets the position of an atom for the link cell setup
- virtual Vector getPositionOfAtomForLinkCells( const unsigned& iatom )=0;
+ virtual Vector getPositionOfAtomForLinkCells( const unsigned& iatom ) const=0;
/// And a virtual function which actually computes the colvar
- virtual double doCalculation( const unsigned& tindex, AtomValuePack& myatoms );
+ virtual double doCalculation( const unsigned& tindex, AtomValuePack& myatoms ) const ;
/// Update the atoms that have derivatives
- virtual void updateActiveAtoms( AtomValuePack& myatoms ){};
+ virtual void updateActiveAtoms( AtomValuePack& myatoms ) const {};
/// This is replaced once we have a function to calculate the cv
- virtual double compute( const unsigned& tindex, AtomValuePack& myatoms )=0;
+ virtual double compute( const unsigned& tindex, AtomValuePack& myatoms ) const=0;
/// Apply the forces from this action
virtual void apply();
/// Get the number of derivatives for this action
@@ -123,11 +121,11 @@ public:
/// Get the index where the central atom is stored
virtual Vector getCentralAtomPos( const unsigned& curr );
/// You can use this to screen contributions that are very small so we can avoid expensive (and pointless) calculations
- virtual void calculateWeight( AtomValuePack& myatoms );
+ virtual void calculateWeight( AtomValuePack& myatoms ) const ;
/// Get the list of indices that have derivatives
// virtual void getIndexList( const unsigned& ntotal, const unsigned& jstore, const unsigned& maxder, std::vector<unsigned>& indices );
/// Is this a density?
- virtual bool isDensity(){ return false; }
+ virtual bool isDensity() const { return false; }
/// Store central atoms so that this can be used in a function
// virtual vesselbase::StoreDataVessel* buildDataStashes( const bool& allow_wcutoff, const double& wtol );
/// Calculate and store getElementValue(uder)/getElementValue(vder) and its derivatives in getElementValue(iout)
diff --git a/src/multicolvar/MultiColvarDensity.cpp b/src/multicolvar/MultiColvarDensity.cpp
index 6edac4cacd89d8417a1c12175c8c2f61236cbafa..cadb11b843d5b704e2ff9ce98c7b29e884efde6a 100644
--- a/src/multicolvar/MultiColvarDensity.cpp
+++ b/src/multicolvar/MultiColvarDensity.cpp
@@ -34,6 +34,7 @@
#include "tools/Grid.h"
#include "tools/KernelFunctions.h"
#include "vesselbase/ActionWithInputVessel.h"
+#include "vesselbase/StoreDataVessel.h"
using namespace std;
@@ -181,16 +182,16 @@ MultiColvarDensity::~MultiColvarDensity(){
void MultiColvarDensity::update(){
+ vesselbase::StoreDataVessel* stash=dynamic_cast<vesselbase::StoreDataVessel*>( getPntrToArgument() );
+ plumed_dbg_assert( stash ); std::vector<double> cvals( mycolv->getNumberOfQuantities() );
Vector origin = getPosition(0); std::vector<double> pp( directions.size() );
- std::vector<double> cvals( mycolv->getNumberOfQuantities()-4 );
for(unsigned i=0;i<mycolv->getFullNumberOfTasks();++i){
- mycolv->getValueForTask( i, cvals );
- Vector apos = pbcDistance( mycolv->getCentralAtomPosition(i), origin );
+ stash->retrieveValue( i, false, cvals );
+ Vector apos = pbcDistance( mycolv->getCentralAtomPos( mycolv->getTaskCode(i) ), origin );
Vector fpos = getPbc().realToScaled( apos );
for(unsigned j=0;j<directions.size();++j) pp[j]=fpos[ directions[j] ];
- KernelFunctions kernel( pp, bw, kerneltype, false, cvals[0], true );
- gg->addKernel( kernel );
- norm += 1.0; // This should be replaced by the proper weight
+ KernelFunctions kernel( pp, bw, kerneltype, false, cvals[0]*cvals[1], true );
+ gg->addKernel( kernel ); norm += cvals[0];
}
// Output and reset the counter if required
diff --git a/src/multicolvar/MultiColvarFilter.cpp b/src/multicolvar/MultiColvarFilter.cpp
index 8e583713309da03f66e53424a8adabdab90a4976..d30976b6ab63df7eedfc8fe1425468a5e44edd7c 100644
--- a/src/multicolvar/MultiColvarFilter.cpp
+++ b/src/multicolvar/MultiColvarFilter.cpp
@@ -43,7 +43,7 @@ void MultiColvarFilter::doJobsRequiredBeforeTaskList(){
ActionWithVessel::doJobsRequiredBeforeTaskList();
}
-void MultiColvarFilter::completeTask( const unsigned& curr, vesselbase::MultiValue& invals, vesselbase::MultiValue& outvals ){
+void MultiColvarFilter::completeTask( const unsigned& curr, MultiValue& invals, MultiValue& outvals ) const {
invals.copyValues( outvals );
if( derivativesAreRequired() ) invals.copyDerivatives( outvals );
diff --git a/src/multicolvar/MultiColvarFilter.h b/src/multicolvar/MultiColvarFilter.h
index 5fa523ede97ba2f521e714288c71d422b1ba874c..51a889d7f3a50b9511c676e2065017266a7a85ca 100644
--- a/src/multicolvar/MultiColvarFilter.h
+++ b/src/multicolvar/MultiColvarFilter.h
@@ -43,9 +43,9 @@ public:
/// Get the number of quantities in the colvar
unsigned getNumberOfQuantities();
/// Actually do what we are asked
- void completeTask( const unsigned& curr, vesselbase::MultiValue& invals, vesselbase::MultiValue& outvals );
+ void completeTask( const unsigned& curr, MultiValue& invals, MultiValue& outvals ) const ;
/// Do the filtering
- virtual double applyFilter( const double& val, double& df )=0;
+ virtual double applyFilter( const double& val, double& df ) const=0;
/// Just checks there are no bridging forces
void addBridgeForces( const std::vector<double>& bb );
};
diff --git a/src/multicolvar/MultiColvarFunction.cpp b/src/multicolvar/MultiColvarFunction.cpp
index a562c64db80fef731b8f6b5bb0fab8b524a5f0f4..130a5467f94685a56d134c34cc5e87a3f7e82929 100644
--- a/src/multicolvar/MultiColvarFunction.cpp
+++ b/src/multicolvar/MultiColvarFunction.cpp
@@ -241,11 +241,11 @@ void MultiColvarFunction::calculateNumericalDerivatives( ActionWithValue* a ){
}
}
-void MultiColvarFunction::updateActiveAtoms( AtomValuePack& myatoms ){
+void MultiColvarFunction::updateActiveAtoms( AtomValuePack& myatoms ) const {
myatoms.updateDynamicList();
}
-void MultiColvarFunction::getVectorDerivatives( const unsigned& ind, const bool& normed, vesselbase::MultiValue& myder ){
+void MultiColvarFunction::getVectorDerivatives( const unsigned& ind, const bool& normed, MultiValue& myder ) const {
plumed_dbg_assert( ind<getFullNumberOfBaseTasks() ); unsigned mmc=colvar_label[ind];
plumed_dbg_assert( mybasedata[mmc]->storedValueIsActive( convertToLocalIndex(ind,mmc) ) );
myder.resize( mybasemulticolvars[mmc]->getNumberOfQuantities(), mybasemulticolvars[mmc]->getNumberOfDerivatives() );
@@ -254,7 +254,7 @@ void MultiColvarFunction::getVectorDerivatives( const unsigned& ind, const bool&
void MultiColvarFunction::mergeVectorDerivatives( const unsigned& ival, const unsigned& start, const unsigned& end,
const unsigned& jatom, const std::vector<double>& der,
- vesselbase::MultiValue& myder, AtomValuePack& myatoms ){
+ MultiValue& myder, AtomValuePack& myatoms ) const {
plumed_dbg_assert( ival<myatoms.getUnderlyingMultiValue().getNumberOfValues() );
plumed_dbg_assert( start<myder.getNumberOfValues() && end<=myder.getNumberOfValues() );
plumed_dbg_assert( der.size()==myder.getNumberOfValues() && jatom<getFullNumberOfBaseTasks() );
@@ -266,7 +266,7 @@ void MultiColvarFunction::mergeVectorDerivatives( const unsigned& ival, const un
// Now get the start of the virial
unsigned virbas = 3*getNumberOfAtoms();
- vesselbase::MultiValue& myvals=myatoms.getUnderlyingMultiValue();
+ MultiValue& myvals=myatoms.getUnderlyingMultiValue();
for(unsigned j=0;j<myder.getNumberActive();++j){
unsigned jder=myder.getActiveIndex(j);
if( jder<3*mybasemulticolvars[mmc]->getNumberOfAtoms() ){
diff --git a/src/multicolvar/MultiColvarFunction.h b/src/multicolvar/MultiColvarFunction.h
index cd8e95ab1b4a5c31d945dff97d8f530b23fb75fb..b7204a014c76dc81b4bee6715161fa7212992e99 100644
--- a/src/multicolvar/MultiColvarFunction.h
+++ b/src/multicolvar/MultiColvarFunction.h
@@ -49,37 +49,17 @@ protected:
/// Get the total number of tasks that this calculation is based on
unsigned getFullNumberOfBaseTasks() const ;
/// Get the derivatives for the central atom with index ind
- CatomPack getCentralAtomPackFromInput( const unsigned& ind );
+ CatomPack getCentralAtomPackFromInput( const unsigned& ind ) const ;
///
- void getVectorForTask( const unsigned& ind, const bool& normed, std::vector<double>& orient0 );
+ void getVectorForTask( const unsigned& ind, const bool& normed, std::vector<double>& orient0 ) const ;
///
- void getVectorDerivatives( const unsigned& ind, const bool& normed, vesselbase::MultiValue& myder0 );
+ void getVectorDerivatives( const unsigned& ind, const bool& normed, MultiValue& myder0 ) const ;
///
void mergeVectorDerivatives( const unsigned& ival, const unsigned& start, const unsigned& end,
const unsigned& jatom, const std::vector<double>& der,
- vesselbase::MultiValue& myder, AtomValuePack& myatoms );
-/// Get the index of the ith colvar we are using
-// unsigned getColvarIndex( const unsigned& ) const;
-/// Get the position of one of the central atoms
-// Vector getPositionOfCentralAtom(const unsigned&) const;
-/// Add derivatives of value wrt to an atomic position
-// void addCentralAtomsDerivatives( const unsigned& , const unsigned& , const Vector& );
-/// Extract the value for the iatom th base task (this ignores current_atoms)
-// void extractValueForBaseTask( const unsigned& iatom, std::vector<double>& vals );
-/// Retrieve the value calculated by the iatom th base task (this uses current_atoms)
-// void getValueForBaseTask( const unsigned& iatom, std::vector<double>& vals );
-/// Retrieve the vector calculated by the iatom th base task (this uses current_atoms)
-// void getVectorForBaseTask( const unsigned& iatom, std::vector<double>& vecs );
-/// Add the derivatives of this quantity (this ignores current_atoms)
-// void extractWeightedAverageAndDerivatives( const unsigned& iatom, const double& weight );
-/// Add the value and derivatives of this quantity (this uses current_atoms)
-// void accumulateWeightedAverageAndDerivatives( const unsigned& iatom, const double& weight );
-/// Add derivative wrt to the position of the central atom
-// void addDerivativeOfCentralAtomPos( const unsigned& iatom, const Tensor& der );
+ MultiValue& myder, AtomValuePack& myatoms ) const ;
/// Convert an index in the global array to an index in the individual base colvars
unsigned convertToLocalIndex( const unsigned& index, const unsigned& mcv_code ) const ;
-/// Add derivatives to the orientation
-// void addOrientationDerivatives( const unsigned& iatom, const std::vector<double>& der );
/// Build sets by taking one multicolvar from each base
void buildSets();
/// Build colvars for atoms as if they were symmetry functions
@@ -89,30 +69,20 @@ protected:
/// Get the number of base multicolvars
unsigned getNumberOfBaseMultiColvars() const ;
/// Get an example of the underlying multicolvar
- MultiColvarBase* getBaseMultiColvar( const unsigned& icolv );
+ MultiColvarBase* getBaseMultiColvar( const unsigned& icolv ) const ;
/// Return the base multicolvar index that this colvar is a part of
unsigned getBaseColvarNumber( const unsigned& iatom ) const ;
public:
MultiColvarFunction(const ActionOptions&);
static void registerKeywords( Keywords& keys );
-/// Active element in atoms_with_derivatives
-// void atomHasDerivative( const unsigned& iatom );
-/// Resize the dynamic arrays
-// void resizeDynamicArrays();
/// Update the atoms that are active
- virtual void updateActiveAtoms( AtomValuePack& myatoms );
+ virtual void updateActiveAtoms( AtomValuePack& myatoms ) const ;
/// Regular calculate
void calculate();
/// Calculate the numerical derivatives for this action
void calculateNumericalDerivatives( ActionWithValue* a=NULL );
-/// Calculate the position of the central atom
-// Vector calculateCentralAtomPosition();
-/// Get the position of the central atom
-// virtual Vector getCentralAtom()=0;
-/// Add derivatives from storage vessels in MultiColvarBase
-// void addStoredDerivative( const unsigned&, const unsigned&, const unsigned&, const double& );
/// This is used in MultiColvarBase only - it is used to setup the link cells
- Vector getPositionOfAtomForLinkCells( const unsigned& iatom );
+ Vector getPositionOfAtomForLinkCells( const unsigned& iatom ) const ;
/// Some things can be inactive in functions
bool isCurrentlyActive( const unsigned& code );
};
@@ -135,7 +105,7 @@ unsigned MultiColvarFunction::getNumberOfBaseMultiColvars() const {
}
inline
-MultiColvarBase* MultiColvarFunction::getBaseMultiColvar( const unsigned& icolv ){
+MultiColvarBase* MultiColvarFunction::getBaseMultiColvar( const unsigned& icolv ) const {
plumed_dbg_assert( icolv<mybasemulticolvars.size() );
return mybasemulticolvars[icolv];
}
@@ -153,13 +123,13 @@ unsigned MultiColvarFunction::getBaseColvarNumber( const unsigned& iatom ) const
}
inline
-Vector MultiColvarFunction::getPositionOfAtomForLinkCells( const unsigned& iatom ){
+Vector MultiColvarFunction::getPositionOfAtomForLinkCells( const unsigned& iatom ) const {
plumed_dbg_assert( iatom<getFullNumberOfBaseTasks() ); unsigned mmc=colvar_label[ iatom ];
return mybasemulticolvars[mmc]->getCentralAtomPos( convertToLocalIndex(iatom,mmc) );
}
inline
-CatomPack MultiColvarFunction::getCentralAtomPackFromInput( const unsigned& ind ){
+CatomPack MultiColvarFunction::getCentralAtomPackFromInput( const unsigned& ind ) const {
plumed_dbg_assert( ind<getFullNumberOfBaseTasks() ); unsigned mmc=colvar_label[ind];
unsigned basen=0;
for(unsigned i=0;i<mmc;++i) basen+=mybasemulticolvars[i]->getNumberOfAtoms();
@@ -167,81 +137,12 @@ CatomPack MultiColvarFunction::getCentralAtomPackFromInput( const unsigned& ind
}
inline
-void MultiColvarFunction::getVectorForTask( const unsigned& ind, const bool& normed, std::vector<double>& orient ){
+void MultiColvarFunction::getVectorForTask( const unsigned& ind, const bool& normed, std::vector<double>& orient ) const {
plumed_dbg_assert( ind<getFullNumberOfBaseTasks() ); unsigned mmc=colvar_label[ind];
plumed_dbg_assert( mybasedata[mmc]->storedValueIsActive( convertToLocalIndex(ind,mmc) ) );
mybasedata[mmc]->retrieveValue( convertToLocalIndex(ind,mmc), normed, orient );
}
-// inline
-// Vector MultiColvarFunction::getPositionOfCentralAtom( const unsigned& iatom ) const {
-// plumed_dbg_assert( iatom<natomsper ); unsigned mmc = colvar_label[ current_atoms[iatom] ];
-// return mybasemulticolvars[mmc]->getCentralAtomPosition( convertToLocalIndex(current_atoms[iatom],mmc) );
-// }
-//
-// inline
-// void MultiColvarFunction::addCentralAtomsDerivatives( const unsigned& iatom, const unsigned& jout, const Vector& der ){
-// if( doNotCalculateDerivatives() ) return ;
-//
-// plumed_dbg_assert( iatom<natomsper ); unsigned mmc = colvar_label[ current_atoms[iatom] ];
-// mybasemulticolvars[mmc]->addCentralAtomDerivativeToFunction( convertToLocalIndex(current_atoms[iatom],mmc), jout, mmc, der, this );
-// }
-
-// inline
-// void MultiColvarFunction::atomHasDerivative( const unsigned& iatom ){
-// plumed_dbg_assert( !doNotCalculateDerivatives() );
-// atoms_with_derivatives.activate( iatom );
-// }
-
-// inline
-// void MultiColvarFunction::addDerivativeOfCentralAtomPos( const unsigned& iatom, const Tensor& der ){
-// if( doNotCalculateDerivatives() ) return;
-//
-// plumed_dbg_assert( iatom<natomsper ); unsigned mmc = colvar_label[ current_atoms[iatom] ]; Vector tmpder;
-// for(unsigned i=0;i<3;++i){
-// for(unsigned j=0;j<3;++j) tmpder[j]=der(i,j);
-// mybasemulticolvars[mmc]->addCentralAtomDerivativeToFunction( convertToLocalIndex(current_atoms[iatom],mmc), (2+i), mmc, tmpder, this );
-// }
-// }
-
-// inline
-// void MultiColvarFunction::extractValueForBaseTask( const unsigned& iatom, std::vector<double>& vals ){
-// unsigned mmc = colvar_label[iatom]; mybasemulticolvars[mmc]->getValueForTask( convertToLocalIndex(iatom,mmc), vals );
-// }
-//
-// inline
-// void MultiColvarFunction::getValueForBaseTask( const unsigned& iatom, std::vector<double>& vals ){
-// plumed_dbg_assert( iatom<natomsper ); extractValueForBaseTask( current_atoms[iatom], vals );
-// }
-
-// inline
-// void MultiColvarFunction::getVectorForBaseTask( const unsigned& iatom, std::vector<double>& vec ){
-// plumed_dbg_assert( vec.size()==mybasemulticolvars[0]->getNumberOfQuantities()-5 && tvals.size()>1 );
-// getValueForBaseTask( iatom, tvals ); for(unsigned i=0;i<vec.size();++i) vec[i]=tvals[i+1];
-// }
-
-// inline
-// void MultiColvarFunction::extractWeightedAverageAndDerivatives( const unsigned& iatom, const double& weight ){
-// if( doNotCalculateDerivatives() ) return;
-// unsigned mmc = colvar_label[iatom];
-// mybasemulticolvars[mmc]->addWeightedValueDerivatives( convertToLocalIndex(iatom, mmc), mmc, weight, this );
-// }
-//
-// inline
-// void MultiColvarFunction::accumulateWeightedAverageAndDerivatives( const unsigned& iatom, const double& weight ){
-// if( doNotCalculateDerivatives() ) return;
-// plumed_dbg_assert( iatom<natomsper ); extractWeightedAverageAndDerivatives( current_atoms[iatom], weight );
-// }
-//
-// inline
-// void MultiColvarFunction::addOrientationDerivatives( const unsigned& iatom , const std::vector<double>& der ){
-// if( doNotCalculateDerivatives() ) return;
-//
-// plumed_dbg_assert( iatom<natomsper ); unsigned mmc = colvar_label[ current_atoms[iatom] ];
-// unsigned jout=2; if( usespecies && iatom==0 ) jout=1;
-// mybasemulticolvars[mmc]->addOrientationDerivativesToBase( convertToLocalIndex(current_atoms[iatom],mmc), jout, mmc, der, this );
-// }
-
}
}
#endif
diff --git a/src/multicolvar/NumberOfLinks.cpp b/src/multicolvar/NumberOfLinks.cpp
index dbc11d8b89556408456f3f7f070da0f58fd92bbb..99d1589854ddff6d487f109c0bf34dda2074325b 100644
--- a/src/multicolvar/NumberOfLinks.cpp
+++ b/src/multicolvar/NumberOfLinks.cpp
@@ -75,9 +75,9 @@ public:
static void registerKeywords( Keywords& keys );
NumberOfLinks(const ActionOptions&);
/// Do the stuff with the switching functions
- void calculateWeight( AtomValuePack& myatoms );
+ void calculateWeight( AtomValuePack& myatoms ) const ;
/// Actually do the calculation
- double compute( const unsigned& tindex, AtomValuePack& myatoms );
+ double compute( const unsigned& tindex, AtomValuePack& myatoms ) const ;
/// Is the variable periodic
bool isPeriodic(){ return false; }
};
@@ -130,7 +130,7 @@ MultiColvarFunction(ao)
readVesselKeywords();
}
-void NumberOfLinks::calculateWeight( AtomValuePack& myatoms ){
+void NumberOfLinks::calculateWeight( AtomValuePack& myatoms ) const {
Vector distance = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
double dfunc, sw = switchingFunction.calculateSqr( distance.modulo2(), dfunc );
myatoms.setValue(0,sw);
@@ -144,7 +144,7 @@ void NumberOfLinks::calculateWeight( AtomValuePack& myatoms ){
}
}
-double NumberOfLinks::compute( const unsigned& tindex, AtomValuePack& myatoms ){
+double NumberOfLinks::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
if( getBaseMultiColvar(0)->getNumberOfQuantities()<3 ) return 1.0;
unsigned ncomp=getBaseMultiColvar(0)->getNumberOfQuantities();
@@ -159,8 +159,8 @@ double NumberOfLinks::compute( const unsigned& tindex, AtomValuePack& myatoms ){
}
if( !doNotCalculateDerivatives() ){
- unsigned nder=getNumberOfDerivatives();
- vesselbase::MultiValue myder0(ncomp,nder), myder1(ncomp,nder);
+ unsigned nder=myatoms.getNumberOfDerivatives(); //getNumberOfDerivatives();
+ MultiValue myder0(ncomp,nder), myder1(ncomp,nder);
getVectorDerivatives( myatoms.getIndex(0), true, myder0 );
mergeVectorDerivatives( 1, 2, orient1.size(), myatoms.getIndex(0), orient1, myder0, myatoms );
getVectorDerivatives( myatoms.getIndex(1), true, myder1 );
diff --git a/src/multicolvar/Sprint.cpp b/src/multicolvar/Sprint.cpp
index c007ef84a4ecafcbef025aad55612eab2dd9a04d..918042f6c935f62f71431a7643f9e75a23b00b84 100644
--- a/src/multicolvar/Sprint.cpp
+++ b/src/multicolvar/Sprint.cpp
@@ -169,7 +169,7 @@ void Sprint::completeCalculation(){
else { rank=comm.Get_rank(); stride=comm.Get_size(); }
// Derivatives
- vesselbase::MultiValue myvals( 2, getNumberOfDerivatives() );
+ MultiValue myvals( 2, getNumberOfDerivatives() );
Matrix<double> mymat_ders( getNumberOfComponents(), getNumberOfDerivatives() );
std::vector<unsigned> catoms(2); unsigned nval = getFullNumberOfBaseTasks(); mymat_ders=0;
for(unsigned i=rank;i<getNumberOfActiveMatrixElements();i+=stride){
diff --git a/src/multicolvar/Torsions.cpp b/src/multicolvar/Torsions.cpp
index 154547a4964784b28e9b37c1070e7077077f98c7..709e7991208fc6d85fb6bd1b90be88dd1e1ee83b 100644
--- a/src/multicolvar/Torsions.cpp
+++ b/src/multicolvar/Torsions.cpp
@@ -75,7 +75,7 @@ class Torsions : public MultiColvar {
public:
static void registerKeywords( Keywords& keys );
Torsions(const ActionOptions&);
- virtual double compute( const unsigned& tindex, AtomValuePack& myatoms );
+ virtual double compute( const unsigned& tindex, AtomValuePack& myatoms ) const ;
bool isPeriodic(){ return true; }
void retrieveDomain( std::string& min, std::string& max ){ min="-pi"; max="pi"; }
};
@@ -101,7 +101,7 @@ PLUMED_MULTICOLVAR_INIT(ao)
checkRead();
}
-double Torsions::compute( const unsigned& tindex, AtomValuePack& myatoms ){
+double Torsions::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
Vector d0,d1,d2;
d0=getSeparation(myatoms.getPosition(1),myatoms.getPosition(0));
d1=getSeparation(myatoms.getPosition(2),myatoms.getPosition(1));
diff --git a/src/multicolvar/VolumeAround.cpp b/src/multicolvar/VolumeAround.cpp
index 1a70cabb8bea7f2f27f6b94b564236a4448adc0e..0883d9186b57aff87e672c839eefec81b335738d 100644
--- a/src/multicolvar/VolumeAround.cpp
+++ b/src/multicolvar/VolumeAround.cpp
@@ -78,7 +78,7 @@ public:
static void registerKeywords( Keywords& keys );
VolumeAround(const ActionOptions& ao);
void setupRegions();
- double calculateNumberInside( const Vector& cpos, HistogramBead& bead, Vector& derivatives, Tensor& vir, std::vector<Vector>& refders );
+ double calculateNumberInside( const Vector& cpos, Vector& derivatives, Tensor& vir, std::vector<Vector>& refders ) const ;
};
PLUMED_REGISTER_ACTION(VolumeAround,"AROUND")
@@ -116,10 +116,12 @@ ActionVolume(ao)
void VolumeAround::setupRegions(){ }
-double VolumeAround::calculateNumberInside( const Vector& cpos, HistogramBead& bead, Vector& derivatives, Tensor& vir, std::vector<Vector>& refders ){
+double VolumeAround::calculateNumberInside( const Vector& cpos, Vector& derivatives, Tensor& vir, std::vector<Vector>& refders ) const {
+ // Setup the histogram bead
+ HistogramBead bead; bead.isNotPeriodic(); bead.setKernelType( getKernelType() );
+
// Calculate position of atom wrt to origin
Vector fpos=pbcDistance( getPosition(0), cpos );
-
double xcontr, ycontr, zcontr, xder, yder, zder;
if( dox ){
bead.set( xlow, xhigh, getSigma() );
diff --git a/src/multicolvar/VolumeCavity.cpp b/src/multicolvar/VolumeCavity.cpp
index d0f5b59a95560c9ccbcd4cf0013800394f5f7acd..78cb96be2fb92689a7df9d2239f3c40011a69daf 100644
--- a/src/multicolvar/VolumeCavity.cpp
+++ b/src/multicolvar/VolumeCavity.cpp
@@ -122,7 +122,7 @@ public:
~VolumeCavity();
void setupRegions();
void update();
- double calculateNumberInside( const Vector& cpos, HistogramBead& bead, Vector& derivatives, Tensor& vir, std::vector<Vector>& refders );
+ double calculateNumberInside( const Vector& cpos, Vector& derivatives, Tensor& vir, std::vector<Vector>& refders ) const ;
};
PLUMED_REGISTER_ACTION(VolumeCavity,"CAVITY")
@@ -340,7 +340,10 @@ void VolumeCavity::update(){
}
}
-double VolumeCavity::calculateNumberInside( const Vector& cpos, HistogramBead& bead, Vector& derivatives, Tensor& vir, std::vector<Vector>& rderiv ){
+double VolumeCavity::calculateNumberInside( const Vector& cpos, Vector& derivatives, Tensor& vir, std::vector<Vector>& rderiv ) const {
+ // Setup the histogram bead
+ HistogramBead bead; bead.isNotPeriodic(); bead.setKernelType( getKernelType() );
+
// Calculate distance of atom from origin of new coordinate frame
Vector datom=pbcDistance( origin, cpos );
double ucontr, uder, vcontr, vder, wcontr, wder;
diff --git a/src/multicolvar/VolumeGradientBase.cpp b/src/multicolvar/VolumeGradientBase.cpp
index 3f08a5e63b9ff76f5276ae480829d26d7e0b65c2..a817c24f72930d91b83673e03a5190a4955a180c 100644
--- a/src/multicolvar/VolumeGradientBase.cpp
+++ b/src/multicolvar/VolumeGradientBase.cpp
@@ -47,7 +47,7 @@ void VolumeGradientBase::doJobsRequiredBeforeTaskList(){
ActionWithVessel::doJobsRequiredBeforeTaskList();
}
-void VolumeGradientBase::completeTask( const unsigned& curr, vesselbase::MultiValue& invals, vesselbase::MultiValue& outvals ){
+void VolumeGradientBase::completeTask( const unsigned& curr, MultiValue& invals, MultiValue& outvals ) const {
if( getPntrToMultiColvar()->isDensity() ){
outvals.setValue(0, 1.0); outvals.setValue(1, 1.0);
} else {
@@ -60,7 +60,7 @@ void VolumeGradientBase::completeTask( const unsigned& curr, vesselbase::MultiVa
void VolumeGradientBase::setNumberInVolume( const unsigned& ivol, const unsigned& curr, const double& weight,
const Vector& wdf, const Tensor& virial, const std::vector<Vector>& refders,
- vesselbase::MultiValue& outvals ){
+ MultiValue& outvals ) const {
MultiColvarBase* mcolv=getPntrToMultiColvar();
if( !mcolv->weightHasDerivatives ){
outvals.setValue(ivol, weight );
diff --git a/src/multicolvar/VolumeGradientBase.h b/src/multicolvar/VolumeGradientBase.h
index 12a9c09344acd345b1118a94c3b79f5a0950fd6a..3e773585444ef5f9b31e6f9c934b902d5f3bc4f1 100644
--- a/src/multicolvar/VolumeGradientBase.h
+++ b/src/multicolvar/VolumeGradientBase.h
@@ -30,8 +30,6 @@ namespace multicolvar {
class VolumeGradientBase : public BridgedMultiColvarFunction {
friend class MultiColvarBase;
private:
-/// This is used for storing positions properly
- Vector tmp_p;
/// This is used to store forces temporarily in apply
std::vector<double> tmpforces;
protected:
@@ -42,20 +40,20 @@ protected:
/// Calculate distance between two points
Vector pbcDistance( const Vector& v1, const Vector& v2) const;
/// Get position of atom
- const Vector & getPosition( int iatom );
+ Vector getPosition( int iatom ) const ;
/// Request the atoms
void requestAtoms( const std::vector<AtomNumber>& atoms );
/// Set the number in the volume
- void setNumberInVolume( const unsigned& , const unsigned& , const double& , const Vector& , const Tensor& , const std::vector<Vector>& , vesselbase::MultiValue& );
+ void setNumberInVolume( const unsigned& , const unsigned& , const double& , const Vector& , const Tensor& , const std::vector<Vector>& , MultiValue& ) const ;
public:
static void registerKeywords( Keywords& keys );
VolumeGradientBase(const ActionOptions&);
/// Do jobs required before tasks are undertaken
void doJobsRequiredBeforeTaskList();
/// Actually do what we are asked
- void completeTask( const unsigned& curr, vesselbase::MultiValue& invals, vesselbase::MultiValue& outvals );
+ void completeTask( const unsigned& curr, MultiValue& invals, MultiValue& outvals ) const ;
/// Calculate what is in the volumes
- virtual void calculateAllVolumes( const unsigned& curr, vesselbase::MultiValue& outvals )=0;
+ virtual void calculateAllVolumes( const unsigned& curr, MultiValue& outvals ) const=0;
/// Setup the regions that this is based on
virtual void setupRegions()=0;
/// Forces here are applied through the bridge
@@ -78,10 +76,10 @@ Vector VolumeGradientBase::pbcDistance( const Vector& v1, const Vector& v2) cons
}
inline
-const Vector & VolumeGradientBase::getPosition( int iatom ){
+Vector VolumeGradientBase::getPosition( int iatom ) const {
if( !checkNumericalDerivatives() ) return ActionAtomistic::getPosition(iatom);
// This is for numerical derivatives of quantity wrt to the local atoms
- tmp_p = ActionAtomistic::getPosition(iatom);
+ Vector tmp_p = ActionAtomistic::getPosition(iatom);
if( bridgeVariable<3*getNumberOfAtoms() ){
if( bridgeVariable>=3*iatom && bridgeVariable<(iatom+1)*3 ) tmp_p[bridgeVariable%3]+=sqrt(epsilon);
}
diff --git a/src/multicolvar/XDistances.cpp b/src/multicolvar/XDistances.cpp
index c6ac26e960b9ad975cb29837b61b8797e947ebae..7612be128337327057f55c0ab00479b886194fb9 100644
--- a/src/multicolvar/XDistances.cpp
+++ b/src/multicolvar/XDistances.cpp
@@ -115,7 +115,7 @@ public:
static void registerKeywords( Keywords& keys );
XDistances(const ActionOptions&);
// active methods:
- virtual double compute( const unsigned& tindex, AtomValuePack& myatoms );
+ virtual double compute( const unsigned& tindex, AtomValuePack& myatoms ) const ;
/// Returns the number of coordinates of the field
bool isPeriodic(){ return false; }
};
@@ -150,7 +150,7 @@ PLUMED_MULTICOLVAR_INIT(ao)
checkRead();
}
-double XDistances::compute( const unsigned& tindex, AtomValuePack& myatoms ){
+double XDistances::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
Vector distance;
distance=getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
const double value=distance[myc];
diff --git a/src/reference/ArgumentOnlyDistance.cpp b/src/reference/ArgumentOnlyDistance.cpp
index f0bef55eeccfcb32aa0e6eedd37d2f44ce9facb9..947ae450d0f50e6fc6c19f68b57df2b1498555d3 100644
--- a/src/reference/ArgumentOnlyDistance.cpp
+++ b/src/reference/ArgumentOnlyDistance.cpp
@@ -30,16 +30,24 @@ ReferenceArguments(ro)
{
}
-double ArgumentOnlyDistance::calculate( const std::vector<Value*>& vals, const bool& squared ){
- clearDerivatives();
- if( tmparg.size()!=vals.size() ) tmparg.resize( vals.size() );
+void ArgumentOnlyDistance::read( const PDB& pdb ){
+ readArgumentsFromPDB( pdb );
+}
+
+double ArgumentOnlyDistance::calculate( const std::vector<Value*>& vals, ReferenceValuePack& myder, const bool& squared ) const {
+ std::vector<double> tmparg( vals.size() );
for(unsigned i=0;i<vals.size();++i) tmparg[i]=vals[i]->get();
- return calc( vals, tmparg, squared );
+ double d=calculateArgumentDistance( vals, tmparg, myder, squared );
+ if( !myder.updateComplete() ) myder.updateDynamicLists();
+ return d;
}
-double ArgumentOnlyDistance::calc( const std::vector<Vector>& pos, const Pbc& pbc, const std::vector<Value*>& vals, const std::vector<double>& arg, const bool& squared ){
+double ArgumentOnlyDistance::calc( const std::vector<Vector>& pos, const Pbc& pbc, const std::vector<Value*>& vals, const std::vector<double>& arg,
+ ReferenceValuePack& myder, const bool& squared ) const {
plumed_dbg_assert( pos.size()==0 );
- return calc( vals, arg, squared );
+ double d=calculateArgumentDistance( vals, arg, myder, squared );
+ if( !myder.updateComplete() ) myder.updateDynamicLists();
+ return d;
}
}
diff --git a/src/reference/ArgumentOnlyDistance.h b/src/reference/ArgumentOnlyDistance.h
index 9a761daa079f8179f5cc4c9262d79199a5985e46..1ff130c8173c01b7abcdd6eb16dd32de8d54bd9c 100644
--- a/src/reference/ArgumentOnlyDistance.h
+++ b/src/reference/ArgumentOnlyDistance.h
@@ -32,13 +32,11 @@ class PDB;
class Pbc;
class ArgumentOnlyDistance : public ReferenceArguments {
-private:
- std::vector<double> tmparg;
public:
ArgumentOnlyDistance( const ReferenceConfigurationOptions& ro );
- double calc( const std::vector<Vector>& pos, const Pbc& pbc, const std::vector<Value*>& vals, const std::vector<double>& arg, const bool& squared );
- double calculate( const std::vector<Value*>& vals, const bool& squared );
- virtual double calc( const std::vector<Value*>& vals, const std::vector<double>& arg, const bool& squared )=0;
+ void read( const PDB& pdb );
+ double calc( const std::vector<Vector>& pos, const Pbc& pbc, const std::vector<Value*>& vals, const std::vector<double>& arg, ReferenceValuePack& myder, const bool& squared ) const ;
+ double calculate( const std::vector<Value*>& vals, ReferenceValuePack& myder, const bool& squared ) const ;
};
}
diff --git a/src/reference/DRMSD.cpp b/src/reference/DRMSD.cpp
index 8ffa7d3d1a16dfcfbfcbad5c7b317b4af69bf76a..3aff7fe58129a8be39c06a66611a3c36fca43b13 100644
--- a/src/reference/DRMSD.cpp
+++ b/src/reference/DRMSD.cpp
@@ -72,7 +72,7 @@ void DRMSD::setup_targets(){
if( targets.size()==0 ) error("drmsd will compare no distances - check upper and lower bounds are sensible");
}
-double DRMSD::calc( const std::vector<Vector>& pos, const Pbc& pbc, const bool& squared ){
+double DRMSD::calc( const std::vector<Vector>& pos, const Pbc& pbc, ReferenceValuePack& myder, const bool& squared ) const {
plumed_dbg_assert( targets.size()>0 );
Vector distance;
@@ -89,9 +89,9 @@ double DRMSD::calc( const std::vector<Vector>& pos, const Pbc& pbc, const bool&
double diff = len - it->second;
drmsd += diff * diff;
- addAtomicDerivatives( i, -( diff / len ) * distance );
- addAtomicDerivatives( j, ( diff / len ) * distance );
- addBoxDerivatives( -( diff / len ) * Tensor(distance,distance) );
+ myder.addAtomDerivatives( i, -( diff / len ) * distance );
+ myder.addAtomDerivatives( j, ( diff / len ) * distance );
+ myder.addBoxDerivatives( -( diff / len ) * Tensor(distance,distance) );
}
double npairs = static_cast<double>(targets.size());
@@ -105,8 +105,9 @@ double DRMSD::calc( const std::vector<Vector>& pos, const Pbc& pbc, const bool&
idrmsd = 1.0/( drmsd * npairs );
}
- virial *= idrmsd;
- for(unsigned i=0;i<getNumberOfAtoms();++i){atom_ders[i] *= idrmsd;}
+ myder.scaleAllDerivatives( idrmsd );
+ // virial *= idrmsd;
+ // for(unsigned i=0;i<getNumberOfAtoms();++i){atom_ders[i] *= idrmsd;}
return drmsd;
}
diff --git a/src/reference/DRMSD.h b/src/reference/DRMSD.h
index a766c86d2982b22ac77264f8ad0818fed2e1e15a..217d484819c8c4223250a90402e620b29528d450 100644
--- a/src/reference/DRMSD.h
+++ b/src/reference/DRMSD.h
@@ -44,7 +44,7 @@ public:
// void check( ReferenceConfiguration* , ReferenceConfiguration* );
void read( const PDB& );
void setReferenceAtoms( const std::vector<Vector>& conf, const std::vector<double>& align_in, const std::vector<double>& displace_in );
- double calc( const std::vector<Vector>& pos, const Pbc& pbc, const bool& squared );
+ double calc( const std::vector<Vector>& pos, const Pbc& pbc, ReferenceValuePack& myder, const bool& squared ) const ;
};
}
diff --git a/src/reference/EuclideanDistance.cpp b/src/reference/EuclideanDistance.cpp
index ecf2acbe6b05cba615a0342e9a05faeadf87d2ae..d035782024217c5f80b0db21a090dafc935118ee 100644
--- a/src/reference/EuclideanDistance.cpp
+++ b/src/reference/EuclideanDistance.cpp
@@ -27,8 +27,6 @@ namespace PLMD {
class EuclideanDistance : public ArgumentOnlyDistance {
public:
EuclideanDistance( const ReferenceConfigurationOptions& ro );
- void read( const PDB& );
- double calc( const std::vector<Value*>& vals, const std::vector<double>& arg, const bool& squared );
};
PLUMED_REGISTER_METRIC(EuclideanDistance,"EUCLIDEAN")
@@ -39,13 +37,4 @@ ArgumentOnlyDistance(ro)
{
}
-void EuclideanDistance::read( const PDB& pdb ){
- readArgumentsFromPDB( pdb );
-}
-
-double EuclideanDistance::calc( const std::vector<Value*>& vals, const std::vector<double>& arg, const bool& squared ){
- return calculateArgumentDistance( vals, arg, squared );
-}
-
-
}
diff --git a/src/reference/FakeFrame.h b/src/reference/FakeFrame.h
index f5bdcd1752fe82cb32c51bb275b26986fc1e1c5d..79e9e5354aa72176fc1c370ce2049b691ced95c3 100644
--- a/src/reference/FakeFrame.h
+++ b/src/reference/FakeFrame.h
@@ -32,7 +32,7 @@ public PLMD::ReferenceConfiguration
public:
FakeFrame( const ReferenceConfigurationOptions& ro ) : ReferenceConfiguration(ro) {}
void read( const PDB& ){ plumed_merror("should not be called"); }
- double calc( const std::vector<Vector>& pos, const Pbc& pbc, const std::vector<Value*>& vals, const std::vector<double>& arg, const bool& squared ){
+ double calc( const std::vector<Vector>& pos, const Pbc& pbc, const std::vector<Value*>& vals, const std::vector<double>& arg, ReferenceValuePack& myder, const bool& squared ) const {
plumed_merror("should not be called"); return 1.0;
}
};
diff --git a/src/reference/MahalanobisDistance.cpp b/src/reference/MahalanobisDistance.cpp
index 21885ddf5af7ef4d881ebf4c995a92362ae495f4..2c6fd1f453711428b67c04963ce597261228c090 100644
--- a/src/reference/MahalanobisDistance.cpp
+++ b/src/reference/MahalanobisDistance.cpp
@@ -27,8 +27,6 @@ namespace PLMD {
class MahalanobisDistance : public ArgumentOnlyDistance {
public:
MahalanobisDistance( const ReferenceConfigurationOptions& ro );
- void read( const PDB& );
- double calc( const std::vector<Value*>& vals, const std::vector<double>& arg, const bool& squared );
};
PLUMED_REGISTER_METRIC(MahalanobisDistance,"MAHALANOBIS")
@@ -40,13 +38,4 @@ ArgumentOnlyDistance(ro)
hasmetric=true;
}
-void MahalanobisDistance::read( const PDB& pdb ){
- readArgumentsFromPDB( pdb );
-}
-
-double MahalanobisDistance::calc( const std::vector<Value*>& vals, const std::vector<double>& arg, const bool& squared ){
- return calculateArgumentDistance( vals, arg, squared );
-}
-
-
}
diff --git a/src/reference/MultiDomainRMSD.cpp b/src/reference/MultiDomainRMSD.cpp
index 017cc4203f9937267d11b38a529e973512661101..2d8c8a768b5991550171fc4103cf6f17e1d41f64 100644
--- a/src/reference/MultiDomainRMSD.cpp
+++ b/src/reference/MultiDomainRMSD.cpp
@@ -60,9 +60,9 @@ void MultiDomainRMSD::read( const PDB& pdb ){
if( !parse("UPPER_CUTTOFF"+num,upper,true) ) upper=std::numeric_limits<double>::max( );
}
domains.push_back( metricRegister().create<SingleDomainRMSD>( ftype ) );
- positions.resize( blocks[i] - blocks[i-1] + 1 );
- align.resize( blocks[i] - blocks[i-1] + 1 );
- displace.resize( blocks[i] - blocks[i-1] + 1 );
+ positions.resize( blocks[i] - blocks[i-1] );
+ align.resize( blocks[i] - blocks[i-1] );
+ displace.resize( blocks[i] - blocks[i-1] );
unsigned n=0;
for(unsigned j=blocks[i-1];j<blocks[i];++j){
positions[n]=pdb.getPositions()[j];
@@ -72,7 +72,7 @@ void MultiDomainRMSD::read( const PDB& pdb ){
}
domains[i-1]->setBoundsOnDistances( true, lower, upper ); // Currently no option for nopbc
domains[i-1]->setReferenceAtoms( positions, align, displace );
- domains[i-1]->setNumberOfAtoms( positions.size() );
+ // domains[i-1]->setNumberOfAtoms( positions.size() );
double ww=0; parse("WEIGHT"+num, ww, true );
if( ww==0 ) weights.push_back( 1.0 );
@@ -86,37 +86,38 @@ void MultiDomainRMSD::setReferenceAtoms( const std::vector<Vector>& conf, const
plumed_error();
}
-double MultiDomainRMSD::calculate( const std::vector<Vector>& pos, const Pbc& pbc, const bool& squared ){
- clearDerivatives(); double totd=0.; Tensor tvirial; std::vector<Vector> mypos;
+double MultiDomainRMSD::calculate( const std::vector<Vector>& pos, const Pbc& pbc, ReferenceValuePack& myder, const bool& squared ) const {
+ //clearDerivatives();
+ double totd=0.; Tensor tvirial; std::vector<Vector> mypos;
+ MultiValue tvals( 2, myder.getNumberOfDerivatives() ); ReferenceValuePack tder( 0, 0, tvals );
for(unsigned i=0;i<domains.size();++i){
// Must extract appropriate positions here
- mypos.resize( blocks[i+1] - blocks[i] + 1 );
- unsigned n=0;
- for(unsigned j=blocks[i];j<blocks[i+1];++j){ mypos[n]=pos[j]; n++; }
+ mypos.resize( blocks[i+1] - blocks[i] );
+ unsigned n=0; tder.resize( 0, mypos.size() );
+ for(unsigned j=blocks[i];j<blocks[i+1];++j){ tder.setAtomIndex(n,j); mypos[n]=pos[j]; n++; }
// This actually does the calculation
- totd += weights[i]*domains[i]->calculate( mypos, pbc, true );
-
- // Must extract derivatives here
- n=0;
- for(unsigned j=blocks[i];j<blocks[i+1];++j){
- addAtomicDerivatives( j, weights[i]*domains[i]->getAtomDerivative(n) ); n++;
- }
- if( domains[i]->getVirial( tvirial ) ){
- addBoxDerivatives( weights[i]*tvirial );
- }
+ totd += weights[i]*domains[i]->calculate( mypos, pbc, tder, true );
+ // Now merge the derivative
+ myder.copyScaledDerivatives( 1, weights[i], tvals );
+ // Make sure virial status is set correctly in output derivative pack
+ // This is only done here so I do this by using class friendship
+ if( tder.virialWasSet() ) myder.boxWasSet=true;
+ // And clear
+ tder.clear();
}
if( !squared ){
totd=sqrt(totd); double xx=0.5/totd;
- for(unsigned iat=0;iat<atom_ders.size();iat++) atom_ders[iat]*=xx;
- if( virialWasSet ) virial*=xx;
+ myder.scaleAllDerivatives( xx );
}
+ if( !myder.updateComplete() ) myder.updateDynamicLists();
return totd;
}
-double MultiDomainRMSD::calc( const std::vector<Vector>& pos, const Pbc& pbc, const std::vector<Value*>& vals, const std::vector<double>& arg, const bool& squared ){
+double MultiDomainRMSD::calc( const std::vector<Vector>& pos, const Pbc& pbc, const std::vector<Value*>& vals, const std::vector<double>& arg,
+ ReferenceValuePack& myder, const bool& squared ) const {
plumed_dbg_assert( vals.size()==0 && pos.size()==getNumberOfAtoms() && arg.size()==0 );
- return calculate( pos, pbc, squared );
+ return calculate( pos, pbc, myder, squared );
}
}
diff --git a/src/reference/MultiDomainRMSD.h b/src/reference/MultiDomainRMSD.h
index 957cfcd8c74365f0fc539cea10ff4a9e3751e219..876bfd46df7873b8efb7fd515eebadd0c83c2ae1 100644
--- a/src/reference/MultiDomainRMSD.h
+++ b/src/reference/MultiDomainRMSD.h
@@ -46,8 +46,8 @@ public:
/// Set the input from an analysis object (don't know how this will work yet so currently just a plumed_error)
void setReferenceAtoms( const std::vector<Vector>& conf, const std::vector<double>& align_in, const std::vector<double>& displace_in );
/// Calculate
- double calc( const std::vector<Vector>& pos, const Pbc& pbc, const std::vector<Value*>& vals, const std::vector<double>& arg, const bool& squared );
- double calculate( const std::vector<Vector>& pos, const Pbc& pbc, const bool& squared );
+ double calc( const std::vector<Vector>& pos, const Pbc& pbc, const std::vector<Value*>& vals, const std::vector<double>& arg, ReferenceValuePack& myder, const bool& squared ) const ;
+ double calculate( const std::vector<Vector>& pos, const Pbc& pbc, ReferenceValuePack& myder, const bool& squared ) const ;
};
}
diff --git a/src/reference/MultiReferenceBase.cpp b/src/reference/MultiReferenceBase.cpp
index bc010684566c1349bcc8b1ac9652187bfd29338a..fcfdd2337ccee5e687be1c70190e3cd3d0f35940 100644
--- a/src/reference/MultiReferenceBase.cpp
+++ b/src/reference/MultiReferenceBase.cpp
@@ -63,12 +63,12 @@ void MultiReferenceBase::getAtomAndArgumentRequirements( std::vector<AtomNumber>
}
}
-void MultiReferenceBase::setNumberOfAtomsAndArguments( const unsigned& natoms, const unsigned& nargs ){
- for(unsigned i=0;i<frames.size();++i){
- frames[i]->setNumberOfAtoms( natoms );
- frames[i]->setNumberOfArguments( nargs );
- }
-}
+// void MultiReferenceBase::setNumberOfAtomsAndArguments( const unsigned& natoms, const unsigned& nargs ){
+// for(unsigned i=0;i<frames.size();++i){
+// frames[i]->setNumberOfAtoms( natoms );
+// frames[i]->setNumberOfArguments( nargs );
+// }
+// }
void MultiReferenceBase::copyFrame( ReferenceConfiguration* frameToCopy ){
// Create a reference configuration of the appropriate type
diff --git a/src/reference/MultiReferenceBase.h b/src/reference/MultiReferenceBase.h
index 031407f12aa4b75d0d47812b491f35da4863d5cd..fddb71530216985fe3edd9fa885ba7b20236f827 100644
--- a/src/reference/MultiReferenceBase.h
+++ b/src/reference/MultiReferenceBase.h
@@ -53,7 +53,7 @@ public:
/// Find what is required of us from the reference frames
void getAtomAndArgumentRequirements( std::vector<AtomNumber>& atoms, std::vector<std::string>& args );
/// Finish setup of frames
- void setNumberOfAtomsAndArguments( const unsigned& natoms, const unsigned& nargs );
+// void setNumberOfAtomsAndArguments( const unsigned& natoms, const unsigned& nargs );
/// Do additional reading required by derived class
virtual void readRestOfFrame(){}
/// Do additional resizing required by derived class
@@ -62,7 +62,7 @@ public:
unsigned getNumberOfReferenceFrames() const ;
/// Calculate the distance from one of the reference points
double calcDistanceFromConfiguration( const unsigned& ifunc, const std::vector<Vector>& pos, const Pbc& pbc,
- const std::vector<Value*>& arg, const bool& squared );
+ const std::vector<Value*>& arg, ReferenceValuePack& myder, const bool& squared ) const ;
/// Return the ith reference frame
ReferenceConfiguration* getFrame( const unsigned& iframe );
/// Copy a reference configuration into the multi reference object
@@ -82,8 +82,8 @@ void MultiReferenceBase::parse(const std::string& key, T& val ){
inline
double MultiReferenceBase::calcDistanceFromConfiguration( const unsigned& ifunc, const std::vector<Vector>& pos, const Pbc& pbc,
- const std::vector<Value*>& arg, const bool& squared ){
- return frames[ifunc]->calculate( pos, pbc, arg, squared );
+ const std::vector<Value*>& arg, ReferenceValuePack& myder, const bool& squared ) const {
+ return frames[ifunc]->calculate( pos, pbc, arg, myder, squared );
}
inline
diff --git a/src/reference/NormalizedEuclideanDistance.cpp b/src/reference/NormalizedEuclideanDistance.cpp
index 4e82bc9419cb2b8d60ed2f38052a0b9b12afa6d4..51ee3be55fe3dea0191f5f6c0fc605c78dd42418 100644
--- a/src/reference/NormalizedEuclideanDistance.cpp
+++ b/src/reference/NormalizedEuclideanDistance.cpp
@@ -27,8 +27,6 @@ namespace PLMD {
class NormalizedEuclideanDistance : public ArgumentOnlyDistance {
public:
NormalizedEuclideanDistance( const ReferenceConfigurationOptions& ro );
- void read( const PDB& );
- double calc( const std::vector<Value*>& vals, const std::vector<double>& arg, const bool& squared );
};
PLUMED_REGISTER_METRIC(NormalizedEuclideanDistance,"NORM-EUCLIDEAN")
@@ -40,13 +38,4 @@ ArgumentOnlyDistance(ro)
hasweights=true;
}
-void NormalizedEuclideanDistance::read( const PDB& pdb ){
- readArgumentsFromPDB( pdb );
-}
-
-double NormalizedEuclideanDistance::calc( const std::vector<Value*>& vals, const std::vector<double>& arg, const bool& squared ){
- return calculateArgumentDistance( vals, arg, squared );
-}
-
-
}
diff --git a/src/reference/OptimalRMSD.cpp b/src/reference/OptimalRMSD.cpp
index 639a0c258a4db9d46cb7232c5d955847758dda64..a158d281d02f5f54dbc67ee5b7265c6d54f0472e 100644
--- a/src/reference/OptimalRMSD.cpp
+++ b/src/reference/OptimalRMSD.cpp
@@ -33,7 +33,7 @@ private:
public:
OptimalRMSD(const ReferenceConfigurationOptions& ro);
void read( const PDB& );
- double calc( const std::vector<Vector>& pos, const bool& squared );
+ double calc( const std::vector<Vector>& pos, ReferenceValuePack& myder, const bool& squared ) const ;
};
PLUMED_REGISTER_METRIC(OptimalRMSD,"OPTIMAL")
@@ -49,14 +49,18 @@ void OptimalRMSD::read( const PDB& pdb ){
readReference( pdb );
}
-double OptimalRMSD::calc( const std::vector<Vector>& pos, const bool& squared ){
+double OptimalRMSD::calc( const std::vector<Vector>& pos, ReferenceValuePack& myder, const bool& squared ) const {
+ double d; std::vector<Vector> atom_ders( pos.size() );
if( fast ){
- if( getAlign()==getDisplace() ) return myrmsd.optimalAlignment<false,true>(getAlign(),getDisplace(),pos,getReferencePositions(),atom_ders,squared);
- return myrmsd.optimalAlignment<false,false>(getAlign(),getDisplace(),pos,getReferencePositions(),atom_ders,squared);
+ if( getAlign()==getDisplace() ) d=myrmsd.optimalAlignment<false,true>(getAlign(),getDisplace(),pos,getReferencePositions(),atom_ders,squared);
+ d=myrmsd.optimalAlignment<false,false>(getAlign(),getDisplace(),pos,getReferencePositions(),atom_ders,squared);
} else {
- if( getAlign()==getDisplace() ) return myrmsd.optimalAlignment<true,true>(getAlign(),getDisplace(),pos,getReferencePositions(),atom_ders,squared);
- return myrmsd.optimalAlignment<true,false>(getAlign(),getDisplace(),pos,getReferencePositions(),atom_ders,squared);
+ if( getAlign()==getDisplace() ) d=myrmsd.optimalAlignment<true,true>(getAlign(),getDisplace(),pos,getReferencePositions(),atom_ders,squared);
+ d=myrmsd.optimalAlignment<true,false>(getAlign(),getDisplace(),pos,getReferencePositions(),atom_ders,squared);
}
+ for(unsigned i=0;i<atom_ders.size();++i) myder.addAtomDerivatives( i, atom_ders[i] );
+ if( !myder.updateComplete() ) myder.updateDynamicLists();
+ return d;
}
}
diff --git a/src/reference/PointWiseMapping.h b/src/reference/PointWiseMapping.h
index 7dc23d41c0318c10b68267a6d887b4cb43a8cc1c..b5358f5133cc70af63f38879635c67536e7ec14f 100644
--- a/src/reference/PointWiseMapping.h
+++ b/src/reference/PointWiseMapping.h
@@ -61,11 +61,11 @@ public:
/// Get the value of the ith property for th jth frame
double getPropertyValue( const unsigned& iframe, const unsigned& jprop ) const ;
/// Get the derivatives wrt to the position of an atom
- Vector getAtomDerivatives( const unsigned& iframe, const unsigned& jatom );
+// Vector getAtomDerivatives( const unsigned& iframe, const unsigned& jatom );
/// Get the derivatives wrt to the box
- bool getVirial( const unsigned& iframe, Tensor& vir );
+// bool getVirial( const unsigned& iframe, Tensor& vir );
/// Ge the derivatives wrt to one of the arguments
- double getArgumentDerivative( const unsigned& iframe, const unsigned& jarg );
+// double getArgumentDerivative( const unsigned& iframe, const unsigned& jarg );
/// Copy derivative information from frame number from to frame number to
void copyFrameDerivatives( const unsigned& from, const unsigned& to );
/// Get a pointer to the matrix of pairwise distances
@@ -112,20 +112,20 @@ double PointWiseMapping::getPropertyValue( const unsigned& iframe, const unsigne
return low_dim[iframe][jprop];
}
-inline
-Vector PointWiseMapping::getAtomDerivatives( const unsigned& iframe, const unsigned& jatom ){
- return frames[iframe]->getAtomDerivative(jatom);
-}
-
-inline
-bool PointWiseMapping::getVirial( const unsigned& iframe, Tensor& vir ){
- return frames[iframe]->getVirial( vir );
-}
-
-inline
-double PointWiseMapping::getArgumentDerivative( const unsigned& iframe, const unsigned& jarg ){
- return frames[iframe]->getArgumentDerivative(jarg);
-}
+// inline
+// Vector PointWiseMapping::getAtomDerivatives( const unsigned& iframe, const unsigned& jatom ){
+// return frames[iframe]->getAtomDerivative(jatom);
+// }
+//
+// inline
+// bool PointWiseMapping::getVirial( const unsigned& iframe, Tensor& vir ){
+// return frames[iframe]->getVirial( vir );
+// }
+
+// inline
+// double PointWiseMapping::getArgumentDerivative( const unsigned& iframe, const unsigned& jarg ){
+// return frames[iframe]->getArgumentDerivative(jarg);
+// }
inline
Matrix<double>& PointWiseMapping::modifyDmat(){
diff --git a/src/reference/RMSDBase.cpp b/src/reference/RMSDBase.cpp
index d8e436962f83b7dc9027e83730fe9f6bd9200bb9..e562fb851160cc3335804311c044414a80f82eb5 100644
--- a/src/reference/RMSDBase.cpp
+++ b/src/reference/RMSDBase.cpp
@@ -29,14 +29,14 @@ SingleDomainRMSD(ro)
{
}
-double RMSDBase::calculate( const std::vector<Vector>& pos, const bool& squared ){
- clearDerivatives();
- return calc( pos, squared );
+double RMSDBase::calculate( const std::vector<Vector>& pos, ReferenceValuePack& myder, const bool& squared ) const {
+// clearDerivatives();
+ return calc( pos, myder, squared );
}
-double RMSDBase::calc( const std::vector<Vector>& pos, const Pbc& pbc, const bool& squared ){
+double RMSDBase::calc( const std::vector<Vector>& pos, const Pbc& pbc, ReferenceValuePack& myder, const bool& squared ) const {
plumed_dbg_assert( pos.size()==getNumberOfAtoms() );
- return calc( pos, squared );
+ return calc( pos, myder, squared );
}
}
diff --git a/src/reference/RMSDBase.h b/src/reference/RMSDBase.h
index f722acd1dee920586dad22f1e2db5269f7c2606d..b055e467f5e3f0568283ecbcbb9f6cd1191011ed 100644
--- a/src/reference/RMSDBase.h
+++ b/src/reference/RMSDBase.h
@@ -33,9 +33,9 @@ class Pbc;
class RMSDBase : public SingleDomainRMSD {
public:
RMSDBase( const ReferenceConfigurationOptions& ro );
- double calc( const std::vector<Vector>& pos, const Pbc& pbc, const bool& squared );
- double calculate( const std::vector<Vector>& pos, const bool& squared );
- virtual double calc( const std::vector<Vector>& pos, const bool& squared )=0;
+ double calc( const std::vector<Vector>& pos, const Pbc& pbc, ReferenceValuePack& myder, const bool& squared ) const;
+ double calculate( const std::vector<Vector>& pos, ReferenceValuePack& myder, const bool& squared ) const ;
+ virtual double calc( const std::vector<Vector>& pos, ReferenceValuePack& myder, const bool& squared ) const=0;
};
}
diff --git a/src/reference/ReferenceArguments.cpp b/src/reference/ReferenceArguments.cpp
index 3788d23938f32fe1441e454e9a6334e8585a8e71..3e30caed5bc5f3cb48e78103b5092be023518d53 100644
--- a/src/reference/ReferenceArguments.cpp
+++ b/src/reference/ReferenceArguments.cpp
@@ -159,8 +159,9 @@ const std::vector<double>& ReferenceArguments::getReferenceMetric(){
return trig_metric;
}
-double ReferenceArguments::calculateArgumentDistance( const std::vector<Value*> & vals, const std::vector<double>& arg, const bool& squared ){
- double r=0;
+double ReferenceArguments::calculateArgumentDistance( const std::vector<Value*> & vals, const std::vector<double>& arg,
+ ReferenceValuePack& myder, const bool& squared ) const {
+ double r=0; std::vector<double> arg_ders( vals.size() );
if( hasmetric ){
double dp_i, dp_j;
for(unsigned i=0;i<reference_args.size();++i){
@@ -185,7 +186,7 @@ double ReferenceArguments::calculateArgumentDistance( const std::vector<Value*>
}
if(!squared){
r=sqrt(r); double ir=1.0/(2.0*r);
- for(unsigned i=0;i<arg_ders.size();++i) arg_ders[i]*=ir;
+ for(unsigned i=0;i<arg_ders.size();++i) myder.addArgumentDerivatives( i, arg_ders[i]*ir );
}
return r;
}
diff --git a/src/reference/ReferenceArguments.h b/src/reference/ReferenceArguments.h
index 92ee822178574f680d2103589460dc87cb5780ac..b2f8bde4b03b2633d394f11d92f6c54ccfedfa95 100644
--- a/src/reference/ReferenceArguments.h
+++ b/src/reference/ReferenceArguments.h
@@ -66,9 +66,11 @@ protected:
void setReferenceArguments();
/// Calculate the euclidean/malanobius distance the atoms have moved from the reference
/// configuration in CV space
- double calculateArgumentDistance( const std::vector<Value*> & vals, const std::vector<double>& arg, const bool& squared );
+ double calculateArgumentDistance( const std::vector<Value*> & vals, const std::vector<double>& arg, ReferenceValuePack& myder, const bool& squared ) const ;
public:
ReferenceArguments( const ReferenceConfigurationOptions& ro );
+/// Get the number of reference arguments
+ unsigned getNumberOfReferenceArguments() const ;
/// Get the arguments required
void getArgumentRequests( std::vector<std::string>&, bool disable_checks=false );
/// Set the names of the arguments
@@ -102,6 +104,11 @@ const std::vector<std::string>& ReferenceArguments::getArgumentNames(){
return arg_names;
}
+inline
+unsigned ReferenceArguments::getNumberOfReferenceArguments() const {
+ return reference_args.size();
+}
+
}
#endif
diff --git a/src/reference/ReferenceAtoms.h b/src/reference/ReferenceAtoms.h
index acbc69f8a755884c66ad67d91ba1888da3a44c2d..69af7004facad6ebc65e89014e4be60a8b2e9fdf 100644
--- a/src/reference/ReferenceAtoms.h
+++ b/src/reference/ReferenceAtoms.h
@@ -75,23 +75,21 @@ protected:
/// Get the position of the ith atom
Vector getReferencePosition( const unsigned& iatom ) const ;
/// Get the reference positions
- const std::vector<Vector> & getReferencePositions();
+ const std::vector<Vector> & getReferencePositions() const ;
/// Add derivatives to iatom th atom in list
- void addAtomicDerivatives( const unsigned& , const Vector& );
+// void addAtomicDerivatives( const unsigned& , const Vector& );
/// Get the atomic derivatives on the ith atom in the list
- Vector retrieveAtomicDerivatives( const unsigned& ) const ;
+// Vector retrieveAtomicDerivatives( const unsigned& ) const ;
/// Add derivatives to the viral
- void addBoxDerivatives( const Tensor& );
+// void addBoxDerivatives( const Tensor& );
/// This does the checks that are always required
void singleDomainRequests( std::vector<AtomNumber>&, bool disable_checks );
+public:
+ ReferenceAtoms( const ReferenceConfigurationOptions& ro );
/// This returns the number of reference atom positions
unsigned getNumberOfReferencePositions() const ;
-/// This returns how many atoms there should be
- unsigned getNumberOfAtoms() const ;
/// This allows us to use a single pos array with RMSD objects using different atom indexes
unsigned getAtomIndex( const unsigned& ) const ;
-public:
- ReferenceAtoms( const ReferenceConfigurationOptions& ro );
/// Get the atoms required (additional checks are required when we have multiple domains)
virtual void getAtomRequests( std::vector<AtomNumber>&, bool disable_checks=false );
/// Set the indices of the reference atoms
@@ -102,6 +100,8 @@ public:
void printAtoms( OFile& ofile ) const ;
/// Return all atom indexes
const std::vector<AtomNumber>& getAbsoluteIndexes();
+/// This returns how many atoms there should be
+ unsigned getNumberOfAtoms() const ;
};
inline
@@ -122,13 +122,13 @@ unsigned ReferenceAtoms::getNumberOfReferencePositions() const {
inline
unsigned ReferenceAtoms::getNumberOfAtoms() const {
- return atom_ders.size();
+ return reference_atoms.size();
}
inline
unsigned ReferenceAtoms::getAtomIndex( const unsigned& iatom ) const {
plumed_dbg_assert( iatom<der_index.size() );
- plumed_dbg_assert( der_index[iatom]<atom_ders.size() );
+ plumed_dbg_assert( der_index[iatom]<reference_atoms.size() );
return der_index[iatom];
}
@@ -139,24 +139,24 @@ Vector ReferenceAtoms::getReferencePosition( const unsigned& iatom ) const {
}
inline
-const std::vector<Vector> & ReferenceAtoms::getReferencePositions(){
+const std::vector<Vector> & ReferenceAtoms::getReferencePositions() const {
return reference_atoms;
}
-inline
-void ReferenceAtoms::addAtomicDerivatives( const unsigned& iatom, const Vector& der ){
- atom_ders[ getAtomIndex(iatom) ]+=der;
-}
+// inline
+// void ReferenceAtoms::addAtomicDerivatives( const unsigned& iatom, const Vector& der ){
+// atom_ders[ getAtomIndex(iatom) ]+=der;
+// }
-inline
-Vector ReferenceAtoms::retrieveAtomicDerivatives( const unsigned& iatom ) const {
- return atom_ders[ getAtomIndex(iatom) ];
-}
+// inline
+// Vector ReferenceAtoms::retrieveAtomicDerivatives( const unsigned& iatom ) const {
+// return atom_ders[ getAtomIndex(iatom) ];
+// }
-inline
-void ReferenceAtoms::addBoxDerivatives( const Tensor& vir ){
- virialWasSet=true; virial+=vir;
-}
+// inline
+// void ReferenceAtoms::addBoxDerivatives( const Tensor& vir ){
+// virialWasSet=true; virial+=vir;
+// }
inline
const std::vector<AtomNumber>& ReferenceAtoms::getAbsoluteIndexes(){
diff --git a/src/reference/ReferenceConfiguration.cpp b/src/reference/ReferenceConfiguration.cpp
index f5ec1663995b137bf3511b67b303528fc6410194..7cb7d510d1a29ef77ee8b2c03afdff677f9777f4 100644
--- a/src/reference/ReferenceConfiguration.cpp
+++ b/src/reference/ReferenceConfiguration.cpp
@@ -44,9 +44,9 @@ std::string ReferenceConfigurationOptions::getMultiRMSDType() const {
}
ReferenceConfiguration::ReferenceConfiguration( const ReferenceConfigurationOptions& ro ):
-name(ro.tt),
-arg_ders(0),
-atom_ders(0)
+name(ro.tt)
+// arg_ders(0),
+// atom_ders(0)
{
}
@@ -68,18 +68,18 @@ void ReferenceConfiguration::set( const PDB& pdb ){
read( pdb );
}
-void ReferenceConfiguration::setNumberOfArguments( const unsigned& n ){
- arg_ders.resize(n); tmparg.resize(n);
-}
+// void ReferenceConfiguration::setNumberOfArguments( const unsigned& n ){
+// arg_ders.resize(n); tmparg.resize(n);
+// }
-void ReferenceConfiguration::setNumberOfAtoms( const unsigned& n ){
- atom_ders.resize(n);
-}
+// void ReferenceConfiguration::setNumberOfAtoms( const unsigned& n ){
+// atom_ders.resize(n);
+// }
-bool ReferenceConfiguration::getVirial( Tensor& virout ) const {
- if(virialWasSet) virout=virial;
- return virialWasSet;
-}
+// bool ReferenceConfiguration::getVirial( Tensor& virout ) const {
+// if(virialWasSet) virout=virial;
+// return virialWasSet;
+// }
void ReferenceConfiguration::parseFlag( const std::string&key, bool&t ){
Tools::parseFlag(line,key,t);
@@ -93,24 +93,24 @@ void ReferenceConfiguration::setNamesAndAtomNumbers( const std::vector<AtomNumbe
ReferenceAtoms* atoms=dynamic_cast<ReferenceAtoms*>( this );
if(!atoms){
plumed_massert( numbers.size()==0, "expecting no atomic positions");
- setNumberOfAtoms( 0 );
+ //setNumberOfAtoms( 0 );
} else {
atoms->setAtomNumbers( numbers );
- setNumberOfAtoms( numbers.size() );
+ // setNumberOfAtoms( numbers.size() );
}
// Copy the arguments to the reference
ReferenceArguments* args=dynamic_cast<ReferenceArguments*>( this );
if(!args){
plumed_massert( arg.size()==0, "expecting no arguments");
- setNumberOfArguments(0);
+ // setNumberOfArguments(0);
} else {
args->setArgumentNames( arg );
- setNumberOfArguments( arg.size() );
+ // setNumberOfArguments( arg.size() );
}
}
void ReferenceConfiguration::setReferenceConfig( const std::vector<Vector>& pos, const std::vector<double>& arg, const std::vector<double>& metric ){
- plumed_dbg_assert( pos.size()==atom_ders.size() && arg.size()==arg_ders.size() );
+// plumed_dbg_assert( pos.size()==atom_ders.size() && arg.size()==arg_ders.size() );
// Copy the atomic positions to the reference
ReferenceAtoms* atoms=dynamic_cast<ReferenceAtoms*>( this );
if(!atoms){
@@ -136,24 +136,26 @@ void ReferenceConfiguration::checkRead(){
}
}
-void ReferenceConfiguration::clearDerivatives(){
- for(unsigned i=0;i<atom_ders.size();++i) atom_ders[i].zero();
- virial.zero(); virialWasSet=false;
- arg_ders.assign(arg_ders.size(),0.0);
-}
+// void ReferenceConfiguration::clearDerivatives(){
+// for(unsigned i=0;i<atom_ders.size();++i) atom_ders[i].zero();
+// virial.zero(); virialWasSet=false;
+// arg_ders.assign(arg_ders.size(),0.0);
+// }
-double ReferenceConfiguration::calculate( const std::vector<Vector>& pos, const Pbc& pbc, const std::vector<Value*>& vals, const bool& squared ){
- clearDerivatives();
+double ReferenceConfiguration::calculate( const std::vector<Vector>& pos, const Pbc& pbc, const std::vector<Value*>& vals,
+ ReferenceValuePack& myder, const bool& squared ) const {
+ // clearDerivatives();
+ std::vector<double> tmparg( vals.size() );
for(unsigned i=0;i<vals.size();++i) tmparg[i]=vals[i]->get();
- return calc( pos, pbc, vals, tmparg, squared );
+ return calc( pos, pbc, vals, tmparg, myder, squared );
}
-void ReferenceConfiguration::copyDerivatives( const ReferenceConfiguration* ref ){
- plumed_dbg_assert( ref->atom_ders.size()==atom_ders.size() && ref->arg_ders.size()==arg_ders.size() );
- for(unsigned i=0;i<atom_ders.size();++i) atom_ders[i]=ref->atom_ders[i];
- for(unsigned i=0;i<arg_ders.size();++i) arg_ders[i]=ref->arg_ders[i];
- virialWasSet=ref->virialWasSet; virial=ref->virial;
-}
+// void ReferenceConfiguration::copyDerivatives( const ReferenceConfiguration* ref ){
+// plumed_dbg_assert( ref->atom_ders.size()==atom_ders.size() && ref->arg_ders.size()==arg_ders.size() );
+// for(unsigned i=0;i<atom_ders.size();++i) atom_ders[i]=ref->atom_ders[i];
+// for(unsigned i=0;i<arg_ders.size();++i) arg_ders[i]=ref->arg_ders[i];
+// virialWasSet=ref->virialWasSet; virial=ref->virial;
+// }
void ReferenceConfiguration::print( OFile& ofile, const double& time, const double& weight, const double& old_norm ){
ofile.printf("REMARK TIME=%f LOG_WEIGHT=%f OLD_NORM=%f\n",time, weight, old_norm );
@@ -169,10 +171,12 @@ void ReferenceConfiguration::print( OFile& ofile, const std::string& fmt ){
}
double distance( const Pbc& pbc, const std::vector<Value*> & vals, ReferenceConfiguration* ref1, ReferenceConfiguration* ref2, const bool& squared ){
- double dist1=ref1->calc( ref2->getReferencePositions(), pbc, vals, ref2->getReferenceArguments(), squared );
+ MultiValue myvals( 2, ref1->getReferenceArguments().size() + 3*ref1->getReferencePositions().size() + 9 );
+ ReferenceValuePack myder( ref1->getReferenceArguments().size() , ref1->getReferencePositions().size() , myvals );
+ double dist1=ref1->calc( ref2->getReferencePositions(), pbc, vals, ref2->getReferenceArguments(), myder, squared );
#ifndef NDEBUG
// Check that A - B = B - A
- double dist2=ref2->calc( ref1->getReferencePositions(), pbc, vals, ref1->getReferenceArguments(), squared );
+ double dist2=ref2->calc( ref1->getReferencePositions(), pbc, vals, ref1->getReferenceArguments(), myder, squared );
plumed_dbg_assert( fabs(dist1-dist2)<epsilon );
#endif
return dist1;
diff --git a/src/reference/ReferenceConfiguration.h b/src/reference/ReferenceConfiguration.h
index 2cf875ba2748e499b30349eb3acfbd02406281ae..a86120d152060951f12e4f20fdb3977212503ed4 100644
--- a/src/reference/ReferenceConfiguration.h
+++ b/src/reference/ReferenceConfiguration.h
@@ -28,6 +28,7 @@
#include "tools/Tensor.h"
#include "tools/Tools.h"
#include "tools/Exception.h"
+#include "ReferenceValuePack.h"
namespace PLMD{
@@ -65,7 +66,7 @@ private:
/// A vector containing all the remarks from the pdb input
std::vector<std::string> line;
/// This is used to store the values of arguments
- std::vector<double> tmparg;
+// std::vector<double> tmparg;
/// These are used to do fake things when we copy frames
std::vector<AtomNumber> fake_atom_numbers;
std::vector<std::string> fake_arg_names;
@@ -75,12 +76,12 @@ private:
std::vector<double> fake_metric;
protected:
/// Derivatives wrt to the arguments
- std::vector<double> arg_ders;
+// std::vector<double> arg_ders;
/// The virial contribution has to be stored
- bool virialWasSet;
- Tensor virial;
+// bool virialWasSet;
+// Tensor virial;
/// Derivatives wrt to the atoms
- std::vector<Vector> atom_ders;
+// std::vector<Vector> atom_ders;
/// Crash with an error
void error(const std::string& msg);
/// Clear the derivatives
@@ -91,14 +92,17 @@ public:
virtual ~ReferenceConfiguration();
/// Return the name of this metric
std::string getName() const ;
+///
+ virtual unsigned getNumberOfReferencePositions() const ;
+ virtual unsigned getNumberOfReferenceArguments() const ;
/// Retrieve the atoms that are required for this guy
virtual void getAtomRequests( std::vector<AtomNumber>&, bool disable_checks=false ){}
/// Retrieve the arguments that are required for this guy
virtual void getArgumentRequests( std::vector<std::string>&, bool disable_checks=false ){}
/// Set the final number of arguments
- virtual void setNumberOfArguments( const unsigned& );
+// virtual void setNumberOfArguments( const unsigned& );
/// Set the final number of atoms
- virtual void setNumberOfAtoms( const unsigned& );
+// virtual void setNumberOfAtoms( const unsigned& );
/// Set the reference configuration using a PDB
virtual void set( const PDB& );
/// Do all local business for setting the configuration
@@ -108,17 +112,18 @@ public:
/// Return the weight for this frame
double getWeight() const ;
/// Calculate the distance from the reference configuration
- double calculate( const std::vector<Vector>& pos, const Pbc& pbc, const std::vector<Value*>& vals, const bool& squared=false );
+ double calculate( const std::vector<Vector>& pos, const Pbc& pbc, const std::vector<Value*>& vals, ReferenceValuePack& myder, const bool& squared=false ) const ;
/// Calculate the distance from the reference configuration
- virtual double calc( const std::vector<Vector>& pos, const Pbc& pbc, const std::vector<Value*>& vals, const std::vector<double>& args, const bool& squared )=0;
-/// Return the derivative wrt to the ith atom
- Vector getAtomDerivative( const unsigned& ) const ;
-/// Return the derivative wrt to the ith argument
- double getArgumentDerivative( const unsigned& ) const ;
+ virtual double calc( const std::vector<Vector>& pos, const Pbc& pbc, const std::vector<Value*>& vals, const std::vector<double>& args,
+ ReferenceValuePack& myder, const bool& squared ) const=0;
+// /// Return the derivative wrt to the ith atom
+// Vector getAtomDerivative( const unsigned& ) const ;
+// /// Return the derivative wrt to the ith argument
+// double getArgumentDerivative( const unsigned& ) const ;
/// Return the derivatives of the distance wrt the cell vectors. This returns false
/// for everything other than DRMSD as these sort of calculations have to be done
/// separately when you use RMSD
- bool getVirial( Tensor& virout ) const ;
+// bool getVirial( Tensor& virout ) const ;
/// Parse something from the pdb remarks
template<class T>
bool parse( const std::string&key, T&t, bool ignore_missing=false );
@@ -142,7 +147,7 @@ public:
virtual double getReferenceArgument( const unsigned& i ){ plumed_error(); return 0.0; }
/// These are overwritten in ReferenceArguments and ReferenceAtoms but are required here
/// to make PLMD::distance work
- virtual const std::vector<Vector>& getReferencePositions();
+ virtual const std::vector<Vector>& getReferencePositions() const ;
virtual const std::vector<double>& getReferenceArguments();
virtual const std::vector<double>& getReferenceMetric();
/// These are overwritten in ReferenceArguments and ReferenceAtoms to make frame copying work
@@ -150,17 +155,17 @@ public:
virtual const std::vector<std::string>& getArgumentNames();
};
-inline
-Vector ReferenceConfiguration::getAtomDerivative( const unsigned& ider ) const {
- plumed_dbg_assert( ider<atom_ders.size() );
- return atom_ders[ider];
-}
+// inline
+// Vector ReferenceConfiguration::getAtomDerivative( const unsigned& ider ) const {
+// plumed_dbg_assert( ider<atom_ders.size() );
+// return atom_ders[ider];
+// }
-inline
-double ReferenceConfiguration::getArgumentDerivative( const unsigned& ider ) const {
- plumed_dbg_assert( ider<arg_ders.size() );
- return arg_ders[ider];
-}
+// inline
+// double ReferenceConfiguration::getArgumentDerivative( const unsigned& ider ) const {
+// plumed_dbg_assert( ider<arg_ders.size() );
+// return arg_ders[ider];
+// }
inline
void ReferenceConfiguration::setWeight( const double& ww ){
@@ -187,7 +192,7 @@ bool ReferenceConfiguration::parseVector(const std::string&key,std::vector<T>&t,
}
inline
-const std::vector<Vector>& ReferenceConfiguration::getReferencePositions(){
+const std::vector<Vector>& ReferenceConfiguration::getReferencePositions() const {
return fake_refatoms;
}
@@ -211,6 +216,16 @@ const std::vector<std::string>& ReferenceConfiguration::getArgumentNames(){
return fake_arg_names;
}
+inline
+unsigned ReferenceConfiguration::getNumberOfReferencePositions() const {
+ return 0;
+}
+
+inline
+unsigned ReferenceConfiguration::getNumberOfReferenceArguments() const {
+ return 0;
+}
+
}
diff --git a/src/reference/ReferenceValuePack.cpp b/src/reference/ReferenceValuePack.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7529674bae594c8745f390a43d4be75ccf398b08
--- /dev/null
+++ b/src/reference/ReferenceValuePack.cpp
@@ -0,0 +1,87 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2013,2014 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 "ReferenceValuePack.h"
+
+namespace PLMD {
+
+ReferenceValuePack::ReferenceValuePack( const unsigned& nargs, const unsigned& natoms, MultiValue& vals ):
+boxWasSet(false),
+numberOfArgs(nargs),
+oind(1),
+myvals(vals)
+{
+ atom_indices.resize( natoms );
+}
+
+void ReferenceValuePack::resize( const unsigned& nargs, const unsigned& natoms ){
+ numberOfArgs=nargs; atom_indices.resize( natoms );
+}
+
+void ReferenceValuePack::updateDynamicLists(){
+ myvals.emptyActiveMembers();
+ for(unsigned i=0;i<numberOfArgs;++i) myvals.updateIndex( i );
+ for(unsigned i=0;i<atom_indices.size();++i){
+ unsigned nbase = numberOfArgs + 3*atom_indices[i];
+ myvals.updateIndex( nbase+0 ); myvals.updateIndex( nbase+1 ); myvals.updateIndex( nbase+2 );
+ }
+ unsigned nbase = myvals.getNumberOfDerivatives() - 9;
+ // zero is added to all virial components to ensure that these are active in the dynamic list
+ // if this is not done there is a problem with secondary structure variables
+ if( atom_indices.size()>0 ){
+ for(unsigned i=0;i<9;++i) myvals.addDerivative( oind, nbase+i, 0.0 );
+ }
+ for(unsigned i=0;i<9;++i) myvals.updateIndex( nbase+i );
+ myvals.sortActiveList();
+}
+
+void ReferenceValuePack::clear(){
+ if( !myvals.updateComplete() ) updateDynamicLists();
+ myvals.clearAll(); boxWasSet=false;
+}
+
+void ReferenceValuePack::scaleAllDerivatives( const double& scalef ){
+ if( !myvals.updateComplete() ) updateDynamicLists();
+
+ for(unsigned i=0;i<myvals.getNumberActive();++i){
+ unsigned ider=myvals.getActiveIndex(i);
+ myvals.setDerivative( oind, ider, scalef*myvals.getDerivative( oind, ider ) );
+ }
+}
+
+void ReferenceValuePack::copyScaledDerivatives( const unsigned& from, const double& scalef, const MultiValue& tvals ){
+ plumed_dbg_assert( tvals.getNumberOfDerivatives()==myvals.getNumberOfDerivatives() );
+ for(unsigned i=0;i<tvals.getNumberActive();++i){
+ unsigned ider=tvals.getActiveIndex(i);
+ myvals.addDerivative( oind, ider, scalef*tvals.getDerivative( from, ider ) );
+ }
+}
+
+void ReferenceValuePack::moveDerivatives( const unsigned& from, const unsigned& to ){
+ if( !myvals.updateComplete() ) updateDynamicLists();
+
+ for(unsigned i=0;i<myvals.getNumberActive();++i){
+ unsigned ider=myvals.getActiveIndex(i);
+ myvals.setDerivative( to, ider, myvals.getDerivative( from, ider ) );
+ }
+}
+
+}
diff --git a/src/reference/ReferenceValuePack.h b/src/reference/ReferenceValuePack.h
new file mode 100644
index 0000000000000000000000000000000000000000..4c39b8bcb8370773428a879f2981d8d6aed8a4fe
--- /dev/null
+++ b/src/reference/ReferenceValuePack.h
@@ -0,0 +1,170 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2013,2014 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_reference_ReferenceValuePack_h
+#define __PLUMED_reference_ReferenceValuePack_h
+
+#include "tools/MultiValue.h"
+
+namespace PLMD {
+
+class ReferenceValuePack {
+friend class MultiDomainRMSD;
+private:
+/// Was the virial set
+ bool boxWasSet;
+///
+ unsigned numberOfArgs;
+///
+ unsigned oind;
+/// Copy of the values that we are adding to
+ MultiValue& myvals;
+/// Ths list of atom indices
+ std::vector<unsigned> atom_indices;
+public:
+ ReferenceValuePack( const unsigned& nargs, const unsigned& natoms, MultiValue& vals );
+///
+ void resize( const unsigned& nargs, const unsigned& natoms );
+///
+ void clear();
+///
+ unsigned getNumberOfDerivatives() const ;
+///
+ unsigned getNumberOfArguments() const ;
+///
+ unsigned getNumberOfAtoms() const ;
+///
+ void setAtomIndex( const unsigned& iatom, const unsigned& jindex );
+///
+ unsigned getAtomIndex( const unsigned& iatom ) const ;
+///
+ void copyScaledDerivatives( const unsigned& from, const double& scalef, const MultiValue& tvals );
+///
+ void addArgumentDerivatives( const unsigned& iarg, const double& der );
+///
+ void addAtomDerivatives( const unsigned& iatom, const Vector& der );
+///
+ void addBoxDerivatives( const Tensor& vir );
+///
+ bool updateComplete() const ;
+///
+ void updateDynamicLists();
+///
+ void scaleAllDerivatives( const double& scalef );
+///
+ void setValIndex( const unsigned& ind );
+///
+ void moveDerivatives( const unsigned& from, const unsigned& to );
+///
+ bool virialWasSet() const ;
+///
+ Vector getAtomDerivative( const unsigned& iatom ) const ;
+///
+ double getArgumentDerivative( const unsigned& ival ) const ;
+///
+ Tensor getBoxDerivatives() const ;
+};
+
+inline
+unsigned ReferenceValuePack::getNumberOfDerivatives() const {
+ return myvals.getNumberOfDerivatives();
+}
+
+inline
+unsigned ReferenceValuePack::getNumberOfArguments() const {
+ return numberOfArgs;
+}
+
+inline
+unsigned ReferenceValuePack::getNumberOfAtoms() const {
+ return atom_indices.size();
+}
+
+inline
+void ReferenceValuePack::setAtomIndex( const unsigned& iatom, const unsigned& jindex ){
+ plumed_dbg_assert( iatom<atom_indices.size() ); atom_indices[iatom]=jindex;
+}
+
+inline
+unsigned ReferenceValuePack::getAtomIndex( const unsigned& iatom ) const {
+ plumed_dbg_assert( iatom<atom_indices.size() );
+ return atom_indices[iatom];
+}
+
+inline
+void ReferenceValuePack::addArgumentDerivatives( const unsigned& iarg, const double& der ){
+ plumed_dbg_assert( iarg<numberOfArgs ); myvals.addDerivative( oind, iarg, der );
+}
+
+inline
+bool ReferenceValuePack::updateComplete() const {
+ return myvals.updateComplete();
+}
+
+inline
+void ReferenceValuePack::addAtomDerivatives( const unsigned& jder, const Vector& der ){
+ plumed_dbg_assert( jder<atom_indices.size() );
+ myvals.addDerivative( oind, numberOfArgs + 3*atom_indices[jder] + 0, der[0] );
+ myvals.addDerivative( oind, numberOfArgs + 3*atom_indices[jder] + 1, der[1] );
+ myvals.addDerivative( oind, numberOfArgs + 3*atom_indices[jder] + 2, der[2] );
+}
+
+inline
+void ReferenceValuePack::addBoxDerivatives( const Tensor& vir ){
+ boxWasSet=true; unsigned nbase = myvals.getNumberOfDerivatives() - 9;
+ for(unsigned i=0;i<3;++i) for(unsigned j=0;j<3;++j) myvals.addDerivative( oind, nbase + 3*i + j , vir(i,j) );
+}
+
+inline
+void ReferenceValuePack::setValIndex( const unsigned& ind ){
+ oind=ind;
+}
+
+inline
+bool ReferenceValuePack::virialWasSet() const {
+ return boxWasSet;
+}
+
+inline
+Vector ReferenceValuePack::getAtomDerivative( const unsigned& iatom ) const {
+ Vector tmp; plumed_dbg_assert( iatom<atom_indices.size() );
+ tmp[0]=myvals.getDerivative( oind, numberOfArgs + 3*atom_indices[iatom] + 0 );
+ tmp[1]=myvals.getDerivative( oind, numberOfArgs + 3*atom_indices[iatom] + 1 );
+ tmp[2]=myvals.getDerivative( oind, numberOfArgs + 3*atom_indices[iatom] + 2 );
+ return tmp;
+}
+
+inline
+double ReferenceValuePack::getArgumentDerivative( const unsigned& ival ) const {
+ plumed_dbg_assert( ival<numberOfArgs );
+ return myvals.getDerivative( oind, ival );
+}
+
+inline
+Tensor ReferenceValuePack::getBoxDerivatives() const {
+ plumed_dbg_assert( boxWasSet ); Tensor tvir; unsigned nbase = myvals.getNumberOfDerivatives() - 9;
+ for(unsigned i=0;i<3;++i) for(unsigned j=0;j<3;++j) tvir(i,j)=myvals.getDerivative( oind, nbase + 3*i + j );
+ return tvir;
+}
+
+}
+
+#endif
diff --git a/src/reference/SimpleRMSD.cpp b/src/reference/SimpleRMSD.cpp
index dfca47a4902cd79e2251019853fbad0661bd20c0..23391e1e4d6017de3dba3c822fe59a5c2085b504 100644
--- a/src/reference/SimpleRMSD.cpp
+++ b/src/reference/SimpleRMSD.cpp
@@ -31,7 +31,7 @@ private:
public:
SimpleRMSD( const ReferenceConfigurationOptions& ro );
void read( const PDB& );
- double calc( const std::vector<Vector>& pos, const bool& squared );
+ double calc( const std::vector<Vector>& pos, ReferenceValuePack& myder, const bool& squared ) const ;
};
PLUMED_REGISTER_METRIC(SimpleRMSD,"SIMPLE")
@@ -46,8 +46,12 @@ void SimpleRMSD::read( const PDB& pdb ){
readReference( pdb );
}
-double SimpleRMSD::calc( const std::vector<Vector>& pos, const bool& squared ){
- return myrmsd.simpleAlignment( getAlign(), getDisplace(), pos, getReferencePositions(), atom_ders, squared );
+double SimpleRMSD::calc( const std::vector<Vector>& pos, ReferenceValuePack& myder, const bool& squared ) const {
+ std::vector<Vector> atom_ders( pos.size() );
+ double d=myrmsd.simpleAlignment( getAlign(), getDisplace(), pos, getReferencePositions(), atom_ders, squared );
+ for(unsigned i=0;i<pos.size();++i) myder.addAtomDerivatives( i, atom_ders[i] );
+ if( !myder.updateComplete() ) myder.updateDynamicLists();
+ return d;
}
}
diff --git a/src/reference/SingleDomainRMSD.cpp b/src/reference/SingleDomainRMSD.cpp
index 6e13ff9148cef6dfc16a8da1fb3a76287ff44649..dbfcef053766a5fb3db58d4a83a70231749c98ae 100644
--- a/src/reference/SingleDomainRMSD.cpp
+++ b/src/reference/SingleDomainRMSD.cpp
@@ -56,17 +56,17 @@ void SingleDomainRMSD::setReferenceAtoms( const std::vector<Vector>& conf, const
center+=conf[i]*align[i]; der_index[i]=i;
}
for(unsigned i=0;i<conf.size();++i) reference_atoms[i]=conf[i]-center;
- setNumberOfAtoms( conf.size() ); setNumberOfArguments( 0 );
+ // setNumberOfAtoms( conf.size() ); setNumberOfArguments( 0 );
}
-double SingleDomainRMSD::calculate( const std::vector<Vector>& pos, const Pbc& pbc, const bool& squared ){
- clearDerivatives();
- return calc( pos, pbc, squared );
+double SingleDomainRMSD::calculate( const std::vector<Vector>& pos, const Pbc& pbc, ReferenceValuePack& myder, const bool& squared ) const {
+ return calc( pos, pbc, myder, squared );
}
-double SingleDomainRMSD::calc( const std::vector<Vector>& pos, const Pbc& pbc, const std::vector<Value*>& vals, const std::vector<double>& arg, const bool& squared ){
+double SingleDomainRMSD::calc( const std::vector<Vector>& pos, const Pbc& pbc, const std::vector<Value*>& vals, const std::vector<double>& arg,
+ ReferenceValuePack& myder, const bool& squared ) const {
plumed_dbg_assert( vals.size()==0 && pos.size()==getNumberOfAtoms() && arg.size()==0 );
- return calc( pos, pbc, squared );
+ return calc( pos, pbc, myder, squared );
}
}
diff --git a/src/reference/SingleDomainRMSD.h b/src/reference/SingleDomainRMSD.h
index 1f7678c09a3146cfe079188cf439247cd9b17d35..bd906cbd2fce78f78f8390287bf4ee65b129dadc 100644
--- a/src/reference/SingleDomainRMSD.h
+++ b/src/reference/SingleDomainRMSD.h
@@ -36,10 +36,10 @@ public:
/// Set the reference structure
virtual void setReferenceAtoms( const std::vector<Vector>& conf, const std::vector<double>& align_in, const std::vector<double>& displace_in );
/// Calculate
- double calc( const std::vector<Vector>& pos, const Pbc& pbc, const std::vector<Value*>& vals, const std::vector<double>& arg, const bool& squared );
- double calculate( const std::vector<Vector>& pos, const Pbc& pbc, const bool& squared );
+ double calc( const std::vector<Vector>& pos, const Pbc& pbc, const std::vector<Value*>& vals, const std::vector<double>& arg, ReferenceValuePack& myder, const bool& squared ) const;
+ double calculate( const std::vector<Vector>& pos, const Pbc& pbc, ReferenceValuePack& myder, const bool& squared ) const;
/// Calculate the distance using the input position
- virtual double calc( const std::vector<Vector>& pos, const Pbc& pbc, const bool& squared )=0;
+ virtual double calc( const std::vector<Vector>& pos, const Pbc& pbc, ReferenceValuePack& myder, const bool& squared ) const=0;
/// This sets upper and lower bounds on distances to be used in DRMSD (here it does nothing)
virtual void setBoundsOnDistances( bool dopbc, double lbound=0.0, double ubound=std::numeric_limits<double>::max( ) ){};
};
diff --git a/src/secondarystructure/SecondaryStructureRMSD.cpp b/src/secondarystructure/SecondaryStructureRMSD.cpp
index 3eaa972e34cf7495484aa93a63db2682a90df4ce..bc4264d8245f3cc194ce1a2abdd33ecbb18e788d 100644
--- a/src/secondarystructure/SecondaryStructureRMSD.cpp
+++ b/src/secondarystructure/SecondaryStructureRMSD.cpp
@@ -161,7 +161,6 @@ void SecondaryStructureRMSD::setSecondaryStructure( std::vector<Vector>& structu
}
if( references.size()==0 ){
- pos.resize( structure.size() );
finishTaskListUpdate();
readVesselKeywords();
@@ -181,7 +180,7 @@ void SecondaryStructureRMSD::setSecondaryStructure( std::vector<Vector>& structu
std::vector<double> align( structure.size(), 1.0 ), displace( structure.size(), 1.0 );
references[nn]->setBoundsOnDistances( true , bondlength ); // We always use pbc
references[nn]->setReferenceAtoms( structure, align, displace );
- references[nn]->setNumberOfAtoms( structure.size() );
+// references[nn]->setNumberOfAtoms( structure.size() );
}
void SecondaryStructureRMSD::prepare(){
@@ -205,8 +204,9 @@ void SecondaryStructureRMSD::calculate(){
runAllTasks();
}
-void SecondaryStructureRMSD::performTask( const unsigned& task_index, const unsigned& current, vesselbase::MultiValue& myvals ){
+void SecondaryStructureRMSD::performTask( const unsigned& task_index, const unsigned& current, MultiValue& myvals ) const {
// Retrieve the positions
+ std::vector<Vector> pos( references[0]->getNumberOfAtoms() );
for(unsigned i=0;i<pos.size();++i) pos[i]=ActionAtomistic::getPosition( getAtomIndex(current,i) );
// This does strands cutoff
@@ -225,42 +225,31 @@ void SecondaryStructureRMSD::performTask( const unsigned& task_index, const unsi
for(unsigned i=15;i<30;++i){
pos[i]+=( origin_new - origin_old );
}
- }
+ }
+ // Create a holder for the derivatives
+ ReferenceValuePack mypack( 0, pos.size(), myvals ); mypack.setValIndex( 1 );
+ for(unsigned i=0;i<pos.size();++i) mypack.setAtomIndex( i, getAtomIndex(current,i) );
// And now calculate the RMSD
double r,nr; const Pbc& pbc=getPbc();
- closest=0; r=references[0]->calculate( pos, pbc, false );
+ unsigned closest=0; r=references[0]->calculate( pos, pbc, mypack, false );
for(unsigned i=1;i<references.size();++i){
- nr=references[i]->calculate( pos, pbc, false );
+ mypack.setValIndex( i+1 );
+ nr=references[i]->calculate( pos, pbc, mypack, false );
if( nr<r ){ closest=i; r=nr; }
}
// Transfer everything to the value
myvals.setValue( 0, 1.0 ); myvals.setValue( 1, r );
- for(unsigned i=0;i<colvar_atoms[current].size();++i){
- unsigned thisatom=3*getAtomIndex(current,i);
- Vector ader=references[closest]->getAtomDerivative(i);
- myvals.addDerivative( 1, thisatom, ader[0] ); thisatom++;
- myvals.addDerivative( 1, thisatom, ader[1] ); thisatom++;
- myvals.addDerivative( 1, thisatom, ader[2] );
- }
- Tensor virial;
- if( !references[closest]->getVirial( virial ) ){
- virial.zero();
- for(unsigned i=0;i<colvar_atoms[current].size();++i){
- virial+=(-1.0*Tensor( pos[i], references[closest]->getAtomDerivative(i) ));
- }
+ if( closest>0 ) mypack.moveDerivatives( closest+1, 1 );
+
+ if( !mypack.virialWasSet() ){
+ Tensor vir; vir.zero();
+ for(unsigned i=0;i<colvar_atoms[current].size();++i){
+ vir+=(-1.0*Tensor( pos[i], mypack.getAtomDerivative(i) ));
+ }
+ mypack.setValIndex(1); mypack.addBoxDerivatives( vir );
}
- unsigned outnat=3*getNumberOfAtoms();
- myvals.addDerivative( 1, outnat, virial(0,0) ); outnat++;
- myvals.addDerivative( 1, outnat, virial(0,1) ); outnat++;
- myvals.addDerivative( 1, outnat, virial(0,2) ); outnat++;
- myvals.addDerivative( 1, outnat, virial(1,0) ); outnat++;
- myvals.addDerivative( 1, outnat, virial(1,1) ); outnat++;
- myvals.addDerivative( 1, outnat, virial(1,2) ); outnat++;
- myvals.addDerivative( 1, outnat, virial(2,0) ); outnat++;
- myvals.addDerivative( 1, outnat, virial(2,1) ); outnat++;
- myvals.addDerivative( 1, outnat, virial(2,2) );
return;
}
diff --git a/src/secondarystructure/SecondaryStructureRMSD.h b/src/secondarystructure/SecondaryStructureRMSD.h
index 141bedba76e5d309e966f272c2de08e4cb54019f..39f441d3a3d12a3e749f10cce32dacfe440a8b0a 100644
--- a/src/secondarystructure/SecondaryStructureRMSD.h
+++ b/src/secondarystructure/SecondaryStructureRMSD.h
@@ -41,8 +41,6 @@ class SecondaryStructureRMSD :
public vesselbase::ActionWithVessel
{
private:
-/// Tempory integer to say which refernce configuration is the closest
- unsigned closest;
/// The type of rmsd we are calculating
std::string alignType;
/// List of all the atoms we require
@@ -60,10 +58,9 @@ private:
unsigned align_atom_1, align_atom_2;
bool verbose_output;
/// Tempory variables for getting positions of atoms and applying forces
- std::vector<Vector> pos;
std::vector<double> forcesToApply;
/// Get the index of an atom
- unsigned getAtomIndex( const unsigned& current, const unsigned& iatom );
+ unsigned getAtomIndex( const unsigned& current, const unsigned& iatom ) const ;
protected:
/// Get the atoms in the backbone
void readBackboneAtoms( const std::string& backnames, std::vector<unsigned>& chain_lengths );
@@ -79,15 +76,22 @@ public:
virtual ~SecondaryStructureRMSD();
unsigned getNumberOfFunctionsInAction();
unsigned getNumberOfDerivatives();
+ unsigned getNumberOfQuantities();
void turnOnDerivatives();
void prepare();
void finishTaskListUpdate();
void calculate();
- void performTask( const unsigned& , const unsigned& , vesselbase::MultiValue& );
+ void performTask( const unsigned& , const unsigned& , MultiValue& ) const ;
void apply();
bool isPeriodic(){ return false; }
};
+inline
+unsigned SecondaryStructureRMSD::getNumberOfQuantities(){
+ return 1 + references.size();
+}
+
+
inline
unsigned SecondaryStructureRMSD::getNumberOfFunctionsInAction(){
return colvar_atoms.size();
@@ -99,7 +103,7 @@ unsigned SecondaryStructureRMSD::getNumberOfDerivatives(){
}
inline
-unsigned SecondaryStructureRMSD::getAtomIndex( const unsigned& current, const unsigned& iatom ){
+unsigned SecondaryStructureRMSD::getAtomIndex( const unsigned& current, const unsigned& iatom ) const {
return all_atoms.linkIndex( colvar_atoms[current][iatom] );
}
diff --git a/src/tools/LinkCells.cpp b/src/tools/LinkCells.cpp
index 35d3339fe38acc1415078f166e37d24a626d4c7b..849f09016a7f3c039ff99bf932fc303b83bf0f0c 100644
--- a/src/tools/LinkCells.cpp
+++ b/src/tools/LinkCells.cpp
@@ -99,7 +99,7 @@ void LinkCells::buildCellLists( const std::vector<Vector>& pos, const std::vecto
#define LINKC_MAX(n) ((n<3)? 1 : 2)
#define LINKC_PBC(n,num) ((n<0)? num-1 : n%num )
-void LinkCells::retrieveNeighboringAtoms( const Vector& pos, unsigned& natomsper, std::vector<unsigned>& atoms ){
+void LinkCells::retrieveNeighboringAtoms( const Vector& pos, unsigned& natomsper, std::vector<unsigned>& atoms ) const {
plumed_assert( natomsper==1 ); // This is really a bug. If you are trying to reuse this ask GAT for help
std::vector<unsigned> celn( findMyCell( pos ) );
diff --git a/src/tools/LinkCells.h b/src/tools/LinkCells.h
index 4c0a2af7095ffbf076558273fbf2bd823490874c..91563482c3bb521621f89b4c516811ed47ac4caf 100644
--- a/src/tools/LinkCells.h
+++ b/src/tools/LinkCells.h
@@ -68,7 +68,7 @@ public:
/// Build the link cell lists
void buildCellLists( const std::vector<Vector>& pos, const std::vector<unsigned>& indices, const Pbc& pbc );
/// Find a list of relevant atoms
- void retrieveNeighboringAtoms( const Vector& pos, unsigned& natomsper, std::vector<unsigned>& atoms );
+ void retrieveNeighboringAtoms( const Vector& pos, unsigned& natomsper, std::vector<unsigned>& atoms ) const ;
};
inline
diff --git a/src/vesselbase/MultiValue.cpp b/src/tools/MultiValue.cpp
similarity index 99%
rename from src/vesselbase/MultiValue.cpp
rename to src/tools/MultiValue.cpp
index 3a7823e3fbf1f1a2e7bbbd859e5f292ec86be746..9bf60a3e8ff83da05a6f169de69d991c976c3919 100644
--- a/src/vesselbase/MultiValue.cpp
+++ b/src/tools/MultiValue.cpp
@@ -22,7 +22,6 @@
#include "MultiValue.h"
namespace PLMD{
-namespace vesselbase{
MultiValue::MultiValue( const unsigned& nvals, const unsigned& nder ):
values(nvals),
@@ -98,4 +97,3 @@ void MultiValue::quotientRule( const unsigned& nder, const unsigned& dder, const
}
}
-}
diff --git a/src/vesselbase/MultiValue.h b/src/tools/MultiValue.h
similarity index 94%
rename from src/vesselbase/MultiValue.h
rename to src/tools/MultiValue.h
index 57203c3231e61cf1987ac0e1c0d3ab99c91f6fab..05021198f28277f9ac33e5ddc70910f89bca9222 100644
--- a/src/vesselbase/MultiValue.h
+++ b/src/tools/MultiValue.h
@@ -19,16 +19,15 @@
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_vesselbase_MultiValue_h
-#define __PLUMED_vesselbase_MultiValue_h
+#ifndef __PLUMED_tools_MultiValue_h
+#define __PLUMED_tools_MultiValue_h
#include <vector>
-#include "tools/Exception.h"
-#include "tools/Matrix.h"
-#include "tools/DynamicList.h"
+#include "Exception.h"
+#include "Matrix.h"
+#include "DynamicList.h"
namespace PLMD{
-namespace vesselbase{
class MultiValue {
private:
@@ -68,7 +67,7 @@ public:
void sortActiveList();
void updateDynamicList();
///
- unsigned getNumberActive();
+ unsigned getNumberActive() const ;
///
unsigned getActiveIndex( const unsigned& ) const ;
/// Transfer derivatives to buffer
@@ -149,7 +148,7 @@ void MultiValue::sortActiveList(){
}
inline
-unsigned MultiValue::getNumberActive(){
+unsigned MultiValue::getNumberActive() const {
return hasDerivatives.getNumberActive();
}
@@ -164,6 +163,5 @@ void MultiValue::updateDynamicList(){
hasDerivatives.updateActiveMembers();
}
-}
}
#endif
diff --git a/src/vesselbase/ActionWithVessel.cpp b/src/vesselbase/ActionWithVessel.cpp
index 9f1639406cdddd7946c4e5d1c5c556b0a47d1951..48854ca6e3fb91865670de85acf88d5702d77b63 100644
--- a/src/vesselbase/ActionWithVessel.cpp
+++ b/src/vesselbase/ActionWithVessel.cpp
@@ -327,7 +327,7 @@ void ActionWithVessel::runAllTasks(){
finishComputations( buffer );
}
-void ActionWithVessel::transformBridgedDerivatives( const unsigned& current, MultiValue& invals, MultiValue& outvals ){
+void ActionWithVessel::transformBridgedDerivatives( const unsigned& current, MultiValue& invals, MultiValue& outvals ) const {
plumed_error();
}
diff --git a/src/vesselbase/ActionWithVessel.h b/src/vesselbase/ActionWithVessel.h
index f8ea4c7987d243a1f6cc2f448470a74deda3f2a2..c5796e539fc26d437bbb9e25937cfcc6d3b5d7c7 100644
--- a/src/vesselbase/ActionWithVessel.h
+++ b/src/vesselbase/ActionWithVessel.h
@@ -26,7 +26,7 @@
#include "core/ActionAtomistic.h"
#include "tools/Exception.h"
#include "tools/DynamicList.h"
-#include "MultiValue.h"
+#include "tools/MultiValue.h"
#include <vector>
namespace PLMD{
@@ -168,9 +168,9 @@ public:
/// Get the code for the ii th task in the list
unsigned getTaskCode( const unsigned& ii ) const ;
/// Calculate one of the functions in the distribution
- virtual void performTask( const unsigned& , const unsigned& , MultiValue& )=0;
+ virtual void performTask( const unsigned& , const unsigned& , MultiValue& ) const=0;
/// Do the task if we have a bridge
- virtual void transformBridgedDerivatives( const unsigned& current, MultiValue& invals, MultiValue& outvals );
+ virtual void transformBridgedDerivatives( const unsigned& current, MultiValue& invals, MultiValue& outvals ) const;
/// Ensure that data required in other vessels is stored
StoreDataVessel* buildDataStashes( const bool& allow_wcutoff, const double& wtol );
/// Apply forces from bridge vessel - this is rarely used - currently only in ActionVolume