From 8435cbf5913ce8fca5fa58a735e18be8a5e8cd68 Mon Sep 17 00:00:00 2001
From: Gareth Tribello <gt@eider.phy.qub.ac.uk>
Date: Fri, 30 Oct 2015 10:01:09 +0000
Subject: [PATCH] Some small changes that allow for opitimisations of
multicolvars
It is now possible to specify a different link cell cutoff for the
third bodies. In addition, the way that connectivity is determined
in adjacency matrix actions works better given that we can now
operate with non-contact matrix contact matrices. Also added some
cutoff stuff to HistogramBead. This is currently not used in many
places
---
src/adjmat/AdjacencyMatrixBase.cpp | 4 ++--
src/adjmat/AdjacencyMatrixBase.h | 9 +++++++
src/adjmat/AdjacencyMatrixVessel.cpp | 33 +++++++++++++++-----------
src/adjmat/ContactAlignedMatrix.cpp | 3 +++
src/adjmat/ContactMatrix.cpp | 2 ++
src/multicolvar/AtomValuePack.cpp | 5 ++--
src/multicolvar/MultiColvarBase.cpp | 10 ++++----
src/multicolvar/MultiColvarBase.h | 9 ++++++-
src/tools/HistogramBead.cpp | 35 +++++++++++++++++++++++++++-
src/tools/HistogramBead.h | 6 +++++
10 files changed, 92 insertions(+), 24 deletions(-)
diff --git a/src/adjmat/AdjacencyMatrixBase.cpp b/src/adjmat/AdjacencyMatrixBase.cpp
index 037ca71e8..1fc3357cd 100644
--- a/src/adjmat/AdjacencyMatrixBase.cpp
+++ b/src/adjmat/AdjacencyMatrixBase.cpp
@@ -143,9 +143,9 @@ void AdjacencyMatrixBase::requestAtoms( const std::vector<AtomNumber>& atoms, co
if( symmetric ){
for(unsigned i=1;i<dims[0];++i){
for(unsigned j=0;j<i;++j){
- bookeeping(i,j).first=getFullNumberOfTasks();
+ bookeeping(j,i).first=bookeeping(i,j).first=getFullNumberOfTasks();
for(unsigned k=0;k<kcount;++k) addTaskToList( i*icoef + j*jcoef + k*kcoef );
- bookeeping(i,j).second=getFullNumberOfTasks();
+ bookeeping(j,i).second=bookeeping(i,j).second=getFullNumberOfTasks();
}
}
} else {
diff --git a/src/adjmat/AdjacencyMatrixBase.h b/src/adjmat/AdjacencyMatrixBase.h
index 019b5934e..3b269b60d 100644
--- a/src/adjmat/AdjacencyMatrixBase.h
+++ b/src/adjmat/AdjacencyMatrixBase.h
@@ -78,6 +78,10 @@ public:
Vector getCentralAtom(){ plumed_merror("cannot find central atoms for adjacency matrix actions"); Vector dum; return dum; }
/// Get the absolute index of an atom
// AtomNumber getAbsoluteIndexOfCentralAtom(const unsigned& i) const ;
+/// Transforms the stored values in whatever way is required
+ virtual double transformStoredValues( const std::vector<double>& myvals, unsigned& vout, double& df ) const ;
+/// Used to check for connections between atoms
+ virtual bool checkForConnection( const std::vector<double>& myvals ) const=0;
};
inline
@@ -98,6 +102,11 @@ void AdjacencyMatrixBase::getOrientationVector( const unsigned& ind, const bool&
mybasedata[mmc]->retrieveValue( convertToLocalIndex(ind,mmc), normed, orient );
}
+inline
+double AdjacencyMatrixBase::transformStoredValues( const std::vector<double>& myvals, unsigned& vout, double& df ) const {
+ plumed_dbg_assert( myvals.size()==2 ); vout=1; df=1; return myvals[1];
+}
+
}
}
diff --git a/src/adjmat/AdjacencyMatrixVessel.cpp b/src/adjmat/AdjacencyMatrixVessel.cpp
index dbfc6989a..180a20a0b 100644
--- a/src/adjmat/AdjacencyMatrixVessel.cpp
+++ b/src/adjmat/AdjacencyMatrixVessel.cpp
@@ -87,7 +87,8 @@ void AdjacencyMatrixVessel::getMatrixIndices( const unsigned& code, unsigned& i,
}
void AdjacencyMatrixVessel::retrieveMatrix( DynamicList<unsigned>& myactive_elements, Matrix<double>& mymatrix ){
- myactive_elements.deactivateAll(); std::vector<double> vals( getNumberOfComponents() );
+ unsigned vin; double df;
+ myactive_elements.deactivateAll(); std::vector<double> vals( getNumberOfComponents() );
for(unsigned i=0;i<getNumberOfStoredValues();++i){
// Ignore any non active members
if( !storedValueIsActive(i) ) continue ;
@@ -95,11 +96,8 @@ void AdjacencyMatrixVessel::retrieveMatrix( DynamicList<unsigned>& myactive_elem
unsigned j, k; getMatrixIndices( i, k, j );
retrieveValue( i, false, vals );
- // Find the maximum element in the stored values
- double max=vals[1]; for(unsigned j=2;j<vals.size();++j){ if( vals[j]>max ) max=vals[j]; }
-
- if( symmetric ) mymatrix(k,j)=mymatrix(j,k)=max / vals[0];
- else mymatrix(k,j) = max / vals[0];
+ if( symmetric ) mymatrix(k,j)=mymatrix(j,k)=function->transformStoredValues( vals, vin, df );
+ else mymatrix(k,j)=function->transformStoredValues( vals, vin, df );
}
myactive_elements.updateActiveMembers();
}
@@ -110,10 +108,17 @@ void AdjacencyMatrixVessel::retrieveAdjacencyLists( std::vector<unsigned>& nneig
for(unsigned i=0;i<nneigh.size();++i) nneigh[i]=0;
// And set up the adjacency list
+ std::vector<double> myvals( getNumberOfComponents() );
for(unsigned i=0;i<getNumberOfStoredValues();++i){
// Ignore any non active members
if( !storedValueIsActive(i) ) continue ;
- unsigned j, k; getMatrixIndices( i, k, j );
+ // Check if atoms are connected
+ retrieveValue( i, false, myvals );
+ unsigned j, k; getMatrixIndices( i, k, j );
+ if( !function->checkForConnection( myvals ) ) continue ;
+
+ // Store if atoms are connected
+ // unsigned j, k; getMatrixIndices( i, k, j );
adj_list(k,nneigh[k])=j; nneigh[k]++;
adj_list(j,nneigh[j])=k; nneigh[j]++;
}
@@ -121,9 +126,14 @@ void AdjacencyMatrixVessel::retrieveAdjacencyLists( std::vector<unsigned>& nneig
void AdjacencyMatrixVessel::retrieveEdgeList( unsigned& nedge, std::vector<std::pair<unsigned,unsigned> >& edge_list ){
plumed_dbg_assert( undirectedGraph() ); nedge=0;
+ std::vector<double> myvals( getNumberOfComponents() );
for(unsigned i=0;i<getNumberOfStoredValues();++i){
// Ignore any non active members
if( !storedValueIsActive(i) ) continue ;
+ // Check if atoms are connected
+ retrieveValue( i, false, myvals );
+ if( !function->checkForConnection( myvals ) ) continue ;
+
getMatrixIndices( i, edge_list[nedge].first, edge_list[nedge].second );
nedge++;
}
@@ -133,16 +143,13 @@ void AdjacencyMatrixVessel::retrieveDerivatives( const unsigned& myelem, const b
StoreDataVessel::retrieveDerivatives( myelem, normed, myvals );
if( !function->weightHasDerivatives ) return ;
- std::vector<double> vals( getNumberOfComponents() ); retrieveValue( myelem, normed, vals );
- unsigned vi=1; double max=vals[1];
- for(unsigned j=2;j<vals.size();++j){
- if( vals[j]>max ){ vi=j; max=vals[j]; }
- }
+ unsigned vi; std::vector<double> vals( getNumberOfComponents() ); retrieveValue( myelem, normed, vals );
+ double df, max=function->transformStoredValues( vals, vi, df );
double pref = max/(vals[0]*vals[0]);
for(unsigned i=0;i<myvals.getNumberActive();++i){
unsigned jder=myvals.getActiveIndex(i);
- myvals.setDerivative( 1, jder, myvals.getDerivative(vi, jder)/vals[0] - pref*myvals.getDerivative(0, jder) );
+ myvals.setDerivative( 1, jder, df*myvals.getDerivative(vi, jder)/vals[0] - pref*myvals.getDerivative(0, jder) );
}
}
diff --git a/src/adjmat/ContactAlignedMatrix.cpp b/src/adjmat/ContactAlignedMatrix.cpp
index 0e4764231..a90c242a2 100644
--- a/src/adjmat/ContactAlignedMatrix.cpp
+++ b/src/adjmat/ContactAlignedMatrix.cpp
@@ -54,6 +54,9 @@ public:
void calculateWeight( const unsigned& taskCode, multicolvar::AtomValuePack& myatoms ) const ;
/// This does nothing
double compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const ;
+///
+ /// Used to check for connections between atoms
+ bool checkForConnection( const std::vector<double>& myvals ) const { return !(myvals[1]>epsilon); }
};
PLUMED_REGISTER_ACTION(ContactAlignedMatrix,"ALIGNED_MATRIX")
diff --git a/src/adjmat/ContactMatrix.cpp b/src/adjmat/ContactMatrix.cpp
index c1284f4bb..111424b22 100644
--- a/src/adjmat/ContactMatrix.cpp
+++ b/src/adjmat/ContactMatrix.cpp
@@ -55,6 +55,8 @@ public:
void calculateWeight( const unsigned& taskCode, multicolvar::AtomValuePack& myatoms ) const ;
/// This does nothing
double compute( const unsigned& tindex, multicolvar::AtomValuePack& myatoms ) const ;
+/// Used to check for connections between atoms
+ virtual bool checkForConnection( const std::vector<double>& myvals ) const { return true; }
};
PLUMED_REGISTER_ACTION(ContactMatrix,"CONTACT_MATRIX")
diff --git a/src/multicolvar/AtomValuePack.cpp b/src/multicolvar/AtomValuePack.cpp
index c9ee792a0..66f7e4cbf 100644
--- a/src/multicolvar/AtomValuePack.cpp
+++ b/src/multicolvar/AtomValuePack.cpp
@@ -43,9 +43,8 @@ myatoms( vals.getAtomVector() )
unsigned AtomValuePack::setupAtomsFromLinkCells( const std::vector<unsigned>& cind, const Vector& cpos, const LinkCells& linkcells ){
natoms=cind.size(); for(unsigned i=0;i<natoms;++i) indices[i]=cind[i];
- linkcells.retrieveNeighboringAtoms( cpos, natoms, indices ); myatoms[0]=Vector(0.0,0.0,0.0);
- for(unsigned i=1;i<cind.size();++i) myatoms[i]=mycolv->getPositionOfAtomForLinkCells( cind[i] ) - cpos;
- for(unsigned i=cind.size();i<natoms;++i) myatoms[i]=mycolv->getPositionOfAtomForLinkCells( indices[i] ) - cpos;
+ linkcells.retrieveNeighboringAtoms( cpos, natoms, indices );
+ for(unsigned i=0;i<natoms;++i) myatoms[i]=mycolv->getPositionOfAtomForLinkCells( indices[i] ) - cpos;
if( mycolv->usesPbc() ) mycolv->applyPbc( myatoms, natoms );
return natoms;
}
diff --git a/src/multicolvar/MultiColvarBase.cpp b/src/multicolvar/MultiColvarBase.cpp
index 422867966..915de80e8 100644
--- a/src/multicolvar/MultiColvarBase.cpp
+++ b/src/multicolvar/MultiColvarBase.cpp
@@ -190,9 +190,11 @@ void MultiColvarBase::turnOnDerivatives(){
forcesToApply.resize( getNumberOfDerivatives() );
}
-void MultiColvarBase::setLinkCellCutoff( const double& lcut ){
+void MultiColvarBase::setLinkCellCutoff( const double& lcut, double tcut ){
plumed_assert( usespecies || ablocks.size()<4 );
- linkcells.setCutoff( lcut ); threecells.setCutoff( lcut );
+ if( tcut<0 ) tcut=lcut;
+ linkcells.setCutoff( lcut );
+ threecells.setCutoff( tcut );
}
void MultiColvarBase::setupLinkCells(){
@@ -395,11 +397,11 @@ bool MultiColvarBase::setupCurrentAtomList( const unsigned& taskCode, AtomValueP
if( usespecies ){
if( isDensity() ) return true;
std::vector<unsigned> task_atoms(1); task_atoms[0]=taskCode;
- unsigned natomsper=myatoms.setupAtomsFromLinkCells( task_atoms, getPositionOfAtomForLinkCells(taskCode), linkcells );
+ unsigned natomsper=myatoms.setupAtomsFromLinkCells( task_atoms, getLinkCellPosition(task_atoms), linkcells );
return natomsper>1;
} else if( allthirdblockintasks ){
plumed_dbg_assert( ablocks.size()==3 ); std::vector<unsigned> atoms(2); decodeIndexToAtoms( taskCode, atoms );
- unsigned natomsper=myatoms.setupAtomsFromLinkCells( atoms, getPositionOfAtomForLinkCells(atoms[0]), threecells );
+ unsigned natomsper=myatoms.setupAtomsFromLinkCells( atoms, getLinkCellPosition(atoms), threecells );
} else if( ablocks.size()<4 ){
std::vector<unsigned> atoms( ablocks.size() );
decodeIndexToAtoms( taskCode, atoms ); myatoms.setNumberOfAtoms( ablocks.size() );
diff --git a/src/multicolvar/MultiColvarBase.h b/src/multicolvar/MultiColvarBase.h
index 55c64f55e..b3eaed4bb 100644
--- a/src/multicolvar/MultiColvarBase.h
+++ b/src/multicolvar/MultiColvarBase.h
@@ -94,7 +94,7 @@ protected:
/// Set which atoms are to be used to calculate the central atom position
void setAtomsForCentralAtom( const std::vector<bool>& catom_ind );
/// Set the value of the cutoff for the link cells
- void setLinkCellCutoff( const double& lcut );
+ void setLinkCellCutoff( const double& lcut, double tcut=-1.0 );
/// Setup link cells in order to make this calculation faster
void setupLinkCells();
/// Get the separation between a pair of vectors
@@ -123,6 +123,8 @@ public:
virtual void updateActiveAtoms( AtomValuePack& myatoms ) const ;
/// This gets the position of an atom for the link cell setup
virtual Vector getPositionOfAtomForLinkCells( const unsigned& iatom ) const ;
+/// Returns the position where we should assume the center is for link cell calculations
+ virtual Vector getLinkCellPosition( const std::vector<unsigned>& atoms ) const ;
/// And a virtual function which actually computes the colvar
virtual double doCalculation( const unsigned& tindex, AtomValuePack& myatoms ) const ;
/// Get the absolute index of the central atom
@@ -187,6 +189,11 @@ Vector MultiColvarBase::getPositionOfAtomForLinkCells( const unsigned& iatom ) c
return ActionAtomistic::getPosition( iatom );
}
+inline
+Vector MultiColvarBase::getLinkCellPosition( const std::vector<unsigned>& atoms ) const {
+ return getPositionOfAtomForLinkCells( atoms[0] );
+}
+
inline
unsigned MultiColvarBase::getNumberOfDerivatives(){
return 3*getNumberOfAtoms()+9;
diff --git a/src/tools/HistogramBead.cpp b/src/tools/HistogramBead.cpp
index 5ca59f9ed..2d7261721 100644
--- a/src/tools/HistogramBead.cpp
+++ b/src/tools/HistogramBead.cpp
@@ -157,7 +157,11 @@ void HistogramBead::set( const std::string& params, std::string& errormsg ){
}
void HistogramBead::set( double l, double h, double w){
- init=true; lowb=l; highb=h; width=w;
+ init=true; lowb=l; highb=h; width=w;
+ const double DP2CUTOFF=6.25;
+ if( type==gaussian ) cutoff=sqrt(2.0*DP2CUTOFF);
+ else if( type==triangular ) cutoff=1.;
+ else plumed_error();
}
void HistogramBead::setKernelType( const std::string& ktype ){
@@ -194,6 +198,35 @@ double HistogramBead::calculate( double x, double& df ) const {
return f;
}
+double HistogramBead::calculateWithCutoff( double x, double& df ) const {
+ plumed_dbg_assert(init && periodicity!=unset );
+
+ double lowB, upperB, f;
+ lowB = difference( x, lowb ) / width ; upperB = difference( x, highb ) / width;
+ if( upperB<=-cutoff || lowB>=cutoff ){ double df=0; return 0; }
+
+ if( type==gaussian ){
+ lowB /= sqrt(2.0); upperB /= sqrt(2.0);
+ df = ( exp( -lowB*lowB ) - exp( -upperB*upperB ) ) / ( sqrt(2*pi)*width );
+ f = 0.5*( erf( upperB ) - erf( lowB ) );
+ } else if( type==triangular ){
+ df=0;
+ if( fabs(lowB)<1. ) df = (1 - fabs(lowB)) / width;
+ if( fabs(upperB)<1. ) df -= (1 - fabs(upperB)) / width;
+ if (upperB<=-1. || lowB >=1.){
+ f=0.;
+ } else {
+ double ia, ib;
+ if( lowB>-1.0 ){ ia=lowB; }else{ ia=-1.0; }
+ if( upperB<1.0 ){ ib=upperB; } else{ ib=1.0; }
+ f = (ib*(2.-fabs(ib))-ia*(2.-fabs(ia)))*0.5;
+ }
+ } else {
+ plumed_merror("function type does not exist");
+ }
+ return f;
+}
+
double HistogramBead::lboundDerivative( const double& x ) const {
double lowB;
if( type==gaussian ){
diff --git a/src/tools/HistogramBead.h b/src/tools/HistogramBead.h
index a8c78b156..3c7b892b9 100644
--- a/src/tools/HistogramBead.h
+++ b/src/tools/HistogramBead.h
@@ -44,6 +44,7 @@ private:
double lowb;
double highb;
double width;
+ double cutoff;
enum {gaussian,triangular} type;
enum {unset,periodic,notperiodic} periodicity;
double min, max, max_minus_min, inv_max_minus_min;
@@ -60,10 +61,12 @@ public:
void set(const std::string& params, std::string& errormsg);
void set(double l, double h, double w);
double calculate(double x, double&df) const;
+ double calculateWithCutoff( double x, double& df ) const;
double lboundDerivative( const double& x ) const;
double uboundDerivative( const double& x ) const;
double getlowb() const ;
double getbigb() const ;
+ double getCutoff() const ;
};
inline
@@ -90,6 +93,9 @@ double HistogramBead::getlowb() const { return lowb; }
inline
double HistogramBead::getbigb() const { return highb; }
+inline
+double HistogramBead::getCutoff() const { return cutoff*width; }
+
inline
double HistogramBead::difference( const double& d1, const double& d2 ) const {
if(periodicity==notperiodic){
--
GitLab