diff --git a/src/multicolvar/ActionVolume.cpp b/src/multicolvar/ActionVolume.cpp
index 7d12425aca7aec1f50db829f193afa2643f0be14..05e26304afbe6175d6998740e5348a5bf5cfea73 100644
--- a/src/multicolvar/ActionVolume.cpp
+++ b/src/multicolvar/ActionVolume.cpp
@@ -86,18 +86,20 @@ bool ActionVolume::performTask( const unsigned& j ){
Vector catom_pos=mycolv->retrieveCentralAtomPos( derivativesOfFractionalCoordinates() );
- double weight; Vector wdf;
+ double weight; Vector wdf;
weight=calculateNumberInside( catom_pos, bead, wdf );
if( not_in ){ weight = 1.0 - weight; wdf *= -1.; }
if( mycolv->isDensity() ){
- setElementValue( 0, weight ); setElementValue( 1, 1.0 );
+ setElementValue( 0, weight ); setElementValue( 1, 1.0 );
for(unsigned i=0;i<mycolv->getNAtoms();++i){
- addElementDerivative( mycolv->getOutputDerivativeIndex(mycolv->current,3*i+0), mycolv->getCentralAtomDerivative(i, 0, wdf ) );
- addElementDerivative( mycolv->getOutputDerivativeIndex(mycolv->current,3*i+1), mycolv->getCentralAtomDerivative(i, 1, wdf ) );
- addElementDerivative( mycolv->getOutputDerivativeIndex(mycolv->current,3*i+2), mycolv->getCentralAtomDerivative(i, 2, wdf ) );
+ unsigned nx=3*linkIndex( i, mycolv->colvar_atoms[mycolv->current], mycolv->all_atoms );
+ addElementDerivative( nx+0, mycolv->getCentralAtomDerivative(i, 0, wdf ) );
+ addElementDerivative( nx+1, mycolv->getCentralAtomDerivative(i, 1, wdf ) );
+ addElementDerivative( nx+2, mycolv->getCentralAtomDerivative(i, 2, wdf ) );
}
} else {
+ // Copy derivatives of the colvar and the value of the colvar
double colv=mycolv->getElementValue(0); setElementValue( 0, colv );
for(unsigned i=0;i<mycolv->nderivatives;++i){
addElementDerivative( mycolv->getOutputDerivativeIndex(mycolv->current,i), mycolv->getElementDerivative(i) );
@@ -106,13 +108,20 @@ bool ActionVolume::performTask( const unsigned& j ){
double ww=mycolv->getElementValue(1);
setElementValue( 1, ww*weight );
- for(unsigned i=0;i<mycolv->getNAtoms();++i){
- addElementDerivative( nderivatives+mycolv->getOutputDerivativeIndex(mycolv->current,3*i+0),
- ww*mycolv->getCentralAtomDerivative(i, 0, wdf ) + weight*mycolv->getElementDerivative(nderivatives+3*i+0) );
- addElementDerivative( nderivatives+mycolv->getOutputDerivativeIndex(mycolv->current,3*i+1),
- ww*mycolv->getCentralAtomDerivative(i, 1, wdf ) + weight*mycolv->getElementDerivative(nderivatives+3*i+1) );
- addElementDerivative( nderivatives+mycolv->getOutputDerivativeIndex(mycolv->current,3*i+2),
- ww*mycolv->getCentralAtomDerivative(i, 2, wdf ) + weight*mycolv->getElementDerivative(nderivatives+3*i+2) );
+ // Add derivatives of weight if we have a weight
+ if( mycolv->weightHasDerivatives ){
+ for(unsigned i=0;i<mycolv->nderivatives;++i){
+ addElementDerivative( nderivatives+mycolv->getOutputDerivativeIndex(mycolv->current,i), weight*mycolv->getElementDerivative(nderivatives+i) );
+ }
+ }
+
+ // Add derivatives of central atoms
+ for(unsigned i=0;i<mycolv->atomsWithCatomDer.getNumberActive();++i){
+ unsigned n=mycolv->atomsWithCatomDer[i];
+ unsigned nx=nderivatives+3*linkIndex( n, mycolv->colvar_atoms[mycolv->current], mycolv->all_atoms);
+ addElementDerivative( nx+0, ww*mycolv->getCentralAtomDerivative(i, 0, wdf ) );
+ addElementDerivative( nx+1, ww*mycolv->getCentralAtomDerivative(i, 1, wdf ) );
+ addElementDerivative( nx+2, ww*mycolv->getCentralAtomDerivative(i, 2, wdf ) );
}
}
diff --git a/src/multicolvar/MultiColvar.cpp b/src/multicolvar/MultiColvar.cpp
index 0a22d60375570501d638c1d9f6f94bf6b8b3a1b2..b73be9095ecd920b19dc98f0dd6578ca04646f3b 100644
--- a/src/multicolvar/MultiColvar.cpp
+++ b/src/multicolvar/MultiColvar.cpp
@@ -103,7 +103,11 @@ void MultiColvar::readAtoms( int& natoms ){
all_atoms.activateAll(); colvar_list.activateAll();
for(unsigned i=0;i<colvar_atoms.size();++i) colvar_atoms[i].activateAll();
ActionAtomistic::requestAtoms( all_atoms.retrieveActiveList() );
- central_derivs.resize( getNumberOfAtoms() );
+ // Set up stuff for central atoms
+ for(unsigned i=0;i<colvar_atoms[0].fullSize();++i) atomsWithCatomDer.addIndexToList( i );
+ atomsWithCatomDer.deactivateAll();
+ central_derivs.resize( colvar_atoms[0].fullSize() );
+ for(unsigned i=0;i<central_derivs.size();++i) central_derivs[i].zero();
}
void MultiColvar::readBackboneAtoms( const std::vector<std::string>& backnames, std::vector<unsigned>& chain_lengths ){
@@ -339,7 +343,7 @@ void MultiColvar::prepare(){
}
if(updatetime){
for(unsigned i=0;i<colvar_list.getNumberActive();++i) activateLinks( colvar_atoms[i], all_atoms );
- all_atoms.updateActiveMembers(); central_derivs.resize( getNumberOfAtoms() );
+ all_atoms.updateActiveMembers();
ActionAtomistic::requestAtoms( all_atoms.retrieveActiveList() ); resizeFunctions();
}
}
@@ -391,13 +395,24 @@ bool MultiColvar::performTask( const unsigned& j ){
Vector MultiColvar::retrieveCentralAtomPos( const bool& frac ){
centralAtomDerivativesAreInFractional=frac;
ibox=getPbc().getInvBox().transpose();
- for(unsigned i=0;i<central_derivs.size();++i) central_derivs[i].zero();
- if(frac) return getPbc().realToScaled( getCentralAtom() );
- return getCentralAtom();
+
+ // Prepare to retrieve central atom
+ for(unsigned i=0;i<atomsWithCatomDer.getNumberActive();++i) central_derivs[ atomsWithCatomDer[i] ].zero();
+ atomsWithCatomDer.deactivateAll();
+
+ // Retrieve the central atom position
+ Vector catom_pos;
+ if(frac) catom_pos=getPbc().realToScaled( getCentralAtom() );
+ else catom_pos=getCentralAtom();
+
+ // Complete the update of the active member list
+ atomsWithCatomDer.updateActiveMembers();
+ return catom_pos;
}
void MultiColvar::addCentralAtomDerivatives( const unsigned& iatom, const Tensor& der ){
plumed_dbg_assert( iatom<colvar_atoms[current].getNumberActive() );
+ atomsWithCatomDer.activate(iatom);
if( centralAtomDerivativesAreInFractional ) central_derivs[iatom] += matmul( ibox, der );
else central_derivs[iatom] += der;
}
@@ -474,6 +489,21 @@ void MultiColvar::apply(){
}
}
+void MultiColvar::buildDerivativeIndexArrays( std::vector< DynamicList<unsigned> >& active_der ) const {
+ active_der.resize( colvar_atoms.size() );
+ for(unsigned i=0;i<colvar_atoms.size();++i){
+ for(unsigned j=0;j<colvar_atoms[i].fullSize();++j){
+ active_der[i].addIndexToList( 3*colvar_atoms[i][j] + 0 );
+ active_der[i].addIndexToList( 3*colvar_atoms[i][j] + 1 );
+ active_der[i].addIndexToList( 3*colvar_atoms[i][j] + 2 );
+ }
+ unsigned virs=3*getNumberOfAtoms();
+ for(unsigned j=0;j<9;++j){
+ active_der[i].addIndexToList( virs ); virs++;
+ }
+ }
+}
+
StoreCentralAtomsVessel* MultiColvar::getCentralAtoms(){
// Look to see if vectors have already been created
StoreCentralAtomsVessel* mycatoms;
@@ -488,6 +518,6 @@ StoreCentralAtomsVessel* MultiColvar::getCentralAtoms(){
addVessel(sv); resizeFunctions(); // This makes sure resizing of vessels is done
return sv;
}
-
+
}
}
diff --git a/src/multicolvar/MultiColvar.h b/src/multicolvar/MultiColvar.h
index 9a67ef7caddbea7e8ed78542b8867dd0b4db28ef..1c3b04e31f3f9ab78ab21682a4335f0ad26dcbad 100644
--- a/src/multicolvar/MultiColvar.h
+++ b/src/multicolvar/MultiColvar.h
@@ -70,6 +70,7 @@ private:
std::vector<Vector> pos;
Tensor ibox;
bool centralAtomDerivativesAreInFractional;
+ DynamicList<unsigned> atomsWithCatomDer;
std::vector<Tensor> central_derivs;
/// Read in the various GROUP keywords
void readGroupsKeyword( int& natoms );
@@ -88,6 +89,8 @@ protected:
void readBackboneAtoms( const std::vector<std::string>& backnames, std::vector<unsigned>& chain_lengths );
/// 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 );
+/// Build lists of indexes of the derivatives from the colvar atoms arrays
+ void buildDerivativeIndexArrays( std::vector< DynamicList<unsigned> >& active_der ) const ;
/// Get the separation between a pair of vectors
Vector getSeparation( const Vector& vec1, const Vector& vec2 ) const ;
/// Find out if it is time to do neighbor list update
diff --git a/src/multicolvar/StoreCentralAtomsVessel.cpp b/src/multicolvar/StoreCentralAtomsVessel.cpp
index 0581240e187cc9429efae31b5fe1998700c4c7d1..631974c5f8509d6b48b8b4e8817dc401025374b7 100644
--- a/src/multicolvar/StoreCentralAtomsVessel.cpp
+++ b/src/multicolvar/StoreCentralAtomsVessel.cpp
@@ -33,6 +33,8 @@ wasforced(false)
{
mycolv=dynamic_cast<MultiColvar*>( getAction() );
plumed_assert( mycolv );
+
+ mycolv->buildDerivativeIndexArrays( active_der );
}
void StoreCentralAtomsVessel::resize(){
@@ -40,7 +42,7 @@ void StoreCentralAtomsVessel::resize(){
unsigned bsize=0; start.resize( nfunc +1 );
for(unsigned i=0;i<nfunc;++i){
start[i] = bsize;
- bsize += 3*( 1 + 3*mycolv->colvar_atoms[i].getNumberActive() );
+ bsize += 3*( 1 + mycolv->getNumberOfDerivatives() );
}
start[nfunc]=bsize;
resizeBuffer( bsize );
@@ -55,25 +57,36 @@ void StoreCentralAtomsVessel::prepare(){
bool StoreCentralAtomsVessel::calculate(){
Vector catom_pos=mycolv->retrieveCentralAtomPos( false );
- unsigned ibuf=start[mycolv->current];
- for(unsigned i=0;i<3;++i){
- addToBufferElement( ibuf, catom_pos[i] ); ibuf++;
- for(unsigned j=0;j<mycolv->getNAtoms();++j){
- for(unsigned k=0;k<3;++k){
- addToBufferElement( ibuf, mycolv->central_derivs[j](i,k) ); ibuf++;
- }
+ // Store the value
+ unsigned ibuf=start[mycolv->current], nspace = 1 + mycolv->getNumberOfDerivatives();
+ for(unsigned i=0;i<3;++i){ addToBufferElement( ibuf, catom_pos[i] ); ibuf+=nspace; }
+ plumed_dbg_assert( ibuf==start[mycolv->current+1] );
+
+ // Store the derivatives
+ active_der[mycolv->current].deactivateAll();
+ for(unsigned j=0;j<mycolv->atomsWithCatomDer.getNumberActive();++j){
+ unsigned n=mycolv->atomsWithCatomDer[j];
+ ibuf = start[mycolv->current] + 1 + 3*mycolv->colvar_atoms[mycolv->current][n];
+ for(unsigned i=0;i<3;++i){
+ active_der[mycolv->current].activate(3*n+i);
+ for(unsigned k=0;k<3;++k) addToBufferElement( ibuf+k, mycolv->central_derivs[n](i,k) );
+ ibuf += nspace;
}
}
- plumed_dbg_assert( ibuf==start[mycolv->current+1] );
+
return true;
}
+void StoreCentralAtomsVessel::finish(){
+ mpi_gatherActiveMembers( comm, active_der );
+}
+
Vector StoreCentralAtomsVessel::getPosition( const unsigned& ivec ) const {
plumed_dbg_assert( ivec<mycolv->getNumberOfFunctionsInAction() );
unsigned pos=start[ivec]; Vector mypos;
for(unsigned i=0;i<3;++i){
mypos[i] = getBufferElement( pos );
- pos+=3*mycolv->colvar_atoms[ivec].getNumberActive()+1;
+ pos+=mycolv->getNumberOfDerivatives() + 1;
}
plumed_dbg_assert( pos==start[ivec+1] );
return mypos;
@@ -96,21 +109,17 @@ bool StoreCentralAtomsVessel::applyForce(std::vector<double>& ff){
void StoreCentralAtomsVessel::chainRuleForCentralAtom( const unsigned& iatom, const unsigned& iderno, const Vector& df, vesselbase::ActionWithVessel* act) const {
plumed_dbg_assert( iatom<mycolv->getNumberOfFunctionsInAction() );
- unsigned nder=3*mycolv->colvar_atoms[iatom].getNumberActive();
- unsigned nder2=mycolv->getNumberOfDerivatives();
- for(unsigned ider=0;ider<nder;++ider){
- unsigned ibuf=start[iatom] + 1 + ider; double thisder=0.0;
+ unsigned nder=mycolv->getNumberOfDerivatives();
+ for(unsigned ider=0;ider<active_der[iatom].getNumberActive();++ider){
+ unsigned ibuf=start[iatom] + 1 + active_der[iatom][ider]; double thisder=0.0;
for(unsigned jcomp=0;jcomp<3;++jcomp){
thisder+=df[jcomp]*getBufferElement(ibuf);
ibuf+=(nder+1);
- //unsigned ibuf=start[iatom] + jcomp*(nder+1) + 1 + ider;
- //act->addElementDerivative( jder, df[jcomp]*getBufferElement(ibuf) );
}
- unsigned jder = iderno*nder2 + mycolv->getOutputDerivativeIndex( iatom, ider );
+ unsigned jder = iderno*nder + active_der[iatom][ider];
act->addElementDerivative( jder, thisder );
}
}
-
}
}
diff --git a/src/multicolvar/StoreCentralAtomsVessel.h b/src/multicolvar/StoreCentralAtomsVessel.h
index c7651b913d0d3d6ec15c4523f73b19c75158efbb..3ddc32d2b83ad400f9c5997bffd302327a9ce88f 100644
--- a/src/multicolvar/StoreCentralAtomsVessel.h
+++ b/src/multicolvar/StoreCentralAtomsVessel.h
@@ -22,6 +22,7 @@
#ifndef __PLUMED_multicolvar_StoreCentralAtomsVessel_h
#define __PLUMED_multicolvar_StoreCentralAtomsVessel_h
+#include "tools/DynamicList.h"
#include "vesselbase/Vessel.h"
#include "MultiColvar.h"
@@ -35,6 +36,7 @@ private:
MultiColvar* mycolv;
std::vector<unsigned> start;
bool wasforced;
+ std::vector< DynamicList<unsigned> > active_der;
std::vector<double> forces;
public:
/// Constructor
@@ -45,8 +47,8 @@ public:
void resize();
/// This does nothing
std::string description(){ return ""; }
-/// This does nothing
- void finish(){}
+/// This should mpi gather the active atoms
+ void finish();
/// Add some force to the atoms
void addForces( const std::vector<double>& );
/// Copy forces to local arrays
diff --git a/src/tools/DynamicList.h b/src/tools/DynamicList.h
index a6b674446cc233da963319b67e11b54207308822..a31ddc406ca5f8e47e37c62410d52d24bb16f295 100644
--- a/src/tools/DynamicList.h
+++ b/src/tools/DynamicList.h
@@ -134,21 +134,22 @@ friend void activateLinks( const DynamicList<unsigned>& , DynamicList<U>& );
template <typename U>
friend void mpi_gatherActiveMembers(Communicator& , std::vector< DynamicList<U> >& );
private:
- bool inactive;
/// This is the list of all the relevent members
std::vector<T> all;
/// This tells us what members of all are on/off at any given time
std::vector<unsigned> onoff;
/// This translates a position from all to active
std::vector<int> translator;
+/// The current number of active members
+ unsigned nactive;
/// This is the list of active members
std::vector<T> active;
public:
/// Constructor
- DynamicList():inactive(true){}
+ DynamicList():nactive(0){}
/// An operator that returns the element from the current active list
inline T operator [] (const unsigned& i) const {
- plumed_dbg_assert(!inactive); plumed_dbg_assert( i<active.size() );
+ plumed_dbg_assert( i<nactive );
return active[i];
}
/// An operator that returns the element from the full list (used in neighbour lists)
@@ -156,7 +157,6 @@ public:
plumed_dbg_assert( i<all.size() );
return all[i];
}
- inline bool get_inactive() const { return inactive; }
/// Clear the list
void clear();
/// Return the total number of elements in the list
@@ -178,17 +178,20 @@ public:
/// Get the list of active members
void updateActiveMembers();
/// Retriee the list of active objects
- std::vector<T>& retrieveActiveList();
+ std::vector<T> retrieveActiveList();
};
template <typename T>
-std::vector<T>& DynamicList<T>::retrieveActiveList(){
- return active;
+std::vector<T> DynamicList<T>::retrieveActiveList(){
+ std::vector<T> this_active(nactive);
+ for(unsigned k=0;k<nactive;++k) this_active[k]=active[k];
+ return this_active;
}
template <typename T>
void DynamicList<T>::clear() {
- all.resize(0); translator.resize(0); onoff.resize(0); active.resize(0);
+ all.resize(0); translator.resize(0);
+ onoff.resize(0); active.resize(0);
}
template <typename T>
@@ -198,39 +201,34 @@ unsigned DynamicList<T>::fullSize() const {
template <typename T>
unsigned DynamicList<T>::getNumberActive() const {
- plumed_dbg_assert( !inactive || active.size()==0 );
- return active.size();
+ return nactive;
}
template <typename T>
void DynamicList<T>::addIndexToList( const T & ii ){
- all.push_back(ii); translator.push_back( all.size()-1 ); onoff.push_back(0);
+ all.push_back(ii); translator.push_back( all.size()-1 );
+ onoff.push_back(0); active.push_back(ii);
}
template <typename T>
void DynamicList<T>::deactivate( const unsigned ii ){
plumed_massert(ii<all.size(),"ii is out of bounds");
- onoff[ii]=0; inactive=true;
- for(unsigned i=0;i<onoff.size();++i){
- if(onoff[i]>0){ inactive=false; break; }
- }
+ onoff[ii]=0;
}
template <typename T>
void DynamicList<T>::deactivateAll(){
- inactive=true;
for(unsigned i=0;i<onoff.size();++i) onoff[i]=0;
}
template <typename T>
void DynamicList<T>::activate( const unsigned ii ){
plumed_massert(ii<all.size(),"ii is out of bounds");
- inactive=false; onoff[ii]=1;
+ onoff[ii]=1;
}
template <typename T>
void DynamicList<T>::activateAll(){
- inactive=false;
for(unsigned i=0;i<onoff.size();++i) onoff[i]=1;
updateActiveMembers();
}
@@ -246,17 +244,16 @@ void DynamicList<T>::mpi_gatherActiveMembers(Communicator& comm){
template <typename T>
void DynamicList<T>::updateActiveMembers(){
- unsigned kk=0; active.resize(0);
+ unsigned kk=0;
for(unsigned i=0;i<all.size();++i){
- if( onoff[i]==1 ){ translator[i]=kk; active.push_back( all[i] ); kk++; }
+ if( onoff[i]==1 ){ translator[i]=kk; active[kk]=all[i]; kk++; }
}
- plumed_assert( kk==active.size() );
- if( kk==0 ) inactive=true;
+ nactive=kk;
}
template <typename U>
unsigned linkIndex( const unsigned ii, const DynamicList<unsigned>& l1, const DynamicList<U>& l2 ){
- plumed_dbg_massert(ii<l1.active.size(),"ii is out of bounds");
+ plumed_dbg_massert(ii<l1.nactive,"ii is out of bounds");
unsigned kk; kk=l1.active[ii];
plumed_dbg_massert(kk<l2.all.size(),"the lists are mismatched");
plumed_dbg_massert( l2.onoff[kk]==1, "This index is not currently in the second list" );
@@ -266,7 +263,7 @@ unsigned linkIndex( const unsigned ii, const DynamicList<unsigned>& l1, const Dy
template <typename U>
void activateLinks( const DynamicList<unsigned>& l1, DynamicList<U>& l2 ){
- for(unsigned i=0;i<l1.active.size();++i) l2.activate( l1.active[i] );
+ for(unsigned i=0;i<l1.nactive;++i) l2.activate( l1.active[i] );
}
template <typename U>
diff --git a/src/vesselbase/Vessel.cpp b/src/vesselbase/Vessel.cpp
index ce4265ee935253326e3c8d43a7bd06f0920a4e76..085464154d4e2e62066e0bf724d72de665306b7e 100644
--- a/src/vesselbase/Vessel.cpp
+++ b/src/vesselbase/Vessel.cpp
@@ -58,8 +58,8 @@ Vessel::Vessel( const VesselOptions& da ):
myname(da.myname),
numlab(da.numlab),
action(da.action),
-comm(da.action->comm),
keywords(da.keywords),
+comm(da.action->comm),
finished_read(false),
log((da.action)->log)
{
diff --git a/src/vesselbase/Vessel.h b/src/vesselbase/Vessel.h
index 1f1689983cae0b8bf6127effe4582754ad44f9eb..3902ae68896ca4dd01ef74fb6640c5a1fb400ed0 100644
--- a/src/vesselbase/Vessel.h
+++ b/src/vesselbase/Vessel.h
@@ -85,8 +85,6 @@ private:
const int numlab;
/// The action that this vessel is created within
ActionWithVessel* action;
-/// A copy of the communicator
- Communicator& comm;
/// Something to store the buffer if this is required
std::vector<double> stash;
/// The start of this Vessel's buffer in buffer in the underlying ActionWithVessel
@@ -102,6 +100,8 @@ private:
/// This just checks we have done checkRead
bool finished_read;
protected:
+/// A copy of the communicator
+ Communicator& comm;
/// Return the numerical label
int getNumericalLabel() const ;
/// Report an error