Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2014-2023 The plumed team
       3             :    (see the PEOPLE file at the root of the distribution for a list of names)
       4             : 
       5             :    See http://www.plumed.org for more information.
       6             : 
       7             :    This file is part of plumed, version 2.
       8             : 
       9             :    plumed is free software: you can redistribute it and/or modify
      10             :    it under the terms of the GNU Lesser General Public License as published by
      11             :    the Free Software Foundation, either version 3 of the License, or
      12             :    (at your option) any later version.
      13             : 
      14             :    plumed is distributed in the hope that it will be useful,
      15             :    but WITHOUT ANY WARRANTY; without even the implied warranty of
      16             :    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      17             :    GNU Lesser General Public License for more details.
      18             : 
      19             :    You should have received a copy of the GNU Lesser General Public License
      20             :    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
      21             : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
      22             : #include "ClusteringBase.h"
      23             : #include "AdjacencyMatrixVessel.h"
      24             : #include "core/ActionRegister.h"
      25             : 
      26             : #ifdef __PLUMED_HAS_BOOST_GRAPH
      27             : #include <boost/graph/adjacency_list.hpp>
      28             : #include <boost/graph/connected_components.hpp>
      29             : #include <boost/graph/graph_utility.hpp>
      30             : #endif
      31             : 
      32             : //+PLUMEDOC MATRIXF DFSCLUSTERING
      33             : /*
      34             : Find the connected components of the matrix using the depth first search clustering algorithm.
      35             : 
      36             : As discussed in the section of the manual on \ref contactmatrix a useful tool for developing complex collective variables is the notion of the
      37             : so called adjacency matrix.  An adjacency matrix is an \f$N \times N\f$ matrix in which the \f$i\f$th, \f$j\f$th element tells you whether
      38             : or not the \f$i\f$th and \f$j\f$th atoms/molecules from a set of \f$N\f$ atoms/molecules are adjacent or not.  As detailed in \cite tribello-clustering
      39             : these matrices provide a representation of a graph and can thus can be analyzed using tools from graph theory.  This particular action performs
      40             : a depth first search clustering to find the connected components of this graph.  You can read more about depth first search here:
      41             : 
      42             : https://en.wikipedia.org/wiki/Depth-first_search
      43             : 
      44             : This action is useful if you are looking at a phenomenon such as nucleation where the aim is to detect the sizes of the crystalline nuclei that have formed
      45             : in your simulation cell.
      46             : 
      47             : \par Examples
      48             : 
      49             : The input below calculates the coordination numbers of atoms 1-100 and then computes the an adjacency
      50             : matrix whose elements measures whether atoms \f$i\f$ and \f$j\f$ are within 0.55 nm of each other.  The action
      51             : labelled dfs then treats the elements of this matrix as zero or ones and thus thinks of the matrix as defining
      52             : a graph.  This dfs action then finds the largest connected component in this graph.  The sum of the coordination
      53             : numbers for the atoms in this largest connected component are then computed and this quantity is output to a colvar
      54             : file.  The way this input can be used is described in detail in \cite tribello-clustering.
      55             : 
      56             : \plumedfile
      57             : lq: COORDINATIONNUMBER SPECIES=1-100 SWITCH={CUBIC D_0=0.45  D_MAX=0.55} LOWMEM
      58             : cm: CONTACT_MATRIX ATOMS=lq  SWITCH={CUBIC D_0=0.45  D_MAX=0.55}
      59             : dfs: DFSCLUSTERING MATRIX=cm
      60             : clust1: CLUSTER_PROPERTIES CLUSTERS=dfs CLUSTER=1 SUM
      61             : PRINT ARG=clust1.* FILE=colvar
      62             : \endplumedfile
      63             : 
      64             : */
      65             : //+ENDPLUMEDOC
      66             : 
      67             : namespace PLMD {
      68             : namespace adjmat {
      69             : 
      70             : class DFSClustering : public ClusteringBase {
      71             : private:
      72             : #ifdef __PLUMED_HAS_BOOST_GRAPH
      73             : /// The list of edges in the graph
      74             :   std::vector<std::pair<unsigned,unsigned> > edge_list;
      75             : #else
      76             : /// The number of neighbors each atom has
      77             :   std::vector<unsigned> nneigh;
      78             : /// The adjacency list
      79             :   Matrix<unsigned> adj_list;
      80             : /// The color that tells us whether a node has been visited
      81             :   std::vector<unsigned> color;
      82             : /// The recursive function at the heart of this method
      83             :   int explore( const unsigned& index );
      84             : #endif
      85             : public:
      86             : /// Create manual
      87             :   static void registerKeywords( Keywords& keys );
      88             : /// Constructor
      89             :   explicit DFSClustering(const ActionOptions&);
      90             : /// Do the clustering
      91             :   void performClustering() override;
      92             : };
      93             : 
      94       13815 : PLUMED_REGISTER_ACTION(DFSClustering,"DFSCLUSTERING")
      95             : 
      96          19 : void DFSClustering::registerKeywords( Keywords& keys ) {
      97          19 :   ClusteringBase::registerKeywords( keys );
      98          38 :   keys.add("compulsory","MAXCONNECT","0","maximum number of connections that can be formed by any given node in the graph. "
      99             :            "By default this is set equal to zero and the number of connections is set equal to the number "
     100             :            "of nodes.  You only really need to set this if you are working with a very large system and "
     101             :            "memory is at a premium");
     102          19 : }
     103             : 
     104          15 : DFSClustering::DFSClustering(const ActionOptions&ao):
     105             :   Action(ao),
     106          15 :   ClusteringBase(ao) {
     107             :   unsigned maxconnections;
     108          15 :   parse("MAXCONNECT",maxconnections);
     109             : #ifdef __PLUMED_HAS_BOOST_GRAPH
     110             :   if( maxconnections>0 ) {
     111             :     edge_list.resize( getNumberOfNodes()*maxconnections );
     112             :   } else {
     113             :     edge_list.resize(0.5*getNumberOfNodes()*(getNumberOfNodes()-1));
     114             :   }
     115             : #else
     116          15 :   nneigh.resize( getNumberOfNodes() );
     117          15 :   color.resize(getNumberOfNodes());
     118          15 :   if( maxconnections>0 ) {
     119           0 :     adj_list.resize(getNumberOfNodes(),maxconnections);
     120             :   } else {
     121          15 :     adj_list.resize(getNumberOfNodes(),getNumberOfNodes());
     122             :   }
     123             : #endif
     124          15 : }
     125             : 
     126          24 : void DFSClustering::performClustering() {
     127             : #ifdef __PLUMED_HAS_BOOST_GRAPH
     128             :   // Get the list of edges
     129             :   unsigned nedges=0;
     130             :   getAdjacencyVessel()->retrieveEdgeList( nedges, edge_list );
     131             : 
     132             :   // Build the graph using boost
     133             :   boost::adjacency_list<boost::vecS,boost::vecS,boost::undirectedS> sg(&edge_list[0],&edge_list[nedges],getNumberOfNodes());
     134             : 
     135             :   // Find the connected components using boost (-1 here for compatibility with non-boost version)
     136             :   number_of_cluster=boost::connected_components(sg,&which_cluster[0]) - 1;
     137             : 
     138             :   // And work out the size of each cluster
     139             :   for(unsigned i=0; i<which_cluster.size(); ++i) {
     140             :     cluster_sizes[which_cluster[i]].first++;
     141             :   }
     142             : #else
     143             :   // Get the adjacency matrix
     144          24 :   getAdjacencyVessel()->retrieveAdjacencyLists( nneigh, adj_list );
     145             : 
     146             :   // Perform clustering
     147          24 :   number_of_cluster=-1;
     148          24 :   color.assign(color.size(),0);
     149        8217 :   for(unsigned i=0; i<getNumberOfNodes(); ++i) {
     150        8193 :     if( color[i]==0 ) {
     151        5814 :       number_of_cluster++;
     152        5814 :       color[i]=explore(i);
     153             :     }
     154             :   }
     155             : #endif
     156          24 : }
     157             : 
     158             : #ifndef __PLUMED_HAS_BOOST_GRAPH
     159        8193 : int DFSClustering::explore( const unsigned& index ) {
     160             : 
     161        8193 :   color[index]=1;
     162       38219 :   for(unsigned i=0; i<nneigh[index]; ++i) {
     163       30026 :     unsigned j=adj_list(index,i);
     164       30026 :     if( color[j]==0 ) {
     165        2379 :       color[j]=explore(j);
     166             :     }
     167             :   }
     168             : 
     169             :   // Count the size of the cluster
     170        8193 :   cluster_sizes[number_of_cluster].first++;
     171        8193 :   which_cluster[index] = number_of_cluster;
     172        8193 :   return color[index];
     173             : }
     174             : #endif
     175             : 
     176             : }
     177             : }