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 13801 : PLUMED_REGISTER_ACTION(OutputCluster,"OUTPUT_CLUSTER")
90 :
91 12 : void OutputCluster::registerKeywords( Keywords& keys ) {
92 12 : Action::registerKeywords( keys );
93 12 : ActionPilot::registerKeywords( keys );
94 24 : keys.add("compulsory","CLUSTERS","the action that performed the clustering");
95 24 : 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 24 : keys.add("compulsory","STRIDE","1","the frequency with which you would like to output the atoms in the cluster");
97 24 : keys.add("compulsory","FILE","the name of the file on which to output the details of the cluster");
98 24 : keys.add("compulsory","MAXDEPTH","6","maximum depth for searches over paths to reconstruct clusters for PBC");
99 24 : keys.add("compulsory","MAXGOES","200","number of times to run searches to reconstruct clusters");
100 24 : keys.addFlag("MAKE_WHOLE",false,"reconstruct the clusters and remove all periodic boundary conditions.");
101 12 : }
102 :
103 8 : OutputCluster::OutputCluster(const ActionOptions& ao):
104 : Action(ao),
105 : ActionPilot(ao),
106 8 : myclusters(NULL) {
107 : // Setup output file
108 8 : ofile.link(*this);
109 : std::string file;
110 16 : parse("FILE",file);
111 8 : if( file.length()==0 ) {
112 0 : error("output file name was not specified");
113 : }
114 : // Search for xyz extension
115 8 : output_xyz=false;
116 8 : if( file.find(".")!=std::string::npos ) {
117 : std::size_t dot=file.find_first_of('.');
118 16 : if( file.substr(dot+1)=="xyz" ) {
119 0 : output_xyz=true;
120 : }
121 : }
122 :
123 8 : ofile.open(file);
124 8 : log.printf(" on file %s \n",file.c_str());
125 8 : parseFlag("MAKE_WHOLE",makewhole);
126 8 : parse("MAXDEPTH",maxdepth);
127 8 : parse("MAXGOES",maxgoes);
128 8 : if( makewhole && !output_xyz) {
129 0 : error("MAKE_WHOLE flag is not compatible with output of non-xyz files");
130 : }
131 :
132 : // Find what action we are taking the clusters from
133 8 : std::vector<std::string> matname(1);
134 8 : parse("CLUSTERS",matname[0]);
135 8 : myclusters = plumed.getActionSet().selectWithLabel<ClusteringBase*>( matname[0] );
136 8 : if( !myclusters ) {
137 0 : error( matname[0] + " does not calculate perform a clustering of the atomic positions");
138 : }
139 : // N.B. the +0.3 is a fudge factor. Reconstrucing PBC doesnt work without this GAT
140 8 : addDependency( myclusters );
141 8 : double rcut=myclusters->getCutoffForConnection() + 0.3;
142 8 : rcut2=rcut*rcut;
143 :
144 : // Read in the cluster we are calculating
145 8 : parse("CLUSTER",clustr);
146 8 : if( clustr<1 ) {
147 0 : error("cannot look for a cluster larger than the largest cluster");
148 : }
149 8 : if( clustr>myclusters->getNumberOfNodes() ) {
150 0 : error("cluster selected is invalid - too few atoms in system");
151 : }
152 8 : log.printf(" outputting atoms in %u th largest cluster found by %s \n",clustr,matname[0].c_str() );
153 16 : }
154 :
155 8 : void OutputCluster::update() {
156 8 : myclusters->retrieveAtomsInCluster( clustr, myatoms );
157 8 : if( output_xyz ) {
158 0 : ofile.printf("%u \n",static_cast<unsigned>(myatoms.size()));
159 0 : ofile.printf("atoms in %u th largest cluster \n",clustr );
160 0 : if( makewhole ) {
161 : // Retrieve the atom positions
162 0 : atomsin.resize( myatoms.size() );
163 0 : for(unsigned i=0; i<myatoms.size(); ++i) {
164 0 : atomsin[i]=myclusters->getPositionOfAtomForLinkCells( myatoms[i] );
165 : }
166 : // Build a connectivity matrix neglecting the pbc
167 0 : nneigh.resize( myatoms.size(), 0 );
168 0 : adj_list.resize( myatoms.size(), myatoms.size() );
169 0 : for(unsigned i=1; i<myatoms.size(); ++i) {
170 0 : for(unsigned j=0; j<i; ++j) {
171 0 : if( delta( atomsin[i], atomsin[j] ).modulo2()<=rcut2 ) {
172 0 : adj_list(i,nneigh[i])=j;
173 0 : adj_list(j,nneigh[j])=i;
174 0 : nneigh[i]++;
175 0 : nneigh[j]++;
176 : }
177 : }
178 : }
179 : // Use DFS to find the largest cluster not broken by PBC
180 0 : number_of_cluster=-1;
181 0 : visited.resize( myatoms.size(), false );
182 0 : cluster_sizes.resize( myatoms.size() );
183 0 : which_cluster.resize( myatoms.size() );
184 0 : for(unsigned i=0; i<cluster_sizes.size(); ++i) {
185 0 : cluster_sizes[i].first=0;
186 0 : cluster_sizes[i].second=i;
187 : }
188 :
189 0 : for(unsigned i=0; i<myatoms.size(); ++i) {
190 0 : if( !visited[i] ) {
191 0 : number_of_cluster++;
192 0 : visited[i]=explore_dfs(i);
193 : }
194 : }
195 0 : std::sort( cluster_sizes.begin(), cluster_sizes.end() );
196 :
197 : // Now set visited so that only those atoms in largest cluster will be start points for PBCing
198 : visited.assign( visited.size(), false );
199 0 : for(unsigned i=0; i<myatoms.size(); ++i) {
200 0 : if( which_cluster[i]==cluster_sizes[cluster_sizes.size()-1].second ) {
201 : visited[i]=true;
202 : }
203 : }
204 :
205 : // Now retrieve the original connectivity matrix (including pbc)
206 0 : nneigh.assign( nneigh.size(), 0 );
207 0 : for(unsigned i=1; i<myatoms.size(); ++i) {
208 0 : for(unsigned j=0; j<i; ++j) {
209 0 : if( myclusters->areConnected( myatoms[i], myatoms[j] ) ) {
210 0 : adj_list(i,nneigh[i])=j;
211 0 : adj_list(j,nneigh[j])=i;
212 0 : nneigh[i]++;
213 0 : nneigh[j]++;
214 : }
215 : }
216 : }
217 :
218 : // Now find broken bonds and run iterative deepening depth first search to reconstruct
219 0 : for(unsigned jj=0; jj<maxgoes; ++jj) {
220 :
221 0 : for(unsigned j=0; j<myatoms.size(); ++j) {
222 0 : if( !visited[j] ) {
223 0 : continue;
224 : }
225 :
226 0 : for(unsigned k=0; k<nneigh[j]; ++k) {
227 0 : if( delta( atomsin[j],atomsin[adj_list(j,k)] ).modulo2()>rcut2 ) {
228 : visited[j]=true;
229 0 : for(unsigned depth=0; depth<=maxdepth; ++depth) {
230 0 : explore( j, depth );
231 : }
232 : }
233 : }
234 : }
235 : }
236 : // And print final positions
237 0 : for(unsigned i=0; i<myatoms.size(); ++i) {
238 0 : ofile.printf( "X %f %f %f \n", atomsin[i][0], atomsin[i][1], atomsin[i][2] );
239 : }
240 : } else {
241 0 : for(unsigned i=0; i<myatoms.size(); ++i) {
242 0 : Vector pos=myclusters->getPositionOfAtomForLinkCells( myatoms[i] );
243 0 : ofile.printf( "X %f %f %f \n", pos[0], pos[1], pos[2] );
244 : }
245 : }
246 : } else {
247 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()) );
248 8 : ofile.printf("INDICES OF ATOMS : ");
249 512 : for(unsigned i=0; i<myatoms.size(); ++i) {
250 504 : ofile.printf("%d ",(myclusters->getAbsoluteIndexOfCentralAtom(myatoms[i])).index());
251 : }
252 8 : ofile.printf("\n");
253 : }
254 8 : }
255 :
256 0 : void OutputCluster::explore( const unsigned& index, const unsigned& depth ) {
257 0 : if( depth==0 ) {
258 : return ;
259 : }
260 :
261 0 : for(unsigned i=0; i<nneigh[index]; ++i) {
262 0 : unsigned j=adj_list(index,i);
263 : visited[j]=true;
264 0 : Vector svec=myclusters->pbcDistance( atomsin[index], atomsin[j] );
265 0 : atomsin[j] = atomsin[index] + svec;
266 0 : explore( j, depth-1 );
267 : }
268 : }
269 :
270 0 : bool OutputCluster::explore_dfs( const unsigned& index ) {
271 0 : visited[index]=true;
272 0 : for(unsigned i=0; i<nneigh[index]; ++i) {
273 0 : unsigned j=adj_list(index,i);
274 0 : if( !visited[j] ) {
275 0 : visited[j]=explore_dfs(j);
276 : }
277 : }
278 :
279 : // Count the size of the cluster
280 0 : cluster_sizes[number_of_cluster].first++;
281 0 : which_cluster[index] = number_of_cluster;
282 0 : return visited[index];
283 : }
284 :
285 : }
286 : }
287 :
288 :
|
