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 unecessary 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 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 : RequiredMatrixElements outmat;
135 : };
136 :
137 : class TorsionsMatrix : public ActionWithMatrix {
138 : public:
139 : using input_type = TorsionsMatrixInput;
140 : using PTM = ParallelTaskManager<TorsionsMatrix>;
141 : private:
142 : PTM taskmanager;
143 : public:
144 : static void registerKeywords( Keywords& keys );
145 : explicit TorsionsMatrix(const ActionOptions&);
146 : unsigned getNumberOfDerivatives() override ;
147 : void prepare() override ;
148 : void calculate() override ;
149 : void applyNonZeroRankForces( std::vector<double>& outforces ) override ;
150 : void getInputData( std::vector<double>& inputdata ) const override ;
151 : static void performTask( std::size_t task_index,
152 : const TorsionsMatrixInput& actiondata,
153 : ParallelActionsInput& input,
154 : ParallelActionsOutput& output );
155 : static int getNumberOfValuesPerTask( std::size_t task_index,
156 : const TorsionsMatrixInput& actiondata );
157 : static void getForceIndices( std::size_t task_index,
158 : std::size_t colno, std::size_t ntotal_force,
159 : const TorsionsMatrixInput& actiondata,
160 : const ParallelActionsInput& input,
161 : ForceIndexHolder force_indices );
162 : };
163 :
164 : PLUMED_REGISTER_ACTION(TorsionsMatrix,"TORSIONS_MATRIX")
165 :
166 14 : void TorsionsMatrix::registerKeywords( Keywords& keys ) {
167 14 : ActionWithMatrix::registerKeywords(keys);
168 28 : keys.addInputKeyword("optional","MASK","matrix","a matrix that is used to used to determine which elements of this matrix to compute");
169 28 : 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");
170 14 : keys.add("atoms","POSITIONS1","the positions to use for the molecules specified using the first argument");
171 14 : keys.add("atoms","POSITIONS2","the positions to use for the molecules specified using the second argument");
172 14 : PTM::registerKeywords( keys );
173 28 : keys.setValueDescription("matrix","the matrix of torsions between the two vectors of input directors");
174 14 : }
175 :
176 7 : TorsionsMatrix::TorsionsMatrix(const ActionOptions&ao):
177 : Action(ao),
178 : ActionWithMatrix(ao),
179 7 : taskmanager(this) {
180 : int nm=getNumberOfMasks();
181 : if( nm<0 ) {
182 : nm = 0;
183 : }
184 7 : if( getNumberOfArguments()-nm!=2 ) {
185 0 : error("should be two arguments to this action, a matrix and a vector");
186 : }
187 7 : if( getPntrToArgument(0)->getRank()!=2 || getPntrToArgument(0)->hasDerivatives() ) {
188 0 : error("first argument to this action should be a matrix");
189 : }
190 7 : if( getPntrToArgument(1)->getRank()!=2 || getPntrToArgument(1)->hasDerivatives() ) {
191 0 : error("second argument to this action should be a matrix");
192 : }
193 7 : if( getPntrToArgument(0)->getShape()[1]!=3 || getPntrToArgument(1)->getShape()[0]!=3 ) {
194 0 : error("number of columns in first matrix and number of rows in second matrix should equal 3");
195 : }
196 :
197 : std::vector<AtomNumber> atoms_a;
198 14 : parseAtomList("POSITIONS1", atoms_a );
199 7 : if( atoms_a.size()!=getPntrToArgument(0)->getShape()[0] ) {
200 0 : error("mismatch between number of atoms specified using POSITIONS1 and number of arguments in vector input");
201 : }
202 7 : log.printf(" using positions of these atoms for vectors in first matrix \n");
203 933 : for(unsigned int i=0; i<atoms_a.size(); ++i) {
204 926 : if ( (i+1) % 25 == 0 ) {
205 36 : log.printf(" \n");
206 : }
207 926 : log.printf(" %d", atoms_a[i].serial());
208 : }
209 7 : log.printf("\n");
210 : std::vector<AtomNumber> atoms_b;
211 14 : parseAtomList("POSITIONS2", atoms_b );
212 7 : if( atoms_b.size()!=getPntrToArgument(1)->getShape()[1] ) {
213 0 : error("mismatch between number of atoms specified using POSITIONS2 and number of arguments in vector input");
214 : }
215 7 : log.printf(" using positions of these atoms for vectors in second matrix \n");
216 1225 : for(unsigned i=0; i<atoms_b.size(); ++i) {
217 1218 : if ( (i+1) % 25 == 0 ) {
218 48 : log.printf(" \n");
219 : }
220 1218 : log.printf(" %d", atoms_b[i].serial());
221 1218 : atoms_a.push_back( atoms_b[i] );
222 : }
223 7 : log.printf("\n");
224 7 : requestAtoms( atoms_a, false );
225 :
226 7 : std::vector<std::size_t> shape(2);
227 7 : shape[0]=getPntrToArgument(0)->getShape()[0];
228 7 : shape[1]=getPntrToArgument(1)->getShape()[1];
229 7 : addValue( shape );
230 14 : setPeriodic("-pi","pi");
231 :
232 7 : if( nm>0 ) {
233 6 : unsigned iarg = getNumberOfArguments()-1;
234 6 : if( getPntrToArgument(iarg)->getRank()!=2
235 6 : || getPntrToArgument(0)->hasDerivatives() ) {
236 0 : error("argument passed to MASK keyword should be a matrix");
237 : }
238 6 : if( getPntrToArgument(iarg)->getShape()[0]!=shape[0]
239 6 : || getPntrToArgument(iarg)->getShape()[1]!=shape[1] ) {
240 0 : error("argument passed to MASK keyword has the wrong shape");
241 : }
242 : }
243 14 : taskmanager.setActionInput( TorsionsMatrixInput() );
244 7 : }
245 :
246 61 : unsigned TorsionsMatrix::getNumberOfDerivatives() {
247 61 : return 3*getNumberOfAtoms()
248 : + 9
249 61 : + getPntrToArgument(0)->getNumberOfStoredValues()
250 61 : + getPntrToArgument(1)->getNumberOfStoredValues();
251 : }
252 :
253 51 : void TorsionsMatrix::prepare() {
254 51 : ActionWithVector::prepare();
255 51 : Value* myval = getPntrToComponent(0);
256 51 : if( myval->getShape()[0]==getPntrToArgument(0)->getShape()[0]
257 51 : && myval->getShape()[1]==getPntrToArgument(1)->getShape()[1] ) {
258 51 : return;
259 : }
260 0 : std::vector<std::size_t> shape(2);
261 0 : shape[0]=getPntrToArgument(0)->getShape()[0];
262 0 : shape[1]=getPntrToArgument(1)->getShape()[1];
263 0 : myval->setShape(shape);
264 0 : myval->reshapeMatrixStore( shape[1] );
265 : }
266 :
267 51 : void TorsionsMatrix::calculate() {
268 51 : updateBookeepingArrays( taskmanager.getActionInput().outmat );
269 51 : taskmanager.setupParallelTaskManager( 21,
270 51 : getNumberOfDerivatives()
271 51 : - getPntrToArgument(0)->getNumberOfStoredValues() );
272 51 : taskmanager.runAllTasks();
273 51 : }
274 :
275 51 : void TorsionsMatrix::getInputData( std::vector<double>& inputdata ) const {
276 51 : std::size_t total_data = getPntrToArgument(0)->getNumberOfStoredValues()
277 51 : + getPntrToArgument(1)->getNumberOfStoredValues()
278 51 : + 3*getNumberOfAtoms();
279 :
280 51 : if( inputdata.size()!=total_data ) {
281 7 : inputdata.resize( total_data );
282 : }
283 :
284 : total_data = 0;
285 : Value* myarg = getPntrToArgument(0);
286 16149 : for(unsigned j=0; j<myarg->getNumberOfStoredValues(); ++j) {
287 16098 : inputdata[total_data] = myarg->get(j,false);
288 16098 : total_data++;
289 : }
290 : myarg = getPntrToArgument(1);
291 25785 : for(unsigned j=0; j<myarg->getNumberOfStoredValues(); ++j) {
292 25734 : inputdata[total_data] = myarg->get(j,false);
293 25734 : total_data++;
294 : }
295 13995 : for(unsigned j=0; j<getNumberOfAtoms(); ++j) {
296 13944 : Vector pos( getPosition(j) );
297 13944 : inputdata[total_data+0] = pos[0];
298 13944 : inputdata[total_data+1] = pos[1];
299 13944 : inputdata[total_data+2] = pos[2];
300 13944 : total_data += 3;
301 : }
302 :
303 51 : }
304 :
305 5544 : void TorsionsMatrix::performTask( std::size_t task_index,
306 : const TorsionsMatrixInput& actiondata,
307 : ParallelActionsInput& input,
308 : ParallelActionsOutput& output ) {
309 5544 : std::size_t fstart = task_index*(1+actiondata.outmat.ncols);
310 : std::size_t nelements = actiondata.outmat[fstart];
311 5544 : auto arg0=ArgumentBookeepingHolder::create( 0, input );
312 5544 : auto arg1=ArgumentBookeepingHolder::create( 1, input );
313 5544 : std::size_t nargdata = arg1.start + arg1.shape[0]*arg1.ncols;
314 :
315 : // Get the position and orientation for the first molecule
316 5544 : std::size_t atbase = nargdata + 3*task_index;
317 5544 : Vector atom1( input.inputdata[atbase],
318 5544 : input.inputdata[atbase+1],
319 5544 : input.inputdata[atbase+2] );
320 5544 : std::size_t agbase = arg0.ncols*task_index;
321 5544 : Vector v1( input.inputdata[agbase],
322 5544 : input.inputdata[agbase+1],
323 5544 : input.inputdata[agbase+2] );
324 :
325 : // Get the distances to all the molecules in the coordination sphere
326 5544 : std::vector<Vector> atom2( nelements );
327 1247276 : for(unsigned i=0; i<nelements; ++i) {
328 : std::size_t at2base = nargdata
329 : + 3*arg0.shape[0]
330 1241732 : + 3*actiondata.outmat[fstart+1+i];
331 1241732 : atom2[i] = Vector( input.inputdata[at2base],
332 1241732 : input.inputdata[at2base+1],
333 1241732 : input.inputdata[at2base+2] ) - atom1;
334 : }
335 5544 : input.pbc->apply( atom2, nelements );
336 :
337 : // Now compute all the torsions
338 : Torsion t;
339 : Vector dv1, dconn, dv2 ;
340 1247276 : for(unsigned i=0; i<nelements; ++i) {
341 1241732 : if( atom2[i].modulo2()<epsilon ) {
342 92 : if( !input.noderiv ) {
343 2 : output.derivatives.subview_n<21>( 21*i ).zero();
344 : }
345 1240498 : continue ;
346 : }
347 :
348 1241640 : std::size_t ag2base = arg1.start + actiondata.outmat[fstart+1+i];
349 1241640 : Vector v2(input.inputdata[ag2base],
350 1241640 : input.inputdata[ag2base+arg1.ncols],
351 1241640 : input.inputdata[ag2base+2*arg1.ncols] );
352 1241640 : output.values[i] = t.compute( v1, atom2[i], v2, dv1, dconn, dv2 );
353 :
354 1241640 : if( input.noderiv ) {
355 1240406 : continue ;
356 : }
357 :
358 1234 : std::size_t base = + 21*i;
359 : output.derivatives.subview_n<3>(base) = dv1;
360 1234 : output.derivatives.subview_n<3>(base+3) = dv2;
361 1234 : output.derivatives.subview_n<3>(base+6) = -dconn;
362 1234 : output.derivatives.subview_n<3>(base+9) = dconn;
363 :
364 1234 : Tensor vir( -extProduct( atom2[i], dconn ) );
365 1234 : View2D<double,3,3> virial( output.derivatives.data() + base + 12 );
366 4936 : for(unsigned j=0; j<3; ++j) {
367 14808 : for(unsigned k=0; k<3; ++k) {
368 11106 : virial[j][k] = vir[j][k];
369 : }
370 : }
371 : }
372 :
373 5544 : }
374 :
375 23 : void TorsionsMatrix::applyNonZeroRankForces( std::vector<double>& outforces ) {
376 23 : taskmanager.applyForces( outforces );
377 23 : }
378 :
379 178 : int TorsionsMatrix::getNumberOfValuesPerTask( std::size_t task_index,
380 : const TorsionsMatrixInput& actiondata ) {
381 178 : std::size_t fstart = task_index*(1+actiondata.outmat.ncols);
382 178 : return actiondata.outmat[fstart];
383 : }
384 :
385 1236 : void TorsionsMatrix::getForceIndices( std::size_t task_index,
386 : std::size_t colno,
387 : std::size_t ntotal_force,
388 : const TorsionsMatrixInput& actiondata,
389 : const ParallelActionsInput& input,
390 : ForceIndexHolder force_indices ) {
391 1236 : auto arg0=ArgumentBookeepingHolder::create( 0, input );
392 1236 : auto arg1=ArgumentBookeepingHolder::create( 1, input );
393 1236 : std::size_t fstart = task_index*(1+actiondata.outmat.ncols);
394 1236 : std::size_t arg1start = task_index*arg0.ncols;
395 1236 : force_indices.indices[0][0] = arg1start;
396 1236 : force_indices.indices[0][1] = arg1start + 1;
397 1236 : force_indices.indices[0][2] = arg1start + 2;
398 1236 : force_indices.threadsafe_derivatives_end[0] = 3;
399 1236 : force_indices.indices[0][3] = arg1.start
400 1236 : + actiondata.outmat[fstart+1+colno];
401 1236 : force_indices.indices[0][4] = arg1.start
402 1236 : + actiondata.outmat[fstart+1+colno]
403 1236 : + arg1.ncols;
404 1236 : force_indices.indices[0][5] = arg1.start
405 1236 : + actiondata.outmat[fstart+1+colno]
406 1236 : + 2*arg1.ncols;
407 1236 : std::size_t atomstart = arg1.start + arg1.shape[0]*arg1.ncols;
408 1236 : std::size_t atom1start = atomstart + 3*task_index;
409 1236 : force_indices.indices[0][6] = atom1start;
410 1236 : force_indices.indices[0][7] = atom1start + 1;
411 1236 : force_indices.indices[0][8] = atom1start + 2;
412 : std::size_t atom2start = atomstart
413 : + 3*arg0.shape[0]
414 1236 : + 3*actiondata.outmat[fstart+1+colno];
415 1236 : force_indices.indices[0][9] = atom2start;
416 1236 : force_indices.indices[0][10] = atom2start + 1;
417 1236 : force_indices.indices[0][11] = atom2start + 2;
418 : std::size_t n=12;
419 12360 : for(unsigned j=ntotal_force-9; j<ntotal_force; ++j) {
420 11124 : force_indices.indices[0][n] = j;
421 11124 : ++n;
422 : }
423 1236 : force_indices.tot_indices[0] = 21;
424 1236 : }
425 :
426 : }
427 : }
|
