From a1fb75b94a10a132932656e79d60634a44913368 Mon Sep 17 00:00:00 2001
From: Gareth Tribello <gareth.tribello@gmail.com>
Date: Sat, 24 Aug 2013 18:57:24 +0100
Subject: [PATCH] Removed colvar_atoms arrays
This modification ensures that multicolvar uses less memory and makes
it a tiny bit faster
---
src/manyrestraints/Sphere.cpp | 4 +-
src/multicolvar/ActionVolume.cpp | 5 +-
src/multicolvar/ActionVolume.h | 2 +-
src/multicolvar/AlphaBeta.cpp | 8 +-
src/multicolvar/Angles.cpp | 20 +--
src/multicolvar/Bridge.cpp | 19 +-
src/multicolvar/CoordinationNumbers.cpp | 16 +-
src/multicolvar/Density.cpp | 6 +-
src/multicolvar/Distances.cpp | 4 +-
src/multicolvar/MultiColvar.cpp | 163 ++++++++++--------
src/multicolvar/MultiColvar.h | 5 +-
src/multicolvar/MultiColvarBase.cpp | 82 ++++++---
src/multicolvar/MultiColvarBase.h | 33 ++--
src/multicolvar/MultiColvarFunction.cpp | 17 +-
src/multicolvar/MultiColvarFunction.h | 30 +---
src/multicolvar/StoreCentralAtomsVessel.cpp | 2 +
.../SecondaryStructureRMSD.cpp | 2 +-
.../SecondaryStructureRMSD.h | 2 +-
src/vesselbase/ActionWithVessel.cpp | 4 +-
src/vesselbase/ActionWithVessel.h | 4 +-
src/vesselbase/BridgeVessel.cpp | 2 +-
src/vesselbase/StoreDataVessel.cpp | 9 +-
src/vesselbase/StoreDataVessel.h | 2 +-
23 files changed, 255 insertions(+), 186 deletions(-)
diff --git a/src/manyrestraints/Sphere.cpp b/src/manyrestraints/Sphere.cpp
index 1d7333455..a094467a9 100644
--- a/src/manyrestraints/Sphere.cpp
+++ b/src/manyrestraints/Sphere.cpp
@@ -39,7 +39,7 @@ public:
static void registerKeywords( Keywords& keys );
Sphere( const ActionOptions& );
bool isPeriodic(){ return false; }
- void performTask( const unsigned& j );
+ void performTask();
void calculate();
};
@@ -90,7 +90,7 @@ void Sphere::calculate(){
runAllTasks();
}
-void Sphere::performTask( const unsigned& j ){
+void Sphere::performTask(){
Vector distance;
if(!nopbc){
diff --git a/src/multicolvar/ActionVolume.cpp b/src/multicolvar/ActionVolume.cpp
index 4f08834e2..4c0e68f80 100644
--- a/src/multicolvar/ActionVolume.cpp
+++ b/src/multicolvar/ActionVolume.cpp
@@ -113,7 +113,8 @@ void ActionVolume::prepare(){
if( updateFreq>0 && (getStep()-lastUpdate)>=updateFreq ){
if( !mycolv->isDensity() ){
mycolv->taskList.activateAll();
- for(unsigned i=0;i<mycolv->taskList.getNumberActive();++i) mycolv->colvar_atoms[i].activateAll();
+ // GAT broken neighbor list
+ // for(unsigned i=0;i<mycolv->taskList.getNumberActive();++i) mycolv->colvar_atoms[i].activateAll();
mycolv->unlockContributors(); mycolv->resizeDynamicArrays();
plumed_dbg_assert( mycolv->getNumberOfVessels()==0 );
} else {
@@ -131,7 +132,7 @@ void ActionVolume::resizeLocalArrays(){
activeAtoms.deactivateAll();
}
-void ActionVolume::performTask( const unsigned& j ){
+void ActionVolume::performTask(){
activeAtoms.deactivateAll(); // Currently no atoms are active so deactivate them all
Vector catom_pos=mycolv->retrieveCentralAtomPos();
diff --git a/src/multicolvar/ActionVolume.h b/src/multicolvar/ActionVolume.h
index d1a8ee86c..5947d0ffb 100644
--- a/src/multicolvar/ActionVolume.h
+++ b/src/multicolvar/ActionVolume.h
@@ -100,7 +100,7 @@ public:
/// Do jobs required before tasks are undertaken
void doJobsRequiredBeforeTaskList();
/// This calculates all the vessels and is called from within a bridge vessel
- void performTask(const unsigned& i );
+ void performTask();
/// Routines that have to be defined so as not to have problems with virtual methods
void deactivate_task();
void calculate(){}
diff --git a/src/multicolvar/AlphaBeta.cpp b/src/multicolvar/AlphaBeta.cpp
index dc2858895..3f9092adc 100644
--- a/src/multicolvar/AlphaBeta.cpp
+++ b/src/multicolvar/AlphaBeta.cpp
@@ -79,7 +79,7 @@ private:
public:
static void registerKeywords( Keywords& keys );
AlphaBeta(const ActionOptions&);
- virtual double compute( const unsigned& j );
+ virtual double compute();
bool isPeriodic(){ return false; }
Vector getCentralAtom();
};
@@ -127,7 +127,7 @@ PLUMED_MULTICOLVAR_INIT(ao)
checkRead();
}
-double AlphaBeta::compute( const unsigned& j ){
+double AlphaBeta::compute(){
Vector d0,d1,d2;
d0=getSeparation(getPosition(1),getPosition(0));
d1=getSeparation(getPosition(2),getPosition(1));
@@ -136,8 +136,8 @@ double AlphaBeta::compute( const unsigned& j ){
Vector dd0,dd1,dd2;
PLMD::Torsion t;
double value = t.compute(d0,d1,d2,dd0,dd1,dd2);
- double svalue = -0.5*sin(value-target[current]);
- double cvalue = 1.+cos(value-target[current]);
+ double svalue = -0.5*sin(value-target[lindex]);
+ double cvalue = 1.+cos(value-target[lindex]);
dd0 *= svalue;
dd1 *= svalue;
diff --git a/src/multicolvar/Angles.cpp b/src/multicolvar/Angles.cpp
index 06e8daa06..f14963b82 100644
--- a/src/multicolvar/Angles.cpp
+++ b/src/multicolvar/Angles.cpp
@@ -93,7 +93,7 @@ public:
static void registerKeywords( Keywords& keys );
Angles(const ActionOptions&);
// active methods:
- virtual double compute( const unsigned& j );
+ virtual double compute();
/// Returns the number of coordinates of the field
void calculateWeight();
bool isPeriodic(){ return false; }
@@ -159,8 +159,8 @@ use_sf(false)
}
void Angles::calculateWeight(){
- dij=getSeparation( getPosition(1), getPosition(2) );
- dik=getSeparation( getPosition(1), getPosition(0) );
+ dij=getSeparation( getPosition(0), getPosition(2) );
+ dik=getSeparation( getPosition(0), getPosition(1) );
if(!use_sf){ setWeight(1.0); return; }
double w1, w2, dw1, dw2, wtot;
@@ -170,19 +170,19 @@ void Angles::calculateWeight(){
setWeight( wtot );
if( wtot<getTolerance() ) return;
- addAtomsDerivativeOfWeight( 0, dw2*dik );
- addAtomsDerivativeOfWeight( 1, -dw1*dij - dw2*dik );
+ addAtomsDerivativeOfWeight( 1, dw2*dik );
+ addAtomsDerivativeOfWeight( 0, -dw1*dij - dw2*dik );
addAtomsDerivativeOfWeight( 2, dw1*dij );
addBoxDerivativesOfWeight( (-dw1)*Tensor(dij,dij) + (-dw2)*Tensor(dik,dik) );
}
-double Angles::compute( const unsigned& j ){
+double Angles::compute(){
Vector ddij,ddik; PLMD::Angle a;
double angle=a.compute(dij,dik,ddij,ddik);
// And finish the calculation
- addAtomsDerivatives( 0, ddik );
- addAtomsDerivatives( 1, - ddik - ddij );
+ addAtomsDerivatives( 1, ddik );
+ addAtomsDerivatives( 0, - ddik - ddij );
addAtomsDerivatives( 2, ddij );
addBoxDerivatives( -(Tensor(dij,ddij)+Tensor(dik,ddik)) );
@@ -190,8 +190,8 @@ double Angles::compute( const unsigned& j ){
}
Vector Angles::getCentralAtom(){
- addCentralAtomDerivatives( 1, Tensor::identity() );
- return getPosition(1);
+ addCentralAtomDerivatives( 0, Tensor::identity() );
+ return getPosition(0);
}
}
diff --git a/src/multicolvar/Bridge.cpp b/src/multicolvar/Bridge.cpp
index 2732a1455..27304d835 100644
--- a/src/multicolvar/Bridge.cpp
+++ b/src/multicolvar/Bridge.cpp
@@ -67,11 +67,10 @@ public:
static void registerKeywords( Keywords& keys );
Bridge(const ActionOptions&);
// active methods:
- virtual double compute( const unsigned& j );
+ virtual double compute();
/// Returns the number of coordinates of the field
void calculateWeight();
bool isPeriodic(){ return false; }
- unsigned getNumberOfAtomsInCentralAtomDerivatives(){ return 1; }
Vector getCentralAtom();
};
@@ -128,30 +127,30 @@ PLUMED_MULTICOLVAR_INIT(ao)
}
void Bridge::calculateWeight(){
- Vector dij=getSeparation( getPosition(1), getPosition(0) );
+ Vector dij=getSeparation( getPosition(0), getPosition(1) );
double dw, w=sf1.calculate( dij.modulo(), dw );
setWeight( w );
if( w<getTolerance() ) return;
- addAtomsDerivativeOfWeight( 1, -dw*dij );
- addAtomsDerivativeOfWeight( 0, dw*dij );
+ addAtomsDerivativeOfWeight( 0, -dw*dij );
+ addAtomsDerivativeOfWeight( 1, dw*dij );
addBoxDerivativesOfWeight( (-dw)*Tensor(dij,dij) );
}
-double Bridge::compute( const unsigned& j ){
- Vector dik=getSeparation( getPosition(1), getPosition(2) );
+double Bridge::compute(){
+ Vector dik=getSeparation( getPosition(0), getPosition(2) );
double dw, w=sf2.calculate( dik.modulo(), dw );
// And finish the calculation
- addAtomsDerivatives( 1, -dw*dik );
+ addAtomsDerivatives( 0, -dw*dik );
addAtomsDerivatives( 2, dw*dik );
addBoxDerivatives( (-dw)*Tensor(dik,dik) );
return w;
}
Vector Bridge::getCentralAtom(){
- addCentralAtomDerivatives( 1, Tensor::identity() );
- return getPosition(1);
+ addCentralAtomDerivatives( 0, Tensor::identity() );
+ return getPosition(0);
}
}
diff --git a/src/multicolvar/CoordinationNumbers.cpp b/src/multicolvar/CoordinationNumbers.cpp
index 7038374b6..7869299a8 100644
--- a/src/multicolvar/CoordinationNumbers.cpp
+++ b/src/multicolvar/CoordinationNumbers.cpp
@@ -70,7 +70,7 @@ public:
static void registerKeywords( Keywords& keys );
CoordinationNumbers(const ActionOptions&);
// active methods:
- virtual double compute( const unsigned& j );
+ virtual double compute();
Vector getCentralAtom();
/// Returns the number of coordinates of the field
bool isPeriodic(){ return false; }
@@ -110,19 +110,13 @@ PLUMED_MULTICOLVAR_INIT(ao)
log.printf(" coordination of central atom and those within %s\n",( switchingFunction.description() ).c_str() );
// Read in the atoms
- int natoms; readAtoms( natoms );
+ int natoms=2; readAtoms( natoms );
// And setup the ActionWithVessel
readVesselKeywords();
-
- // Create the groups for the neighbor list
- std::vector<AtomNumber> ga_lista, gb_lista; AtomNumber aa;
- aa.setIndex(0); ga_lista.push_back(aa);
- for(unsigned i=1;i<natoms;++i){ aa.setIndex(i); gb_lista.push_back(aa); }
- // And check everything has been read in correctly
checkRead();
}
-double CoordinationNumbers::compute( const unsigned& j ){
+double CoordinationNumbers::compute(){
double value=0, dfunc; Vector distance;
// Calculate the coordination number
@@ -130,8 +124,8 @@ double CoordinationNumbers::compute( const unsigned& j ){
for(unsigned i=1;i<getNAtoms();++i){
distance=getSeparation( getPosition(0), getPosition(i) );
sw = switchingFunction.calculate( distance.modulo(), dfunc );
- if( sw>=getTolerance() ){ // nl_cut<0 ){
- value += sw; // switchingFunction.calculate( distance.modulo(), dfunc );
+ if( sw>=getTolerance() ){
+ value += sw;
addAtomsDerivatives( 0, (-dfunc)*distance );
addAtomsDerivatives( i, (dfunc)*distance );
addBoxDerivatives( (-dfunc)*Tensor(distance,distance) );
diff --git a/src/multicolvar/Density.cpp b/src/multicolvar/Density.cpp
index 0794f598c..365da74c6 100644
--- a/src/multicolvar/Density.cpp
+++ b/src/multicolvar/Density.cpp
@@ -54,7 +54,7 @@ public:
static void registerKeywords( Keywords& keys );
Density(const ActionOptions&);
// active methods:
- virtual double compute( const unsigned& j );
+ virtual double compute();
Vector getCentralAtom();
/// Returns the number of coordinates of the field
bool isPeriodic(){ return false; }
@@ -71,13 +71,13 @@ void Density::registerKeywords( Keywords& keys ){
Density::Density(const ActionOptions&ao):
PLUMED_MULTICOLVAR_INIT(ao)
{
- int nat; readAtoms( nat );
+ int nat=1; readAtoms( nat );
readVesselKeywords();
// And check everything has been read in correctly
checkRead();
}
-double Density::compute( const unsigned& j ){
+double Density::compute(){
return 1.0;
}
diff --git a/src/multicolvar/Distances.cpp b/src/multicolvar/Distances.cpp
index bbe066e2b..970c06b27 100644
--- a/src/multicolvar/Distances.cpp
+++ b/src/multicolvar/Distances.cpp
@@ -80,7 +80,7 @@ public:
static void registerKeywords( Keywords& keys );
Distances(const ActionOptions&);
// active methods:
- virtual double compute( const unsigned& j );
+ virtual double compute();
/// Returns the number of coordinates of the field
bool isPeriodic(){ return false; }
Vector getCentralAtom();
@@ -111,7 +111,7 @@ PLUMED_MULTICOLVAR_INIT(ao)
checkRead();
}
-double Distances::compute( const unsigned& j ){
+double Distances::compute(){
Vector distance;
distance=getSeparation( getPosition(0), getPosition(1) );
const double value=distance.modulo();
diff --git a/src/multicolvar/MultiColvar.cpp b/src/multicolvar/MultiColvar.cpp
index 0f9ad26ea..204a5cfbb 100644
--- a/src/multicolvar/MultiColvar.cpp
+++ b/src/multicolvar/MultiColvar.cpp
@@ -59,15 +59,6 @@ verbose_output(false)
parseFlag("VERBOSE",verbose_output);
}
-void MultiColvar::addColvar( const std::vector<unsigned>& newatoms ){
- if( verbose_output ){
- log.printf(" Colvar %d is calculated from atoms : ", colvar_atoms.size()+1);
- for(unsigned i=0;i<newatoms.size();++i) log.printf("%d ",all_atoms(newatoms[i]).serial() );
- log.printf("\n");
- }
- MultiColvarBase::addColvar( newatoms );
-}
-
void MultiColvar::readAtoms( int& natoms ){
if( keywords.exists("ATOMS") ) readAtomsLikeKeyword( "ATOMS", natoms );
if( keywords.exists("GROUP") ) readGroupsKeyword( natoms );
@@ -79,9 +70,10 @@ void MultiColvar::readAtoms( int& natoms ){
}
void MultiColvar::readAtomsLikeKeyword( const std::string & key, int& natoms ){
+ plumed_assert( !usespecies );
if( all_atoms.fullSize()>0 ) return;
- std::vector<AtomNumber> t; std::vector<unsigned> newlist;
+ std::vector<AtomNumber> t;
for(int i=1;;++i ){
parseAtomList(key, i, t );
if( t.empty() ) break;
@@ -92,21 +84,29 @@ void MultiColvar::readAtomsLikeKeyword( const std::string & key, int& natoms ){
log.printf("\n");
}
- if( i==1 && natoms<0 ) natoms=t.size();
+ if( i==1 && natoms<0 ){ natoms=t.size(); ablocks.resize(natoms); }
+ else if( i==1 ) ablocks.resize(natoms);
if( t.size()!=natoms ){
- std::string ss; Tools::convert(i,ss);
- error(key + ss + " keyword has the wrong number of atoms");
+ std::string ss; Tools::convert(i,ss);
+ error(key + ss + " keyword has the wrong number of atoms");
}
for(unsigned j=0;j<natoms;++j){
- newlist.push_back( natoms*(i-1)+j );
- all_atoms.addIndexToList( t[j] );
+ ablocks[j].push_back( natoms*(i-1)+j ); all_atoms.addIndexToList( t[j] );
+ }
+ t.resize(0);
+ }
+ if( all_atoms.fullSize()>0 ){
+ current_atoms.resize( natoms ); nblock=ablocks[0].size();
+ for(unsigned i=0;i<nblock;++i){
+ unsigned cvcode=0, tmpc=1;
+ for(unsigned j=0;j<natoms;++j){ cvcode += i*tmpc; tmpc *= nblock; }
+ taskList.addIndexToList( cvcode );
}
- t.resize(0); addColvar( newlist );
- newlist.clear();
}
}
void MultiColvar::readGroupsKeyword( int& natoms ){
+ plumed_assert( !usespecies );
if( all_atoms.fullSize()>0 ) return;
if( natoms==2 ){
@@ -123,22 +123,20 @@ void MultiColvar::readGroupsKeyword( int& natoms ){
std::vector<AtomNumber> t;
parseAtomList("GROUP",t);
if( !t.empty() ){
+ ablocks.resize( natoms ); current_atoms.resize( natoms );
for(unsigned i=0;i<t.size();++i) all_atoms.addIndexToList( t[i] );
- std::vector<unsigned> newlist;
if(natoms==2){
+ nblock=t.size(); for(unsigned i=0;i<2;++i) ablocks[i].resize(nblock);
+ for(unsigned i=0;i<t.size();++i){ ablocks[0][i]=i; ablocks[1][i]=i; }
for(unsigned i=1;i<t.size();++i){
- for(unsigned j=0;j<i;++j){
- newlist.push_back(i); newlist.push_back(j); addColvar( newlist ); newlist.clear();
- }
+ for(unsigned j=0;j<i;++j) taskList.addIndexToList( i*nblock + j );
}
} else if(natoms==3){
+ nblock=t.size(); for(unsigned i=0;i<3;++i) ablocks[i].resize(nblock);
+ for(unsigned i=0;i<t.size();++i){ ablocks[0][i]=i; ablocks[1][i]=i; ablocks[2][i]=i; }
for(unsigned i=2;i<t.size();++i){
for(unsigned j=1;j<i;++j){
- for(unsigned k=0;k<j;++k){
- newlist.push_back(i); newlist.push_back(j);
- newlist.push_back(k); addColvar( newlist );
- newlist.clear();
- }
+ for(unsigned k=0;k<j;++k) taskList.addIndexToList( i*nblock*nblock + j*nblock + k );
}
}
}
@@ -160,17 +158,26 @@ void MultiColvar::readGroupsKeyword( int& natoms ){
void MultiColvar::readTwoGroups( const std::string& key1, const std::string& key2 ){
plumed_assert( all_atoms.fullSize()==0 );
+ ablocks.resize( 2 ); current_atoms.resize( 2 );
std::vector<AtomNumber> t1, t2; std::vector<unsigned> newlist;
parseAtomList(key1,t1); parseAtomList(key2,t2);
if( t1.empty() ) error("missing atom specification " + key1);
if ( t2.empty() ) error("missing atom specification " + key2);
- for(unsigned i=0;i<t1.size();++i) all_atoms.addIndexToList( t1[i] );
- for(unsigned i=0;i<t2.size();++i) all_atoms.addIndexToList( t2[i] );
+
+ if( t1.size()>t2.size() ) nblock = t1.size();
+ else nblock=t2.size();
+
+ ablocks[0].resize( t1.size() );
for(unsigned i=0;i<t1.size();++i){
- for(unsigned j=0;j<t2.size();++j){
- newlist.push_back(i); newlist.push_back( t1.size() + j ); addColvar( newlist ); newlist.clear();
- }
+ all_atoms.addIndexToList( t1[i] ); ablocks[0][i]=i;
+ }
+ ablocks[1].resize( t2.size() );
+ for(unsigned i=0;i<t2.size();++i){
+ all_atoms.addIndexToList( t2[i] ); ablocks[1][i]=t1.size() + i;
+ }
+ for(unsigned i=0;i<t1.size();++i){
+ for(unsigned j=0;j<t2.size();++j) taskList.addIndexToList( i*nblock + j );
}
if( !verbose_output ){
log.printf(" constructing colvars from two groups containing %d and %d atoms respectively\n",t1.size(),t2.size() );
@@ -185,25 +192,32 @@ void MultiColvar::readTwoGroups( const std::string& key1, const std::string& key
void MultiColvar::readThreeGroups( const std::string& key1, const std::string& key2, const std::string& key3, const bool& allow2 ){
plumed_assert( all_atoms.fullSize()==0 );
+ ablocks.resize( 3 ); current_atoms.resize( 3 );
std::vector<AtomNumber> t1, t2, t3; std::vector<unsigned> newlist;
parseAtomList(key1,t1); parseAtomList(key2,t2);
if( t1.empty() ) error("missing atom specification " + key1);
if( t2.empty() ) error("missing atom specification " + key2);
- for(unsigned i=0;i<t1.size();++i) all_atoms.addIndexToList( t1[i] );
- for(unsigned i=0;i<t2.size();++i) all_atoms.addIndexToList( t2[i] );
+ ablocks[0].resize( t1.size() );
+ for(unsigned i=0;i<t1.size();++i){
+ all_atoms.addIndexToList( t1[i] ); ablocks[0][i]=i;
+ }
+ ablocks[1].resize( t2.size() );
+ for(unsigned i=0;i<t2.size();++i){
+ all_atoms.addIndexToList( t2[i] ); ablocks[1][i] = t1.size() + i;
+ }
parseAtomList(key3,t3);
if( t3.empty() && !allow2 ){
error("missing atom specification " + key3);
} else if( t3.empty() ){
+ if( t2.size()>t1.size() ) nblock=t2.size();
+ else nblock=t1.size();
+
+ ablocks[2].resize( t2.size() );
+ for(unsigned i=0;i<t2.size();++i) ablocks[2][i]=t1.size() + i;
for(unsigned i=0;i<t1.size();++i){
for(unsigned j=1;j<t2.size();++j){
- for(unsigned k=0;k<j;++k){
- newlist.push_back( t1.size() + j );
- newlist.push_back(i);
- newlist.push_back( t1.size() + k );
- addColvar( newlist ); newlist.clear();
- }
+ for(unsigned k=0;k<j;++k) taskList.addIndexToList( nblock*nblock*i + nblock*j + k );
}
}
if( !verbose_output ){
@@ -216,15 +230,17 @@ void MultiColvar::readThreeGroups( const std::string& key1, const std::string& k
log.printf("\n");
}
} else {
- for(unsigned i=0;i<t3.size();++i) all_atoms.addIndexToList( t3[i] );
+ if( t2.size()>t1.size() ) nblock=t2.size();
+ else nblock=t1.size();
+ if( t3.size()>nblock ) nblock=t3.size();
+
+ ablocks[2].resize( t3.size() );
+ for(unsigned i=0;i<t3.size();++i){
+ all_atoms.addIndexToList( t3[i] ); ablocks[2][i] = t1.size() + t2.size() + i;
+ }
for(unsigned i=0;i<t1.size();++i){
for(unsigned j=0;j<t2.size();++j){
- for(unsigned k=0;k<t3.size();++k){
- newlist.push_back( t1.size() + j );
- newlist.push_back(i);
- newlist.push_back( t1.size() + t2.size() + k );
- addColvar( newlist ); newlist.clear();
- }
+ for(unsigned k=0;k<t3.size();++k) taskList.addIndexToList( nblock*nblock*i + nblock*j + k );
}
}
if( !verbose_output ){
@@ -243,22 +259,25 @@ void MultiColvar::readThreeGroups( const std::string& key1, const std::string& k
}
void MultiColvar::readSpeciesKeyword( int& natoms ){
+ plumed_assert( usespecies );
if( all_atoms.fullSize()>0 ) return ;
+ ablocks.resize( natoms-1 );
std::vector<AtomNumber> t;
parseAtomList("SPECIES",t);
if( !t.empty() ){
- natoms=t.size();
for(unsigned i=0;i<t.size();++i) all_atoms.addIndexToList( t[i] );
- std::vector<unsigned> newlist;
if( keywords.exists("SPECIESA") && keywords.exists("SPECIESB") ){
- for(unsigned i=0;i<t.size();++i){
- newlist.push_back(i);
- for(unsigned j=0;j<t.size();++j){
- if(i!=j) newlist.push_back(j);
- }
- addColvar( newlist ); newlist.clear();
- }
+ plumed_assert( natoms==2 ); current_atoms.resize( t.size() );
+ for(unsigned i=0;i<t.size();++i) taskList.addIndexToList(i);
+ ablocks[0].resize( t.size() ); for(unsigned i=0;i<t.size();++i) ablocks[0][i]=i;
+ // for(unsigned i=0;i<t.size();++i){
+ // newlist.push_back(i);
+ // for(unsigned j=0;j<t.size();++j){
+ // if(i!=j) newlist.push_back(j);
+ // }
+ // addColvar( newlist ); newlist.clear();
+ // }
if( !verbose_output ){
log.printf(" generating colvars from %d atoms of a particular type\n",t.size() );
log.printf(" atoms involved : ");
@@ -268,9 +287,9 @@ void MultiColvar::readSpeciesKeyword( int& natoms ){
} else if( !( keywords.exists("SPECIESA") && keywords.exists("SPECIESB") ) ){
std::vector<unsigned> newlist; verbose_output=false; // Make sure we don't do verbose output
log.printf(" involving atoms : ");
+ current_atoms.resize(1);
for(unsigned i=0;i<t.size();++i){
- newlist.push_back(i); addColvar( newlist ); newlist.clear();
- log.printf(" %d",t[i].serial() );
+ taskList.addIndexToList(i); log.printf(" %d",t[i].serial() );
}
log.printf("\n");
} else {
@@ -282,17 +301,16 @@ void MultiColvar::readSpeciesKeyword( int& natoms ){
if( !t1.empty() ){
parseAtomList("SPECIESB",t2);
if ( t2.empty() ) error("SPECIESB keyword defines no atoms or is missing. Use either SPECIESA and SPECIESB or just SPECIES");
- natoms=1+t2.size();
- for(unsigned i=0;i<t1.size();++i) all_atoms.addIndexToList( t1[i] );
- for(unsigned i=0;i<t2.size();++i) all_atoms.addIndexToList( t2[i] );
- std::vector<unsigned> newlist;
- for(unsigned i=0;i<t1.size();++i){
- newlist.push_back(i);
- for(unsigned j=0;j<t2.size();++j){
- if( t1[i]!=t2[j] ) newlist.push_back( t1.size() + j );
- }
- addColvar( newlist ); newlist.clear();
- }
+ current_atoms.resize( 1+ t2.size() );
+ for(unsigned i=0;i<t1.size();++i){ all_atoms.addIndexToList( t1[i] ); taskList.addIndexToList(i); }
+ ablocks[0].resize( t2.size() ); for(unsigned i=0;i<t2.size();++i){ all_atoms.addIndexToList( t2[i] ); ablocks[0][i]=t1.size() + i; }
+ // for(unsigned i=0;i<t1.size();++i){
+ // newlist.push_back(i);
+ // for(unsigned j=0;j<t2.size();++j){
+ // if( t1[i]!=t2[j] ) newlist.push_back( t1.size() + j );
+ // }
+ // addColvar( newlist ); newlist.clear();
+ // }
if( !verbose_output ){
log.printf(" generating colvars from a group of %d central atoms and %d other atoms\n",t1.size(), t2.size() );
log.printf(" central atoms are : ");
@@ -309,10 +327,11 @@ void MultiColvar::readSpeciesKeyword( int& natoms ){
void MultiColvar::resizeDynamicArrays(){
all_atoms.deactivateAll();
for(unsigned i=0;i<taskList.getNumberActive();++i){
- unsigned n=taskList[i];
- for(unsigned j=0;j<colvar_atoms[n].getNumberActive();++j){
- all_atoms.activate( colvar_atoms[n][j] );
- }
+ lindex=taskList.linkIndex(i);
+ current=taskList[i];
+ bool check=setupCurrentAtomList();
+ plumed_dbg_assert( check );
+ for(unsigned j=0;j<natomsper;++j) all_atoms.activate( current_atoms[j] );
}
all_atoms.updateActiveMembers();
// Request the atoms
diff --git a/src/multicolvar/MultiColvar.h b/src/multicolvar/MultiColvar.h
index 37e66714f..4455837d0 100644
--- a/src/multicolvar/MultiColvar.h
+++ b/src/multicolvar/MultiColvar.h
@@ -92,8 +92,9 @@ public:
inline
unsigned MultiColvar::getAtomIndex( const unsigned& iatom ) const {
- plumed_dbg_assert( iatom<colvar_atoms[current].getNumberActive() );
- return all_atoms.linkIndex( colvar_atoms[current][iatom] );
+ plumed_dbg_assert( iatom<natomsper );
+ return current_atoms[iatom];
+// return all_atoms.linkIndex( colvar_atoms[current][iatom] );
}
inline
diff --git a/src/multicolvar/MultiColvarBase.cpp b/src/multicolvar/MultiColvarBase.cpp
index 17e591704..e88f245ca 100644
--- a/src/multicolvar/MultiColvarBase.cpp
+++ b/src/multicolvar/MultiColvarBase.cpp
@@ -48,24 +48,19 @@ ActionWithValue(ao),
ActionWithVessel(ao),
usepbc(false),
updateFreq(0),
-lastUpdate(0)
+lastUpdate(0),
+usespecies(false)
{
if( keywords.exists("NOPBC") ){
bool nopbc=!usepbc; parseFlag("NOPBC",nopbc);
usepbc=!nopbc;
}
+ if( keywords.exists("SPECIES") ) usespecies=true;
if( keywords.exists("NL_STRIDE") ) parse("NL_STRIDE",updateFreq);
if(updateFreq>0) log.printf(" Updating contributors every %d steps.\n",updateFreq);
else log.printf(" Updating contributors every step.\n");
}
-void MultiColvarBase::addColvar( const std::vector<unsigned>& newatoms ){
- DynamicList<unsigned> newlist; newlist.setupMPICommunication( comm );
- for(unsigned i=0;i<newatoms.size();++i) newlist.addIndexToList( newatoms[i] );
- taskList.addIndexToList( colvar_atoms.size() );
- colvar_atoms.push_back( newlist );
-}
-
void MultiColvarBase::copyAtomListToFunction( MultiColvarBase* myfunction ){
for(unsigned i=0;i<all_atoms.fullSize();++i) myfunction->all_atoms.addIndexToList( all_atoms(i) );
}
@@ -83,7 +78,22 @@ void MultiColvarBase::copyActiveAtomsToFunction( MultiColvarBase* myfunction ){
void MultiColvarBase::setupMultiColvarBase(){
// Activate everything
taskList.activateAll();
- for(unsigned i=0;i<colvar_atoms.size();++i) colvar_atoms[i].activateAll();
+ // Setup decoder array
+ if( !usespecies ){
+ decoder.resize( ablocks.size() ); unsigned code=1;
+ for(unsigned i=0;i<ablocks.size();++i){ decoder[ablocks.size()-1-i]=code; code *= nblock; }
+ } else if( ablocks.size()>0 ) {
+ plumed_assert( ablocks.size()==1 );
+ // Setup coordination sphere
+ csphere_atoms.resize( taskList.fullSize() );
+ for(unsigned i=0;i<taskList.fullSize();++i){
+ for(unsigned j=0;j<ablocks[0].size();++j){
+ if( taskList(i)!=ablocks[0][j] ) csphere_atoms[i].addIndexToList( ablocks[0][j] );
+ }
+ csphere_atoms[i].activateAll();
+ }
+ }
+
// Resize stuff in derived classes
resizeDynamicArrays();
// Resize local arrays
@@ -100,12 +110,14 @@ void MultiColvarBase::prepare(){
bool updatetime=false;
if( contributorsAreUnlocked ){
taskList.mpi_gatherActiveMembers( comm );
- mpi_gatherActiveMembers( comm, colvar_atoms );
+ if( usespecies ) mpi_gatherActiveMembers( comm, csphere_atoms );
lockContributors(); updatetime=true;
}
if( updateFreq>0 && (getStep()-lastUpdate)>=updateFreq ){
taskList.activateAll();
- for(unsigned i=0;i<taskList.getNumberActive();++i) colvar_atoms[i].activateAll();
+ if(usespecies){
+ for(unsigned i=0;i<taskList.getNumberActive();++i) csphere_atoms[i].activateAll();
+ }
unlockContributors(); updatetime=true; lastUpdate=getStep();
}
if(updatetime){
@@ -130,17 +142,43 @@ void MultiColvarBase::resizeLocalArrays(){
forcesToApply.resize( getNumberOfDerivatives() );
}
-void MultiColvarBase::performTask( const unsigned& j ){
+bool MultiColvarBase::setupCurrentAtomList(){
+ if( usespecies ){
+ natomsper=1;
+ current_atoms[0]=all_atoms.linkIndex( current );
+ for(unsigned j=0;j<ablocks.size();++j){
+ for(unsigned i=0;i<csphere_atoms[current].getNumberActive();++i){
+ current_atoms[natomsper]=all_atoms.linkIndex( csphere_atoms[current][i] );
+ natomsper++;
+ }
+ }
+ if( natomsper==1 ) return isDensity();
+ } else {
+ natomsper=current_atoms.size();
+ unsigned scode = current;
+ for(unsigned i=0;i<ablocks.size();++i){
+ unsigned ind=std::floor( scode / decoder[i] );
+ current_atoms[i]=all_atoms.linkIndex( ablocks[i][ind] );
+ scode -= ind*decoder[i];
+ }
+ }
+ return true;
+}
+
+void MultiColvarBase::performTask(){
// Currently no atoms have derivatives so deactivate those that are active
atoms_with_derivatives.deactivateAll();
// Currently no central atoms have derivatives so deactive them all
atomsWithCatomDer.deactivateAll();
+ // Retrieve the atom list
+ if( !setupCurrentAtomList() ) return;
// Do nothing if there are no active atoms in the colvar
- if( colvar_atoms[current].getNumberActive()==0 ){
- setElementValue(1,0.0);
- return;
- }
+// if( colvar_atoms[current].getNumberActive()==0 ){
+// setElementValue(1,0.0);
+// return;
+// } Add retrieve atoms here
+
// Do a quick check on the size of this contribution
calculateWeight();
if( getElementValue(1)<getTolerance() ){
@@ -149,14 +187,14 @@ void MultiColvarBase::performTask( const unsigned& j ){
}
// Compute everything
- double vv=doCalculation( j );
+ double vv=doCalculation();
// Set the value of this element in ActionWithVessel
setElementValue( 0, vv );
return;
}
-double MultiColvarBase::doCalculation( const unsigned& j ){
- double val=compute(j); updateActiveAtoms();
+double MultiColvarBase::doCalculation(){
+ double val=compute(); updateActiveAtoms();
return val;
}
@@ -198,10 +236,10 @@ Vector MultiColvarBase::getSeparation( const Vector& vec1, const Vector& vec2 )
}
unsigned MultiColvarBase::getInternalIndex( const AtomNumber& iatom ) const {
+ plumed_massert( usespecies && ablocks.size()==1, "This should only be used to interogate atom centered multicolvars");
unsigned katom; bool found=false;
- for(unsigned i=0;i<colvar_atoms.size();++i){
- unsigned jatom = colvar_atoms[i][0];
- if( all_atoms[ all_atoms.linkIndex(jatom) ]==iatom ){
+ for(unsigned i=0;i<ablocks[0].size();++i){
+ if( all_atoms[ all_atoms.linkIndex(ablocks[0][i]) ]==iatom ){
katom=i; found=true;
}
}
diff --git a/src/multicolvar/MultiColvarBase.h b/src/multicolvar/MultiColvarBase.h
index 728e8c4f1..648cb8cdf 100644
--- a/src/multicolvar/MultiColvarBase.h
+++ b/src/multicolvar/MultiColvarBase.h
@@ -56,15 +56,25 @@ private:
DynamicList<unsigned> atomsWithCatomDer;
/// The forces we are going to apply to things
std::vector<double> forcesToApply;
+/// Neighbor lists for coordination numbers
+ std::vector<DynamicList<unsigned> > csphere_atoms;
/// This resizes the local arrays after neighbor list updates and during initialization
void resizeLocalArrays();
protected:
/// A dynamic list containing those atoms with derivatives
DynamicList<unsigned> atoms_with_derivatives;
-/// The lists of the atoms involved in each of the individual colvars
- std::vector< DynamicList<unsigned> > colvar_atoms;
-/// Add a colvar to the set of colvars we are calculating (in practise just a list of atoms)
- void addColvar( const std::vector<unsigned>& newatoms );
+/// Using the species keyword to read in atoms
+ bool usespecies;
+/// Number of atoms in each block
+ unsigned nblock;
+/// This is used when turning cvcodes into atom numbers
+ std::vector<unsigned> decoder;
+/// Blocks of atom numbers
+ std::vector< std::vector<unsigned> > ablocks;
+/// Number of atoms in the cv - set at start of calculation
+ unsigned natomsper;
+/// Vector containing the indices of the current atoms
+ std::vector<unsigned> current_atoms;
/// Finish setting up the multicolvar base
void setupMultiColvarBase();
/// Request the atoms from action atomistic
@@ -95,6 +105,8 @@ protected:
void getCentralAtomIndexList( const unsigned& ntotal, const unsigned& jstore, const unsigned& maxder, std::vector<unsigned>& indices ) const ;
/// Calculate and store getElementValue(uder)/getElementValue(vder) and its derivatives in getElementValue(iout)
void quotientRule( const unsigned& uder, const unsigned& vder, const unsigned& iout );
+/// This sets up the list of atoms that are involved in this colvar
+ bool setupCurrentAtomList();
public:
MultiColvarBase(const ActionOptions&);
~MultiColvarBase(){}
@@ -103,13 +115,13 @@ public:
void prepare();
virtual void resizeDynamicArrays()=0;
/// Perform one of the tasks
- void performTask( const unsigned& j );
+ void performTask();
/// And a virtual function which actually computes the colvar
- virtual double doCalculation( const unsigned& j );
+ virtual double doCalculation();
/// Update the atoms that have derivatives
virtual void updateActiveAtoms()=0;
/// This is replaced once we have a function to calculate the cv
- virtual double compute( const unsigned& j )=0;
+ virtual double compute()=0;
/// These replace the functions in ActionWithVessel to make the code faster
void mergeDerivatives( const unsigned& ider, const double& df );
void clearDerivativesAfterTask( const unsigned& ider );
@@ -152,15 +164,16 @@ unsigned MultiColvarBase::getNumberOfDerivatives(){
inline
void MultiColvarBase::removeAtomRequest( const unsigned& i, const double& weight ){
+ plumed_dbg_assert( usespecies );
if( !contributorsAreUnlocked ) return;
plumed_dbg_assert( weight<getTolerance() );
- if( weight<getNLTolerance() ) colvar_atoms[current].deactivate( i );
+// GAT neighbor list
+ if( weight<getNLTolerance() ) csphere_atoms[current].deactivate( i );
}
inline
void MultiColvarBase::deactivate_task(){
if( !contributorsAreUnlocked ) return; // Deactivating tasks only possible during neighbor list update
- colvar_atoms[current].deactivateAll(); // Deactivate all atom requests for this colvar
ActionWithVessel::deactivate_task(); // Deactivate the colvar from the list
}
@@ -176,7 +189,7 @@ unsigned MultiColvarBase::getNumberOfQuantities(){
inline
unsigned MultiColvarBase::getNAtoms() const {
- return colvar_atoms[current].getNumberActive();
+ return natomsper; // colvar_atoms[current].getNumberActive();
}
inline
diff --git a/src/multicolvar/MultiColvarFunction.cpp b/src/multicolvar/MultiColvarFunction.cpp
index 970eda643..614c422e8 100644
--- a/src/multicolvar/MultiColvarFunction.cpp
+++ b/src/multicolvar/MultiColvarFunction.cpp
@@ -59,13 +59,24 @@ MultiColvarBase(ao)
log.printf(" calculating function for atoms ");
for(unsigned i=0;i<atoms.size();++i){
log.printf("%d ",atoms[i].serial() );
- subatomlist.push_back( mycolv->getInternalIndex(atoms[i]) );
+ taskList.addIndexToList( mycolv->getInternalIndex(atoms[i]) );
}
log.printf("\n");
} else {
- for(unsigned i=0;i<mycolv->colvar_atoms.size();++i) subatomlist.push_back( i );
+ for(unsigned i=0;i<mycolv->taskList.fullSize();++i) taskList.addIndexToList( i );
}
- }
+ usespecies=true; ablocks.resize(1); ablocks[0].resize( getNumberOfBaseFunctions() );
+ for(unsigned i=0;i<getNumberOfBaseFunctions();++i) ablocks[0][i]=i;
+ current_atoms.resize( 1 + ablocks[0].size() );
+ } else {
+ usespecies=false; nblock=getNumberOfBaseFunctions(); ablocks.resize(2);
+ for(unsigned i=0;i<2;++i) ablocks[i].resize(nblock);
+ for(unsigned i=0;i<getNumberOfBaseFunctions();++i){ ablocks[0][i]=i; ablocks[1][i]=i; }
+ for(unsigned i=1;i<getNumberOfBaseFunctions();++i){
+ for(unsigned j=0;j<i;++j) taskList.addIndexToList( i*nblock + j );
+ }
+ current_atoms.resize( 2 );
+ }
}
void MultiColvarFunction::completeSetup(){
diff --git a/src/multicolvar/MultiColvarFunction.h b/src/multicolvar/MultiColvarFunction.h
index eb57c303f..a4392d13d 100644
--- a/src/multicolvar/MultiColvarFunction.h
+++ b/src/multicolvar/MultiColvarFunction.h
@@ -30,15 +30,11 @@ namespace multicolvar {
class MultiColvarFunction : public MultiColvarBase {
private:
-/// The list of atoms that are involved in this colvar
- std::vector<unsigned> subatomlist;
/// The multicolvar from which we construct these quantities
multicolvar::MultiColvarBase* mycolv;
/// The central atom positions
multicolvar::StoreCentralAtomsVessel* catoms;
protected:
-/// Get the number of colvar functions we are calculating
- unsigned getNumberOfColvarFunctions() const;
/// Get the index of the ith colvar we are using
unsigned getColvarIndex( const unsigned& ) const;
/// Get the number of functions in the multicolvar we are operating on
@@ -83,25 +79,13 @@ public:
inline
unsigned MultiColvarFunction::getNumberOfBaseFunctions() const {
- return mycolv->colvar_atoms.size();
-}
-
-inline
-unsigned MultiColvarFunction::getNumberOfColvarFunctions() const {
- plumed_dbg_assert( subatomlist.size()>0 );
- return subatomlist.size();
-}
-
-inline
-unsigned MultiColvarFunction::getColvarIndex( const unsigned& jindex ) const {
- plumed_dbg_assert( subatomlist.size()>0 && jindex<subatomlist.size() );
- return subatomlist[jindex];
+ return mycolv->taskList.fullSize();
}
inline
unsigned MultiColvarFunction::getAtomIndex( const unsigned& iatom ) const {
- plumed_dbg_assert( iatom<colvar_atoms[current].getNumberActive() );
- return mycolv->getIndexForTask( colvar_atoms[current][iatom] );
+ plumed_dbg_assert( iatom<natomsper );
+ return mycolv->getIndexForTask( current_atoms[iatom] );
}
inline
@@ -111,13 +95,13 @@ MultiColvarBase* MultiColvarFunction::getPntrToMultiColvar(){
inline
Vector MultiColvarFunction::getPositionOfCentralAtom( const unsigned& iatom ) const {
- plumed_dbg_assert( iatom<colvar_atoms[current].getNumberActive() );
+ plumed_dbg_assert( iatom<natomsper );
return catoms->getPosition( getAtomIndex(iatom) );
}
inline
void MultiColvarFunction::addCentralAtomsDerivatives( const unsigned& iatom, const Vector& der ){
- plumed_dbg_assert( iatom<colvar_atoms[current].getNumberActive() );
+ plumed_dbg_assert( iatom<natomsper );
if( usingLowMem() ){
catoms->recompute( getAtomIndex(iatom), 0 );
catoms->addAtomsDerivatives( 0, 0, der, this );
@@ -128,7 +112,7 @@ void MultiColvarFunction::addCentralAtomsDerivatives( const unsigned& iatom, con
inline
void MultiColvarFunction::addCentralAtomsDerivativesOfWeight( const unsigned& iatom, const Vector& der ){
- plumed_dbg_assert( iatom<colvar_atoms[current].getNumberActive() );
+ plumed_dbg_assert( iatom<natomsper );
if( usingLowMem() ){
catoms->recompute( getAtomIndex(iatom), 0 );
catoms->addAtomsDerivatives( 0, 1, der, this );
@@ -144,7 +128,7 @@ void MultiColvarFunction::atomHasDerivative( const unsigned& iatom ){
inline
void MultiColvarFunction::addDerivativeOfCentralAtomPos( const unsigned& iatom, const Tensor& der ){
- plumed_dbg_assert( iatom<colvar_atoms[current].getNumberActive() );
+ plumed_dbg_assert( iatom<natomsper );
if( usingLowMem() ){
catoms->recompute( getAtomIndex(iatom), 0 ); Vector tmpder;
for(unsigned i=0;i<3;++i){
diff --git a/src/multicolvar/StoreCentralAtomsVessel.cpp b/src/multicolvar/StoreCentralAtomsVessel.cpp
index 6577c3104..cda1acafc 100644
--- a/src/multicolvar/StoreCentralAtomsVessel.cpp
+++ b/src/multicolvar/StoreCentralAtomsVessel.cpp
@@ -54,6 +54,8 @@ Vector StoreCentralAtomsVessel::getPosition( const unsigned& iatom ){
void StoreCentralAtomsVessel::performTask( const unsigned& itask ){
mycolv->atomsWithCatomDer.deactivateAll();
+ bool check=mycolv->setupCurrentAtomList();
+ plumed_dbg_assert( !check );
Vector ignore = mycolv->retrieveCentralAtomPos();
}
diff --git a/src/secondarystructure/SecondaryStructureRMSD.cpp b/src/secondarystructure/SecondaryStructureRMSD.cpp
index ed9835bb9..61c32650d 100644
--- a/src/secondarystructure/SecondaryStructureRMSD.cpp
+++ b/src/secondarystructure/SecondaryStructureRMSD.cpp
@@ -217,7 +217,7 @@ void SecondaryStructureRMSD::calculate(){
runAllTasks();
}
-void SecondaryStructureRMSD::performTask( const unsigned& j ){
+void SecondaryStructureRMSD::performTask(){
// Retrieve the positions
for(unsigned i=0;i<pos.size();++i) pos[i]=ActionAtomistic::getPosition( getAtomIndex(i) );
diff --git a/src/secondarystructure/SecondaryStructureRMSD.h b/src/secondarystructure/SecondaryStructureRMSD.h
index 915fdd26e..2a72dc7fd 100644
--- a/src/secondarystructure/SecondaryStructureRMSD.h
+++ b/src/secondarystructure/SecondaryStructureRMSD.h
@@ -90,7 +90,7 @@ public:
unsigned getNumberOfDerivatives();
void prepare();
void calculate();
- void performTask( const unsigned& j );
+ void performTask();
void clearDerivativesAfterTask( const unsigned& );
void apply();
void mergeDerivatives( const unsigned& , const double& );
diff --git a/src/vesselbase/ActionWithVessel.cpp b/src/vesselbase/ActionWithVessel.cpp
index d1d495005..69d4ad0f9 100644
--- a/src/vesselbase/ActionWithVessel.cpp
+++ b/src/vesselbase/ActionWithVessel.cpp
@@ -182,10 +182,12 @@ void ActionWithVessel::runAllTasks(){
doJobsRequiredBeforeTaskList();
for(unsigned i=rank;i<taskList.getNumberActive();i+=stride){
+ // This is the position of the task in the dynamic list
+ lindex=taskList.linkIndex(i);
// Store the task we are currently working on
current=taskList[i];
// Calculate the stuff in the loop for this action
- performTask(i);
+ performTask();
// Weight should be between zero and one
plumed_dbg_assert( thisval[1]>=0 && thisval[1]<=1.0 );
diff --git a/src/vesselbase/ActionWithVessel.h b/src/vesselbase/ActionWithVessel.h
index af098eb4d..fc814c3b3 100644
--- a/src/vesselbase/ActionWithVessel.h
+++ b/src/vesselbase/ActionWithVessel.h
@@ -78,6 +78,8 @@ protected:
bool contributorsAreUnlocked;
/// Does the weight have derivatives
bool weightHasDerivatives;
+/// The position of the current task in the taskList
+ unsigned lindex;
/// The numerical index of the task we are curently performing
unsigned current;
/// This is used for numerical derivatives of bridge variables
@@ -153,7 +155,7 @@ public:
/// Get the index for a particular numbered task
unsigned getIndexForTask( const unsigned& itask ) const ;
/// Calculate one of the functions in the distribution
- virtual void performTask( const unsigned& j )=0;
+ virtual void performTask()=0;
/// Return a pointer to the field
Vessel* getVessel( const std::string& name );
/// Set the derivative of the jth element wrt to a numbered element
diff --git a/src/vesselbase/BridgeVessel.cpp b/src/vesselbase/BridgeVessel.cpp
index c206c15dc..72b018e38 100644
--- a/src/vesselbase/BridgeVessel.cpp
+++ b/src/vesselbase/BridgeVessel.cpp
@@ -57,7 +57,7 @@ void BridgeVessel::prepare(){
}
bool BridgeVessel::calculate(){
- myOutputAction->performTask( getAction()->current );
+ myOutputAction->performTask();
if( myOutputAction->thisval[1]<myOutputAction->getTolerance() ){
myOutputAction->clearAfterTask();
return ( !myOutputAction->contributorsAreUnlocked || myOutputAction->thisval[1]>=myOutputAction->getNLTolerance() );
diff --git a/src/vesselbase/StoreDataVessel.cpp b/src/vesselbase/StoreDataVessel.cpp
index 86e7cb1b1..e41314546 100644
--- a/src/vesselbase/StoreDataVessel.cpp
+++ b/src/vesselbase/StoreDataVessel.cpp
@@ -69,7 +69,8 @@ void StoreDataVessel::recompute( const unsigned& ivec, const unsigned& jstore ){
ActionWithVessel* act = getAction();
// Set the task we want to reperform
- act->current = ivec;
+ act->current = act->taskList(ivec);
+ act->lindex = act->taskList.linkIndex(ivec);
// Reperform the task
performTask( jstore );
// Store the derivatives
@@ -88,6 +89,7 @@ void StoreDataVessel::storeValues( const unsigned& myelem ){
}
void StoreDataVessel::storeDerivativesHighMem( const unsigned& myelem ){
+ plumed_dbg_assert( myelem<getAction()->getNumberOfTasks() );
ActionWithVessel* act=getAction();
getIndexList( act->getNumberOfTasks(), myelem, nspace-1, active_der );
@@ -121,7 +123,7 @@ void StoreDataVessel::storeDerivativesLowMem( const unsigned& jstore ){
}
bool StoreDataVessel::calculate(){
- unsigned myelem = getAction()->current;
+ unsigned myelem = getAction()->lindex;
// Normalize vector if it is required
finishTask( myelem );
@@ -245,7 +247,7 @@ void StoreDataVessel::chainRule( const unsigned& ival, const std::vector<double>
void StoreDataVessel::transformComponents( const unsigned& jstore, const double& weight, double& wdf, const std::vector<double>& dfvec ){
plumed_dbg_assert( dfvec.size()==vecsize );
ActionWithVessel* act = getAction();
- unsigned myelem = act->current;
+ unsigned myelem = act->lindex;
unsigned ibuf = myelem * vecsize * nspace;
for(unsigned icomp=0;icomp<vecsize;++icomp){
@@ -254,6 +256,7 @@ void StoreDataVessel::transformComponents( const unsigned& jstore, const double&
}
if( !act->lowmem ) {
+ plumed_dbg_assert( jstore<getAction()->getNumberOfTasks() );
for(unsigned ider=0;ider<active_der[myelem];++ider){
double comp2=0.0; unsigned ibuf = myelem * vecsize * nspace + 1 + ider;
for(unsigned jcomp=0;jcomp<vecsize;++jcomp){
diff --git a/src/vesselbase/StoreDataVessel.h b/src/vesselbase/StoreDataVessel.h
index cf8c9c74c..96b1aa1e9 100644
--- a/src/vesselbase/StoreDataVessel.h
+++ b/src/vesselbase/StoreDataVessel.h
@@ -132,7 +132,7 @@ bool StoreDataVessel::usingLowMem(){
inline
void StoreDataVessel::performTask( const unsigned& ivec ){
- if( usingLowMem() ) getAction()->performTask( ivec );
+ if( usingLowMem() ) getAction()->performTask();
}
inline
--
GitLab