From 06734ff18d590d117d681dc1e272b568e90e0643 Mon Sep 17 00:00:00 2001
From: Gareth Tribello <gareth.tribello@gmail.com>
Date: Mon, 18 Feb 2013 15:31:24 +0100
Subject: [PATCH] Some small changes to core of MultiColvar and Vessels
I want to reuse stuff from these classes to do some other things and
I needed to make a few small changes to make the new things work.
---
src/multicolvar/MultiColvar.cpp | 23 +++-
src/multicolvar/MultiColvar.h | 11 +-
src/multicolvar/Region.cpp | 3 +-
src/multicolvar/StoreCentralAtomsVessel.cpp | 121 ++++++++++++++++++++
src/multicolvar/StoreCentralAtomsVessel.h | 69 +++++++++++
src/vesselbase/ActionWithVessel.h | 7 +-
src/vesselbase/ShortcutVessel.h | 2 +-
src/vesselbase/Sum.cpp | 6 +-
src/vesselbase/Vessel.cpp | 5 +-
src/vesselbase/Vessel.h | 7 ++
10 files changed, 239 insertions(+), 15 deletions(-)
create mode 100644 src/multicolvar/StoreCentralAtomsVessel.cpp
create mode 100644 src/multicolvar/StoreCentralAtomsVessel.h
diff --git a/src/multicolvar/MultiColvar.cpp b/src/multicolvar/MultiColvar.cpp
index dd5117505..557120462 100644
--- a/src/multicolvar/MultiColvar.cpp
+++ b/src/multicolvar/MultiColvar.cpp
@@ -326,7 +326,7 @@ void MultiColvar::readSpeciesKeyword( int& natoms ){
}
void MultiColvar::prepare(){
- bool updatetime=false;
+ updatetime=false;
if( reduceAtNextStep ){
colvar_list.mpi_gatherActiveMembers( comm );
mpi_gatherActiveMembers( comm, colvar_atoms );
@@ -387,7 +387,8 @@ bool MultiColvar::performTask( const unsigned& j ){
return false;
}
-Vector MultiColvar::retrieveCentralAtomPos(){
+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();
return getPbc().realToScaled( getCentralAtom() );
@@ -470,5 +471,23 @@ void MultiColvar::apply(){
}
}
}
+
+StoreCentralAtomsVessel* MultiColvar::getCentralAtoms(){
+ // Look to see if vectors have already been created
+ StoreCentralAtomsVessel* mycatoms;
+ for(unsigned i=0;i<getNumberOfVessels();++i){
+ mycatoms=dynamic_cast<StoreCentralAtomsVessel*>( getPntrToVessel(i) );
+ if( mycatoms ) return mycatoms;
+ }
+
+ // Create the vessel
+ std::string input; addVessel( "CATOM_STASH", input );
+ resizeFunctions(); // This makes sure resizing of vessels is done
+
+ // And create a pointer to it here
+ mycatoms = dynamic_cast<StoreCentralAtomsVessel*>( getPntrToVessel( getNumberOfVessels() - 1 ) );
+ return mycatoms;
+}
+
}
}
diff --git a/src/multicolvar/MultiColvar.h b/src/multicolvar/MultiColvar.h
index f32741878..473021cb5 100644
--- a/src/multicolvar/MultiColvar.h
+++ b/src/multicolvar/MultiColvar.h
@@ -26,6 +26,7 @@
#include "core/ActionWithValue.h"
#include "tools/DynamicList.h"
#include "vesselbase/ActionWithVessel.h"
+#include "StoreCentralAtomsVessel.h"
#include <vector>
#define PLUMED_MULTICOLVAR_INIT(ao) Action(ao),MultiColvar(ao)
@@ -45,6 +46,7 @@ class MultiColvar :
public vesselbase::ActionWithVessel
{
friend class Region;
+friend class StoreCentralAtomsVessel;
private:
bool usepbc;
bool readatoms;
@@ -77,6 +79,8 @@ private:
/// Update the atoms request
void requestAtoms();
protected:
+/// Are we on an update step
+ bool updatetime;
/// Read in all the keywords that can be used to define atoms
void readAtoms( int& natoms );
/// Read in the atoms that form the backbone of a polymeric chain
@@ -118,7 +122,7 @@ public:
/// Calculate the multicolvar
void calculate();
/// Prepare for the calculation
- void prepare();
+ virtual void prepare();
/// Apply the forces on the values
void apply();
/// Get the number of derivatives for this action
@@ -128,7 +132,7 @@ public:
/// Return the number of derivatives for a given colvar
unsigned getNumberOfDerivatives( const unsigned& j );
/// Retrieve the position of the central atom
- Vector retrieveCentralAtomPos();
+ Vector retrieveCentralAtomPos( const bool& frac );
/// Make sure we calculate the position of the central atom
void useCentralAtom();
/// Merge the derivatives
@@ -159,6 +163,9 @@ public:
double getCharge(unsigned) const ;
/// Get the absolute index of atom iatom
AtomNumber getAbsoluteIndex(unsigned) const ;
+/// Return a pointer to the vessel that stores the positions of
+/// all the central atoms
+ StoreCentralAtomsVessel* getCentralAtoms();
};
inline
diff --git a/src/multicolvar/Region.cpp b/src/multicolvar/Region.cpp
index cada68362..254a3c056 100644
--- a/src/multicolvar/Region.cpp
+++ b/src/multicolvar/Region.cpp
@@ -135,7 +135,6 @@ not_in(false)
if( volname.substr(0,1)=="!" ){ volname=volname.substr(1); not_in=true; }
myvol = (mycolv->plumed).getActionSet().selectWithLabel<ActionVolume*>(volname);
if(!myvol){ error( "in REGION " + volname + " is not a valid volume element"); return; }
- mycolv->centralAtomDerivativesAreInFractional=myvol->derivativesOfFractionalCoordinates();
mycolv->addDependency( myvol );
parse("SIGMA",sigma); myvol->setSigma( sigma );
@@ -152,7 +151,7 @@ std::string Region::function_description(){
}
bool Region::calculate(){
- Vector catom_pos=mycolv->retrieveCentralAtomPos();
+ Vector catom_pos=mycolv->retrieveCentralAtomPos( myvol->derivativesOfFractionalCoordinates() );
double weight; Vector wdf;
weight=myvol->calculateNumberInside( catom_pos, bead, wdf );
diff --git a/src/multicolvar/StoreCentralAtomsVessel.cpp b/src/multicolvar/StoreCentralAtomsVessel.cpp
new file mode 100644
index 000000000..565c9420e
--- /dev/null
+++ b/src/multicolvar/StoreCentralAtomsVessel.cpp
@@ -0,0 +1,121 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2012 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.0.
+
+ 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 "vesselbase/VesselRegister.h"
+#include "vesselbase/ActionWithVessel.h"
+#include "MultiColvar.h"
+#include "StoreCentralAtomsVessel.h"
+
+namespace PLMD {
+namespace multicolvar {
+
+PLUMED_REGISTER_VESSEL(StoreCentralAtomsVessel,"CATOM_STASH")
+
+void StoreCentralAtomsVessel::registerKeywords( Keywords& keys ){
+ Vessel::registerKeywords( keys );
+ plumed_assert( keys.size()==0 );
+}
+
+void StoreCentralAtomsVessel::reserveKeyword( Keywords& keys ){
+ keys.reserve("optional","CATOM_STASH","This documentation shouldn't appear");
+}
+
+StoreCentralAtomsVessel::StoreCentralAtomsVessel( const vesselbase::VesselOptions& da ):
+Vessel(da),
+wasforced(false)
+{
+ mycolv=dynamic_cast<MultiColvar*>( getAction() );
+ plumed_assert( mycolv );
+}
+
+void StoreCentralAtomsVessel::resize(){
+ unsigned nfunc=mycolv->getNumberOfFunctionsInAction();
+ unsigned bsize=0; start.resize( nfunc +1 );
+ for(unsigned i=0;i<nfunc;++i){
+ start[i] = bsize;
+ bsize += 3*( 1 + mycolv->getNumberOfDerivatives(i) );
+ }
+ start[nfunc]=bsize;
+ resizeBuffer( bsize );
+ forces.resize( mycolv->getNumberOfDerivatives() );
+}
+
+void StoreCentralAtomsVessel::prepare(){
+ wasforced=false;
+ forces.assign( forces.size(), 0.0 );
+}
+
+bool StoreCentralAtomsVessel::calculate(){
+ Vector catom_pos=mycolv->retrieveCentralAtomPos( false );
+ Vector wdf; wdf.zero();
+
+ unsigned ibuf=start[mycolv->current];
+ for(unsigned i=0;i<3;++i){
+ addToBufferElement( ibuf, catom_pos[i] ); ibuf++; wdf[i]=1.0;
+ for(unsigned j=0;j<mycolv->getNAtoms();++j){
+ addToBufferElement( ibuf, mycolv->getCentralAtomDerivative( j, i, wdf ) ); ibuf++;
+ }
+ wdf.zero();
+ }
+ plumed_dbg_assert( ibuf==start[mycolv->current+1] );
+ return true;
+}
+
+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+=mycolv->getNumberOfDerivatives(ivec)+1;
+ }
+ plumed_dbg_assert( pos==start[ivec+1] );
+ return mypos;
+}
+
+void StoreCentralAtomsVessel::addForces( const std::vector<double>& ff){
+ plumed_dbg_assert( ff.size()==forces.size() );
+ wasforced=true;
+ for(unsigned i=0;i<forces.size();++i) forces[i]+=ff[i];
+}
+
+bool StoreCentralAtomsVessel::applyForce(std::vector<double>& ff){
+ plumed_dbg_assert( ff.size()==forces.size() );
+ if(wasforced){
+ for(unsigned i=0;i<forces.size();++i) ff[i]=forces[i];
+ }
+ return wasforced;
+}
+
+void StoreCentralAtomsVessel::chainRuleForCentralAtom( const unsigned& iatom, const unsigned& iderno, const Vector& df, vesselbase::ActionWithVessel* act) const {
+ plumed_dbg_assert( iatom<mycolv->getNumberOfFunctionsInAction() );
+
+ unsigned nder=mycolv->getNumberOfDerivatives(iatom);
+ unsigned nder2=mycolv->getNumberOfDerivatives();
+ for(unsigned ider=0;ider<nder;++ider){
+ for(unsigned jcomp=0;jcomp<3;++jcomp){
+ unsigned ibuf=start[iatom] + jcomp*(nder+1) + ider;
+ act->addElementDerivative( iderno*nder2 + mycolv->getOutputDerivativeIndex(iatom, ider), df[jcomp]*getBufferElement(ibuf) );
+ }
+ }
+}
+
+
+}
+}
diff --git a/src/multicolvar/StoreCentralAtomsVessel.h b/src/multicolvar/StoreCentralAtomsVessel.h
new file mode 100644
index 000000000..eaa113e95
--- /dev/null
+++ b/src/multicolvar/StoreCentralAtomsVessel.h
@@ -0,0 +1,69 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2012 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.0.
+
+ 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_multicolvar_StoreCentralAtomsVessel_h
+#define __PLUMED_multicolvar_StoreCentralAtomsVessel_h
+
+#include "vesselbase/Vessel.h"
+#include "MultiColvar.h"
+
+namespace PLMD {
+namespace multicolvar {
+
+class MultiColvar;
+
+class StoreCentralAtomsVessel : public vesselbase::Vessel {
+private:
+ MultiColvar* mycolv;
+ std::vector<unsigned> start;
+ bool wasforced;
+ std::vector<double> forces;
+public:
+ static void reserveKeyword( Keywords& keys );
+ static void registerKeywords( Keywords& keys );
+/// Constructor
+ StoreCentralAtomsVessel( const vesselbase::VesselOptions& );
+/// Return the number of terms
+ unsigned getNumberOfTerms(){ return 2; }
+/// This does the resizing of the buffer
+ void resize();
+/// This does nothing
+ std::string description(){ return""; }
+/// This does nothing
+ void finish(){}
+/// Add some force to the atoms
+ void addForces( const std::vector<double>& );
+/// Copy forces to local arrays
+ bool applyForce(std::vector<double>&);
+/// Set all forces equal to zero before main loop
+ void prepare();
+/// This makes sure all vectors are stored
+ bool calculate();
+/// Get the orientation of the ith vector
+ Vector getPosition( const unsigned& ) const ;
+/// Chain rule for central atom
+ void chainRuleForCentralAtom( const unsigned& iatom, const unsigned& iderno, const Vector& df, vesselbase::ActionWithVessel* act ) const ;
+};
+
+}
+}
+#endif
+
diff --git a/src/vesselbase/ActionWithVessel.h b/src/vesselbase/ActionWithVessel.h
index 616be2f70..220cd8538 100644
--- a/src/vesselbase/ActionWithVessel.h
+++ b/src/vesselbase/ActionWithVessel.h
@@ -85,8 +85,6 @@ protected:
void runAllTasks( const unsigned& ntasks );
/// Resize all the functions when the number of derivatives change
void resizeFunctions();
-/// Add some derivative of the quantity in the sum wrt to a numbered element
- void addElementDerivative( const unsigned&, const double& );
/// Set the derivative of the jth element wrt to a numbered element
void setElementDerivative( const unsigned&, const double& );
public:
@@ -116,6 +114,8 @@ public:
virtual bool performTask( const unsigned& j )=0;
/// Return a pointer to the field
Vessel* getVessel( const std::string& name );
+/// Add some derivative of the quantity in the sum wrt to a numbered element
+ void addElementDerivative( const unsigned&, const double& );
/// Set the value of the element
void setElementValue( const unsigned& , const double& );
/// Get the value of this element
@@ -183,7 +183,8 @@ void ActionWithVessel::addElementDerivative( const unsigned& ider, const double&
inline
void ActionWithVessel::setElementDerivative( const unsigned& ider, const double& der ){
#ifndef NDEBUG
- if(ider==1) plumed_dbg_assert( weightHasDerivatives );
+ unsigned ndertmp=getNumberOfDerivatives();
+ if( ider>=ndertmp && ider<2*ndertmp ) plumed_dbg_assert( weightHasDerivatives );
#endif
plumed_dbg_assert( ider<derivatives.size() );
derivatives[ider] = der;
diff --git a/src/vesselbase/ShortcutVessel.h b/src/vesselbase/ShortcutVessel.h
index 374ec80ee..2ef876cb4 100644
--- a/src/vesselbase/ShortcutVessel.h
+++ b/src/vesselbase/ShortcutVessel.h
@@ -37,7 +37,7 @@ protected:
public:
static void registerKeywords( Keywords& keys );
ShortcutVessel( const VesselOptions& );
- std::string description(){ plumed_error(); }
+ std::string description(){ return ""; }
unsigned getNumberOfTerms(){ plumed_error(); return 0; }
void resize(){ plumed_error(); }
bool calculate(){ plumed_error(); }
diff --git a/src/vesselbase/Sum.cpp b/src/vesselbase/Sum.cpp
index 973b13d72..b3f8283cc 100644
--- a/src/vesselbase/Sum.cpp
+++ b/src/vesselbase/Sum.cpp
@@ -49,7 +49,7 @@ void Sum::reserveKeyword( Keywords& keys ){
Sum::Sum( const VesselOptions& da ) :
FunctionVessel(da),
-df(1,1.0)
+df(2)
{
}
@@ -69,7 +69,9 @@ bool Sum::calculate(){
}
void Sum::finish(){
- setOutputValue( getFinalValue(0) ); mergeFinalDerivatives( df );
+ setOutputValue( getFinalValue(0) );
+ df[0]=1.0; df[1]=0.0;
+ mergeFinalDerivatives( df );
}
diff --git a/src/vesselbase/Vessel.cpp b/src/vesselbase/Vessel.cpp
index 51bf34e0b..ce4265ee9 100644
--- a/src/vesselbase/Vessel.cpp
+++ b/src/vesselbase/Vessel.cpp
@@ -21,7 +21,6 @@
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
#include "ActionWithVessel.h"
#include "Vessel.h"
-#include "ShortcutVessel.h"
#include "tools/Exception.h"
#include "tools/Communicator.h"
#include "tools/Log.h"
@@ -116,8 +115,8 @@ void Vessel::checkRead(){
error(msg);
}
finished_read=true;
- ShortcutVessel* sv=dynamic_cast<ShortcutVessel*>(this);
- if(!sv) log.printf(" %s\n", description().c_str() );
+ std::string describe=description();
+ if( describe.length()>0 ) log.printf(" %s\n", describe.c_str() );
}
void Vessel::error( const std::string& msg ){
diff --git a/src/vesselbase/Vessel.h b/src/vesselbase/Vessel.h
index 464825d19..1f1689983 100644
--- a/src/vesselbase/Vessel.h
+++ b/src/vesselbase/Vessel.h
@@ -145,6 +145,8 @@ public:
std::string getLabel() const ;
/// Check that readin was fine
void checkRead();
+/// Set all the buffer elements to zero
+ void zero();
/// Add something to the ith element in the buffer
void addToBufferElement( const unsigned& i, const double& val);
/// Return a description of the vessel contents
@@ -230,6 +232,11 @@ ActionWithVessel* Vessel::getAction(){
return action;
}
+inline
+void Vessel::zero(){
+ for(unsigned i=0;i<bufsize;++i) setBufferElement( i, 0.0 );
+}
+
inline
void Vessel::setBufferElement( const unsigned& i, const double& val){
plumed_dbg_assert( i<bufsize );
--
GitLab