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 "AdjacencyMatrixVessel.h"
24 : #include "core/ActionRegister.h"
25 : #include "tools/SwitchingFunction.h"
26 :
27 : //+PLUMEDOC CONCOMP CLUSTER_DISTRIBUTION
28 : /*
29 : Calculate functions of the distribution of properties in your connected components.
30 :
31 : This collective variable was developed for looking at nucleation phenomena, where you are
32 : interested in using studying the behavior of atoms in small aggregates or nuclei. In these sorts of
33 : problems you might be interested in the distribution of the sizes of the clusters in your system.
34 : A detailed description of this CV can be found in \cite tribello-clustering.
35 :
36 : \par Examples
37 :
38 : The input provided below calculates the local q6 Steinhardt parameter on each atom. The coordination number
39 : that atoms with a high value for the local q6 Steinhardt parameter have with other atoms that have a high
40 : value for the local q6 Steinhardt parameter is then computed. A contact matrix is then computed that measures
41 : whether atoms atoms \f$i\f$ and \f$j\f$ have a high value for this coordination number and if they are within
42 : 3.6 nm of each other. The connected components of this matrix are then found using a depth first clustering
43 : algorithm on the corresponding graph. The number of componets in this graph that contain more than 27 atoms is then computed.
44 : As discussed in \cite tribello-clustering this input was used to analyse the formation of a polycrystal of GeTe from amorphous
45 : GeTe.
46 :
47 : \verbatim
48 : q6: Q6 SPECIES=1-32768 SWITCH={GAUSSIAN D_0=5.29 R_0=0.01 D_MAX=5.3} LOWMEM
49 : lq6: LOCAL_Q6 SPECIES=q6 SWITCH={GAUSSIAN D_0=5.29 R_0=0.01 D_MAX=5.3} LOWMEM
50 : flq6: MFILTER_MORE DATA=lq6 SWITCH={GAUSSIAN D_0=0.19 R_0=0.01 D_MAX=0.2}
51 : cc: COORDINATIONNUMBER SPECIES=flq6 SWITCH={GAUSSIAN D_0=3.59 R_0=0.01 D_MAX=3.6}
52 : fcc: MFILTER_MORE DATA=cc SWITCH={GAUSSIAN D_0=5.99 R_0=0.01 D_MAX=6.0}
53 : mat: CONTACT_MATRIX ATOMS=fcc SWITCH={GAUSSIAN D_0=3.59 R_0=0.01 D_MAX=3.6}
54 : dfs: DFSCLUSTERING MATRIX=mat
55 : nclust: CLUSTER_DISTRIBUTION CLUSTERS=dfs TRANSFORM={GAUSSIAN D_0=5.99 R_0=0.01 D_MAX=6.0} MORE_THAN={GAUSSIAN D_0=26.99 R_0=0.01 D_MAX=27}
56 : PRINT ARG=nclust.* FILE=colvar
57 : \endverbatim
58 :
59 : */
60 : //+ENDPLUMEDOC
61 :
62 : namespace PLMD {
63 : namespace adjmat {
64 :
65 2 : class ClusterDistribution : public ClusterAnalysisBase {
66 : private:
67 : unsigned nderivatives;
68 : ///
69 : bool use_switch, inverse;
70 : //
71 : SwitchingFunction sf;
72 : public:
73 : /// Create manual
74 : static void registerKeywords( Keywords& keys );
75 : /// Constructor
76 : explicit ClusterDistribution(const ActionOptions&);
77 : /// Do the calculation
78 : void calculate();
79 : /// We can use ActionWithVessel to run all the calculation
80 : void performTask( const unsigned&, const unsigned&, MultiValue& ) const ;
81 : };
82 :
83 2524 : PLUMED_REGISTER_ACTION(ClusterDistribution,"CLUSTER_DISTRIBUTION")
84 :
85 2 : void ClusterDistribution::registerKeywords( Keywords& keys ) {
86 2 : ClusterAnalysisBase::registerKeywords( keys );
87 2 : keys.add("compulsory","TRANSFORM","none","the switching function to use to convert the crystallinity parameter to a number between zero and one");
88 : keys.addFlag("INVERSE_TRANSFORM",false,"when TRANSFORM appears alone the input symmetry functions, \\f$x\\f$ are transformed used \\f$1-s(x)\\f$ "
89 2 : "where \\f$s(x)\\f$ is a switching function. When this option is used you instead transform using \\f$s(x)\\f$ only.");
90 2 : keys.use("MORE_THAN"); keys.use("LESS_THAN"); keys.use("BETWEEN");
91 2 : keys.use("HISTOGRAM"); keys.use("ALT_MIN"); keys.use("MIN"); keys.use("MAX");
92 2 : }
93 :
94 1 : ClusterDistribution::ClusterDistribution(const ActionOptions&ao):
95 : Action(ao),
96 : ClusterAnalysisBase(ao),
97 1 : nderivatives(0)
98 : {
99 1 : use_switch=false;
100 2 : std::string input, errors; parse("TRANSFORM",input);
101 1 : if( input!="none" ) {
102 0 : use_switch=true; sf.set( input, errors );
103 0 : if( errors.length()!=0 ) error("problem reading SWITCH keyword : " + errors );
104 : }
105 1 : parseFlag("INVERSE_TRANSFORM",inverse);
106 1 : if( inverse && !use_switch ) error("INVERSE_TRANSFORM option was specified but no TRANSOFRM switching function was given");
107 :
108 : // Create all tasks by copying those from underlying DFS object (which is actually MultiColvar)
109 1 : for(unsigned i=0; i<getNumberOfNodes(); ++i) addTaskToList(i);
110 :
111 : // And now finish the setup of everything in the base
112 2 : std::vector<AtomNumber> fake_atoms; setupMultiColvarBase( fake_atoms );
113 1 : }
114 :
115 1 : void ClusterDistribution::calculate() {
116 : // Activate the relevant tasks
117 1 : nderivatives = getNumberOfDerivatives();
118 1 : deactivateAllTasks();
119 1 : for(unsigned i=0; i<getNumberOfClusters(); ++i) taskFlags[i]=1;
120 1 : lockContributors();
121 : // Now do the calculation
122 1 : runAllTasks();
123 1 : }
124 :
125 3 : void ClusterDistribution::performTask( const unsigned& task_index, const unsigned& current, MultiValue& myvals ) const {
126 3 : std::vector<unsigned> myatoms; retrieveAtomsInCluster( current+1, myatoms );
127 : // This deals with filters
128 6 : if( myatoms.size()==1 && !nodeIsActive(myatoms[0]) ) return ;
129 :
130 6 : std::vector<double> vals( getNumberOfQuantities() );
131 6 : MultiValue tvals( getNumberOfQuantities(), nderivatives );
132 :
133 : // And this builds everything for this particular atom
134 3 : double vv, df, tval=0;
135 9 : for(unsigned j=0; j<myatoms.size(); ++j) {
136 6 : unsigned i=myatoms[j];
137 6 : getPropertiesOfNode( i, vals );
138 6 : if( use_switch && !inverse ) {
139 0 : vv = 1.0 - sf.calculate( vals[1], df );
140 0 : tval += vals[0]*vv; df=-df*vals[1];
141 6 : } else if( use_switch ) {
142 0 : vv = sf.calculate( vals[1], df );
143 0 : tval += vals[0]*vv; df=df*vals[1];
144 : } else {
145 6 : tval += vals[0]*vals[1]; df=1.; vv=vals[1];
146 : }
147 6 : if( !doNotCalculateDerivatives() ) {
148 6 : getNodePropertyDerivatives( i, tvals );
149 168 : for(unsigned k=0; k<tvals.getNumberActive(); ++k) {
150 162 : unsigned kat=tvals.getActiveIndex(k);
151 162 : myvals.addDerivative( 1, kat, vals[0]*df*tvals.getDerivative(1,kat) + vv*tvals.getDerivative(0,kat) );
152 : }
153 6 : tvals.clearAll();
154 : }
155 : }
156 6 : myvals.setValue( 0, 1.0 ); myvals.setValue( 1, tval );
157 : }
158 :
159 : }
160 2523 : }
|
