Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2019 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 ARG=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 neighbouring atoms in a system containg 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 ARG=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 12 : 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 ;
80 : /// Actually do the calculation
81 : double compute( const unsigned& tindex, AtomValuePack& myatoms ) const ;
82 : /// Is the variable periodic
83 0 : bool isPeriodic() { return false; }
84 : };
85 :
86 6456 : PLUMED_REGISTER_ACTION(NumberOfLinks,"NLINKS")
87 :
88 5 : void NumberOfLinks::registerKeywords( Keywords& keys ) {
89 5 : MultiColvarBase::registerKeywords( keys );
90 20 : keys.add("atoms","GROUP","");
91 20 : keys.add("atoms-1","GROUPA","");
92 20 : keys.add("atoms-1","GROUPB","");
93 25 : keys.add("compulsory","NN","6","The n parameter of the switching function ");
94 25 : keys.add("compulsory","MM","0","The m parameter of the switching function; 0 implies 2*NN");
95 25 : keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
96 20 : keys.add("compulsory","R_0","The r_0 parameter of the switching function");
97 20 : keys.add("optional","SWITCH","This keyword is used if you want to employ an alternative to the continuous swiching function defined above. "
98 : "The following provides information on the \\ref switchingfunction that are available. "
99 : "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
100 : // Use actionWithDistributionKeywords
101 10 : keys.remove("LOWMEM");
102 15 : keys.addFlag("LOWMEM",false,"lower the memory requirements");
103 5 : }
104 :
105 4 : NumberOfLinks::NumberOfLinks(const ActionOptions& ao):
106 : Action(ao),
107 8 : MultiColvarBase(ao)
108 : {
109 : // The weight of this does have derivatives
110 4 : weightHasDerivatives=true;
111 :
112 : // Read in the switching function
113 8 : std::string sw, errors; parse("SWITCH",sw);
114 4 : if(sw.length()>0) {
115 4 : switchingFunction.set(sw,errors);
116 : } else {
117 0 : double r_0=-1.0, d_0; int nn, mm;
118 0 : parse("NN",nn); parse("MM",mm);
119 0 : parse("R_0",r_0); parse("D_0",d_0);
120 0 : if( r_0<0.0 ) error("you must set a value for R_0");
121 0 : switchingFunction.set(nn,mm,r_0,d_0);
122 : }
123 12 : log.printf(" calculating number of links with atoms separation of %s\n",( switchingFunction.description() ).c_str() );
124 16 : std::vector<AtomNumber> all_atoms; readTwoGroups( "GROUP", "GROUPA", "GROUPB", all_atoms );
125 4 : setupMultiColvarBase( all_atoms ); setLinkCellCutoff( switchingFunction.get_dmax() );
126 :
127 12 : for(unsigned i=0; i<getNumberOfBaseMultiColvars(); ++i) {
128 4 : if( !getBaseMultiColvar(i)->hasDifferentiableOrientation() ) error("cannot use multicolvar of type " + getBaseMultiColvar(i)->getName() );
129 : }
130 :
131 : // Create holders for the collective variable
132 4 : readVesselKeywords();
133 4 : plumed_assert( getNumberOfVessels()==0 );
134 8 : std::string input; addVessel( "SUM", input, -1 );
135 4 : readVesselKeywords();
136 4 : }
137 :
138 8064 : double NumberOfLinks::calculateWeight( const unsigned& taskCode, const double& weight, AtomValuePack& myatoms ) const {
139 16128 : Vector distance = getSeparation( myatoms.getPosition(0), myatoms.getPosition(1) );
140 8064 : double dfunc, sw = switchingFunction.calculateSqr( distance.modulo2(), dfunc );
141 :
142 8064 : if( !doNotCalculateDerivatives() ) {
143 8064 : addAtomDerivatives( 0, 0, (-dfunc)*weight*distance, myatoms );
144 8064 : addAtomDerivatives( 0, 1, (dfunc)*weight*distance, myatoms );
145 8064 : myatoms.addBoxDerivatives( 0, (-dfunc)*weight*Tensor(distance,distance) );
146 : }
147 8064 : return sw;
148 : }
149 :
150 8064 : double NumberOfLinks::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
151 8064 : if( getBaseMultiColvar(0)->getNumberOfQuantities()<3 ) return 1.0;
152 :
153 8064 : unsigned ncomp=getBaseMultiColvar(0)->getNumberOfQuantities();
154 :
155 8064 : std::vector<double> orient0( ncomp ), orient1( ncomp );
156 8064 : getInputData( 0, true, myatoms, orient0 );
157 8064 : getInputData( 1, true, myatoms, orient1 );
158 :
159 : double dot=0;
160 548352 : for(unsigned k=2; k<orient0.size(); ++k) {
161 354816 : dot+=orient0[k]*orient1[k];
162 : }
163 :
164 8064 : if( !doNotCalculateDerivatives() ) {
165 8064 : MultiValue& myder0=getInputDerivatives( 0, true, myatoms );
166 16128 : mergeInputDerivatives( 1, 2, orient1.size(), 0, orient1, myder0, myatoms );
167 8064 : MultiValue& myder1=getInputDerivatives( 1, true, myatoms );
168 16128 : mergeInputDerivatives( 1, 2, orient0.size(), 1, orient0, myder1, myatoms );
169 : }
170 :
171 : return dot;
172 : }
173 :
174 : }
175 4839 : }
|
