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/ParallelTaskManager.h"
24 : #include "core/ActionRegister.h"
25 : #include "tools/LeptonCall.h"
26 : #include "tools/Angle.h"
27 :
28 : namespace PLMD {
29 : namespace symfunc {
30 :
31 : //+PLUMEDOC COLVAR GSYMFUNC_THREEBODY
32 : /*
33 : Calculate functions of the coordinates of the coordinates of all pairs of bonds in the first coordination sphere of an atom
34 :
35 : This shortcut can be used to calculate [symmetry function](https://www.plumed-tutorials.org/lessons/23/001/data/SymmetryFunction.html) that
36 : are like those defined by Behler in the paper that is cited in the bibliography below. The particular symmetry functions that are computed
37 : 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
38 : the angular symmetry functions that Behler introduces is:
39 :
40 : $$
41 : 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})
42 : $$
43 :
44 : 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
45 : $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
46 : 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.
47 :
48 : ```plumed
49 : # Calculate the contact matrix and the x,y and z components of the bond vectors
50 : # This action calculates 4 100x100 matrices
51 : cmat: CONTACT_MATRIX GROUP=1-100 SWITCH={CUSTOM R_0=4.5 D_MAX=4.5 FUNC=0.5*(cos(pi*x)+1)} COMPONENTS
52 :
53 : # Compute the symmetry function for the 100 atoms from the 4 100x100 matrices output
54 : # by cmat. The output from this action is a vector with 100 elements
55 : beh3: GSYMFUNC_THREEBODY ...
56 : WEIGHT=cmat.w ARG=cmat.x,cmat.y,cmat.z
57 : FUNCTION1={FUNC=0.25*exp(-0.1*(rij+rik))*(1+3*cos(ajik))^3 LABEL=g5}
58 : ...
59 :
60 : # Print the 100 symmetry function values to a file
61 : PRINT ARG=beh3.g5 FILE=colvar
62 : ```
63 :
64 : 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
65 : function of the following four quantities:
66 :
67 : * `rij` - the distance between atom $i$ and atom $j$
68 : * `rik` - the distance between atom $i$ and atom $k$
69 : * `rjk` - the distance between atom $j$ and atom $k$
70 : * `ajik` - the angle between the vector connecting atom $i$ to atom $j$ and the vector connecting atom $i$ to atom $k$.
71 :
72 : Furthermore we can calculate more than one function of these four quantities at a time as illustrated by the input below:
73 :
74 : ```plumed
75 : # Calculate the contact matrix and the x,y and z components of the bond vectors
76 : # This action calculates 4 100x100 matrices
77 : cmat: CONTACT_MATRIX GROUP=1-100 SWITCH={CUSTOM R_0=4.5 D_MAX=4.5 FUNC=0.5*(cos(pi*x)+1)} COMPONENTS
78 :
79 : # Compute the 4 symmetry function below for the 100 atoms from the 4 100x100 matrices output
80 : # by cmat. The output from this action is a vector with 100 elements
81 : beh3: GSYMFUNC_THREEBODY ...
82 : WEIGHT=cmat.w ARG=cmat.x,cmat.y,cmat.z
83 : FUNCTION1={FUNC=0.25*(cos(pi*sqrt(rjk)/4.5)+1)*exp(-0.1*(rij+rik+rjk))*(1+2*cos(ajik))^2 LABEL=g4}
84 : FUNCTION2={FUNC=0.25*exp(-0.1*(rij+rik))*(1+3.5*cos(ajik))^3 LABEL=g5}
85 : FUNCTION3={FUNC=0.125*(1+6.6*cos(ajik))^4 LABEL=g6}
86 : FUNCTION4={FUNC=sin(3.0*(ajik-1)) LABEL=g7}
87 : ...
88 :
89 : # Print the 4 sets of 100 symmetry function values to a file
90 : PRINT ARG=beh3.g4,beh3.g5,beh3.g6,beh3.g7 FILE=colvar
91 : ```
92 :
93 : You can even use this action in tandem with the features that are in the [volumes module](module_volumes.md) as shown below:
94 :
95 : ```plumed
96 : # The atoms that are of interest
97 : ow: GROUP ATOMS=1-16500
98 : # Fixed virtual atom which serves as the probe volume's center (pos. in nm)
99 : center: FIXEDATOM AT=2.5,2.5,2.5
100 : # Vector in which element i is one if atom i is in sphere of interest and zero otherwise
101 : sphere: INSPHERE ATOMS=ow CENTER=center RADIUS={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52}
102 : # The distance matrix
103 : dmap: DISTANCE_MATRIX COMPONENTS GROUP=ow CUTOFF=1.0 MASK=sphere
104 : # Find the four nearest neighbors
105 : acv_neigh: NEIGHBORS ARG=dmap.w NLOWEST=4 MASK=sphere
106 : # Compute a function for the atoms that are in the first coordination sphere
107 : acv_g8: GSYMFUNC_THREEBODY ...
108 : WEIGHT=acv_neigh ARG=dmap.x,dmap.y,dmap.z
109 : FUNCTION1={FUNC=(cos(ajik)+1/3)^2 LABEL=g8}
110 : MASK=sphere
111 : ...
112 : # Now compute the value of the function above for those atoms that are in the
113 : # sphere of interest
114 : acv: CUSTOM ARG=acv_g8.g8,sphere FUNC=y*(1-(3*x/8)) PERIODIC=NO
115 : # And now compute the final average
116 : acv_sum: SUM ARG=acv PERIODIC=NO
117 : acv_norm: SUM ARG=sphere PERIODIC=NO
118 : mean: CUSTOM ARG=acv_sum,acv_norm FUNC=x/y PERIODIC=NO
119 : PRINT ARG=mean FILE=colvar
120 : ```
121 :
122 : 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).
123 :
124 : */
125 : //+ENDPLUMEDOC
126 :
127 6 : class ThreeBodyGFunctionsInput {
128 : public:
129 : bool multi_action_input;
130 : std::vector<std::string> funcstr;
131 : std::vector<LeptonCall> functions;
132 6 : ThreeBodyGFunctionsInput& operator=( const ThreeBodyGFunctionsInput& m ) {
133 6 : multi_action_input = m.multi_action_input;
134 17 : for(unsigned i=0; i<m.funcstr.size(); ++i) {
135 22 : addFunction( m.funcstr[i], NULL );
136 : }
137 6 : return *this;
138 : }
139 22 : void addFunction( std::string myfunc, ActionWithVector* action ) {
140 22 : funcstr.push_back( myfunc );
141 22 : functions.push_back( LeptonCall() );
142 22 : std::vector<std::string> argnames(1);
143 : argnames[0]="ajik";
144 22 : if( myfunc.find("rij")!=std::string::npos ) {
145 16 : argnames.push_back("rij");
146 : }
147 22 : if( myfunc.find("rik")!=std::string::npos ) {
148 8 : if( action && argnames.size()<2 ) {
149 0 : action->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");
150 : }
151 16 : argnames.push_back("rik");
152 : }
153 22 : if( myfunc.find("rjk")!=std::string::npos ) {
154 6 : if( action && argnames.size()<2 ) {
155 0 : action->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");
156 : }
157 12 : argnames.push_back("rjk");
158 : }
159 22 : functions[functions.size()-1].set( myfunc, argnames, action, true );
160 22 : }
161 : };
162 :
163 : class ThreeBodyGFunctions : public ActionWithVector {
164 : public:
165 : using input_type = ThreeBodyGFunctionsInput;
166 : using PTM = ParallelTaskManager<ThreeBodyGFunctions>;
167 : private:
168 : PTM taskmanager;
169 : public:
170 : static void registerKeywords( Keywords& keys );
171 : explicit ThreeBodyGFunctions(const ActionOptions&);
172 : std::string getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const override ;
173 : void calculate() override ;
174 : void applyNonZeroRankForces( std::vector<double>& outforces ) override ;
175 : unsigned getNumberOfDerivatives() override;
176 : static std::size_t getIndex( std::size_t irow, std::size_t jcol, const ArgumentBookeepingHolder& mat );
177 : static void performTask( std::size_t task_index, const ThreeBodyGFunctionsInput& actiondata, ParallelActionsInput& input, ParallelActionsOutput& output );
178 : static int getNumberOfValuesPerTask( std::size_t task_index, const ThreeBodyGFunctionsInput& actiondata );
179 : static void getForceIndices( std::size_t task_index, std::size_t colno, std::size_t ntotal_force, const ThreeBodyGFunctionsInput& actiondata, const ParallelActionsInput& input, ForceIndexHolder force_indices );
180 : };
181 :
182 : PLUMED_REGISTER_ACTION(ThreeBodyGFunctions,"GSYMFUNC_THREEBODY")
183 :
184 11 : void ThreeBodyGFunctions::registerKeywords( Keywords& keys ) {
185 11 : ActionWithVector::registerKeywords( keys );
186 22 : keys.addInputKeyword("optional","MASK","vector","a vector that is used to used to determine which symmetry functions should be calculated");
187 22 : keys.addInputKeyword("compulsory","ARG","matrix","three matrices containing the bond vectors of interest");
188 22 : keys.addInputKeyword("compulsory","WEIGHT","matrix","the matrix that contains the weights that should be used for each connection");
189 11 : keys.add("numbered","FUNCTION","the parameters of the function you would like to compute");
190 11 : PTM::registerKeywords( keys );
191 11 : ActionWithValue::useCustomisableComponents( keys );
192 11 : keys.addDOI("10.1063/1.3553717");
193 11 : }
194 :
195 6 : ThreeBodyGFunctions::ThreeBodyGFunctions(const ActionOptions&ao):
196 : Action(ao),
197 : ActionWithVector(ao),
198 6 : taskmanager(this) {
199 : unsigned nargs = getNumberOfArguments();
200 6 : if( getNumberOfMasks()>0 ) {
201 0 : nargs = nargs - getNumberOfMasks();
202 : }
203 6 : if( nargs!=3 ) {
204 0 : error("found wrong number of arguments in input");
205 : }
206 : std::vector<Value*> wval;
207 12 : parseArgumentList("WEIGHT",wval);
208 6 : if( wval.size()!=1 ) {
209 0 : error("keyword WEIGHT should be provided with the label of a single action");
210 : }
211 :
212 24 : for(unsigned i=0; i<3; ++i) {
213 18 : if( getPntrToArgument(i)->getRank()!=2 ) {
214 0 : error("input argument should be a matrix");
215 : }
216 18 : if( wval[0]->getShape()[0]!=getPntrToArgument(i)->getShape()[0] || wval[0]->getShape()[1]!=getPntrToArgument(i)->getShape()[1] ) {
217 0 : error("mismatched shapes of matrices in input");
218 : }
219 : }
220 6 : log.printf(" using bond weights from matrix labelled %s \n",wval[0]->getName().c_str() );
221 : // Rerequest the arguments
222 6 : std::vector<Value*> myargs( getArguments() );
223 6 : myargs.push_back( wval[0] );
224 6 : requestArguments( myargs );
225 6 : std::vector<std::size_t> shape(1);
226 6 : shape[0] = getPntrToArgument(0)->getShape()[0];
227 :
228 : // And now read the functions to compute
229 : ThreeBodyGFunctionsInput input;
230 6 : for(int i=1;; i++) {
231 : std::string myfunc, mystr, lab, num;
232 17 : Tools::convert(i,num);
233 34 : if( !parseNumbered("FUNCTION",i,mystr ) ) {
234 : break;
235 : }
236 11 : std::vector<std::string> data=Tools::getWords(mystr);
237 22 : if( !Tools::parse(data,"LABEL",lab ) ) {
238 0 : error("found no LABEL in FUNCTION" + num + " specification");
239 : }
240 11 : addComponent( lab, shape );
241 11 : componentIsNotPeriodic( lab );
242 22 : if( !Tools::parse(data,"FUNC",myfunc) ) {
243 0 : error("found no FUNC in FUNCTION" + num + " specification");
244 : }
245 11 : log.printf(" component labelled %s is computed using %s \n",lab.c_str(), myfunc.c_str() );
246 11 : input.addFunction( myfunc, this );
247 22 : }
248 6 : checkRead();
249 6 : input.multi_action_input = getPntrToArgument(3)->getPntrToAction()!=getPntrToArgument(0)->getPntrToAction();
250 6 : taskmanager.setupParallelTaskManager( 4*wval[0]->getShape()[1], 0 );
251 6 : taskmanager.setActionInput( input );
252 6 : }
253 :
254 2 : std::string ThreeBodyGFunctions::getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const {
255 3 : for(unsigned i=0; i<getNumberOfComponents(); ++i) {
256 6 : if( getConstPntrToComponent(i)->getName().find(cname)!=std::string::npos ) {
257 : std::string num;
258 2 : Tools::convert( i+1, num );
259 4 : return "the function defined by the FUNCTION" + num + " keyword";
260 : }
261 : }
262 0 : plumed_error();
263 : return "";
264 : }
265 :
266 36 : unsigned ThreeBodyGFunctions::getNumberOfDerivatives() {
267 36 : return 0;
268 : }
269 :
270 19 : void ThreeBodyGFunctions::calculate() {
271 19 : taskmanager.runAllTasks();
272 19 : }
273 :
274 3072 : std::size_t ThreeBodyGFunctions::getIndex( std::size_t irow, std::size_t jcol, const ArgumentBookeepingHolder& mat ) {
275 3072 : if( mat.shape[1]==mat.ncols ) {
276 0 : return irow*mat.ncols + jcol;
277 : }
278 99105 : for(unsigned i=0; i<mat.bookeeping[(1+mat.ncols)*irow]; ++i) {
279 99105 : if( mat.bookeeping[(1+mat.ncols)*irow+1+i]==jcol ) {
280 3072 : return irow*mat.ncols+i;
281 : }
282 : }
283 0 : plumed_merror("could not find index");
284 : return 0;
285 : }
286 :
287 600 : void ThreeBodyGFunctions::performTask( std::size_t task_index, const ThreeBodyGFunctionsInput& actiondata, ParallelActionsInput& input, ParallelActionsOutput& output ) {
288 : const double* xpntr=NULL;
289 : const double* ypntr=NULL;
290 : const double* zpntr=NULL;
291 : /// This function uses lepton so it is unlikely that we will GPU it. that is why I have allowed myself to create vectors here
292 : std::vector<double> xvals, yvals, zvals;
293 600 : auto arg0 = ArgumentBookeepingHolder::create( 0, input );
294 600 : auto arg1 = ArgumentBookeepingHolder::create( 1, input );
295 600 : auto arg2 = ArgumentBookeepingHolder::create( 2, input);
296 600 : auto arg3 = ArgumentBookeepingHolder::create( 3, input );
297 :
298 600 : std::size_t rowlen = arg3.bookeeping[(1+arg3.ncols)*task_index];
299 600 : View<const double> wval( input.inputdata + arg3.start + arg3.ncols*task_index, rowlen );
300 600 : if( actiondata.multi_action_input ) {
301 256 : xvals.resize( rowlen );
302 256 : yvals.resize( rowlen );
303 256 : zvals.resize( rowlen );
304 256 : const auto wbooks = arg3.bookeeping.subview((1+arg3.ncols)*task_index+1, rowlen);
305 1280 : for(unsigned i=0; i<rowlen; ++i) {
306 1024 : xvals[i] = input.inputdata[arg0.start + getIndex( task_index, wbooks[i], arg0) ];
307 1024 : yvals[i] = input.inputdata[arg1.start + getIndex( task_index, wbooks[i], arg1) ];
308 1024 : zvals[i] = input.inputdata[arg2.start + getIndex( task_index, wbooks[i], arg2) ];
309 : }
310 : xpntr = xvals.data();
311 : ypntr = yvals.data();
312 : zpntr = zvals.data();
313 : } else {
314 344 : xpntr = input.inputdata + arg0.start + arg0.ncols*task_index;
315 344 : ypntr = input.inputdata + arg1.start + arg1.ncols*task_index;
316 344 : zpntr = input.inputdata + arg2.start + arg2.ncols*task_index;
317 : }
318 : View<const double> xval( xpntr, rowlen );
319 : View<const double> yval( ypntr, rowlen );
320 : View<const double> zval( zpntr, rowlen );
321 279512 : for(unsigned i=0; i<output.derivatives.size(); ++i) {
322 278912 : output.derivatives[i] = 0;
323 : }
324 600 : View2D<double> derivatives( output.derivatives.data(), actiondata.functions.size(), 4*arg3.shape[1] );
325 :
326 : Angle angle;
327 : Vector disti, distj;
328 600 : std::vector<double> values(4);
329 600 : std::vector<Vector> der_i(4), der_j(4);
330 21878 : for(unsigned i=0; i<rowlen; ++i) {
331 21278 : if( wval[i]<epsilon ) {
332 9584 : continue;
333 : }
334 11694 : disti[0] = xval[i];
335 11694 : disti[1] = yval[i];
336 11694 : disti[2] = zval[i];
337 11694 : values[1] = disti.modulo2();
338 11694 : der_i[1]=2*disti;
339 : der_i[2].zero();
340 371037 : for(unsigned j=0; j<i; ++j) {
341 359343 : if( wval[j]<epsilon) {
342 92168 : continue;
343 : }
344 267175 : distj[0] = xval[j];
345 267175 : distj[1] = yval[j];
346 267175 : distj[2] = zval[j];
347 267175 : values[2] = distj.modulo2();
348 : der_j[1].zero();
349 267175 : der_j[2]=2*distj;
350 267175 : der_i[3] = ( disti - distj );
351 267175 : values[3] = der_i[3].modulo2();
352 267175 : der_i[3] = 2*der_i[3];
353 267175 : der_j[3] = -der_i[3];
354 : // Compute angle between bonds
355 267175 : values[0] = angle.compute( disti, distj, der_i[0], der_j[0] );
356 : // Compute product of weights
357 267175 : double weightij = wval[i]*wval[j];
358 815378 : for(unsigned n=0; n<actiondata.functions.size(); ++n) {
359 548203 : double nonweight = actiondata.functions[n].evaluate( values );
360 548203 : output.values[n] += nonweight*weightij;
361 :
362 548203 : if( input.noderiv ) {
363 277523 : continue;
364 : }
365 :
366 567230 : for(unsigned m=0; m<actiondata.functions[n].getNumberOfArguments(); ++m) {
367 296550 : double der = weightij*actiondata.functions[n].evaluateDeriv( m, values );
368 296550 : derivatives[n][i] += der*der_i[m][0];
369 296550 : derivatives[n][rowlen+i] += der*der_i[m][1];
370 296550 : derivatives[n][2*rowlen+i] += der*der_i[m][2];
371 296550 : derivatives[n][j] += der*der_j[m][0];
372 296550 : derivatives[n][rowlen+j] += der*der_j[m][1];
373 296550 : derivatives[n][2*rowlen+j] += der*der_j[m][2];
374 : }
375 270680 : derivatives[n][3*rowlen+i] += nonweight*wval[j];
376 270680 : derivatives[n][3*rowlen+j] += nonweight*wval[i];
377 : }
378 : }
379 : }
380 600 : }
381 :
382 2 : void ThreeBodyGFunctions::applyNonZeroRankForces( std::vector<double>& outforces ) {
383 2 : taskmanager.applyForces( outforces );
384 2 : }
385 :
386 128 : int ThreeBodyGFunctions::getNumberOfValuesPerTask( std::size_t task_index, const ThreeBodyGFunctionsInput& actiondata ) {
387 128 : return 1;
388 : }
389 :
390 128 : void ThreeBodyGFunctions::getForceIndices( std::size_t task_index,
391 : std::size_t colno,
392 : std::size_t ntotal_force,
393 : const ThreeBodyGFunctionsInput& actiondata,
394 : const ParallelActionsInput& input,
395 : ForceIndexHolder force_indices ) {
396 128 : auto arg0 = ArgumentBookeepingHolder::create( 0, input );
397 128 : auto arg1 = ArgumentBookeepingHolder::create( 1, input );
398 128 : auto arg2 = ArgumentBookeepingHolder::create( 2, input);
399 128 : auto arg3 = ArgumentBookeepingHolder::create( 3, input );
400 128 : std::size_t rowlen = arg3.bookeeping[(1+arg3.ncols)*task_index];
401 128 : if( actiondata.multi_action_input ) {
402 0 : View<const std::size_t> wbooks( arg3.bookeeping.data()+(1+arg3.ncols)*task_index+1, rowlen);
403 0 : for(unsigned j=0; j<rowlen; ++j) {
404 0 : std::size_t matpos = task_index*arg3.ncols + j;
405 0 : std::size_t xpos = getIndex( task_index, wbooks[j], arg0 );
406 0 : std::size_t ypos = getIndex( task_index, wbooks[j], arg1 );
407 0 : std::size_t zpos = getIndex( task_index, wbooks[j], arg2 );
408 0 : for(unsigned i=0; i<input.ncomponents; ++i) {
409 0 : force_indices.indices[i][j] = arg0.start + xpos;
410 0 : force_indices.indices[i][rowlen+j] = arg1.start + ypos;
411 0 : force_indices.indices[i][2*rowlen+j] = arg2.start + zpos;
412 0 : force_indices.indices[i][3*rowlen+j] = arg3.start + matpos;
413 : }
414 : }
415 : } else {
416 8192 : for(unsigned j=0; j<rowlen; ++j) {
417 8064 : std::size_t matpos = task_index*arg3.ncols + j;
418 32256 : for(unsigned i=0; i<input.ncomponents; ++i) {
419 24192 : force_indices.indices[i][j] = arg0.start + matpos;
420 24192 : force_indices.indices[i][rowlen+j] = arg1.start + matpos;
421 24192 : force_indices.indices[i][2*rowlen+j] = arg2.start + matpos;
422 24192 : force_indices.indices[i][3*rowlen+j] = arg3.start + matpos;
423 : }
424 : }
425 : }
426 512 : for(unsigned i=0; i<input.ncomponents; ++i) {
427 384 : force_indices.threadsafe_derivatives_end[i] = force_indices.tot_indices[i] = 4*rowlen;
428 : }
429 128 : }
430 :
431 : }
432 : }
|
