Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-2018 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 "ClusterAnalysisBase.h"
23 : #include "core/ActionRegister.h"
24 :
25 : //+PLUMEDOC CONCOMP CLUSTER_DIAMETER
26 : /*
27 : Print out the diameter of one of the connected components
28 :
29 : As discussed in the section of the manual on \ref contactmatrix a useful tool for developing complex collective variables is the notion of the
30 : 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
31 : 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
32 : 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
33 : to output the largest of the distances between the paris of atoms that are connected together in a particular connected component. It is important to
34 : note that the quantity that is output by this action is not differentiable. As such it cannot be used as a collective variable in a biased simulation.
35 :
36 : \par Examples
37 :
38 : 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 greater
39 : than 2.0 and if they are within 6.0 nm of each other. Depth first search clustering is used to find the connected components in this matrix. The distance
40 : between every pair of atoms that are within the largest of the clusters found is then calculated and the largest of these distances is output to a file named
41 : colvar.
42 :
43 : \verbatim
44 : # Calculate coordination numbers
45 : c1: COORDINATIONNUMBER SPECIES=1-512 SWITCH={EXP D_0=4.0 R_0=0.5 D_MAX=6.0}
46 : # Select coordination numbers that are more than 2.0
47 : cf: MFILTER_MORE DATA=c1 SWITCH={RATIONAL D_0=2.0 R_0=0.1} LOWMEM
48 : # Build a contact matrix
49 : mat: CONTACT_MATRIX ATOMS=cf SWITCH={EXP D_0=4.0 R_0=0.5 D_MAX=6.0}
50 : # Find largest cluster
51 : dfs: DFSCLUSTERING MATRIX=mat
52 : clust1: CLUSTER_PROPERTIES CLUSTERS=dfs CLUSTER=1
53 : dia: CLUSTER_DIAMETER CLUSTERS=dfs CLUSTER=1
54 : PRINT ARG=dia FILE=colvar
55 : \endverbatim
56 :
57 : */
58 : //+ENDPLUMEDOC
59 :
60 : namespace PLMD {
61 : namespace adjmat {
62 :
63 4 : class ClusterDiameter : public ClusterAnalysisBase {
64 : private:
65 : /// The cluster we are looking for
66 : unsigned clustr;
67 : public:
68 : /// Create manual
69 : static void registerKeywords( Keywords& keys );
70 : /// Constructor
71 : explicit ClusterDiameter(const ActionOptions&);
72 : ///
73 : void calculate();
74 : ///
75 : void performTask( const unsigned& task_index, const unsigned& current, MultiValue& myvals ) const ;
76 : ///
77 : void turnOnDerivatives();
78 : };
79 :
80 2525 : PLUMED_REGISTER_ACTION(ClusterDiameter,"CLUSTER_DIAMETER")
81 :
82 3 : void ClusterDiameter::registerKeywords( Keywords& keys ) {
83 3 : ClusterAnalysisBase::registerKeywords( keys );
84 3 : 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.");
85 3 : }
86 :
87 2 : ClusterDiameter::ClusterDiameter(const ActionOptions&ao):
88 : Action(ao),
89 2 : ClusterAnalysisBase(ao)
90 : {
91 : // Find out which cluster we want
92 2 : parse("CLUSTER",clustr);
93 :
94 2 : if( clustr<1 ) error("cannot look for a cluster larger than the largest cluster");
95 2 : if( clustr>getNumberOfNodes() ) error("cluster selected is invalid - too few atoms in system");
96 :
97 : // Create the task list
98 3994 : for(unsigned i=0; i<getNumberOfNodes(); ++i) {
99 3992 : for(unsigned j=0; j<getNumberOfNodes(); ++j) addTaskToList( i*getNumberOfNodes() + j );
100 : }
101 : // Now create a higest vessel
102 2 : addVessel("HIGHEST", "", -1); std::vector<AtomNumber> fake_atoms; setupMultiColvarBase( fake_atoms );
103 2 : }
104 :
105 0 : void ClusterDiameter::turnOnDerivatives() {
106 0 : error("cannot calculate derivatives of cluster radius. This quantity is not differentiable");
107 0 : }
108 :
109 2 : void ClusterDiameter::calculate() {
110 : // Retrieve the atoms in the largest cluster
111 2 : std::vector<unsigned> myatoms; retrieveAtomsInCluster( clustr, myatoms );
112 : // Activate the relevant tasks
113 2 : deactivateAllTasks();
114 128 : for(unsigned i=1; i<myatoms.size(); ++i) {
115 126 : for(unsigned j=0; j<i; ++j) taskFlags[ myatoms[i]*getNumberOfNodes() + myatoms[j] ] = 1;
116 : }
117 2 : lockContributors();
118 : // Now do the calculation
119 2 : runAllTasks();
120 2 : }
121 :
122 2016 : void ClusterDiameter::performTask( const unsigned& task_index, const unsigned& current, MultiValue& myvals ) const {
123 2016 : unsigned iatom=std::floor(current/getNumberOfNodes()), jatom = current - iatom*getNumberOfNodes();
124 2016 : Vector distance=getSeparation( getPosition(iatom), getPosition(jatom) );
125 2016 : double dd = distance.modulo();
126 2016 : myvals.setValue( 0, 1.0 ); myvals.setValue( 1, dd );
127 2016 : }
128 :
129 : }
130 2523 : }
|
