Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2017 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 "core/ActionWithVector.h"
23 : #include "core/ActionRegister.h"
24 : #include "tools/LeptonCall.h"
25 : #include "tools/Angle.h"
26 :
27 : namespace PLMD {
28 : namespace symfunc {
29 :
30 : //+PLUMEDOC COLVAR GSYMFUNC_THREEBODY
31 : /*
32 : Calculate functions of the coordinates of the coordinates of all pairs of bonds in the first coordination sphere of an atom
33 :
34 : This shortcut can be used to calculate [symmetry function](https://www.plumed-tutorials.org/lessons/23/001/data/SymmetryFunction.html) that
35 : are like those defined by Behler in the paper that is cited in the bibliography below. The particular symmetry functions that are computed
36 : by this shortcut are the angular ones that are functions of the set pairs of atoms in the coordination sphere of the central atom. One of
37 : the angular symmetry functions that Behler introduces is:
38 :
39 : $$
40 : G^5_i = 2^{1-\zeta} \sum_{j,k\ne i} (1 + \lambda\cos\theta_{ijk})^\zeta e^{-\nu(R_{ij}^2 + R_{ik}^2)} f_c(R_{ij}) f_c(R_{ik})
41 : $$
42 :
43 : In this expression $\zeta$, $\nu$ and $\lambda$ are all parameters. $f_c$ is a switching function which acts upon $R_{ij}$, the distance between atom $i$ and atom $j$, and
44 : $R_{ik}$, the distance between atom $i$ and atom $k$. $\theta_{ijk}$ is then the angle between the vector that points from atom $i$ to atom $j$ and the vector that points from
45 : atom $i$ to atom $k$. THe input below can be used to get PLUMED to calculate the 100 values for this symmetry function for the 100 atoms in a system.
46 :
47 : ```plumed
48 : # Calculate the contact matrix and the x,y and z components of the bond vectors
49 : # This action calculates 4 100x100 matrices
50 : cmat: CONTACT_MATRIX GROUP=1-100 SWITCH={CUSTOM R_0=4.5 D_MAX=4.5 FUNC=0.5*(cos(pi*x)+1)} COMPONENTS
51 :
52 : # Compute the symmetry function for the 100 atoms from the 4 100x100 matrices output
53 : # by cmat. The output from this action is a vector with 100 elements
54 : beh3: GSYMFUNC_THREEBODY ...
55 : WEIGHT=cmat.w ARG=cmat.x,cmat.y,cmat.z
56 : FUNCTION1={FUNC=0.25*exp(-0.1*(rij+rik))*(1+3*cos(ajik))^3 LABEL=g5}
57 : ...
58 :
59 : # Print the 100 symmetry function values to a file
60 : PRINT ARG=beh3.g5 FILE=colvar
61 : ```
62 :
63 : The GSYMFUNC_THREEBODY action sums over all the distinct triples of atoms that are identified in the contact matrix. This action uses the same functionality as [CUSTOM](CUSTOM.md) and can thus compute any
64 : function of the following four quantities:
65 :
66 : * `rij` - the distance between atom $i$ and atom $j$
67 : * `rik` - the distance between atom $i$ and atom $k$
68 : * `rjk` - the distance between atom $j$ and atom $k$
69 : * `ajik` - the angle between the vector connecting atom $i$ to atom $j$ and the vector connecting atom $i$ to atom $k$.
70 :
71 : Furthermore we can calculate more than one function of these four quantities at a time as illustrated by the input below:
72 :
73 : ```plumed
74 : # Calculate the contact matrix and the x,y and z components of the bond vectors
75 : # This action calculates 4 100x100 matrices
76 : cmat: CONTACT_MATRIX GROUP=1-100 SWITCH={CUSTOM R_0=4.5 D_MAX=4.5 FUNC=0.5*(cos(pi*x)+1)} COMPONENTS
77 :
78 : # Compute the 4 symmetry function below for the 100 atoms from the 4 100x100 matrices output
79 : # by cmat. The output from this action is a vector with 100 elements
80 : beh3: GSYMFUNC_THREEBODY ...
81 : WEIGHT=cmat.w ARG=cmat.x,cmat.y,cmat.z
82 : FUNCTION1={FUNC=0.25*(cos(pi*sqrt(rjk)/4.5)+1)*exp(-0.1*(rij+rik+rjk))*(1+2*cos(ajik))^2 LABEL=g4}
83 : FUNCTION2={FUNC=0.25*exp(-0.1*(rij+rik))*(1+3.5*cos(ajik))^3 LABEL=g5}
84 : FUNCTION3={FUNC=0.125*(1+6.6*cos(ajik))^4 LABEL=g6}
85 : FUNCTION4={FUNC=sin(3.0*(ajik-1)) LABEL=g7}
86 : ...
87 :
88 : # Print the 4 sets of 100 symmetry function values to a file
89 : PRINT ARG=beh3.g4,beh3.g5,beh3.g6,beh3.g7 FILE=colvar
90 : ```
91 :
92 : You can read more about how to calculate more Behler-type symmetry functions [here](https://www.plumed-tutorials.org/lessons/23/001/data/Behler.html).
93 :
94 : */
95 : //+ENDPLUMEDOC
96 :
97 : class ThreeBodyGFunctions : public ActionWithVector {
98 : private:
99 : std::vector<LeptonCall> functions;
100 : public:
101 : static void registerKeywords( Keywords& keys );
102 : explicit ThreeBodyGFunctions(const ActionOptions&);
103 : std::string getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const override ;
104 : void calculate() override ;
105 : unsigned getNumberOfDerivatives() override;
106 : void performTask( const unsigned& task_index, MultiValue& myvals ) const override ;
107 : };
108 :
109 : PLUMED_REGISTER_ACTION(ThreeBodyGFunctions,"GSYMFUNC_THREEBODY")
110 :
111 11 : void ThreeBodyGFunctions::registerKeywords( Keywords& keys ) {
112 11 : ActionWithVector::registerKeywords( keys );
113 22 : keys.addInputKeyword("compulsory","ARG","matrix","three matrices containing the bond vectors of interest");
114 22 : keys.addInputKeyword("compulsory","WEIGHT","matrix","the matrix that contains the weights that should be used for each connection");
115 11 : keys.add("numbered","FUNCTION","the parameters of the function you would like to compute");
116 11 : ActionWithValue::useCustomisableComponents( keys );
117 11 : keys.addDOI("10.1063/1.3553717");
118 11 : }
119 :
120 6 : ThreeBodyGFunctions::ThreeBodyGFunctions(const ActionOptions&ao):
121 : Action(ao),
122 6 : ActionWithVector(ao) {
123 6 : if( getNumberOfArguments()!=3 ) {
124 0 : error("found wrong number of arguments in input");
125 : }
126 : std::vector<Value*> wval;
127 12 : parseArgumentList("WEIGHT",wval);
128 6 : if( wval.size()!=1 ) {
129 0 : error("keyword WEIGHT should be provided with the label of a single action");
130 : }
131 :
132 24 : for(unsigned i=0; i<3; ++i) {
133 18 : if( getPntrToArgument(i)->getRank()!=2 ) {
134 0 : error("input argument should be a matrix");
135 : }
136 18 : if( wval[0]->getShape()[0]!=getPntrToArgument(i)->getShape()[0] || wval[0]->getShape()[1]!=getPntrToArgument(i)->getShape()[1] ) {
137 0 : error("mismatched shapes of matrices in input");
138 : }
139 : }
140 6 : log.printf(" using bond weights from matrix labelled %s \n",wval[0]->getName().c_str() );
141 : // Rerequest the arguments
142 6 : std::vector<Value*> myargs( getArguments() );
143 6 : myargs.push_back( wval[0] );
144 6 : requestArguments( myargs );
145 30 : for(unsigned i=0; i<myargs.size(); ++i) {
146 24 : myargs[i]->buildDataStore();
147 : }
148 6 : std::vector<unsigned> shape(1);
149 6 : shape[0] = getPntrToArgument(0)->getShape()[0];
150 :
151 : // And now read the functions to compute
152 6 : for(int i=1;; i++) {
153 : std::string myfunc, mystr, lab, num;
154 17 : Tools::convert(i,num);
155 34 : if( !parseNumbered("FUNCTION",i,mystr ) ) {
156 : break;
157 : }
158 11 : std::vector<std::string> data=Tools::getWords(mystr);
159 22 : if( !Tools::parse(data,"LABEL",lab ) ) {
160 0 : error("found no LABEL in FUNCTION" + num + " specification");
161 : }
162 11 : addComponent( lab, shape );
163 11 : componentIsNotPeriodic( lab );
164 22 : if( !Tools::parse(data,"FUNC",myfunc) ) {
165 0 : error("found no FUNC in FUNCTION" + num + " specification");
166 : }
167 11 : log.printf(" component labelled %s is computed using %s \n",lab.c_str(), myfunc.c_str() );
168 11 : functions.push_back( LeptonCall() );
169 11 : std::vector<std::string> argnames(1);
170 : argnames[0]="ajik";
171 11 : if( myfunc.find("rij")!=std::string::npos ) {
172 8 : argnames.push_back("rij");
173 : }
174 11 : if( myfunc.find("rik")!=std::string::npos ) {
175 4 : if( argnames.size()<2 ) {
176 0 : error("if you have a function of rik it must also be a function of rij -- email gareth.tribello@gmail.com if this is a problem");
177 : }
178 8 : argnames.push_back("rik");
179 : }
180 11 : if( myfunc.find("rjk")!=std::string::npos ) {
181 3 : if( argnames.size()<2 ) {
182 0 : error("if you have a function of rjk it must also be a function of rij and rik -- email gareth.tribello@gmail.com if this is a problem");
183 : }
184 6 : argnames.push_back("rjk");
185 : }
186 11 : functions[i-1].set( myfunc, argnames, this, true );
187 22 : }
188 6 : checkRead();
189 6 : }
190 :
191 2 : std::string ThreeBodyGFunctions::getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const {
192 3 : for(unsigned i=0; i<getNumberOfComponents(); ++i) {
193 6 : if( getConstPntrToComponent(i)->getName().find(cname)!=std::string::npos ) {
194 : std::string num;
195 2 : Tools::convert( i+1, num );
196 4 : return "the function defined by the FUNCTION" + num + " keyword";
197 : }
198 : }
199 0 : plumed_error();
200 : return "";
201 : }
202 :
203 36 : unsigned ThreeBodyGFunctions::getNumberOfDerivatives() {
204 36 : return 0;
205 : }
206 :
207 19 : void ThreeBodyGFunctions::calculate() {
208 19 : runAllTasks();
209 19 : }
210 :
211 472 : void ThreeBodyGFunctions::performTask( const unsigned& task_index, MultiValue& myvals ) const {
212 : const Value* wval = getPntrToArgument(3);
213 : const Value* xval = getPntrToArgument(0);
214 : const Value* yval = getPntrToArgument(1);
215 : const Value* zval = getPntrToArgument(2);
216 : Angle angle;
217 472 : Vector disti, distj;
218 472 : unsigned matsize = wval->getNumberOfValues();
219 472 : std::vector<double> values(4);
220 472 : std::vector<Vector> der_i(4), der_j(4);
221 472 : unsigned nbonds = wval->getRowLength( task_index ), ncols = wval->getShape()[1];
222 13880 : for(unsigned i=0; i<nbonds; ++i) {
223 13408 : unsigned ipos = ncols*task_index + wval->getRowIndex( task_index, i );
224 13408 : double weighti = wval->get( ipos );
225 13408 : if( weighti<epsilon ) {
226 6584 : continue ;
227 : }
228 6824 : disti[0] = xval->get( ipos );
229 6824 : disti[1] = yval->get( ipos );
230 6824 : disti[2] = zval->get( ipos );
231 6824 : values[1] = disti.modulo2();
232 6824 : der_i[1]=2*disti;
233 6824 : der_i[2].zero();
234 189327 : for(unsigned j=0; j<i; ++j) {
235 182503 : unsigned jpos = ncols*task_index + wval->getRowIndex( task_index, j );
236 182503 : double weightj = wval->get( jpos );
237 182503 : if( weightj<epsilon ) {
238 45494 : continue ;
239 : }
240 137009 : distj[0] = xval->get( jpos );
241 137009 : distj[1] = yval->get( jpos );
242 137009 : distj[2] = zval->get( jpos );
243 137009 : values[2] = distj.modulo2();
244 137009 : der_j[1].zero();
245 137009 : der_j[2]=2*distj;
246 137009 : der_i[3] = ( disti - distj );
247 137009 : values[3] = der_i[3].modulo2();
248 137009 : der_i[3] = 2*der_i[3];
249 137009 : der_j[3] = -der_i[3];
250 : // Compute angle between bonds
251 137009 : values[0] = angle.compute( disti, distj, der_i[0], der_j[0] );
252 : // Compute product of weights
253 137009 : double weightij = weighti*weightj;
254 : // Now compute all symmetry functions
255 414532 : for(unsigned n=0; n<functions.size(); ++n) {
256 277523 : unsigned ostrn = getConstPntrToComponent(n)->getPositionInStream();
257 277523 : double nonweight = functions[n].evaluate( values );
258 277523 : myvals.addValue( ostrn, nonweight*weightij );
259 277523 : if( doNotCalculateDerivatives() ) {
260 6843 : continue;
261 : }
262 :
263 567230 : for(unsigned m=0; m<functions[n].getNumberOfArguments(); ++m) {
264 296550 : double der = weightij*functions[n].evaluateDeriv( m, values );
265 296550 : myvals.addDerivative( ostrn, ipos, der*der_i[m][0] );
266 296550 : myvals.addDerivative( ostrn, matsize+ipos, der*der_i[m][1] );
267 296550 : myvals.addDerivative( ostrn, 2*matsize+ipos, der*der_i[m][2] );
268 296550 : myvals.addDerivative( ostrn, jpos, der*der_j[m][0] );
269 296550 : myvals.addDerivative( ostrn, matsize+jpos, der*der_j[m][1] );
270 296550 : myvals.addDerivative( ostrn, 2*matsize+jpos, der*der_j[m][2] );
271 : }
272 270680 : myvals.addDerivative( ostrn, 3*matsize+ipos, nonweight*weightj );
273 270680 : myvals.addDerivative( ostrn, 3*matsize+jpos, nonweight*weighti );
274 : }
275 : }
276 : }
277 472 : if( doNotCalculateDerivatives() ) {
278 : return ;
279 : }
280 :
281 : // And update the elements that have derivatives
282 : // Needs a separate loop here as there may be forces from j
283 8320 : for(unsigned i=0; i<nbonds; ++i) {
284 8192 : unsigned ipos = ncols*task_index + wval->getRowIndex( task_index, i );
285 8192 : double weighti = wval->get( ipos );
286 8192 : if( weighti<epsilon ) {
287 3322 : continue ;
288 : }
289 :
290 16286 : for(unsigned n=0; n<functions.size(); ++n) {
291 11416 : unsigned ostrn = getConstPntrToComponent(n)->getPositionInStream();
292 11416 : myvals.updateIndex( ostrn, ipos );
293 11416 : myvals.updateIndex( ostrn, matsize+ipos );
294 11416 : myvals.updateIndex( ostrn, 2*matsize+ipos );
295 11416 : myvals.updateIndex( ostrn, 3*matsize+ipos );
296 : }
297 : }
298 : }
299 :
300 : }
301 : }
|
