Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-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 "MultiColvarBase.h"
23 : #include "AtomValuePack.h"
24 : #include "core/ActionRegister.h"
25 : #include "tools/SwitchingFunction.h"
26 :
27 : #include <string>
28 : #include <cmath>
29 :
30 : //+PLUMEDOC MCOLVARF NLINKS
31 : /*
32 : Calculate number of pairs of atoms/molecules that are linked
33 :
34 : In its simplest guise this coordinate calculates a coordination number. Each pair
35 : of atoms is assumed "linked" if they are within some cutoff of each other. In more
36 : complex applications each entity is a vector and this quantity measures whether
37 : pairs of vectors are (a) within a certain cutoff and (b) if the two vectors have
38 : similar orientations. The vectors on individual atoms could be Steinhardt parameters
39 : (see \ref Q3, \ref Q4 and \ref Q6) or they could describe some internal vector in a molecule.
40 :
41 : \par Examples
42 :
43 : The following calculates how many bonds there are in a system containing 64 atoms and outputs
44 : this quantity to a file.
45 :
46 : \plumedfile
47 : DENSITY SPECIES=1-64 LABEL=d1
48 : NLINKS GROUP=d1 SWITCH={RATIONAL D_0=1.3 R_0=0.2} LABEL=dd
49 : PRINT ARG=dd FILE=colvar
50 : \endplumedfile
51 :
52 : The following calculates how many pairs of neighboring atoms in a system containing 64 atoms have
53 : similar dispositions for the atoms in their coordination sphere. This calculation uses the
54 : dot product of the Q6 vectors on adjacent atoms to measure whether or not two atoms have the same
55 : ``orientation"
56 :
57 : \plumedfile
58 : Q6 SPECIES=1-64 SWITCH={RATIONAL D_0=1.3 R_0=0.2} LABEL=q6
59 : NLINKS GROUP=q6 SWITCH={RATIONAL D_0=1.3 R_0=0.2} LABEL=dd
60 : PRINT ARG=dd FILE=colvar
61 : \endplumedfile
62 :
63 : */
64 : //+ENDPLUMEDOC
65 :
66 : namespace PLMD {
67 : namespace multicolvar {
68 :
69 : class NumberOfLinks : public MultiColvarBase {
70 : private:
71 : /// The values of the quantities in the dot products
72 : std::vector<double> orient0, orient1;
73 : /// The switching function that tells us if atoms are close enough together
74 : SwitchingFunction switchingFunction;
75 : public:
76 : static void registerKeywords( Keywords& keys );
77 : explicit NumberOfLinks(const ActionOptions&);
78 : /// Do the stuff with the switching functions
79 : double calculateWeight( const unsigned& taskCode, const double& weight, AtomValuePack& myatoms ) const override;
80 : /// Actually do the calculation
81 : double compute( const unsigned& tindex, AtomValuePack& myatoms ) const override;
82 : /// Is the variable periodic
83 0 : bool isPeriodic() override {
84 0 : return false;
85 : }
86 : };
87 :
88 13793 : PLUMED_REGISTER_ACTION(NumberOfLinks,"NLINKS")
89 :
90 8 : void NumberOfLinks::registerKeywords( Keywords& keys ) {
91 8 : MultiColvarBase::registerKeywords( keys );
92 16 : keys.add("atoms","GROUP","");
93 16 : keys.add("atoms-1","GROUPA","");
94 16 : keys.add("atoms-1","GROUPB","");
95 16 : keys.add("compulsory","NN","6","The n parameter of the switching function ");
96 16 : keys.add("compulsory","MM","0","The m parameter of the switching function; 0 implies 2*NN");
97 16 : keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
98 16 : keys.add("compulsory","R_0","The r_0 parameter of the switching function");
99 16 : keys.add("optional","SWITCH","This keyword is used if you want to employ an alternative to the continuous switching function defined above. "
100 : "The following provides information on the \\ref switchingfunction that are available. "
101 : "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
102 : // Use actionWithDistributionKeywords
103 8 : keys.remove("LOWMEM");
104 16 : keys.addFlag("LOWMEM",false,"lower the memory requirements");
105 8 : }
106 :
107 4 : NumberOfLinks::NumberOfLinks(const ActionOptions& ao):
108 : Action(ao),
109 4 : MultiColvarBase(ao) {
110 : // The weight of this does have derivatives
111 4 : weightHasDerivatives=true;
112 :
113 : // Read in the switching function
114 : std::string sw, errors;
115 8 : parse("SWITCH",sw);
116 4 : if(sw.length()>0) {
117 4 : switchingFunction.set(sw,errors);
118 : } else {
119 0 : double r_0=-1.0, d_0;
120 : int nn, mm;
121 0 : parse("NN",nn);
122 0 : parse("MM",mm);
123 0 : parse("R_0",r_0);
124 0 : parse("D_0",d_0);
125 0 : if( r_0<0.0 ) {
126 0 : error("you must set a value for R_0");
127 : }
128 0 : switchingFunction.set(nn,mm,r_0,d_0);
129 : }
130 8 : log.printf(" calculating number of links with atoms separation of %s\n",( switchingFunction.description() ).c_str() );
131 : std::vector<AtomNumber> all_atoms;
132 8 : readTwoGroups( "GROUP", "GROUPA", "GROUPB", all_atoms );
133 4 : setupMultiColvarBase( all_atoms );
134 4 : setLinkCellCutoff( switchingFunction.get_dmax() );
135 :
136 8 : for(unsigned i=0; i<getNumberOfBaseMultiColvars(); ++i) {
137 4 : if( !getBaseMultiColvar(i)->hasDifferentiableOrientation() ) {
138 0 : error("cannot use multicolvar of type " + getBaseMultiColvar(i)->getName() );
139 : }
140 : }
141 :
142 : // Create holders for the collective variable
143 4 : readVesselKeywords();
144 4 : plumed_assert( getNumberOfVessels()==0 );
145 : std::string input;
146 4 : addVessel( "SUM", input, -1 );
147 4 : readVesselKeywords();
148 4 : }
149 :
150 8064 : double NumberOfLinks::calculateWeight( const unsigned& taskCode, const double& weight, AtomValuePack& myatoms ) const {
151 8064 : Vector distance = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
152 8064 : double dfunc, sw = switchingFunction.calculateSqr( distance.modulo2(), dfunc );
153 :
154 8064 : if( !doNotCalculateDerivatives() ) {
155 8064 : addAtomDerivatives( 0, 0, (-dfunc)*weight*distance, myatoms );
156 8064 : addAtomDerivatives( 0, 1, (dfunc)*weight*distance, myatoms );
157 8064 : myatoms.addBoxDerivatives( 0, (-dfunc)*weight*Tensor(distance,distance) );
158 : }
159 8064 : return sw;
160 : }
161 :
162 8064 : double NumberOfLinks::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
163 8064 : if( getBaseMultiColvar(0)->getNumberOfQuantities()<3 ) {
164 : return 1.0;
165 : }
166 :
167 8064 : unsigned ncomp=getBaseMultiColvar(0)->getNumberOfQuantities();
168 :
169 8064 : std::vector<double> orient0( ncomp ), orient1( ncomp );
170 8064 : getInputData( 0, true, myatoms, orient0 );
171 8064 : getInputData( 1, true, myatoms, orient1 );
172 :
173 : double dot=0;
174 185472 : for(unsigned k=2; k<orient0.size(); ++k) {
175 177408 : dot+=orient0[k]*orient1[k];
176 : }
177 :
178 8064 : if( !doNotCalculateDerivatives() ) {
179 8064 : MultiValue& myder0=getInputDerivatives( 0, true, myatoms );
180 8064 : mergeInputDerivatives( 1, 2, orient1.size(), 0, orient1, myder0, myatoms );
181 8064 : MultiValue& myder1=getInputDerivatives( 1, true, myatoms );
182 8064 : mergeInputDerivatives( 1, 2, orient0.size(), 1, orient0, myder1, myatoms );
183 : }
184 :
185 : return dot;
186 : }
187 :
188 : }
189 : }
|
