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 : // #ifdef __PLUMED_HAS_OPENACC
23 : // #define __PLUMED_USE_OPENACC 1
24 : // #endif //__PLUMED_HAS_OPENACC
25 : #include "core/ActionWithVector.h"
26 : #include "core/ParallelTaskManager.h"
27 : #include "core/ActionRegister.h"
28 : #include "tools/Torsion.h"
29 :
30 :
31 : #include <iostream>
32 :
33 : namespace PLMD {
34 : namespace crystdistrib {
35 :
36 : //+PLUMEDOC MCOLVAR QUATERNION_BOND_PRODUCT_MATRIX
37 : /*
38 : Calculate the product between a matrix of quaternions and the bonds connecting molecules
39 :
40 : Calculate the product of a quaternion, times a vector, times the conjugate of the quaternion. Geometrically, this is the given vector projected onto the coordinate system defined by the quaternion. For context, the QUATERNION action defines an orthogonal coordinate frame given 3 atoms (i.e. one can define a coordinate system based on the structure of a given molecule). This action can then project a vector from the given molecular frame, toward another molecule, essentially pointing toward molecule 2, from the perspective of molecule 1. See QUATERNION for information about the molecular coordinate frame. Given a quaternion centered on molecule 1 $\mathbf{q1}$, and the vector connecting molecule 1, and some other molecule 2, $\mathbf{r_{21}}$, the following calculation is performed:
41 :
42 : $$
43 : \mathbf{r} = \overline{\mathbf{q_1}} \mathbf{r_{21}} \mathbf{q_1}
44 : $$
45 :
46 : where the overline denotes the quaternion conjugate. Internally, the vector $\mathbf{r_{21}}$ is treated as a quaternion with zero real part. Such a multiplication will always yield another quaternion with zero real part, and the results can be interpreted as an ordinary vector in $\mathbb{R}^3$ Nevertheless, this action has four components, the first of which, w, will always be entirely zeros. Finally, the resulting vector is normalized within the action, and the real length can be returned by multiplying each component by the norm of the vector given to the action. The quaternion should be a vector, and the distances a matrix.
47 :
48 : In this example, 3 quaternion frames are calculated, and multiplied element-wise onto a distance matrix, yielding 9 vectors overall.
49 :
50 : ```plumed
51 : #calculate the quaternion frames for 3 molecules
52 : quat: QUATERNION ATOMS1=1,2,3 ATOMS2=4,5,6 ATOMS3=7,8,9
53 : #also find the distance between the 3 origins of the molecule frames
54 : c1: DISTANCE_MATRIX GROUP=1,4,7 CUTOFF=100.0 COMPONENTS
55 : qp: QUATERNION_BOND_PRODUCT_MATRIX ARG=quat.*,c1.*
56 : #this is now a matrix showing how each molecule is oriented in 3D space
57 : #relative to eachother's origins
58 : #now use in a simple collective variable
59 : ops: CUSTOM ARG=qp.* VAR=w,i,j,k FUNC=w+i+j+k PERIODIC=NO
60 : #w could have been ignored because it is always zero.
61 :
62 : ```
63 :
64 : */
65 : //+ENDPLUMEDOC
66 :
67 : struct QuatBondProdMatInput {
68 : void toACCDevice() const {
69 : #pragma acc enter data copyin(this[0:1])
70 : }
71 : void removeFromACCDevice() const {
72 : #pragma acc exit data delete(this[0:1])
73 : }
74 :
75 : };
76 :
77 : class QuaternionBondProductMatrix : public ActionWithVector {
78 : public:
79 : using input_type = QuatBondProdMatInput;
80 : using PTM = ParallelTaskManager<QuaternionBondProductMatrix>;
81 : private:
82 : PTM taskmanager;
83 : std::vector<unsigned> active_tasks;
84 : // const Vector4d& rightMultiply(Tensor4d&, Vector4d&);
85 : public:
86 : static void registerKeywords( Keywords& keys );
87 : explicit QuaternionBondProductMatrix(const ActionOptions&);
88 : unsigned getNumberOfDerivatives();
89 : void prepare() override ;
90 : void calculate() override ;
91 : void applyNonZeroRankForces( std::vector<double>& outforces ) override ;
92 : void getNumberOfTasks( unsigned& ntasks ) override ;
93 : std::vector<unsigned>& getListOfActiveTasks( ActionWithVector* action ) override ;
94 : static void performTask( std::size_t task_index,
95 : const QuatBondProdMatInput& actiondata,
96 : const ParallelActionsInput& input,
97 : ParallelActionsOutput& output );
98 : static int getNumberOfValuesPerTask( std::size_t task_index, const QuatBondProdMatInput& actiondata );
99 : static void getForceIndices( std::size_t task_index, std::size_t colno, std::size_t ntotal_force, const QuatBondProdMatInput& actiondata, const ParallelActionsInput& input, ForceIndexHolder force_indices );
100 : };
101 :
102 : PLUMED_REGISTER_ACTION(QuaternionBondProductMatrix,"QUATERNION_BOND_PRODUCT_MATRIX")
103 :
104 :
105 : //const Vector4d& QuaternionBondMatrix::rightMultiply(Tensor4d& pref, Vector4d& quat) {
106 : // Vector4d temp;
107 : // int sumTemp;
108 : // for (int i=0; i<4; i++){ //rows
109 : // sumTemp=0;
110 : // for (int j=0; j<4; j++){ //cols
111 : // sumTemp+=pref(i,j)*quat[j];
112 : // }
113 : // temp[i]=sumTemp;
114 : // }
115 : // return temp;
116 : //}
117 :
118 :
119 :
120 :
121 8 : void QuaternionBondProductMatrix::registerKeywords( Keywords& keys ) {
122 8 : ActionWithVector::registerKeywords(keys);
123 16 : keys.addInputKeyword("compulsory","ARG","vector/matrix","this action takes 8 arguments. The first four should be the w,i,j and k components of a quaternion vector. The second four should be contact matrix and the matrices should be the x, y and z components of the bond vectors");
124 8 : PTM::registerKeywords( keys );
125 16 : keys.addOutputComponent("w","default","matrix","the real component of quaternion");
126 16 : keys.addOutputComponent("i","default","matrix","the i component of the quaternion");
127 16 : keys.addOutputComponent("j","default","matrix","the j component of the quaternion");
128 16 : keys.addOutputComponent("k","default","matrix","the k component of the quaternion");
129 8 : }
130 :
131 5 : QuaternionBondProductMatrix::QuaternionBondProductMatrix(const ActionOptions&ao):
132 : Action(ao),
133 : ActionWithVector(ao),
134 5 : taskmanager(this) {
135 5 : if( getNumberOfArguments()!=8 ) {
136 0 : error("should be eight arguments to this action, 4 quaternion components and 4 matrices");
137 : }
138 5 : unsigned nquat = getPntrToArgument(0)->getNumberOfValues();
139 25 : for(unsigned i=0; i<4; ++i) {
140 : Value* myarg=getPntrToArgument(i);
141 20 : if( myarg->getRank()!=1 ) {
142 0 : error("first four arguments to this action should be vectors");
143 : }
144 20 : if( (myarg->getPntrToAction())->getName()!="QUATERNION_VECTOR" ) {
145 0 : error("first four arguments to this action should be quaternions");
146 : }
147 20 : std::string mylab=getPntrToArgument(i)->getName();
148 20 : std::size_t dot=mylab.find_first_of(".");
149 30 : if( i==0 && mylab.substr(dot+1)!="w" ) {
150 0 : error("quaternion arguments are in wrong order");
151 : }
152 30 : if( i==1 && mylab.substr(dot+1)!="i" ) {
153 0 : error("quaternion arguments are in wrong order");
154 : }
155 30 : if( i==2 && mylab.substr(dot+1)!="j" ) {
156 0 : error("quaternion arguments are in wrong order");
157 : }
158 30 : if( i==3 && mylab.substr(dot+1)!="k" ) {
159 0 : error("quaternion arguments are in wrong order");
160 : }
161 : }
162 5 : std::vector<std::size_t> shape( getPntrToArgument(4)->getShape() );
163 25 : for(unsigned i=4; i<8; ++i) {
164 : Value* myarg=getPntrToArgument(i);
165 20 : if( myarg->getRank()!=2 ) {
166 0 : error("second four arguments to this action should be matrices");
167 : }
168 20 : if( myarg->getShape()[0]!=shape[0] || myarg->getShape()[1]!=shape[1] ) {
169 0 : error("matrices should all have the same shape");
170 : }
171 20 : if( myarg->getShape()[0]!=nquat ) {
172 0 : error("number of rows in matrix should equal number of input quaternions");
173 : }
174 20 : std::string mylab=getPntrToArgument(i)->getName();
175 20 : std::size_t dot=mylab.find_first_of(".");
176 30 : if( i==5 && mylab.substr(dot+1)!="x" ) {
177 0 : error("quaternion arguments are in wrong order");
178 : }
179 30 : if( i==6 && mylab.substr(dot+1)!="y" ) {
180 0 : error("quaternion arguments are in wrong order");
181 : }
182 30 : if( i==7 && mylab.substr(dot+1)!="z" ) {
183 0 : error("quaternion arguments are in wrong order");
184 : }
185 : }
186 : // We now swap round the order of the arguments to make gathering of forces easier for parallel task manager
187 5 : std::vector<Value*> args( getArguments() ), newargs;
188 25 : for(unsigned i=0; i<4; ++i) {
189 20 : newargs.push_back(args[4+i]);
190 : }
191 25 : for(unsigned i=0; i<4; ++i) {
192 20 : newargs.push_back(args[i]);
193 : }
194 5 : requestArguments( newargs );
195 :
196 5 : addComponent( "w", shape );
197 5 : componentIsNotPeriodic("w");
198 5 : addComponent( "i", shape );
199 5 : componentIsNotPeriodic("i");
200 5 : addComponent( "j", shape );
201 5 : componentIsNotPeriodic("j");
202 5 : addComponent( "k", shape );
203 5 : componentIsNotPeriodic("k");
204 : std::size_t nonthreadsafe = 0;
205 25 : for(unsigned i=4; i<8; ++i) {
206 20 : nonthreadsafe += getPntrToArgument(i)->getNumberOfStoredValues();
207 : }
208 5 : taskmanager.setupParallelTaskManager( 8, nonthreadsafe );
209 : taskmanager.setActionInput( QuatBondProdMatInput() );
210 5 : }
211 :
212 64 : unsigned QuaternionBondProductMatrix::getNumberOfDerivatives() {
213 : unsigned nder=0;
214 576 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
215 512 : nder += getPntrToArgument(i)->getNumberOfStoredValues();
216 : }
217 64 : return nder;
218 : }
219 :
220 27 : void QuaternionBondProductMatrix::prepare() {
221 27 : ActionWithVector::prepare();
222 27 : active_tasks.resize(0);
223 27 : }
224 :
225 5 : void QuaternionBondProductMatrix::getNumberOfTasks( unsigned& ntasks ) {
226 5 : ntasks=getPntrToComponent(0)->getNumberOfStoredValues();
227 5 : }
228 :
229 54 : std::vector<unsigned>& QuaternionBondProductMatrix::getListOfActiveTasks( ActionWithVector* action ) {
230 54 : if( active_tasks.size()>0 ) {
231 27 : return active_tasks;
232 : }
233 :
234 : unsigned atsize = 0;
235 : Value* myarg = getPntrToArgument(0);
236 27 : unsigned nrows = myarg->getShape()[0];
237 1668 : for(unsigned i=0; i<nrows; ++i) {
238 3282 : atsize += myarg->getRowLength(i);
239 : }
240 27 : active_tasks.resize( atsize );
241 :
242 : unsigned base=0, k=0;
243 1668 : for(unsigned i=0; i<nrows; ++i) {
244 1641 : unsigned ncols = myarg->getRowLength(i);
245 383061 : for(unsigned j=0; j<ncols; ++j) {
246 381420 : active_tasks[k] = base+j;
247 381420 : ++k;
248 : }
249 1641 : base += myarg->getNumberOfColumns();
250 : }
251 : return active_tasks;
252 : }
253 :
254 762840 : void QuaternionBondProductMatrix::performTask( std::size_t task_index,
255 : const QuatBondProdMatInput& /*actiondata*/,
256 : const ParallelActionsInput& input,
257 : ParallelActionsOutput& output ) {
258 :
259 : View2D<double, 4, 8> derivatives( output.derivatives.data() );
260 762840 : std::size_t index1 = std::floor( task_index / input.ncols[0] );
261 :
262 762840 : std::array<Tensor4d,2> dqt; //dqt[0] -> derivs w.r.t quat [dwt/dw1 dwt/di1 dwt/dj1 dwt/dk1]
263 : //[dit/dw1 dit/di1 dit/dj1 dit/dk1] etc, and dqt[1] is w.r.t the vector-turned-quaternion called bond
264 :
265 : // Retrieve the quaternion
266 : const std::array<double,4> quat{
267 762840 : input.inputdata[ input.argstarts[4+0] + index1 ],
268 762840 : input.inputdata[ input.argstarts[4+1] + index1 ],
269 762840 : input.inputdata[ input.argstarts[4+2] + index1 ],
270 762840 : input.inputdata[ input.argstarts[4+3] + index1 ]
271 762840 : };
272 : // Retrieve the components of the matrix
273 762840 : const std::array<double,4> bond{
274 : 0.0,
275 762840 : input.inputdata[ input.argstarts[1] + task_index ],
276 762840 : input.inputdata[ input.argstarts[2] + task_index ],
277 762840 : input.inputdata[ input.argstarts[3] + task_index ]
278 762840 : };
279 :
280 : //make a conjugate of q1 my own sanity
281 762840 : Vector4d quat_conj{ quat[0],-quat[1],-quat[2],-quat[3]};
282 : //declaring some signs to simplify the loop
283 762840 : constexpr std::array<double,4> pref_i{1,1,1,-1};
284 762840 : constexpr std::array<double,4> pref2_i{1,1,-1,1};
285 :
286 762840 : constexpr std::array<double,4> pref_j{1,-1,1,1};
287 762840 : constexpr std::array<double,4> pref2_j{1,1,1,-1};
288 :
289 762840 : constexpr std::array<double,4> pref_k{1,1,-1,1};
290 762840 : constexpr std::array<double,4> pref2_k{1,-1,1,1};
291 :
292 762840 : constexpr std::array<double,4> conj{1,-1,-1,-1};
293 762840 : std::array<double,4> quatTemp{0.0,0.0,0.0,0.0};
294 : //q1_conj * r first, while keep track of derivs
295 3814200 : for(unsigned i=0; i<4; ++i) {
296 : //real part of q1*q2
297 3051360 : quatTemp[0]+= quat[i]*bond[i];//conj*pref 1*1 for 0, -1 * -1 for 1++
298 3051360 : dqt[0](0,i) = bond[i];//conj*pref 1*1 for 0, -1 * -1 for 1++
299 3051360 : dqt[1](0,i) = quat[i];//conj*pref 1*1 for 0, -1 * -1 for 1++
300 : //i component
301 3051360 : quatTemp[1]+= pref_i[i]*quat_conj[i]*bond[(5-i)%4];
302 3051360 : dqt[0](1,i) = pref_i[i]*conj[i]*bond[(5-i)%4];
303 3051360 : dqt[1](1,i) = pref2_i[i]*quat_conj[(5-i)%4];
304 : //j component
305 3051360 : quatTemp[2]+= pref_j[i]*quat_conj[i]*bond[(i+2)%4];
306 3051360 : dqt[0](2,i) = pref_j[i]*conj[i]*bond[(i+2)%4];
307 3051360 : dqt[1](2,i) = pref2_j[i]*quat_conj[(i+2)%4];
308 : //k component
309 3051360 : quatTemp[3]+= pref_k[i]*quat_conj[i]*bond[3-i];
310 3051360 : dqt[0](3,i) = pref_k[i]*conj[i]*bond[3-i];
311 3051360 : dqt[1](3,i) = pref2_k[i]*quat_conj[3-i];
312 : }
313 : //now previous ^ product times quat again, not conjugated
314 :
315 : // calculate normalization factor
316 762840 : const double normFac = (bond[1] == 0.0 && bond[2]==0.0 && bond[3]==0) ?
317 : //just for the case where im comparing a quat to itself, itll be 0 at the end anyway
318 762840 : 1: 1/sqrt(bond[1]*bond[1] + bond[2]*bond[2] + bond[3]*bond[3]);
319 : //I hold off on normalizing because this can be done at the very end, and it
320 : // makes the derivatives with respect to 'bond' more simple
321 :
322 : double wf=0,xf=0,yf=0,zf=0;
323 :
324 3814200 : for(unsigned i=0; i<4; ++i) {
325 : //real part of q1*q2
326 3051360 : output.values[0] += normFac*conj[i]*quatTemp[i]*quat[i];
327 : wf+=normFac*conj[i]*quatTemp[i]*quat[i];
328 : //i component
329 3051360 : output.values[1] += normFac*pref_i[i]*quatTemp[i]*quat[(5-i)%4];
330 3051360 : xf+=normFac*pref_i[i]*quatTemp[i]*quat[(5-i)%4];
331 : //j component
332 3051360 : output.values[2] += normFac*pref_j[i]*quatTemp[i]*quat[(i+2)%4];
333 3051360 : yf+=normFac*pref_j[i]*quatTemp[i]*quat[(i+2)%4];
334 : //k component
335 3051360 : output.values[3] += normFac*pref_k[i]*quatTemp[i]*quat[(3-i)];
336 3051360 : zf+=normFac*pref_k[i]*quatTemp[i]*quat[(3-i)];
337 : }
338 : //had to split because bond's derivatives depend on the value of the overall quaternion component
339 762840 : if( input.noderiv ) {
340 381420 : return ;
341 : }
342 381420 : const auto normFacSQ=normFac*normFac;
343 381420 : xf*=normFacSQ;
344 381420 : yf*=normFacSQ;
345 381420 : zf*=normFacSQ;
346 : {
347 : //I unrolled the loop to opt out the ifs
348 : //I copy pasted the first component and constexpr'd the index
349 : constexpr int i = 0;
350 : //real part of q1*q2
351 381420 : derivatives[0][4+i] = (dotProduct(quat_conj, dqt[0].getCol(i)) + conj[i]*quatTemp[i])*normFac;
352 : //i component
353 381420 : derivatives[1][4+i] = (dotProduct(Vector4d(quat[1],quat[0],quat[3],-quat[2]),
354 381420 : dqt[0].getCol(i)) + pref2_i[i]*quatTemp[(5-i)%4])*normFac;
355 : //j component
356 381420 : derivatives[2][4+i] = (dotProduct(Vector4d(quat[2],-quat[3],quat[0],quat[1]),
357 381420 : dqt[0].getCol(i)) + pref2_j[i]*quatTemp[(i+2)%4])*normFac;
358 : //k component
359 381420 : derivatives[3][4+i] = (dotProduct(Vector4d(quat[3],quat[2],-quat[1],quat[0]),
360 381420 : dqt[0].getCol(i)) + pref2_k[i]*quatTemp[(3-i)])*normFac;
361 : }
362 1525680 : for(unsigned i=1; i<4; ++i) {
363 : //real part of q1*q2
364 1144260 : derivatives[0][4+i] = (dotProduct(quat_conj, dqt[0].getCol(i)) + conj[i]*quatTemp[i])*normFac;
365 1144260 : derivatives[0][i] = dotProduct(quat_conj, dqt[1].getCol(i))*normFac;
366 : //i component
367 1144260 : derivatives[1][4+i] =(dotProduct(Vector4d(quat[1],quat[0],quat[3],-quat[2]),
368 1144260 : dqt[0].getCol(i)) + pref2_i[i]*quatTemp[(5-i)%4])*normFac;
369 1144260 : derivatives[1][i] = dotProduct(Vector4d(quat[1],quat[0],quat[3],-quat[2]),
370 1144260 : dqt[1].getCol(i))*normFac
371 1144260 : +(-bond[i]*xf);
372 : //j component
373 1144260 : derivatives[2][4+i] = (dotProduct(Vector4d(quat[2],-quat[3],quat[0],quat[1]),
374 1144260 : dqt[0].getCol(i)) + pref2_j[i]*quatTemp[(i+2)%4])*normFac;
375 1144260 : derivatives[2][i] = dotProduct(Vector4d(quat[2],-quat[3],quat[0],quat[1]),
376 1144260 : dqt[1].getCol(i))*normFac+(-bond[i]*yf);
377 : //k component
378 1144260 : derivatives[3][4+i] = (dotProduct(Vector4d(quat[3],quat[2],-quat[1],quat[0]),
379 1144260 : dqt[0].getCol(i)) + pref2_k[i]*quatTemp[(3-i)])*normFac;
380 1144260 : derivatives[3][i] = dotProduct(Vector4d(quat[3],quat[2],-quat[1],quat[0]),
381 1144260 : dqt[1].getCol(i))*normFac+(-bond[i]*zf);
382 : }
383 : }
384 :
385 27 : void QuaternionBondProductMatrix::applyNonZeroRankForces( std::vector<double>& outforces ) {
386 27 : taskmanager.applyForces( outforces );
387 27 : }
388 :
389 381420 : int QuaternionBondProductMatrix::getNumberOfValuesPerTask( std::size_t task_index, const QuatBondProdMatInput& actiondata ) {
390 381420 : return 1;
391 : }
392 :
393 381420 : void QuaternionBondProductMatrix::getForceIndices( std::size_t task_index,
394 : std::size_t colno,
395 : std::size_t ntotal_force,
396 : const QuatBondProdMatInput& actiondata,
397 : const ParallelActionsInput& input,
398 : ForceIndexHolder force_indices ) {
399 381420 : std::size_t nquat = input.shapedata[input.shapestarts[4]];
400 381420 : std::size_t index1 = std::floor( task_index / input.ncols[0] );
401 1907100 : for(unsigned i=0; i<4; ++i) {
402 7628400 : for(unsigned j=0; j<4; ++j) {
403 6102720 : force_indices.indices[i][j] = input.argstarts[j] + task_index;
404 : }
405 1525680 : force_indices.threadsafe_derivatives_end[i] = 4;
406 7628400 : for(unsigned j=0; j<4; ++j) {
407 6102720 : force_indices.indices[i][j+4] = j*nquat + index1;
408 : }
409 1525680 : force_indices.tot_indices[i] = 8;
410 : }
411 381420 : }
412 :
413 27 : void QuaternionBondProductMatrix::calculate() {
414 : // Copy bookeeping arrays from input matrices to output matrices
415 135 : for(unsigned i=0; i<4; ++i) {
416 108 : getPntrToComponent(i)->copyBookeepingArrayFromArgument( getPntrToArgument(i) );
417 : }
418 27 : taskmanager.runAllTasks();
419 27 : }
420 :
421 : }
422 : }
|
