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 LINKCELL_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=false;
130 : std::vector<std::string> funcstr;
131 : std::vector<LeptonCall> functions;
132 : // the assignment is customized (see below): copy construction is the default one
133 6 : ThreeBodyGFunctionsInput() = default;
134 : ThreeBodyGFunctionsInput(const ThreeBodyGFunctionsInput&) = default;
135 6 : ThreeBodyGFunctionsInput& operator=( const ThreeBodyGFunctionsInput& m ) {
136 6 : multi_action_input = m.multi_action_input;
137 17 : for(unsigned i=0; i<m.funcstr.size(); ++i) {
138 22 : addFunction( m.funcstr[i], NULL );
139 : }
140 6 : return *this;
141 : }
142 22 : void addFunction( std::string myfunc, ActionWithVector* action ) {
143 22 : funcstr.push_back( myfunc );
144 44 : functions.push_back( LeptonCall() );
145 22 : std::vector<std::string> argnames(1);
146 : argnames[0]="ajik";
147 22 : if( myfunc.find("rij")!=std::string::npos ) {
148 16 : argnames.push_back("rij");
149 : }
150 22 : if( myfunc.find("rik")!=std::string::npos ) {
151 8 : if( action && argnames.size()<2 ) {
152 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");
153 : }
154 16 : argnames.push_back("rik");
155 : }
156 22 : if( myfunc.find("rjk")!=std::string::npos ) {
157 6 : if( action && argnames.size()<2 ) {
158 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");
159 : }
160 12 : argnames.push_back("rjk");
161 : }
162 22 : functions[functions.size()-1].set( myfunc, argnames, action, true );
163 22 : }
164 : };
165 :
166 : class ThreeBodyGFunctions : public ActionWithVector {
167 : public:
168 : using input_type = ThreeBodyGFunctionsInput;
169 : using PTM = ParallelTaskManager<ThreeBodyGFunctions>;
170 : typedef typename PTM::ParallelActionsInput ParallelActionsInput;
171 : typedef typename PTM::ParallelActionsOutput ParallelActionsOutput;
172 : private:
173 : PTM taskmanager;
174 : public:
175 : static void registerKeywords( Keywords& keys );
176 : explicit ThreeBodyGFunctions(const ActionOptions&);
177 : std::string getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const override ;
178 : void calculate() override ;
179 : void applyNonZeroRankForces( std::vector<double>& outforces ) override ;
180 : unsigned getNumberOfDerivatives() override;
181 : static std::size_t getIndex( std::size_t irow, std::size_t jcol, const ArgumentBookeepingHolder& mat );
182 : static void performTask( std::size_t task_index, const ThreeBodyGFunctionsInput& actiondata, ParallelActionsInput& input, ParallelActionsOutput& output );
183 : static int getNumberOfValuesPerTask( std::size_t task_index, const ThreeBodyGFunctionsInput& actiondata );
184 : 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 );
185 : };
186 :
187 : PLUMED_REGISTER_ACTION(ThreeBodyGFunctions,"GSYMFUNC_THREEBODY")
188 :
189 13 : void ThreeBodyGFunctions::registerKeywords( Keywords& keys ) {
190 13 : ActionWithVector::registerKeywords( keys );
191 26 : keys.addInputKeyword("optional","MASK","vector","a vector that is used to used to determine which symmetry functions should be calculated");
192 26 : keys.addInputKeyword("compulsory","ARG","matrix","three matrices containing the bond vectors of interest");
193 26 : keys.addInputKeyword("compulsory","WEIGHT","matrix","the matrix that contains the weights that should be used for each connection");
194 13 : keys.add("numbered","FUNCTION","the parameters of the function you would like to compute");
195 13 : PTM::registerKeywords( keys );
196 13 : ActionWithValue::useCustomisableComponents( keys );
197 13 : keys.addDOI("10.1063/1.3553717");
198 13 : }
199 :
200 6 : ThreeBodyGFunctions::ThreeBodyGFunctions(const ActionOptions&ao):
201 : Action(ao),
202 : ActionWithVector(ao),
203 6 : taskmanager(this) {
204 : unsigned nargs = getNumberOfArguments();
205 6 : if( getNumberOfMasks()>0 ) {
206 0 : nargs = nargs - getNumberOfMasks();
207 : }
208 6 : if( nargs!=3 ) {
209 0 : error("found wrong number of arguments in input");
210 : }
211 6 : std::vector<Value*> wval;
212 12 : parseArgumentList("WEIGHT",wval);
213 6 : if( wval.size()!=1 ) {
214 0 : error("keyword WEIGHT should be provided with the label of a single action");
215 : }
216 :
217 24 : for(unsigned i=0; i<3; ++i) {
218 18 : if( getPntrToArgument(i)->getRank()!=2 ) {
219 0 : error("input argument should be a matrix");
220 : }
221 18 : if( wval[0]->getShape()[0]!=getPntrToArgument(i)->getShape()[0] || wval[0]->getShape()[1]!=getPntrToArgument(i)->getShape()[1] ) {
222 0 : error("mismatched shapes of matrices in input");
223 : }
224 : }
225 6 : log.printf(" using bond weights from matrix labelled %s \n",wval[0]->getName().c_str() );
226 : // Rerequest the arguments
227 6 : std::vector<Value*> myargs( getArguments() );
228 6 : myargs.push_back( wval[0] );
229 6 : requestArguments( myargs );
230 6 : std::vector<std::size_t> shape(1);
231 6 : shape[0] = getPntrToArgument(0)->getShape()[0];
232 :
233 : // And now read the functions to compute
234 6 : ThreeBodyGFunctionsInput input;
235 6 : for(int i=1;; i++) {
236 : std::string myfunc, mystr, lab, num;
237 17 : Tools::convert(i,num);
238 34 : if( !parseNumbered("FUNCTION",i,mystr ) ) {
239 : break;
240 : }
241 11 : std::vector<std::string> data=Tools::getWords(mystr);
242 22 : if( !Tools::parse(data,"LABEL",lab ) ) {
243 0 : error("found no LABEL in FUNCTION" + num + " specification");
244 : }
245 11 : addComponent( lab, shape );
246 11 : componentIsNotPeriodic( lab );
247 22 : if( !Tools::parse(data,"FUNC",myfunc) ) {
248 0 : error("found no FUNC in FUNCTION" + num + " specification");
249 : }
250 11 : log.printf(" component labelled %s is computed using %s \n",lab.c_str(), myfunc.c_str() );
251 11 : input.addFunction( myfunc, this );
252 22 : }
253 6 : checkRead();
254 6 : input.multi_action_input = getPntrToArgument(3)->getPntrToAction()!=getPntrToArgument(0)->getPntrToAction();
255 6 : taskmanager.setupParallelTaskManager( 4*wval[0]->getShape()[1], 0 );
256 6 : taskmanager.setActionInput( input );
257 12 : }
258 :
259 2 : std::string ThreeBodyGFunctions::getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const {
260 3 : for(unsigned i=0; i<getNumberOfComponents(); ++i) {
261 6 : if( getConstPntrToComponent(i)->getName().find(cname)!=std::string::npos ) {
262 : std::string num;
263 2 : Tools::convert( i+1, num );
264 4 : return "the function defined by the FUNCTION" + num + " keyword";
265 : }
266 : }
267 0 : plumed_error();
268 : return "";
269 : }
270 :
271 36 : unsigned ThreeBodyGFunctions::getNumberOfDerivatives() {
272 36 : return 0;
273 : }
274 :
275 19 : void ThreeBodyGFunctions::calculate() {
276 19 : taskmanager.runAllTasks();
277 19 : }
278 :
279 3072 : std::size_t ThreeBodyGFunctions::getIndex( std::size_t irow, std::size_t jcol, const ArgumentBookeepingHolder& mat ) {
280 3072 : if( mat.shape[1]==mat.ncols ) {
281 0 : return irow*mat.ncols + jcol;
282 : }
283 99105 : for(unsigned i=0; i<mat.bookeeping[(1+mat.ncols)*irow]; ++i) {
284 99105 : if( mat.bookeeping[(1+mat.ncols)*irow+1+i]==jcol ) {
285 3072 : return irow*mat.ncols+i;
286 : }
287 : }
288 0 : plumed_merror("could not find index");
289 : return 0;
290 : }
291 :
292 600 : void ThreeBodyGFunctions::performTask( std::size_t task_index, const ThreeBodyGFunctionsInput& actiondata, ParallelActionsInput& input, ParallelActionsOutput& output ) {
293 : const double* xpntr=NULL;
294 : const double* ypntr=NULL;
295 : const double* zpntr=NULL;
296 : /// This function uses lepton so it is unlikely that we will GPU it. that is why I have allowed myself to create vectors here
297 600 : std::vector<double> xvals, yvals, zvals;
298 600 : auto arg0 = ArgumentBookeepingHolder::create( 0, input );
299 600 : auto arg1 = ArgumentBookeepingHolder::create( 1, input );
300 600 : auto arg2 = ArgumentBookeepingHolder::create( 2, input);
301 600 : auto arg3 = ArgumentBookeepingHolder::create( 3, input );
302 :
303 600 : std::size_t rowlen = arg3.bookeeping[(1+arg3.ncols)*task_index];
304 600 : View<const double> wval( input.inputdata + arg3.start + arg3.ncols*task_index, rowlen );
305 600 : if( actiondata.multi_action_input ) {
306 256 : xvals.resize( rowlen );
307 256 : yvals.resize( rowlen );
308 256 : zvals.resize( rowlen );
309 256 : const auto wbooks = arg3.bookeeping.subview((1+arg3.ncols)*task_index+1, rowlen);
310 1280 : for(unsigned i=0; i<rowlen; ++i) {
311 1024 : xvals[i] = input.inputdata[arg0.start + getIndex( task_index, wbooks[i], arg0) ];
312 1024 : yvals[i] = input.inputdata[arg1.start + getIndex( task_index, wbooks[i], arg1) ];
313 1024 : zvals[i] = input.inputdata[arg2.start + getIndex( task_index, wbooks[i], arg2) ];
314 : }
315 : xpntr = xvals.data();
316 : ypntr = yvals.data();
317 : zpntr = zvals.data();
318 : } else {
319 344 : xpntr = input.inputdata + arg0.start + arg0.ncols*task_index;
320 344 : ypntr = input.inputdata + arg1.start + arg1.ncols*task_index;
321 344 : zpntr = input.inputdata + arg2.start + arg2.ncols*task_index;
322 : }
323 : View<const double> xval( xpntr, rowlen );
324 : View<const double> yval( ypntr, rowlen );
325 : View<const double> zval( zpntr, rowlen );
326 279512 : for(unsigned i=0; i<output.derivatives.size(); ++i) {
327 278912 : output.derivatives[i] = 0;
328 : }
329 600 : View2D<double> derivatives( output.derivatives.data(), actiondata.functions.size(), 4*arg3.shape[1] );
330 :
331 : Angle angle;
332 : Vector disti, distj;
333 600 : std::vector<double> values(4);
334 600 : std::vector<Vector> der_i(4), der_j(4);
335 12294 : for(unsigned i=0; i<rowlen; ++i) {
336 11694 : if( wval[i]<epsilon ) {
337 0 : continue;
338 : }
339 11694 : disti[0] = xval[i];
340 11694 : disti[1] = yval[i];
341 11694 : disti[2] = zval[i];
342 11694 : values[1] = disti.modulo2();
343 11694 : der_i[1]=2*disti;
344 : der_i[2].zero();
345 278869 : for(unsigned j=0; j<i; ++j) {
346 267175 : if( wval[j]<epsilon) {
347 0 : continue;
348 : }
349 267175 : distj[0] = xval[j];
350 267175 : distj[1] = yval[j];
351 267175 : distj[2] = zval[j];
352 267175 : values[2] = distj.modulo2();
353 : der_j[1].zero();
354 267175 : der_j[2]=2*distj;
355 267175 : der_i[3] = ( disti - distj );
356 267175 : values[3] = der_i[3].modulo2();
357 267175 : der_i[3] = 2*der_i[3];
358 267175 : der_j[3] = -der_i[3];
359 : // Compute angle between bonds
360 267175 : values[0] = angle.compute( disti, distj, der_i[0], der_j[0] );
361 : // Compute product of weights
362 267175 : double weightij = wval[i]*wval[j];
363 815378 : for(unsigned n=0; n<actiondata.functions.size(); ++n) {
364 548203 : double nonweight = actiondata.functions[n].evaluate( values );
365 548203 : output.values[n] += nonweight*weightij;
366 :
367 548203 : if( input.noderiv ) {
368 277523 : continue;
369 : }
370 :
371 567230 : for(unsigned m=0; m<actiondata.functions[n].getNumberOfArguments(); ++m) {
372 296550 : double der = weightij*actiondata.functions[n].evaluateDeriv( m, values );
373 296550 : derivatives[n][i] += der*der_i[m][0];
374 296550 : derivatives[n][rowlen+i] += der*der_i[m][1];
375 296550 : derivatives[n][2*rowlen+i] += der*der_i[m][2];
376 296550 : derivatives[n][j] += der*der_j[m][0];
377 296550 : derivatives[n][rowlen+j] += der*der_j[m][1];
378 296550 : derivatives[n][2*rowlen+j] += der*der_j[m][2];
379 : }
380 270680 : derivatives[n][3*rowlen+i] += nonweight*wval[j];
381 270680 : derivatives[n][3*rowlen+j] += nonweight*wval[i];
382 : }
383 : }
384 : }
385 1200 : }
386 :
387 2 : void ThreeBodyGFunctions::applyNonZeroRankForces( std::vector<double>& outforces ) {
388 2 : taskmanager.applyForces( outforces );
389 2 : }
390 :
391 128 : int ThreeBodyGFunctions::getNumberOfValuesPerTask( std::size_t task_index, const ThreeBodyGFunctionsInput& actiondata ) {
392 128 : return 1;
393 : }
394 :
395 128 : void ThreeBodyGFunctions::getForceIndices( std::size_t task_index,
396 : std::size_t colno,
397 : std::size_t ntotal_force,
398 : const ThreeBodyGFunctionsInput& actiondata,
399 : const ParallelActionsInput& input,
400 : ForceIndexHolder force_indices ) {
401 128 : auto arg0 = ArgumentBookeepingHolder::create( 0, input );
402 128 : auto arg1 = ArgumentBookeepingHolder::create( 1, input );
403 128 : auto arg2 = ArgumentBookeepingHolder::create( 2, input);
404 128 : auto arg3 = ArgumentBookeepingHolder::create( 3, input );
405 128 : std::size_t rowlen = arg3.bookeeping[(1+arg3.ncols)*task_index];
406 128 : if( actiondata.multi_action_input ) {
407 0 : View<const std::size_t> wbooks( arg3.bookeeping.data()+(1+arg3.ncols)*task_index+1, rowlen);
408 0 : for(unsigned j=0; j<rowlen; ++j) {
409 0 : std::size_t matpos = task_index*arg3.ncols + j;
410 0 : std::size_t xpos = getIndex( task_index, wbooks[j], arg0 );
411 0 : std::size_t ypos = getIndex( task_index, wbooks[j], arg1 );
412 0 : std::size_t zpos = getIndex( task_index, wbooks[j], arg2 );
413 0 : for(unsigned i=0; i<input.ncomponents; ++i) {
414 0 : force_indices.indices[i][j] = arg0.start + xpos;
415 0 : force_indices.indices[i][rowlen+j] = arg1.start + ypos;
416 0 : force_indices.indices[i][2*rowlen+j] = arg2.start + zpos;
417 0 : force_indices.indices[i][3*rowlen+j] = arg3.start + matpos;
418 : }
419 : }
420 : } else {
421 4998 : for(unsigned j=0; j<rowlen; ++j) {
422 4870 : std::size_t matpos = task_index*arg3.ncols + j;
423 16286 : for(unsigned i=0; i<input.ncomponents; ++i) {
424 11416 : force_indices.indices[i][j] = arg0.start + matpos;
425 11416 : force_indices.indices[i][rowlen+j] = arg1.start + matpos;
426 11416 : force_indices.indices[i][2*rowlen+j] = arg2.start + matpos;
427 11416 : force_indices.indices[i][3*rowlen+j] = arg3.start + matpos;
428 : }
429 : }
430 : }
431 512 : for(unsigned i=0; i<input.ncomponents; ++i) {
432 384 : force_indices.threadsafe_derivatives_end[i] = force_indices.tot_indices[i] = 4*rowlen;
433 : }
434 128 : }
435 :
436 : }
437 : }
|
