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