Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2014-2019 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 analysed 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 45 : 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();
92 : };
93 :
94 6467 : PLUMED_REGISTER_ACTION(DFSClustering,"DFSCLUSTERING")
95 :
96 16 : void DFSClustering::registerKeywords( Keywords& keys ) {
97 16 : ClusteringBase::registerKeywords( keys );
98 80 : 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 16 : }
103 :
104 15 : DFSClustering::DFSClustering(const ActionOptions&ao):
105 : Action(ao),
106 30 : ClusteringBase(ao)
107 : {
108 30 : unsigned maxconnections; parse("MAXCONNECT",maxconnections);
109 : #ifdef __PLUMED_HAS_BOOST_GRAPH
110 : if( maxconnections>0 ) edge_list.resize( getNumberOfNodes()*maxconnections );
111 : else edge_list.resize(0.5*getNumberOfNodes()*(getNumberOfNodes()-1));
112 : #else
113 15 : nneigh.resize( getNumberOfNodes() ); color.resize(getNumberOfNodes());
114 15 : if( maxconnections>0 ) adj_list.resize(getNumberOfNodes(),maxconnections);
115 15 : else adj_list.resize(getNumberOfNodes(),getNumberOfNodes());
116 : #endif
117 15 : }
118 :
119 24 : void DFSClustering::performClustering() {
120 : #ifdef __PLUMED_HAS_BOOST_GRAPH
121 : // Get the list of edges
122 : unsigned nedges=0; getAdjacencyVessel()->retrieveEdgeList( nedges, edge_list );
123 :
124 : // Build the graph using boost
125 : boost::adjacency_list<boost::vecS,boost::vecS,boost::undirectedS> sg(&edge_list[0],&edge_list[nedges],getNumberOfNodes());
126 :
127 : // Find the connected components using boost (-1 here for compatibility with non-boost version)
128 : number_of_cluster=boost::connected_components(sg,&which_cluster[0]) - 1;
129 :
130 : // And work out the size of each cluster
131 : for(unsigned i=0; i<which_cluster.size(); ++i) cluster_sizes[which_cluster[i]].first++;
132 : #else
133 : // Get the adjacency matrix
134 24 : getAdjacencyVessel()->retrieveAdjacencyLists( nneigh, adj_list );
135 :
136 : // Perform clustering
137 48 : number_of_cluster=-1; color.assign(color.size(),0);
138 8217 : for(unsigned i=0; i<getNumberOfNodes(); ++i) {
139 22200 : if( color[i]==0 ) { number_of_cluster++; color[i]=explore(i); }
140 : }
141 : #endif
142 24 : }
143 :
144 : #ifndef __PLUMED_HAS_BOOST_GRAPH
145 8193 : int DFSClustering::explore( const unsigned& index ) {
146 :
147 16386 : color[index]=1;
148 106464 : for(unsigned i=0; i<nneigh[index]; ++i) {
149 30026 : unsigned j=adj_list(index,i);
150 62431 : if( color[j]==0 ) color[j]=explore(j);
151 : }
152 :
153 : // Count the size of the cluster
154 16386 : cluster_sizes[number_of_cluster].first++;
155 16386 : which_cluster[index] = number_of_cluster;
156 16386 : return color[index];
157 : }
158 : #endif
159 :
160 : }
161 4839 : }
|
