From 1d8567cb47218a9e6e5720b307975cfe8d3e21bd Mon Sep 17 00:00:00 2001
From: Gareth Tribello <gareth.tribello@gmail.com>
Date: Thu, 11 May 2017 20:33:51 +0100
Subject: [PATCH] Deleted non-working link cell stuff for topology matrix
---
src/adjmat/TopologyMatrix.cpp | 30 +-------------
src/multicolvar/AtomValuePack.cpp | 2 +-
src/multicolvar/MultiColvarBase.cpp | 9 +----
src/multicolvar/MultiColvarBase.h | 11 -----
src/tools/LinkCells.cpp | 63 -----------------------------
src/tools/LinkCells.h | 6 ---
6 files changed, 5 insertions(+), 116 deletions(-)
diff --git a/src/adjmat/TopologyMatrix.cpp b/src/adjmat/TopologyMatrix.cpp
index 92f56f839..72227515f 100644
--- a/src/adjmat/TopologyMatrix.cpp
+++ b/src/adjmat/TopologyMatrix.cpp
@@ -32,8 +32,6 @@ Adjacency matrix in which two atoms are adjacent if they are connected topologic
\par Examples
-\bug Link cells not working
-\bug Virial contribution not calculated correctly
*/
//+ENDPLUMEDOC
@@ -68,8 +66,6 @@ public:
unsigned getNumberOfQuantities() const ;
/// Create the ith, ith switching function
void setupConnector( const unsigned& id, const unsigned& i, const unsigned& j, const std::vector<std::string>& desc );
-/// Get the position that we should use as the center of our link cells
- Vector getLinkCellPosition( const std::vector<unsigned>& atoms ) const ;
/// This actually calculates the value of the contact function
double calculateWeight( const unsigned& taskCode, const double& weight, multicolvar::AtomValuePack& myatoms ) const ;
/// This does nothing
@@ -77,9 +73,6 @@ public:
/// Calculate the contribution from one of the atoms in the third element of the pack
void calculateForThreeAtoms( const unsigned& iat, const Vector& d1, const double& d1_len,
HistogramBead& bead, multicolvar::AtomValuePack& myatoms ) const ;
-///
- void buildListOfLinkCells( const std::vector<unsigned>& cind, const LinkCells& linkc,
- unsigned& ncells_required, std::vector<unsigned>& cells_required ) const ;
};
PLUMED_REGISTER_ACTION(TopologyMatrix,"TOPOLOGY_MATRIX")
@@ -151,8 +144,8 @@ maxbins(0)
beadrad = bead.getCutoff();
// Set the link cell cutoff
- log.printf(" setting cutoffs %f %f \n", sfmax, rfmax );
- setLinkCellCutoff( sfmax, std::numeric_limits<double>::max() ); // rfmax ); // std::numeric_limits<double>::max() );
+ log.printf(" setting cutoff %f \n", sfmax );
+ setLinkCellCutoff( sfmax, std::numeric_limits<double>::max() );
double maxsize=0;
for(unsigned i=0;i<getNumberOfNodeTypes();++i){
@@ -195,25 +188,6 @@ void TopologyMatrix::setupConnector( const unsigned& id, const unsigned& i, cons
}
}
-void TopologyMatrix::buildListOfLinkCells( const std::vector<unsigned>& cind, const LinkCells& linkc,
- unsigned& ncells_required, std::vector<unsigned>& cells_required ) const {
- // Shift ends of cylinder from positions of atoms to generate line of interest
- Vector distance=getSeparation( getPositionOfAtomForLinkCells( cind[0] ), getPositionOfAtomForLinkCells( cind[1] ) );
- double len = distance.modulo(); distance /= len;
- std::vector<Vector> p(2); p[0] = getPositionOfAtomForLinkCells( cind[0] ) - beadrad*distance;
- double binw = binw_mat( getBaseColvarNumber( cind[0] ), getBaseColvarNumber( cind[1] ) );
- double lcylinder = (std::floor( len / binw ) + 1)*binw;
- p[1] = p[0] + (lcylinder + 2*beadrad)*distance; pbcApply( p, 2 );
- // And get the appropriate atoms
- unsigned ncells=0; ncells_required = 0;
- linkc.getCellsThatLinePassesThrough( p[0], p[1], ncells_required, cells_required );
-}
-
-Vector TopologyMatrix::getLinkCellPosition( const std::vector<unsigned>& atoms ) const {
- Vector myatom = getPositionOfAtomForLinkCells( atoms[0] );
- return myatom + 0.5*pbcDistance( getPositionOfAtomForLinkCells( atoms[1] ), myatom );
-}
-
double TopologyMatrix::calculateWeight( const unsigned& taskCode, const double& weight, multicolvar::AtomValuePack& myatoms ) const {
Vector distance = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
if( distance.modulo2()<switchingFunction( getBaseColvarNumber( myatoms.getIndex(0) ), getBaseColvarNumber( myatoms.getIndex(1) ) ).get_dmax2() ) return 1.0;
diff --git a/src/multicolvar/AtomValuePack.cpp b/src/multicolvar/AtomValuePack.cpp
index c35883d7c..b9bf9793e 100644
--- a/src/multicolvar/AtomValuePack.cpp
+++ b/src/multicolvar/AtomValuePack.cpp
@@ -44,7 +44,7 @@ myatoms( vals.getAtomVector() )
unsigned AtomValuePack::setupAtomsFromLinkCells( const std::vector<unsigned>& cind, const Vector& cpos, const LinkCells& linkcells ){
if( cells_required.size()!=linkcells.getNumberOfCells() ) cells_required.resize( linkcells.getNumberOfCells() );
// Build the list of cells that we need
- unsigned ncells_required=0; mycolv->buildListOfLinkCells( cind, linkcells, ncells_required, cells_required );
+ unsigned ncells_required=0; linkcells.addRequiredCells( linkcells.findMyCell( cpos ), ncells_required, cells_required );
// Now build the list of atoms we need
natoms=cind.size(); for(unsigned i=0;i<natoms;++i) indices[i]=cind[i];
linkcells.retrieveAtomsInCells( ncells_required, cells_required, natoms, indices );
diff --git a/src/multicolvar/MultiColvarBase.cpp b/src/multicolvar/MultiColvarBase.cpp
index 66e8fb1fb..139ff16c7 100644
--- a/src/multicolvar/MultiColvarBase.cpp
+++ b/src/multicolvar/MultiColvarBase.cpp
@@ -643,14 +643,14 @@ bool MultiColvarBase::setupCurrentAtomList( const unsigned& taskCode, AtomValueP
myatoms.setNumberOfAtoms( 1 ); myatoms.setAtom( 0, taskCode ); return true;
} else if( usespecies ){
std::vector<unsigned> task_atoms(1); task_atoms[0]=taskCode;
- unsigned natomsper=myatoms.setupAtomsFromLinkCells( task_atoms, getLinkCellPosition(task_atoms), linkcells );
+ unsigned natomsper=myatoms.setupAtomsFromLinkCells( task_atoms, getPositionOfAtomForLinkCells( taskCode ), linkcells );
return natomsper>1;
} else if( matsums ){
myatoms.setNumberOfAtoms( getNumberOfAtoms() );
for(unsigned i=0;i<getNumberOfAtoms();++i) myatoms.setAtom( i, i );
} else if( allthirdblockintasks ){
plumed_dbg_assert( ablocks.size()==3 ); std::vector<unsigned> atoms(2); decodeIndexToAtoms( taskCode, atoms );
- myatoms.setupAtomsFromLinkCells( atoms, getLinkCellPosition(atoms), threecells );
+ myatoms.setupAtomsFromLinkCells( atoms, getPositionOfAtomForLinkCells( atoms[0] ), threecells );
} else if( nblock>0 ){
std::vector<unsigned> atoms( ablocks.size() );
decodeIndexToAtoms( taskCode, atoms ); myatoms.setNumberOfAtoms( ablocks.size() );
@@ -662,11 +662,6 @@ bool MultiColvarBase::setupCurrentAtomList( const unsigned& taskCode, AtomValueP
return true;
}
-void MultiColvarBase::buildListOfLinkCells( const std::vector<unsigned>& cind, const LinkCells& linkc,
- unsigned& ncells_required, std::vector<unsigned>& cells_required ) const {
- linkc.addRequiredCells( linkc.findMyCell( getPositionOfAtomForLinkCells(cind[0]) ), ncells_required, cells_required );
-}
-
void MultiColvarBase::setupActiveTaskSet( std::vector<unsigned>& active_tasks, const std::string& input_label ){
if( !setup_completed ){
bool justVolumes=false;
diff --git a/src/multicolvar/MultiColvarBase.h b/src/multicolvar/MultiColvarBase.h
index 0d5d8a664..0cd946996 100644
--- a/src/multicolvar/MultiColvarBase.h
+++ b/src/multicolvar/MultiColvarBase.h
@@ -172,14 +172,8 @@ public:
virtual void performTask( const unsigned& , const unsigned& , MultiValue& ) const ;
/// Update the active atoms
virtual void updateActiveAtoms( AtomValuePack& myatoms ) const ;
-/// This builds the list of link cells that we need during this calculation. Its a virtual function
-/// because it has to be redefined in TopologyMatrix
- virtual void buildListOfLinkCells( const std::vector<unsigned>& cind, const LinkCells& linkc,
- unsigned& ncells_required, std::vector<unsigned>& cells_required ) 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 ;
/// Get the absolute index of the central atom
virtual AtomNumber getAbsoluteIndexOfCentralAtom( const unsigned& i ) const ;
/// This is replaced once we have a function to calculate the cv
@@ -246,11 +240,6 @@ Vector MultiColvarBase::getPositionOfAtomForLinkCells( const unsigned& iatom ) c
return ActionAtomistic::getPosition( atom_lab[iatom].second );
}
-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/LinkCells.cpp b/src/tools/LinkCells.cpp
index 31aa4a54e..8805ff113 100644
--- a/src/tools/LinkCells.cpp
+++ b/src/tools/LinkCells.cpp
@@ -102,69 +102,6 @@ void LinkCells::buildCellLists( const std::vector<Vector>& pos, const std::vecto
}
}
-bool LinkCells::checkLineBox( const double& dist1, const double& dist2, const Vector& fpos1, const Vector& fdir,
- const Vector& plow, const Vector& phigh, const unsigned& axis ) const {
- if( dist1*dist2>=0.0 || dist1==dist2 ) return false;
- Vector hit = fpos1 + fdir * ( -dist1/(dist2-dist1) );
- if( axis==0 && hit[2]>plow[2] && hit[2]<phigh[2] && hit[1]>plow[1] && hit[1]<phigh[1] ) return true;
- if( axis==1 && hit[2]>plow[2] && hit[2]<phigh[2] && hit[0]>plow[0] && hit[0]<phigh[0] ) return true;
- if( axis==2 && hit[0]>plow[0] && hit[0]<phigh[0] && hit[1]>plow[1] && hit[1]<phigh[1] ) return true;
- return false;
-}
-
-void LinkCells::getCellsThatLinePassesThrough( const Vector& pos1, const Vector& pos2, unsigned& ncells_required,
- std::vector<unsigned>& cells_required ) const {
- // Retrieve the cell indices of the extrema for the line segment
- std::vector<double> delx(3); std::vector<int> celln(3); std::vector<unsigned> celli(3);
- for(unsigned i=0;i<3;++i) delx[i] = 1.0 / static_cast<double>(ncells[i]);
-
- // Get the vector connecting the two points
- Vector dir = mypbc.distance( pos1, pos2 );
- // Get the indices of link cell containing the first point
- std::vector<int> cd(3), c2(3); std::vector<unsigned> c1 = findMyCell( pos1 );
- // Now find the position of the second point in nearest cell
- Vector wpos2 = pos1 + dir; Vector fpos2 = mypbc.realToScaled( wpos2 );
- // And the link cell the second point is in
- for(unsigned j=0;j<3;++j){
- c2[j] = std::floor( ( fpos2[j] + 0.5 ) * ncells[j] );
- cd[j] = (c2[j]<static_cast<int>(c1[j]))? -1 : +1;
- c2[j] += cd[j];
- }
- // Now loop over relevant cells
- Vector plow, phigh, cplow, cphigh;
- for(celln[0]=c1[0];celln[0]!=c2[0];celln[0]+=cd[0]){
- plow[0] = -0.5 + celln[0] * delx[0]; phigh[0] = plow[0] + delx[0];
- for(celln[1]=c1[1];celln[1]!=c2[1];celln[1]+=cd[1]){
- plow[1] = -0.5 + celln[1] * delx[1]; phigh[1] = plow[1] + delx[1];
- for(celln[2]=c1[2];celln[2]!=c2[2];celln[2]+=cd[2]){
- plow[2] = -0.5 + celln[2] * delx[2]; phigh[2] = plow[2] + delx[2];
- cplow = mypbc.scaledToReal( plow ); cphigh = mypbc.scaledToReal( phigh );
- for(unsigned j=0;j<3;++j){
- celli[j] = (celln[j]<0)? (ncells[j]+celln[j])%ncells[j] : celln[j]%ncells[j];
- plumed_assert( celli[j]>=0 && celli[j]<ncells[j] );
- }
- addRequiredCells( celli, ncells_required, cells_required );
- // if( pos1[0]>cplow[0] && pos1[0]<cphigh[0] &&
- // pos1[1]>cplow[1] && pos1[1]<cphigh[1] &&
- // pos1[2]>cplow[2] && pos1[2]<cphigh[2] ){
- // for(unsigned j=0;j<3;++j) celli[j] = (celln[j]<0)? (ncells[j]+celln[j])%ncells[j] : celln[j]%ncells[j];
- // addRequiredCells( celli, ncells_required, cells_required );
- // continue;
- // }
- // if( checkLineBox( pos1[0] - cplow[0], wpos2[0] - cplow[0], pos1, dir, cplow, cphigh, 0 ) ||
- // checkLineBox( pos1[1] - cplow[1], wpos2[1] - cplow[1], pos1, dir, cplow, cphigh, 1 ) ||
- // checkLineBox( pos1[2] - cplow[2], wpos2[2] - cplow[2], pos1, dir, cplow, cphigh, 2 ) ||
- // checkLineBox( pos1[0] - cphigh[0], wpos2[0] - cphigh[0], pos1, dir, cplow, cphigh, 0 ) ||
- // checkLineBox( pos1[1] - cphigh[1], wpos2[1] - cphigh[1], pos1, dir, cplow, cphigh, 1 ) ||
- // checkLineBox( pos1[2] - cphigh[2], wpos2[2] - cphigh[2], pos1, dir, cplow, cphigh, 2 ) ){
- // for(unsigned j=0;j<3;++j) celli[j] = (celln[j]<0)? (ncells[j]+celln[j])%ncells[j] : celln[j]%ncells[j];
- // addRequiredCells( celli, ncells_required, cells_required );
- // }
- }
- }
- }
-}
-
#define LINKC_MIN(n) ((n<2)? 0 : -1)
#define LINKC_MAX(n) ((n<3)? 1 : 2)
#define LINKC_PBC(n,num) ((n<0)? num-1 : n%num )
diff --git a/src/tools/LinkCells.h b/src/tools/LinkCells.h
index f336d612c..e5659c2d1 100644
--- a/src/tools/LinkCells.h
+++ b/src/tools/LinkCells.h
@@ -54,9 +54,6 @@ private:
std::vector<unsigned> lcell_tots;
/// The atoms ordered by link cells
std::vector<unsigned> lcell_lists;
-/// Checks if a line passes through a link cell
- bool checkLineBox( const double& dist1, const double& dist2, const Vector& fpos1, const Vector& fpos2,
- const Vector& plow, const Vector& phigh, const unsigned& axis ) const ;
public:
///
explicit LinkCells( Communicator& comm );
@@ -76,9 +73,6 @@ public:
unsigned findCell( const Vector& pos ) const ;
/// Find the cell in which this position is contained
std::vector<unsigned> findMyCell( const Vector& pos ) const ;
-/// Get the cells that a line passes through so we can build a list of atoms
- void getCellsThatLinePassesThrough( const Vector& pos1, const Vector& pos2, unsigned& ncells_required,
- std::vector<unsigned>& cells_required ) const ;
/// Get the list of cells we need to surround the a particular cell
void addRequiredCells( const std::vector<unsigned>& celn, unsigned& ncells_required,
std::vector<unsigned>& cells_required ) const ;
--
GitLab