Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-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 "tools/OFile.h"
24 : #include "core/PlumedMain.h"
25 : #include "core/ActionSet.h"
26 : #include "core/ActionPilot.h"
27 : #include "core/ActionRegister.h"
28 :
29 : //+PLUMEDOC CONCOMP OUTPUT_CLUSTER
30 : /*
31 : Output the indices of the atoms in one of the clusters identified by a clustering object
32 :
33 : This action provides one way of getting output from a \ref DFSCLUSTERING calculation.
34 : The output in question here is either
35 :
36 : - a file that contains a list of the atom indices that form part of one of the clusters that was identified using \ref DFSCLUSTERING
37 : - an xyz file containing the positions of the atoms in one of the the clusters that was identified using \ref DFSCLUSTERING
38 :
39 : Notice also that if you choose to output an xyz file you can ask PLUMED to try to reconstruct the cluster
40 : taking the periodic boundary conditions into account by using the MAKE_WHOLE flag.
41 :
42 : \par Examples
43 :
44 : The input shown below identifies those atoms with a coordination number less than 13
45 : and then constructs a contact matrix that describes the connectivity between the atoms
46 : that satisfy this criteria. The DFS algorithm is then used to find the connected components
47 : in this matrix and the indices of the atoms in the largest connected component are then output
48 : to a file.
49 :
50 : \plumedfile
51 : c1: COORDINATIONNUMBER SPECIES=1-1996 SWITCH={CUBIC D_0=0.34 D_MAX=0.38}
52 : cf: MFILTER_LESS DATA=c1 SWITCH={CUBIC D_0=13 D_MAX=13.5}
53 : mat: CONTACT_MATRIX ATOMS=cf SWITCH={CUBIC D_0=0.34 D_MAX=0.38}
54 : dfs: DFSCLUSTERING MATRIX=mat
55 : OUTPUT_CLUSTER CLUSTERS=dfs CLUSTER=1 FILE=dfs.dat
56 : \endplumedfile
57 :
58 : */
59 : //+ENDPLUMEDOC
60 :
61 : namespace PLMD {
62 : namespace adjmat {
63 :
64 : class OutputCluster : public ActionPilot {
65 : private:
66 : bool makewhole, output_xyz;
67 : OFile ofile;
68 : ClusteringBase* myclusters;
69 : double rcut2;
70 : unsigned clustr, maxdepth, maxgoes;
71 : std::vector<bool> visited;
72 : std::vector<unsigned> myatoms;
73 : std::vector<Vector> atomsin;
74 : std::vector<unsigned> nneigh;
75 : Matrix<unsigned> adj_list;
76 : int number_of_cluster;
77 : std::vector< std::pair<unsigned,unsigned> > cluster_sizes;
78 : std::vector<unsigned> which_cluster;
79 : bool explore_dfs( const unsigned& index );
80 : void explore( const unsigned& index, const unsigned& depth );
81 : public:
82 : static void registerKeywords( Keywords& keys );
83 : explicit OutputCluster(const ActionOptions&);
84 8 : void calculate() override {}
85 8 : void apply() override {}
86 : void update() override;
87 : };
88 :
89 10435 : PLUMED_REGISTER_ACTION(OutputCluster,"OUTPUT_CLUSTER")
90 :
91 9 : void OutputCluster::registerKeywords( Keywords& keys ) {
92 9 : Action::registerKeywords( keys );
93 9 : ActionPilot::registerKeywords( keys );
94 18 : keys.add("compulsory","CLUSTERS","the action that performed the clustering");
95 18 : keys.add("compulsory","CLUSTER","1","which cluster would you like to look at 1 is the largest cluster, 2 is the second largest, 3 is the the third largest and so on");
96 18 : keys.add("compulsory","STRIDE","1","the frequency with which you would like to output the atoms in the cluster");
97 18 : keys.add("compulsory","FILE","the name of the file on which to output the details of the cluster");
98 18 : keys.add("compulsory","MAXDEPTH","6","maximum depth for searches over paths to reconstruct clusters for PBC");
99 18 : keys.add("compulsory","MAXGOES","200","number of times to run searches to reconstruct clusters");
100 18 : keys.addFlag("MAKE_WHOLE",false,"reconstruct the clusters and remove all periodic boundary conditions.");
101 9 : }
102 :
103 8 : OutputCluster::OutputCluster(const ActionOptions& ao):
104 : Action(ao),
105 : ActionPilot(ao),
106 8 : myclusters(NULL)
107 : {
108 : // Setup output file
109 16 : ofile.link(*this); std::string file; parse("FILE",file);
110 8 : if( file.length()==0 ) error("output file name was not specified");
111 : // Search for xyz extension
112 8 : output_xyz=false;
113 8 : if( file.find(".")!=std::string::npos ) {
114 : std::size_t dot=file.find_first_of('.');
115 16 : if( file.substr(dot+1)=="xyz" ) output_xyz=true;
116 : }
117 :
118 8 : ofile.open(file); log.printf(" on file %s \n",file.c_str());
119 24 : parseFlag("MAKE_WHOLE",makewhole); parse("MAXDEPTH",maxdepth); parse("MAXGOES",maxgoes);
120 8 : if( makewhole && !output_xyz) error("MAKE_WHOLE flag is not compatible with output of non-xyz files");
121 :
122 : // Find what action we are taking the clusters from
123 8 : std::vector<std::string> matname(1); parse("CLUSTERS",matname[0]);
124 8 : myclusters = plumed.getActionSet().selectWithLabel<ClusteringBase*>( matname[0] );
125 8 : if( !myclusters ) error( matname[0] + " does not calculate perform a clustering of the atomic positions");
126 : // N.B. the +0.3 is a fudge factor. Reconstrucing PBC doesnt work without this GAT
127 8 : addDependency( myclusters ); double rcut=myclusters->getCutoffForConnection() + 0.3; rcut2=rcut*rcut;
128 :
129 : // Read in the cluster we are calculating
130 8 : parse("CLUSTER",clustr);
131 8 : if( clustr<1 ) error("cannot look for a cluster larger than the largest cluster");
132 8 : if( clustr>myclusters->getNumberOfNodes() ) error("cluster selected is invalid - too few atoms in system");
133 8 : log.printf(" outputting atoms in %u th largest cluster found by %s \n",clustr,matname[0].c_str() );
134 16 : }
135 :
136 8 : void OutputCluster::update() {
137 8 : myclusters->retrieveAtomsInCluster( clustr, myatoms );
138 8 : if( output_xyz ) {
139 0 : ofile.printf("%u \n",static_cast<unsigned>(myatoms.size()));
140 0 : ofile.printf("atoms in %u th largest cluster \n",clustr );
141 0 : if( makewhole ) {
142 : // Retrieve the atom positions
143 0 : atomsin.resize( myatoms.size() );
144 0 : for(unsigned i=0; i<myatoms.size(); ++i) atomsin[i]=myclusters->getPositionOfAtomForLinkCells( myatoms[i] );
145 : // Build a connectivity matrix neglecting the pbc
146 0 : nneigh.resize( myatoms.size(), 0 ); adj_list.resize( myatoms.size(), myatoms.size() );
147 0 : for(unsigned i=1; i<myatoms.size(); ++i) {
148 0 : for(unsigned j=0; j<i; ++j) {
149 0 : if( delta( atomsin[i], atomsin[j] ).modulo2()<=rcut2 ) { adj_list(i,nneigh[i])=j; adj_list(j,nneigh[j])=i; nneigh[i]++; nneigh[j]++; }
150 : }
151 : }
152 : // Use DFS to find the largest cluster not broken by PBC
153 0 : number_of_cluster=-1; visited.resize( myatoms.size(), false );
154 0 : cluster_sizes.resize( myatoms.size() ); which_cluster.resize( myatoms.size() );
155 0 : for(unsigned i=0; i<cluster_sizes.size(); ++i) { cluster_sizes[i].first=0; cluster_sizes[i].second=i; }
156 :
157 0 : for(unsigned i=0; i<myatoms.size(); ++i) {
158 0 : if( !visited[i] ) { number_of_cluster++; visited[i]=explore_dfs(i); }
159 : }
160 0 : std::sort( cluster_sizes.begin(), cluster_sizes.end() );
161 :
162 : // Now set visited so that only those atoms in largest cluster will be start points for PBCing
163 : visited.assign( visited.size(), false );
164 0 : for(unsigned i=0; i<myatoms.size(); ++i) {
165 0 : if( which_cluster[i]==cluster_sizes[cluster_sizes.size()-1].second ) visited[i]=true;
166 : }
167 :
168 : // Now retrieve the original connectivity matrix (including pbc)
169 0 : nneigh.assign( nneigh.size(), 0 );
170 0 : for(unsigned i=1; i<myatoms.size(); ++i) {
171 0 : for(unsigned j=0; j<i; ++j) {
172 0 : if( myclusters->areConnected( myatoms[i], myatoms[j] ) ) { adj_list(i,nneigh[i])=j; adj_list(j,nneigh[j])=i; nneigh[i]++; nneigh[j]++; }
173 : }
174 : }
175 :
176 : // Now find broken bonds and run iterative deepening depth first search to reconstruct
177 0 : for(unsigned jj=0; jj<maxgoes; ++jj) {
178 :
179 0 : for(unsigned j=0; j<myatoms.size(); ++j) {
180 0 : if( !visited[j] ) continue;
181 :
182 0 : for(unsigned k=0; k<nneigh[j]; ++k) {
183 0 : if( delta( atomsin[j],atomsin[adj_list(j,k)] ).modulo2()>rcut2 ) {
184 : visited[j]=true;
185 0 : for(unsigned depth=0; depth<=maxdepth; ++depth) explore( j, depth );
186 : }
187 : }
188 : }
189 : }
190 : // And print final positions
191 0 : for(unsigned i=0; i<myatoms.size(); ++i) ofile.printf( "X %f %f %f \n", atomsin[i][0], atomsin[i][1], atomsin[i][2] );
192 : } else {
193 0 : for(unsigned i=0; i<myatoms.size(); ++i) {
194 0 : Vector pos=myclusters->getPositionOfAtomForLinkCells( myatoms[i] );
195 0 : ofile.printf( "X %f %f %f \n", pos[0], pos[1], pos[2] );
196 : }
197 : }
198 : } else {
199 8 : ofile.printf("CLUSTERING RESULTS AT TIME %f : NUMBER OF ATOMS IN %u TH LARGEST CLUSTER EQUALS %u \n",getTime(),clustr,static_cast<unsigned>(myatoms.size()) );
200 8 : ofile.printf("INDICES OF ATOMS : ");
201 512 : for(unsigned i=0; i<myatoms.size(); ++i) ofile.printf("%d ",(myclusters->getAbsoluteIndexOfCentralAtom(myatoms[i])).index());
202 8 : ofile.printf("\n");
203 : }
204 8 : }
205 :
206 0 : void OutputCluster::explore( const unsigned& index, const unsigned& depth ) {
207 0 : if( depth==0 ) return ;
208 :
209 0 : for(unsigned i=0; i<nneigh[index]; ++i) {
210 0 : unsigned j=adj_list(index,i); visited[j]=true;
211 0 : Vector svec=myclusters->pbcDistance( atomsin[index], atomsin[j] );
212 0 : atomsin[j] = atomsin[index] + svec;
213 0 : explore( j, depth-1 );
214 : }
215 : }
216 :
217 0 : bool OutputCluster::explore_dfs( const unsigned& index ) {
218 0 : visited[index]=true;
219 0 : for(unsigned i=0; i<nneigh[index]; ++i) {
220 0 : unsigned j=adj_list(index,i);
221 0 : if( !visited[j] ) visited[j]=explore_dfs(j);
222 : }
223 :
224 : // Count the size of the cluster
225 0 : cluster_sizes[number_of_cluster].first++;
226 0 : which_cluster[index] = number_of_cluster;
227 0 : return visited[index];
228 : }
229 :
230 : }
231 : }
232 :
233 :
|
