Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2023 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/ActionWithMatrix.h"
23 : #include "core/ActionRegister.h"
24 : #include "core/ParallelTaskManager.h"
25 : #include "tools/Torsion.h"
26 :
27 : //+PLUMEDOC MCOLVAR TORSIONS_MATRIX
28 : /*
29 : Calculate the matrix of torsions between two vectors of molecules
30 :
31 : This action was implemented to ensure that we can calculate the [SMAC](SMAC.md) collective variable that is discussed in
32 : [this paper](https://www.sciencedirect.com/science/article/abs/pii/S0009250914004503?via%3Dihub). This particular action
33 : tracks the relative orientations for all the pairs of molecules in a set much like the variables described in the crystdistrib module.
34 :
35 : The orientations of molecules can be specified using either [PLANE](PLANE.md) or [DISTANCE](DISTANCE.md). The example below shows how you can use
36 : internal vectors connecting two atoms in the molecules to define the orientation of that molecule. Three of these internal
37 : vectors are calculated using a DISTANCE command in the input below. The matrix of torsional angles between these various
38 : vectors is then computed:
39 :
40 : ```plumed
41 : d1: DISTANCE ATOMS1=1,5 ATOMS2=11,15 ATOMS3=21,25 COMPONENTS
42 : s: VSTACK ARG=d1.x,d1.y,d1.z
43 : sT: TRANSPOSE ARG=s
44 : m: TORSIONS_MATRIX ARG=s,sT POSITIONS1=1,11,21 POSITIONS2=1,11,21
45 : PRINT ARG=m FILE=matrix
46 : ```
47 :
48 : In this example, the torsional angle in element $(1,2)$ of the matrix with label `m` is the angle between the plane containing atoms 1,5 and 10 and the plane
49 : connecting atoms 1,10 and 15. In other words, the elements in this matrix are the torsional angles between the vectors in the input matrices
50 : around the vector connecting the corresponding atomic positions that are specified using the `POSTIONS` keyword.
51 :
52 : You can also calculate a matrix of torsional angles between two different groups of molecules by using an input like the one below:
53 :
54 : ```plumed
55 : pA: PLANE ATOMS1=1,2,3 ATOMS2=11,12,13
56 : sA: VSTACK ARG=pA.x,pA.y,pA.z
57 : pB: PLANE ATOMS1=21,22,23 ATOMS2=31,32,33 ATOMS3=41,42,43
58 : sB: VSTACK ARG=pB.x,pB.y,pB.z
59 : sBT: TRANSPOSE ARG=sB
60 : m: TORSIONS_MATRIX ARG=sA,sBT POSITIONS1=1,11 POSITIONS2=21,31,41
61 : PRINT ARG=m FILE=matrix
62 : ```
63 :
64 : In this example, the orientations of the molecules are specified using the [PLANE](PLANE.md) action and is given by a normal to the plane containing the three atoms from the molecule
65 : that was specified. The final output is $2 \times 3$ matrix that contains all the torsional angles between the molecules defined by the two PLANE actions.
66 :
67 : ## Performance considerations
68 :
69 : Suppose that you are using an input like the one shown below to calculate the average torsion angle between neighboring molecules:
70 :
71 : ```plumed
72 : # Notice that in a realistic version of this calculation you would likely
73 : # by calculating the orientations of many more molecules using this command
74 : d: DISTANCE COMPONENTS ATOMS1=1,5 ATOMS2=11,15 ATOMS3=21,25 ATOMS4=31,35
75 : s: VSTACK ARG=d.x,d.y,d.z
76 : sT: TRANSPOSE ARG=s
77 : # Calculate a contact matrix in which element i,j is 1 if molecules
78 : # i and j are neighbors.
79 : c: CONTACT_MATRIX GROUP=1,11,21,35 SWITCH={RATIONAL R_0=0.1 D_MAX=0.3}
80 : # Now calculate all the torsions
81 : t: TORSIONS_MATRIX ARG=s,sT POSITIONS1=1,11,21,31 POSITIONS2=1,11,21,31
82 : # And the product between the contact matrix and the torsions
83 : tc: CUSTOM ARG=c,t FUNC=x*y PERIODIC=NO
84 : # Total of all the torsional angles in the first coordination sphere
85 : tsum: SUM ARG=tc PERIODIC=NO
86 : # Total number of neighbouring atoms
87 : bsum: SUM ARG=c PERIODIC=NO
88 : # And finally the average torsion angle
89 : avt: CUSTOM ARG=tc,c FUNC=x/y PERIODIC=NO
90 : ```
91 :
92 : If you have a large number of molecules the most expensive part of this calculation will be the evalulation of the TORSIONS_MATRIX as you need to evaluate one torsion anlge for every pair of molecules.
93 : Furthermore, this computation is unnecessary as most pairs of molecules will __not__ be neighbors. In other words, for the majority of the molecular pairs element the corresponding element of the
94 : [CONTACT_MATRIX](CONTACT_MATRIX.md) will be zero. Consequently, when you compute the corresponding element `tc` by multiplying the torsion $i,j$ by the crresponding $i,j$ element of the
95 : [CONTACT_MATRIX](CONTACT_MATRIX.md) you will get zero.
96 :
97 : Thankfully PLUMED allows you to exploit this fact through the MASK keyword as illustrated below:
98 :
99 : ```plumed
100 : # Notice that in a realistic version of this calculation you would likely
101 : # by calculating the orientations of many more molecules using this command
102 : d: DISTANCE COMPONENTS ATOMS1=1,5 ATOMS2=11,15 ATOMS3=21,25 ATOMS4=31,35
103 : s: VSTACK ARG=d.x,d.y,d.z
104 : sT: TRANSPOSE ARG=s
105 : # Calculate a contact matrix in which element i,j is 1 if molecules
106 : # i and j are neighbors.
107 : c: CONTACT_MATRIX GROUP=1,11,21,35 SWITCH={RATIONAL R_0=0.1 D_MAX=0.3}
108 : # Now calculate all the torsions
109 : t: TORSIONS_MATRIX ...
110 : MASK=c ARG=s,sT
111 : POSITIONS1=1,11,21,31 POSITIONS2=1,11,21,31
112 : ...
113 : # And the product between the contact matrix and the torsions
114 : tc: CUSTOM ARG=c,t FUNC=x*y PERIODIC=NO
115 : # Total of all the torsional angles in the first coordination sphere
116 : tsum: SUM ARG=tc PERIODIC=NO
117 : # Total number of neighbouring atoms
118 : bsum: SUM ARG=c PERIODIC=NO
119 : # And finally the average torsion angle
120 : avt: CUSTOM ARG=tc,c FUNC=x/y PERIODIC=NO
121 : ```
122 :
123 : Adding the instruction `MASK=c` to the TORSIONS_MATRIX command here ensures that element $i,j$ of the matrix `t` is only computed if the corresponding element of `c` is non-zero. By using this command
124 : you thus avoid the computational expense associated with evaluating the full set of pairwise torsions.
125 :
126 : */
127 : //+ENDPLUMEDOC
128 :
129 : namespace PLMD {
130 : namespace adjmat {
131 :
132 0 : class TorsionsMatrixInput {
133 : public:
134 : bool no_thread_gather;
135 : bool gatherForceOnColumns;
136 : RequiredMatrixElements outmat;
137 : };
138 :
139 : class TorsionsMatrix : public ActionWithMatrix {
140 : public:
141 : using input_type = TorsionsMatrixInput;
142 : using PTM = ParallelTaskManager<TorsionsMatrix>;
143 : typedef typename PTM::ParallelActionsInput ParallelActionsInput;
144 : typedef typename PTM::ParallelActionsOutput ParallelActionsOutput;
145 : private:
146 : PTM taskmanager;
147 : public:
148 : static void registerKeywords( Keywords& keys );
149 : explicit TorsionsMatrix(const ActionOptions&);
150 : unsigned getNumberOfDerivatives() override ;
151 : void prepare() override ;
152 : void calculate() override ;
153 : void applyNonZeroRankForces( std::vector<double>& outforces ) override ;
154 : void getAtomicInputData( std::vector<double>& inputdata ) const override ;
155 : void getAtomicInputData( std::vector<float>& inputdata ) const override ;
156 : void getArgumentInputData( std::vector<double>& inputdata ) const override ;
157 : void getArgumentInputData( std::vector<float>& inputdata ) const override ;
158 : static void performTask( std::size_t task_index,
159 : const TorsionsMatrixInput& actiondata,
160 : ParallelActionsInput& input,
161 : ParallelActionsOutput& output );
162 : static int getNumberOfValuesPerTask( std::size_t task_index,
163 : const TorsionsMatrixInput& actiondata );
164 : static void getForceIndices( std::size_t task_index,
165 : std::size_t colno, std::size_t ntotal_force,
166 : const TorsionsMatrixInput& actiondata,
167 : const ParallelActionsInput& input,
168 : ForceIndexHolder force_indices );
169 : };
170 :
171 : PLUMED_REGISTER_ACTION(TorsionsMatrix,"TORSIONS_MATRIX")
172 :
173 18 : void TorsionsMatrix::registerKeywords( Keywords& keys ) {
174 18 : ActionWithMatrix::registerKeywords(keys);
175 36 : keys.addInputKeyword("optional","MASK","matrix","a matrix that is used to used to determine which elements of this matrix to compute");
176 36 : keys.addInputKeyword("compulsory","ARG","matrix","an Nx3 and a 3xN matrix that contain the bond vectors that you would like to determine the torsion angles between");
177 18 : keys.add("atoms","POSITIONS1","the positions to use for the molecules specified using the first argument");
178 18 : keys.add("atoms","POSITIONS2","the positions to use for the molecules specified using the second argument");
179 18 : PTM::registerKeywords( keys );
180 36 : keys.setValueDescription("matrix","the matrix of torsions between the two vectors of input directors");
181 18 : }
182 :
183 9 : TorsionsMatrix::TorsionsMatrix(const ActionOptions&ao):
184 : Action(ao),
185 : ActionWithMatrix(ao),
186 9 : taskmanager(this) {
187 : int nm=getNumberOfMasks();
188 9 : if( nm<0 ) {
189 : nm = 0;
190 : }
191 9 : if( getNumberOfArguments()-nm!=2 ) {
192 0 : error("should be two arguments to this action, a matrix and a vector");
193 : }
194 9 : if( getPntrToArgument(0)->getRank()!=2 || getPntrToArgument(0)->hasDerivatives() ) {
195 0 : error("first argument to this action should be a matrix");
196 : }
197 9 : if( getPntrToArgument(1)->getRank()!=2 || getPntrToArgument(1)->hasDerivatives() ) {
198 0 : error("second argument to this action should be a matrix");
199 : }
200 9 : if( getPntrToArgument(0)->getShape()[1]!=3 || getPntrToArgument(1)->getShape()[0]!=3 ) {
201 0 : error("number of columns in first matrix and number of rows in second matrix should equal 3");
202 : }
203 :
204 9 : std::vector<AtomNumber> atoms_a;
205 18 : parseAtomList("POSITIONS1", atoms_a );
206 9 : if( atoms_a.size()!=getPntrToArgument(0)->getShape()[0] ) {
207 0 : error("mismatch between number of atoms specified using POSITIONS1 and number of arguments in vector input");
208 : }
209 9 : log.printf(" using positions of these atoms for vectors in first matrix \n");
210 939 : for(unsigned int i=0; i<atoms_a.size(); ++i) {
211 930 : if ( (i+1) % 25 == 0 ) {
212 36 : log.printf(" \n");
213 : }
214 930 : log.printf(" %d", atoms_a[i].serial());
215 : }
216 9 : log.printf("\n");
217 9 : std::vector<AtomNumber> atoms_b;
218 18 : parseAtomList("POSITIONS2", atoms_b );
219 9 : if( atoms_b.size()!=getPntrToArgument(1)->getShape()[1] ) {
220 0 : error("mismatch between number of atoms specified using POSITIONS2 and number of arguments in vector input");
221 : }
222 9 : log.printf(" using positions of these atoms for vectors in second matrix \n");
223 1232 : for(unsigned i=0; i<atoms_b.size(); ++i) {
224 1223 : if ( (i+1) % 25 == 0 ) {
225 48 : log.printf(" \n");
226 : }
227 1223 : log.printf(" %d", atoms_b[i].serial());
228 1223 : atoms_a.push_back( atoms_b[i] );
229 : }
230 9 : log.printf("\n");
231 9 : requestAtoms( atoms_a, false );
232 :
233 9 : std::vector<std::size_t> shape(2);
234 9 : shape[0]=getPntrToArgument(0)->getShape()[0];
235 9 : shape[1]=getPntrToArgument(1)->getShape()[1];
236 9 : addValue( shape );
237 18 : setPeriodic("-pi","pi");
238 :
239 9 : if( nm>0 ) {
240 6 : unsigned iarg = getNumberOfArguments()-1;
241 6 : if( getPntrToArgument(iarg)->getRank()!=2
242 6 : || getPntrToArgument(0)->hasDerivatives() ) {
243 0 : error("argument passed to MASK keyword should be a matrix");
244 : }
245 6 : if( getPntrToArgument(iarg)->getShape()[0]!=shape[0]
246 6 : || getPntrToArgument(iarg)->getShape()[1]!=shape[1] ) {
247 0 : error("argument passed to MASK keyword has the wrong shape");
248 : }
249 : }
250 : TorsionsMatrixInput torinput;
251 9 : torinput.gatherForceOnColumns = false;
252 9 : torinput.no_thread_gather = no_thread_gather;
253 9 : taskmanager.setActionInput( torinput );
254 9 : }
255 :
256 65 : unsigned TorsionsMatrix::getNumberOfDerivatives() {
257 65 : return 3*getNumberOfAtoms()
258 : + 9
259 65 : + getPntrToArgument(0)->getNumberOfStoredValues()
260 65 : + getPntrToArgument(1)->getNumberOfStoredValues();
261 : }
262 :
263 55 : void TorsionsMatrix::prepare() {
264 55 : ActionWithVector::prepare();
265 55 : Value* myval = getPntrToComponent(0);
266 55 : if( myval->getShape()[0]==getPntrToArgument(0)->getShape()[0]
267 55 : && myval->getShape()[1]==getPntrToArgument(1)->getShape()[1] ) {
268 55 : return;
269 : }
270 0 : std::vector<std::size_t> shape(2);
271 0 : shape[0]=getPntrToArgument(0)->getShape()[0];
272 0 : shape[1]=getPntrToArgument(1)->getShape()[1];
273 0 : myval->setShape(shape);
274 0 : myval->reshapeMatrixStore( shape[1] );
275 0 : }
276 :
277 55 : void TorsionsMatrix::calculate() {
278 55 : updateBookeepingArrays( taskmanager.getActionInput().outmat );
279 55 : if( no_thread_gather ) {
280 4 : taskmanager.setupParallelTaskManager( 15, 9 );
281 : } else {
282 51 : taskmanager.setupParallelTaskManager( 21,
283 51 : getNumberOfDerivatives()
284 51 : - getPntrToArgument(0)->getNumberOfStoredValues() );
285 : }
286 55 : taskmanager.runAllTasks();
287 55 : }
288 :
289 55 : void TorsionsMatrix::getAtomicInputData( std::vector<double>& inputdata ) const {
290 55 : std::size_t total_data = getPntrToArgument(0)->getNumberOfStoredValues()
291 55 : + getPntrToArgument(1)->getNumberOfStoredValues()
292 55 : + 3*getNumberOfAtoms();
293 :
294 55 : if( inputdata.size()!=total_data ) {
295 9 : inputdata.resize( total_data );
296 : }
297 :
298 55 : total_data = total_data - 3*getNumberOfAtoms();
299 14018 : for(unsigned j=0; j<getNumberOfAtoms(); ++j) {
300 13963 : Vector pos( getPosition(j) );
301 13963 : inputdata[total_data+0] = pos[0];
302 13963 : inputdata[total_data+1] = pos[1];
303 13963 : inputdata[total_data+2] = pos[2];
304 13963 : total_data += 3;
305 : }
306 55 : }
307 :
308 0 : void TorsionsMatrix::getAtomicInputData( std::vector<float>& inputdata ) const {
309 0 : std::size_t total_data = getPntrToArgument(0)->getNumberOfStoredValues()
310 0 : + getPntrToArgument(1)->getNumberOfStoredValues()
311 0 : + 3*getNumberOfAtoms();
312 :
313 0 : if( inputdata.size()!=total_data ) {
314 0 : inputdata.resize( total_data );
315 : }
316 :
317 0 : total_data = total_data - 3*getNumberOfAtoms();
318 0 : for(unsigned j=0; j<getNumberOfAtoms(); ++j) {
319 0 : Vector pos( getPosition(j) );
320 0 : inputdata[total_data+0] = pos[0];
321 0 : inputdata[total_data+1] = pos[1];
322 0 : inputdata[total_data+2] = pos[2];
323 0 : total_data += 3;
324 : }
325 0 : }
326 :
327 55 : void TorsionsMatrix::getArgumentInputData( std::vector<double>& inputdata ) const {
328 : unsigned total_data = 0;
329 : Value* myarg = getPntrToArgument(0);
330 16177 : for(unsigned j=0; j<myarg->getNumberOfStoredValues(); ++j) {
331 16122 : inputdata[total_data] = myarg->get(j,false);
332 16122 : total_data++;
333 : }
334 : myarg = getPntrToArgument(1);
335 25822 : for(unsigned j=0; j<myarg->getNumberOfStoredValues(); ++j) {
336 25767 : inputdata[total_data] = myarg->get(j,false);
337 25767 : total_data++;
338 : }
339 55 : }
340 :
341 0 : void TorsionsMatrix::getArgumentInputData( std::vector<float>& inputdata ) const {
342 : unsigned total_data = 0;
343 : Value* myarg = getPntrToArgument(0);
344 0 : for(unsigned j=0; j<myarg->getNumberOfStoredValues(); ++j) {
345 0 : inputdata[total_data] = myarg->get(j,false);
346 0 : total_data++;
347 : }
348 : myarg = getPntrToArgument(1);
349 0 : for(unsigned j=0; j<myarg->getNumberOfStoredValues(); ++j) {
350 0 : inputdata[total_data] = myarg->get(j,false);
351 0 : total_data++;
352 : }
353 0 : }
354 :
355 : // Perform one parallel task for TORSIONS_MATRIX.
356 : // The task index identifies a compact row in `outmat` that stores:
357 : // - the number of active matrix elements for this task (at offset fstart)
358 : // - then the column indices of those elements.
359 : // For each active element, this routine reads orientation vectors from ARG0/ARG1,
360 : // reads the corresponding atomic positions, computes the torsion, and writes
361 : // values/forces into the task-local output buffers.
362 : //
363 : // Data layout notes:
364 : // - `arg0` and `arg1` describe where the two argument matrices are stored in `input`.
365 : // - `nargdata` marks the end of argument data and the start of packed atomic coordinates.
366 : // - `gatherForceOnColumns` controls whether force contributions are accumulated by
367 : // iterating primarily over columns (true) or rows (false).
368 5571 : void TorsionsMatrix::performTask( std::size_t task_index,
369 : const TorsionsMatrixInput& actiondata,
370 : ParallelActionsInput& input,
371 : ParallelActionsOutput& output ) {
372 : // Locate the list of matrix elements assigned to this task.
373 5571 : std::size_t fstart = task_index*(1+actiondata.outmat.ncols);
374 : std::size_t nelements = actiondata.outmat[fstart];
375 :
376 : // Build bookkeeping descriptors for the two input argument matrices.
377 5571 : auto arg0=ArgumentBookeepingHolder::create( 0, input );
378 5571 : auto arg1=ArgumentBookeepingHolder::create( 1, input );
379 :
380 : // Offset of packed atomic position data in the task input buffer.
381 5571 : std::size_t nargdata = arg1.start + arg1.shape[0]*arg1.ncols;
382 :
383 5571 : Vector v1, v2;
384 5571 : std::vector<Vector> atom2( nelements );
385 : // Two execution modes: gather force contributions on columns or on rows.
386 5571 : if( actiondata.gatherForceOnColumns ) {
387 11 : std::size_t at2base = nargdata + 3*arg0.shape[0] + 3*task_index;
388 11 : Vector atom1( input.inputdata[at2base],
389 11 : input.inputdata[at2base+1],
390 11 : input.inputdata[at2base+2] );
391 11 : std::size_t ag2base = arg1.start + task_index;
392 11 : v2 = Vector( input.inputdata[ag2base],
393 11 : input.inputdata[ag2base+arg1.ncols],
394 11 : input.inputdata[ag2base+2*arg1.ncols] );
395 33 : for(unsigned i=0; i<nelements; ++i) {
396 22 : std::size_t atbase = nargdata + 3*i;
397 22 : atom2[i] = atom1 - Vector( input.inputdata[atbase],
398 22 : input.inputdata[atbase+1],
399 22 : input.inputdata[atbase+2] );
400 : }
401 : } else {
402 5560 : std::size_t atbase = nargdata + 3*task_index;
403 5560 : Vector atom1( input.inputdata[atbase],
404 5560 : input.inputdata[atbase+1],
405 5560 : input.inputdata[atbase+2] );
406 5560 : std::size_t agbase = arg0.ncols*task_index;
407 5560 : v1 = Vector( input.inputdata[agbase],
408 5560 : input.inputdata[agbase+1],
409 5560 : input.inputdata[agbase+2] );
410 :
411 1044400 : for(unsigned i=0; i<nelements; ++i) {
412 : std::size_t at2base = nargdata
413 : + 3*arg0.shape[0]
414 1038840 : + 3*actiondata.outmat[fstart+1+i];
415 1038840 : atom2[i] = Vector( input.inputdata[at2base],
416 1038840 : input.inputdata[at2base+1],
417 1038840 : input.inputdata[at2base+2] ) - atom1;
418 : }
419 : }
420 5571 : input.pbc->apply( atom2, nelements );
421 :
422 : // Now compute all the torsions
423 : Torsion t;
424 5571 : Vector dv1, dconn, dv2 ;
425 : unsigned nderper = 21;
426 5571 : if( actiondata.no_thread_gather ) {
427 : nderper = 15;
428 : }
429 1044433 : for(unsigned i=0; i<nelements; ++i) {
430 1038862 : if( atom2[i].modulo2()<epsilon ) {
431 98 : if( !input.noderiv ) {
432 6 : output.derivatives.subview( nderper, static_cast<std::size_t>(nderper) * static_cast<std::size_t>(i) ).zero();
433 : }
434 98 : continue ;
435 : }
436 :
437 1038764 : if( actiondata.gatherForceOnColumns ) {
438 20 : std::size_t ag1base = arg0.ncols*actiondata.outmat[fstart+1+i];
439 20 : v1 = Vector( input.inputdata[ag1base],
440 20 : input.inputdata[ag1base+1],
441 20 : input.inputdata[ag1base+2] );
442 : } else {
443 1038744 : std::size_t ag2base = arg1.start + actiondata.outmat[fstart+1+i];
444 1038744 : v2 = Vector( input.inputdata[ag2base],
445 1038744 : input.inputdata[ag2base+arg1.ncols],
446 1038744 : input.inputdata[ag2base+2*arg1.ncols] );
447 : }
448 1038764 : output.values[i] = t.compute( v1, atom2[i], v2, dv1, dconn, dv2 );
449 :
450 1038764 : if( input.noderiv ) {
451 1037490 : continue ;
452 : }
453 :
454 1274 : std::size_t base = static_cast<std::size_t>(nderper) * static_cast<std::size_t>(i);
455 1274 : if( !actiondata.no_thread_gather ) {
456 : output.derivatives.subview_n<3>(base) = dv1;
457 1234 : output.derivatives.subview_n<3>(base+3) = -dconn;
458 1234 : output.derivatives.subview_n<3>(base+6) = dv2;
459 1234 : output.derivatives.subview_n<3>(base+9) = dconn;
460 1234 : Tensor vir( -extProduct( atom2[i], dconn ) );
461 1234 : View2D<double,3,3> virial( output.derivatives.data() + base + 12 );
462 4936 : for(unsigned j=0; j<3; ++j) {
463 14808 : for(unsigned k=0; k<3; ++k) {
464 11106 : virial[j][k] = vir[j][k];
465 : }
466 : }
467 40 : } else if( actiondata.gatherForceOnColumns ) {
468 : output.derivatives.subview_n<3>(base) = dv2;
469 20 : output.derivatives.subview_n<3>(base+3) = dconn;
470 : } else {
471 : output.derivatives.subview_n<3>(base) = dv1;
472 20 : output.derivatives.subview_n<3>(base+3) = -dconn;
473 20 : Tensor vir( -extProduct( atom2[i], dconn ) );
474 20 : View2D<double,3,3> virial( output.derivatives.data() + base + 6 );
475 80 : for(unsigned j=0; j<3; ++j) {
476 240 : for(unsigned k=0; k<3; ++k) {
477 180 : virial[j][k] = vir[j][k];
478 : }
479 : }
480 : }
481 : }
482 :
483 5571 : }
484 :
485 27 : void TorsionsMatrix::applyNonZeroRankForces( std::vector<double>& outforces ) {
486 27 : taskmanager.applyForces( outforces );
487 27 : if( no_thread_gather ) {
488 4 : getColumnBookeepingArrays( taskmanager.getActionInput().outmat );
489 4 : taskmanager.getActionInput().gatherForceOnColumns = gatherForceOnColumns = true;
490 4 : taskmanager.setNForceScalars( getNumberOfComponents()*maxcolsize );
491 4 : taskmanager.applyForces( outforces, false );
492 4 : taskmanager.setNForceScalars( getNumberOfComponents()*getPntrToComponent(0)->getNumberOfColumns() );
493 4 : taskmanager.getActionInput().gatherForceOnColumns = gatherForceOnColumns = false;
494 : }
495 27 : }
496 :
497 197 : int TorsionsMatrix::getNumberOfValuesPerTask( std::size_t task_index,
498 : const TorsionsMatrixInput& actiondata ) {
499 197 : std::size_t fstart = task_index*(1+actiondata.outmat.ncols);
500 197 : return actiondata.outmat[fstart];
501 : }
502 :
503 1280 : void TorsionsMatrix::getForceIndices( std::size_t task_index,
504 : std::size_t colno,
505 : std::size_t ntotal_force,
506 : const TorsionsMatrixInput& actiondata,
507 : const ParallelActionsInput& input,
508 : ForceIndexHolder force_indices ) {
509 1280 : auto arg0=ArgumentBookeepingHolder::create( 0, input );
510 1280 : auto arg1=ArgumentBookeepingHolder::create( 1, input );
511 1280 : if( !actiondata.no_thread_gather ) {
512 1236 : std::size_t fstart = task_index*(1+actiondata.outmat.ncols);
513 1236 : std::size_t arg1start = task_index*arg0.ncols;
514 1236 : force_indices.indices[0][0] = arg1start;
515 1236 : force_indices.indices[0][1] = arg1start + 1;
516 1236 : force_indices.indices[0][2] = arg1start + 2;
517 1236 : std::size_t atomstart = arg1.start + arg1.shape[0]*arg1.ncols;
518 1236 : std::size_t atom1start = atomstart + 3*task_index;
519 1236 : force_indices.indices[0][3] = atom1start;
520 1236 : force_indices.indices[0][4] = atom1start + 1;
521 1236 : force_indices.indices[0][5] = atom1start + 2;
522 1236 : force_indices.threadsafe_derivatives_end[0] = 6;
523 1236 : force_indices.indices[0][6] = arg1.start
524 1236 : + actiondata.outmat[fstart+1+colno];
525 1236 : force_indices.indices[0][7] = arg1.start
526 1236 : + actiondata.outmat[fstart+1+colno]
527 1236 : + arg1.ncols;
528 1236 : force_indices.indices[0][8] = arg1.start
529 1236 : + actiondata.outmat[fstart+1+colno]
530 1236 : + 2*arg1.ncols;
531 : std::size_t atom2start = atomstart
532 : + 3*arg0.shape[0]
533 1236 : + 3*actiondata.outmat[fstart+1+colno];
534 1236 : force_indices.indices[0][9] = atom2start;
535 1236 : force_indices.indices[0][10] = atom2start + 1;
536 1236 : force_indices.indices[0][11] = atom2start + 2;
537 : std::size_t n=12;
538 12360 : for(unsigned j=ntotal_force-9; j<ntotal_force; ++j) {
539 11124 : force_indices.indices[0][n] = j;
540 11124 : ++n;
541 : }
542 1236 : force_indices.tot_indices[0] = 21;
543 44 : } else if( actiondata.gatherForceOnColumns ) {
544 22 : std::size_t fstart = colno*(1+actiondata.outmat.ncols);
545 22 : force_indices.indices[0][0] = arg1.start
546 22 : + actiondata.outmat[fstart+1+task_index];
547 22 : force_indices.indices[0][1] = arg1.start
548 22 : + actiondata.outmat[fstart+1+task_index]
549 22 : + arg1.ncols;
550 22 : force_indices.indices[0][2] = arg1.start
551 22 : + actiondata.outmat[fstart+1+task_index]
552 22 : + 2*arg1.ncols;
553 22 : std::size_t atomstart = arg1.start + arg1.shape[0]*arg1.ncols;
554 : std::size_t atom2start = atomstart
555 : + 3*arg0.shape[0]
556 22 : + 3*actiondata.outmat[fstart+1+task_index];
557 22 : force_indices.indices[0][3] = atom2start;
558 22 : force_indices.indices[0][4] = atom2start + 1;
559 22 : force_indices.indices[0][5] = atom2start + 2;
560 22 : force_indices.threadsafe_derivatives_end[0] = 6;
561 22 : force_indices.tot_indices[0] = 6;
562 : } else {
563 22 : std::size_t arg1start = task_index*arg0.ncols;
564 22 : force_indices.indices[0][0] = arg1start;
565 22 : force_indices.indices[0][1] = arg1start + 1;
566 22 : force_indices.indices[0][2] = arg1start + 2;
567 22 : std::size_t atomstart = arg1.start + arg1.shape[0]*arg1.ncols;
568 22 : std::size_t atom1start = atomstart + 3*task_index;
569 22 : force_indices.indices[0][3] = atom1start;
570 22 : force_indices.indices[0][4] = atom1start + 1;
571 22 : force_indices.indices[0][5] = atom1start + 2;
572 22 : force_indices.threadsafe_derivatives_end[0] = 6;
573 : std::size_t n=6;
574 220 : for(unsigned j=ntotal_force-9; j<ntotal_force; ++j) {
575 198 : force_indices.indices[0][n] = j;
576 198 : ++n;
577 : }
578 22 : force_indices.tot_indices[0] = 15;
579 : }
580 1280 : }
581 :
582 : }
583 : }
|
