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