Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-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 "AdjacencyMatrixBase.h"
25 : #include "core/ActionRegister.h"
26 :
27 : //+PLUMEDOC MATRIXF CLUSTER_WITHSURFACE
28 : /*
29 : Take a connected component that was found using a clustering algorithm and create a new cluster that contains those atoms that are in the cluster together with those atoms that are within a certain cutoff of the cluster.
30 :
31 : As discussed in the section of the manual on \ref contactmatrix a useful tool for developing complex collective variables is the notion of the
32 : 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
33 : 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. When analysing these matrix
34 : we can treat them as a graph and find connected components using some clustering algorithm. This action is used in tandem with this form of analysis
35 : and takes one of the connected components that was found during this analysis and creates a new cluster that includes all the atoms within the
36 : connected component that was found together that were within a certain cutoff distance of the atoms in the connected component. This form of analysis
37 : has been used sucessfully in the forward flux sampling simulations described in this paper \cite gab-ice-kaolinite
38 :
39 : \par Examples
40 :
41 : The following input uses PLUMED to calculate a adjacency matrix that connects a pair of atoms if they both have a coordination number that is less
42 : than 13.5 and if they are within 0.38 nm of each other. Depth first search clustering is used to find the connected components in this matrix. The
43 : number of atoms with indices that are between 1 and 1996 and that are either in the second largest cluster or that are within within 0.3 nm of one of the
44 : atoms within the the second largest cluster are then counted and this number of atoms is output to a file called size. In addition the indices of the atoms
45 : that were counted are output to a file called dfs2.dat.
46 :
47 : \plumedfile
48 : c1: COORDINATIONNUMBER SPECIES=1-1996 SWITCH={CUBIC D_0=0.34 D_MAX=0.38}
49 : cf: MFILTER_LESS DATA=c1 SWITCH={CUBIC D_0=13 D_MAX=13.5}
50 : mat: CONTACT_MATRIX ATOMS=cf SWITCH={CUBIC D_0=0.34 D_MAX=0.38}
51 : dfs: DFSCLUSTERING MATRIX=mat
52 : clust2a: CLUSTER_WITHSURFACE CLUSTERS=dfs RCUT_SURF=0.3
53 : size2a: CLUSTER_NATOMS CLUSTERS=clust2a CLUSTER=2
54 : PRINT ARG=size2a FILE=size FMT=%8.4f
55 : OUTPUT_CLUSTER CLUSTERS=clust2a CLUSTER=2 FILE=dfs2.dat
56 : \endplumedfile
57 :
58 :
59 : */
60 : //+ENDPLUMEDOC
61 :
62 : namespace PLMD {
63 : namespace adjmat {
64 :
65 4 : class ClusterWithSurface : public ClusteringBase {
66 : private:
67 : /// The clusters that we are adding surface atoms to
68 : ClusteringBase* myclusters;
69 : /// The cutoff for surface atoms
70 : double rcut_surf2;
71 : public:
72 : /// Create manual
73 : static void registerKeywords( Keywords& keys );
74 : /// Constructor
75 : explicit ClusterWithSurface(const ActionOptions&);
76 : ///
77 : unsigned getNumberOfDerivatives();
78 : ///
79 : unsigned getNumberOfNodes() const ;
80 : ///
81 : AtomNumber getAbsoluteIndexOfCentralAtom(const unsigned& i) const ;
82 : ///
83 : void retrieveAtomsInCluster( const unsigned& clust, std::vector<unsigned>& myatoms ) const ;
84 : ///
85 : void getInputData( const unsigned& ind, const bool& normed, const multicolvar::AtomValuePack& myatoms, std::vector<double>& orient0 ) const ;
86 : ///
87 : MultiValue& getInputDerivatives( const unsigned& ind, const bool& normed, const multicolvar::AtomValuePack& myatoms ) const ;
88 : ///
89 : unsigned getNumberOfQuantities() const ;
90 : /// Do the calculation
91 2 : void performClustering() {};
92 : ///
93 : double getCutoffForConnection() const ;
94 : ///
95 : Vector getPositionOfAtomForLinkCells( const unsigned& taskIndex ) const ;
96 : };
97 :
98 6454 : PLUMED_REGISTER_ACTION(ClusterWithSurface,"CLUSTER_WITHSURFACE")
99 :
100 3 : void ClusterWithSurface::registerKeywords( Keywords& keys ) {
101 3 : ClusteringBase::registerKeywords( keys );
102 6 : keys.remove("MATRIX");
103 12 : keys.add("compulsory","CLUSTERS","the label of the action that does the clustering");
104 12 : keys.add("compulsory","RCUT_SURF","you also have the option to find the atoms on the surface of the cluster. An atom must be within this distance of one of the atoms "
105 : "of the cluster in order to be considered a surface atom");
106 3 : }
107 :
108 2 : ClusterWithSurface::ClusterWithSurface(const ActionOptions&ao):
109 : Action(ao),
110 2 : ClusteringBase(ao)
111 : {
112 : std::vector<AtomNumber> fake_atoms;
113 4 : if( !parseMultiColvarAtomList("CLUSTERS",-1,fake_atoms ) ) error("unable to find CLUSTERS input");
114 2 : if( mybasemulticolvars.size()!=1 ) error("should be exactly one multicolvar input");
115 :
116 : // Retrieve the adjacency matrix of interest
117 4 : atom_lab.resize(0); myclusters = dynamic_cast<ClusteringBase*>( mybasemulticolvars[0] );
118 2 : if( !myclusters ) error( mybasemulticolvars[0]->getLabel() + " does not calculate clusters");
119 :
120 : // Setup switching function for surface atoms
121 4 : double rcut_surf; parse("RCUT_SURF",rcut_surf);
122 2 : if( rcut_surf>0 ) log.printf(" counting surface atoms that are within %f of the cluster atoms \n",rcut_surf);
123 2 : rcut_surf2=rcut_surf*rcut_surf;
124 :
125 : // And now finish the setup of everything in the base
126 2 : setupMultiColvarBase( fake_atoms );
127 2 : }
128 :
129 4 : unsigned ClusterWithSurface::getNumberOfDerivatives() {
130 4 : return myclusters->getNumberOfDerivatives();
131 : }
132 :
133 16155834 : unsigned ClusterWithSurface::getNumberOfNodes() const {
134 16155834 : return myclusters->getNumberOfNodes();
135 : }
136 :
137 128 : AtomNumber ClusterWithSurface::getAbsoluteIndexOfCentralAtom(const unsigned& i) const {
138 128 : return myclusters->getAbsoluteIndexOfCentralAtom(i);
139 : }
140 :
141 0 : void ClusterWithSurface::getInputData( const unsigned& ind, const bool& normed, const multicolvar::AtomValuePack& myatoms, std::vector<double>& orient0 ) const {
142 0 : myclusters->getInputData( ind, normed, myatoms, orient0 );
143 0 : }
144 :
145 0 : MultiValue& ClusterWithSurface::getInputDerivatives( const unsigned& ind, const bool& normed, const multicolvar::AtomValuePack& myatoms ) const {
146 0 : return myclusters->getInputDerivatives( ind, normed, myatoms );
147 : }
148 :
149 14 : unsigned ClusterWithSurface::getNumberOfQuantities() const {
150 14 : return myclusters->getNumberOfQuantities();
151 : }
152 :
153 2 : double ClusterWithSurface::getCutoffForConnection() const {
154 2 : double tcut = myclusters->getCutoffForConnection();
155 2 : if( tcut>sqrt(rcut_surf2) ) return tcut;
156 0 : return sqrt(rcut_surf2);
157 : }
158 :
159 6 : void ClusterWithSurface::retrieveAtomsInCluster( const unsigned& clust, std::vector<unsigned>& myatoms ) const {
160 6 : std::vector<unsigned> tmpat; myclusters->retrieveAtomsInCluster( clust, tmpat );
161 :
162 : // Prevent double counting
163 6 : std::vector<bool> incluster( getNumberOfNodes(), false );
164 348 : for(unsigned i=0; i<tmpat.size(); ++i) incluster[tmpat[i]]=true;
165 :
166 : // Find the atoms in the the clusters
167 6 : std::vector<bool> surface_atom( getNumberOfNodes(), false );
168 264 : for(unsigned i=0; i<tmpat.size(); ++i) {
169 335412 : for(unsigned j=0; j<getNumberOfNodes(); ++j) {
170 335328 : if( incluster[j] ) continue;
171 499464 : double dist2=getSeparation( getPosition(tmpat[i]), getPosition(j) ).modulo2();
172 166488 : if( dist2<rcut_surf2 ) { surface_atom[j]=true; }
173 : }
174 : }
175 : unsigned nsurf_at=0;
176 23958 : for(unsigned j=0; j<getNumberOfNodes(); ++j) {
177 23952 : if( surface_atom[j] ) nsurf_at++;
178 : }
179 12 : myatoms.resize( nsurf_at + tmpat.size() );
180 264 : for(unsigned i=0; i<tmpat.size(); ++i) myatoms[i]=tmpat[i];
181 6 : unsigned nn=tmpat.size();
182 23958 : for(unsigned j=0; j<getNumberOfNodes(); ++j) {
183 24252 : if( surface_atom[j] ) { myatoms[nn]=j; nn++; }
184 : }
185 12 : plumed_assert( nn==myatoms.size() );
186 6 : }
187 :
188 0 : Vector ClusterWithSurface::getPositionOfAtomForLinkCells( const unsigned& iatom ) const {
189 0 : return myclusters->getPositionOfAtomForLinkCells( iatom );
190 : }
191 :
192 : }
193 4839 : }
|
