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 "tools/Torsion.h"
25 :
26 :
27 : #include <iostream>
28 :
29 : namespace PLMD {
30 : namespace crystdistrib {
31 :
32 : //+PLUMEDOC MCOLVAR QUATERNION_BOND_PRODUCT_MATRIX
33 : /*
34 : Calculate the product between a matrix of quaternions and the bonds connecting molecules
35 :
36 : 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:
37 :
38 : $$
39 : \mathbf{r} = \overline{\mathbf{q_1}} \mathbf{r_{21}} \mathbf{q_1}
40 : $$
41 :
42 : 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.
43 :
44 : In this example, 3 quaternion frames are calculated, and multiplied element-wise onto a distance matrix, yielding 9 vectors overall.
45 :
46 : ```plumed
47 : #calculate the quaternion frames for 3 molecules
48 : quat: QUATERNION ATOMS1=1,2,3 ATOMS2=4,5,6 ATOMS3=7,8,9
49 : #also find the distance between the 3 origins of the molecule frames
50 : c1: DISTANCE_MATRIX GROUP=1,4,7 CUTOFF=100.0 COMPONENTS
51 : qp: QUATERNION_BOND_PRODUCT_MATRIX ARG=quat.*,c1.*
52 : #this is now a matrix showing how each molecule is oriented in 3D space
53 : #relative to eachother's origins
54 : #now use in a simple collective variable
55 : ops: CUSTOM ARG=qp.* VAR=w,i,j,k FUNC=w+i+j+k PERIODIC=NO
56 : #w could have been ignored because it is always zero.
57 :
58 : ```
59 :
60 : */
61 : //+ENDPLUMEDOC
62 :
63 : class QuaternionBondProductMatrix : public ActionWithMatrix {
64 : private:
65 : unsigned nderivatives;
66 : std::vector<bool> stored;
67 : // const Vector4d& rightMultiply(Tensor4d&, Vector4d&);
68 : public:
69 : static void registerKeywords( Keywords& keys );
70 : explicit QuaternionBondProductMatrix(const ActionOptions&);
71 : unsigned getNumberOfDerivatives();
72 : unsigned getNumberOfColumns() const override ;
73 : void setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const ;
74 : void performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const override;
75 : void runEndOfRowJobs( const unsigned& ival, const std::vector<unsigned> & indices, MultiValue& myvals ) const override ;
76 : };
77 :
78 : PLUMED_REGISTER_ACTION(QuaternionBondProductMatrix,"QUATERNION_BOND_PRODUCT_MATRIX")
79 :
80 :
81 : //const Vector4d& QuaternionBondMatrix::rightMultiply(Tensor4d& pref, Vector4d& quat) {
82 : // Vector4d temp;
83 : // int sumTemp;
84 : // for (int i=0; i<4; i++){ //rows
85 : // sumTemp=0;
86 : // for (int j=0; j<4; j++){ //cols
87 : // sumTemp+=pref(i,j)*quat[j];
88 : // }
89 : // temp[i]=sumTemp;
90 : // }
91 : // return temp;
92 : //}
93 :
94 :
95 :
96 :
97 7 : void QuaternionBondProductMatrix::registerKeywords( Keywords& keys ) {
98 7 : ActionWithMatrix::registerKeywords(keys);
99 14 : 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");
100 14 : keys.addOutputComponent("w","default","matrix","the real component of quaternion");
101 14 : keys.addOutputComponent("i","default","matrix","the i component of the quaternion");
102 14 : keys.addOutputComponent("j","default","matrix","the j component of the quaternion");
103 14 : keys.addOutputComponent("k","default","matrix","the k component of the quaternion");
104 7 : }
105 :
106 4 : QuaternionBondProductMatrix::QuaternionBondProductMatrix(const ActionOptions&ao):
107 : Action(ao),
108 4 : ActionWithMatrix(ao) {
109 4 : if( getNumberOfArguments()!=8 ) {
110 0 : error("should be eight arguments to this action, 4 quaternion components and 4 matrices");
111 : }
112 4 : unsigned nquat = getPntrToArgument(0)->getNumberOfValues();
113 20 : for(unsigned i=0; i<4; ++i) {
114 : Value* myarg=getPntrToArgument(i);
115 16 : myarg->buildDataStore();
116 16 : if( myarg->getRank()!=1 ) {
117 0 : error("first four arguments to this action should be vectors");
118 : }
119 16 : if( (myarg->getPntrToAction())->getName()!="QUATERNION_VECTOR" ) {
120 0 : error("first four arguments to this action should be quaternions");
121 : }
122 16 : std::string mylab=getPntrToArgument(i)->getName();
123 16 : std::size_t dot=mylab.find_first_of(".");
124 24 : if( i==0 && mylab.substr(dot+1)!="w" ) {
125 0 : error("quaternion arguments are in wrong order");
126 : }
127 24 : if( i==1 && mylab.substr(dot+1)!="i" ) {
128 0 : error("quaternion arguments are in wrong order");
129 : }
130 24 : if( i==2 && mylab.substr(dot+1)!="j" ) {
131 0 : error("quaternion arguments are in wrong order");
132 : }
133 24 : if( i==3 && mylab.substr(dot+1)!="k" ) {
134 0 : error("quaternion arguments are in wrong order");
135 : }
136 : }
137 4 : std::vector<unsigned> shape( getPntrToArgument(4)->getShape() );
138 20 : for(unsigned i=4; i<8; ++i) {
139 : Value* myarg=getPntrToArgument(i);
140 16 : if( myarg->getRank()!=2 ) {
141 0 : error("second four arguments to this action should be matrices");
142 : }
143 16 : if( myarg->getShape()[0]!=shape[0] || myarg->getShape()[1]!=shape[1] ) {
144 0 : error("matrices should all have the same shape");
145 : }
146 16 : if( myarg->getShape()[0]!=nquat ) {
147 0 : error("number of rows in matrix should equal number of input quaternions");
148 : }
149 16 : std::string mylab=getPntrToArgument(i)->getName();
150 16 : std::size_t dot=mylab.find_first_of(".");
151 24 : if( i==5 && mylab.substr(dot+1)!="x" ) {
152 0 : error("quaternion arguments are in wrong order");
153 : }
154 24 : if( i==6 && mylab.substr(dot+1)!="y" ) {
155 0 : error("quaternion arguments are in wrong order");
156 : }
157 24 : if( i==7 && mylab.substr(dot+1)!="z" ) {
158 0 : error("quaternion arguments are in wrong order");
159 : }
160 : }
161 4 : addComponent( "w", shape );
162 4 : componentIsNotPeriodic("w");
163 4 : addComponent( "i", shape );
164 4 : componentIsNotPeriodic("i");
165 4 : addComponent( "j", shape );
166 4 : componentIsNotPeriodic("j");
167 4 : addComponent( "k", shape );
168 4 : componentIsNotPeriodic("k");
169 4 : done_in_chain=true;
170 4 : nderivatives = buildArgumentStore(0);
171 :
172 4 : std::string headstr=getFirstActionInChain()->getLabel();
173 4 : stored.resize( getNumberOfArguments() );
174 36 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
175 32 : stored[i] = getPntrToArgument(i)->ignoreStoredValue( headstr );
176 : }
177 4 : }
178 :
179 32 : unsigned QuaternionBondProductMatrix::getNumberOfDerivatives() {
180 32 : return nderivatives;
181 : }
182 :
183 96 : unsigned QuaternionBondProductMatrix::getNumberOfColumns() const {
184 96 : const ActionWithMatrix* am=dynamic_cast<const ActionWithMatrix*>( getPntrToArgument(4)->getPntrToAction() );
185 96 : plumed_assert( am );
186 96 : return am->getNumberOfColumns();
187 : }
188 :
189 0 : void QuaternionBondProductMatrix::setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const {
190 0 : unsigned start_n = getPntrToArgument(4)->getShape()[0], size_v = getPntrToArgument(4)->getShape()[1];
191 0 : if( indices.size()!=size_v+1 ) {
192 0 : indices.resize( size_v+1 );
193 : }
194 0 : for(unsigned i=0; i<size_v; ++i) {
195 0 : indices[i+1] = start_n + i;
196 : }
197 : myvals.setSplitIndex( size_v + 1 );
198 0 : }
199 :
200 381366 : void QuaternionBondProductMatrix::performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const {
201 381366 : unsigned ind2=index2;
202 381366 : if( index2>=getPntrToArgument(0)->getShape()[0] ) {
203 0 : ind2 = index2 - getPntrToArgument(0)->getShape()[0];
204 : }
205 :
206 381366 : std::vector<double> quat(4), bond(4), quatTemp(4);
207 381366 : std::vector<Tensor4d> dqt(2); //dqt[0] -> derivs w.r.t quat [dwt/dw1 dwt/di1 dwt/dj1 dwt/dk1]
208 : //[dit/dw1 dit/di1 dit/dj1 dit/dk1] etc, and dqt[1] is w.r.t the vector-turned-quaternion called bond
209 :
210 : // Retrieve the quaternion
211 1906830 : for(unsigned i=0; i<4; ++i) {
212 1525464 : quat[i] = getArgumentElement( i, index1, myvals );
213 : }
214 :
215 : // Retrieve the components of the matrix
216 381366 : double weight = getElementOfMatrixArgument( 4, index1, ind2, myvals );
217 1525464 : for(unsigned i=1; i<4; ++i) {
218 1144098 : bond[i] = getElementOfMatrixArgument( 4+i, index1, ind2, myvals );
219 : }
220 :
221 : // calculate normalization factor
222 381366 : bond[0]=0.0;
223 381366 : double normFac = 1/sqrt(bond[1]*bond[1] + bond[2]*bond[2] + bond[3]*bond[3]);
224 381366 : if (bond[1] == 0.0 && bond[2]==0.0 && bond[3]==0) {
225 : normFac=1; //just for the case where im comparing a quat to itself, itll be 0 at the end anyway
226 : }
227 : double normFac3 = normFac*normFac*normFac;
228 : //I hold off on normalizing because this can be done at the very end, and it makes the derivatives with respect to 'bond' more simple
229 :
230 :
231 :
232 381366 : std::vector<double> quat_conj(4);
233 381366 : quat_conj[0] = quat[0];
234 381366 : quat_conj[1] = -1*quat[1];
235 381366 : quat_conj[2] = -1*quat[2];
236 381366 : quat_conj[3] = -1*quat[3];
237 : //make a conjugate of q1 my own sanity
238 :
239 :
240 :
241 :
242 : //q1_conj * r first, while keep track of derivs
243 : double pref=1;
244 : double conj=1;
245 : double pref2=1;
246 : //real part of q1*q2
247 :
248 1906830 : for(unsigned i=0; i<4; ++i) {
249 1525464 : if( i>0 ) {
250 : pref=-1;
251 : conj=-1;
252 : pref2=-1;
253 : }
254 1525464 : quatTemp[0]+=pref*quat_conj[i]*bond[i];
255 1525464 : dqt[0](0,i) = conj*pref*bond[i];
256 1525464 : dqt[1](0,i) = pref2*quat_conj[i];
257 : //addDerivativeOnVectorArgument( false, 0, i, index1, conj*pref*bond[i], myvals );
258 : //addDerivativeOnVectorArgument( false, 0, 4+i, ind2, conj*pref*quat[i], myvals );
259 : }
260 : //i component
261 : pref=1;
262 : conj=1;
263 : pref2=1;
264 :
265 1906830 : for (unsigned i=0; i<4; i++) {
266 1525464 : if(i==3) {
267 : pref=-1;
268 : } else {
269 : pref=1;
270 : }
271 1525464 : if(i==2) {
272 : pref2=-1;
273 : } else {
274 : pref2=1;
275 : }
276 1525464 : if (i>0) {
277 : conj=-1;
278 : }
279 :
280 1525464 : quatTemp[1]+=pref*quat_conj[i]*bond[(5-i)%4];
281 1525464 : dqt[0](1,i) =conj*pref*bond[(5-i)%4];
282 1525464 : dqt[1](1,i) = pref2*quat_conj[(5-i)%4];
283 : //addDerivativeOnVectorArgument( false, 1, i, index1, conj*pref*bond[(5-i)%4], myvals );
284 : //addDerivativeOnVectorArgument( false, 1, 4+i, ind2, conj*pref*quat[i], myvals );
285 : }
286 :
287 : //j component
288 : pref=1;
289 : pref2=1;
290 : conj=1;
291 :
292 1906830 : for (unsigned i=0; i<4; i++) {
293 1525464 : if(i==1) {
294 : pref=-1;
295 : } else {
296 : pref=1;
297 : }
298 1525464 : if (i==3) {
299 : pref2=-1;
300 : } else {
301 : pref2=1;
302 : }
303 1525464 : if (i>0) {
304 : conj=-1;
305 : }
306 :
307 1525464 : quatTemp[2]+=pref*quat_conj[i]*bond[(i+2)%4];
308 1525464 : dqt[0](2,i)=conj*pref*bond[(i+2)%4];
309 1525464 : dqt[1](2,i)=pref2*quat_conj[(i+2)%4];
310 : //addDerivativeOnVectorArgument( false, 2, i, index1, conj*pref*bond[(i+2)%4], myvals );
311 : //addDerivativeOnVectorArgument( false, 2, 4+i, ind2, conj*pref*quat[i], myvals );
312 : }
313 :
314 : //k component
315 : pref=1;
316 : pref2=1;
317 : conj=1;
318 :
319 1906830 : for (unsigned i=0; i<4; i++) {
320 1525464 : if(i==2) {
321 : pref=-1;
322 : } else {
323 : pref=1;
324 : }
325 1525464 : if(i==1) {
326 : pref2=-1;
327 : } else {
328 : pref2=1;
329 : }
330 1525464 : if(i>0) {
331 : conj=-1;
332 : }
333 1525464 : quatTemp[3]+=pref*quat_conj[i]*bond[(3-i)];
334 1525464 : dqt[0](3,i)=conj*pref*bond[3-i];
335 1525464 : dqt[1](3,i)= pref2*quat_conj[3-i];
336 : //addDerivativeOnVectorArgument( false, 3, i, index1, conj*pref*bond[3-i], myvals );
337 : //addDerivativeOnVectorArgument( false, 3, 4+i, ind2, conj*pref*quat[i], myvals );
338 :
339 : }
340 :
341 :
342 : //now previous ^ product times quat again, not conjugated
343 : //real part of q1*q2
344 381366 : double tempDot=0,wf=0,xf=0,yf=0,zf=0;
345 : pref=1;
346 : pref2=1;
347 1906830 : for(unsigned i=0; i<4; ++i) {
348 1525464 : if( i>0 ) {
349 : pref=-1;
350 : pref2=-1;
351 : }
352 1525464 : myvals.addValue( getConstPntrToComponent(0)->getPositionInStream(), normFac*pref*quatTemp[i]*quat[i] );
353 : wf+=normFac*pref*quatTemp[i]*quat[i];
354 1525464 : if( doNotCalculateDerivatives() ) {
355 0 : continue ;
356 : }
357 1525464 : tempDot=(dotProduct(Vector4d(quat[0],-quat[1],-quat[2],-quat[3]), dqt[0].getCol(i)) + pref2*quatTemp[i])*normFac;
358 1525464 : addDerivativeOnVectorArgument( stored[i], 0, i, index1, tempDot, myvals);
359 : }
360 : //had to split because bond's derivatives depend on the value of the overall quaternion component
361 : //addDerivativeOnMatrixArgument( false, 0, 4, index1, ind2, 0.0, myvals );
362 1906830 : for(unsigned i=0; i<4; ++i) {
363 1525464 : tempDot=dotProduct(Vector4d(quat[0],-quat[1],-quat[2],-quat[3]), dqt[1].getCol(i))*normFac;
364 1525464 : if (i!=0 ) {
365 1144098 : addDerivativeOnMatrixArgument( stored[4+i], 0, 4+i, index1, ind2, tempDot, myvals );
366 : } else {
367 381366 : addDerivativeOnMatrixArgument( stored[4+i], 0, 4+i, index1, ind2, 0.0, myvals );
368 : }
369 : }
370 : // for (unsigned i=0; i<4; ++i) {
371 : //myvals.addValue( getConstPntrToComponent(0)->getPositionInStream(), 0.0 );
372 : //if( doNotCalculateDerivatives() ) continue ;
373 : //addDerivativeOnVectorArgument( false, 0, i, index1, 0.0, myvals);
374 : //addDerivativeOnVectorArgument( false, 0, 4+i, ind2, 0.0 , myvals);
375 : // }
376 : //the w component should always be zero, barring some catastrophe, but we calculate it out anyway
377 :
378 : //i component
379 : pref=1;
380 : pref2=1;
381 1906830 : for (unsigned i=0; i<4; i++) {
382 1525464 : if(i==3) {
383 : pref=-1;
384 : } else {
385 : pref=1;
386 : }
387 1525464 : myvals.addValue( getConstPntrToComponent(1)->getPositionInStream(), normFac*pref*quatTemp[i]*quat[(5-i)%4]);
388 1525464 : xf+=normFac*pref*quatTemp[i]*quat[(5-i)%4];
389 1525464 : if(i==2) {
390 : pref2=-1;
391 : } else {
392 : pref2=1;
393 : }
394 1525464 : if( doNotCalculateDerivatives() ) {
395 0 : continue ;
396 : }
397 1525464 : tempDot=(dotProduct(Vector4d(quat[1],quat[0],quat[3],-quat[2]), dqt[0].getCol(i)) + pref2*quatTemp[(5-i)%4])*normFac;
398 1525464 : addDerivativeOnVectorArgument( stored[i], 1, i, index1, tempDot, myvals);
399 : }
400 : //addDerivativeOnMatrixArgument( false, 1, 4, index1, ind2, 0.0, myvals );
401 :
402 1906830 : for(unsigned i=0; i<4; ++i) {
403 1525464 : tempDot=dotProduct(Vector4d(quat[1],quat[0],quat[3],-quat[2]), dqt[1].getCol(i))*normFac;
404 1525464 : if (i!=0) {
405 1144098 : addDerivativeOnMatrixArgument( stored[4+i], 1, 4+i, index1, ind2, tempDot+(-bond[i]*normFac*normFac*xf), myvals );
406 : } else {
407 381366 : addDerivativeOnMatrixArgument( stored[4+i], 1, 4+i, index1, ind2, 0.0, myvals );
408 : }
409 :
410 : }
411 :
412 :
413 : //j component
414 : pref=1;
415 : pref2=1;
416 1906830 : for (unsigned i=0; i<4; i++) {
417 1525464 : if(i==1) {
418 : pref=-1;
419 : } else {
420 : pref=1;
421 : }
422 1525464 : if (i==3) {
423 : pref2=-1;
424 : } else {
425 : pref2=1;
426 : }
427 :
428 1525464 : myvals.addValue( getConstPntrToComponent(2)->getPositionInStream(), normFac*pref*quatTemp[i]*quat[(i+2)%4]);
429 1525464 : yf+=normFac*pref*quatTemp[i]*quat[(i+2)%4];
430 1525464 : if( doNotCalculateDerivatives() ) {
431 0 : continue ;
432 : }
433 1525464 : tempDot=(dotProduct(Vector4d(quat[2],-quat[3],quat[0],quat[1]), dqt[0].getCol(i)) + pref2*quatTemp[(i+2)%4])*normFac;
434 1525464 : addDerivativeOnVectorArgument( stored[i], 2, i, index1, tempDot, myvals);
435 : }
436 : // addDerivativeOnMatrixArgument( false, 2, 4, index1, ind2,0.0 , myvals );
437 :
438 1906830 : for(unsigned i=0; i<4; ++i) {
439 1525464 : tempDot=dotProduct(Vector4d(quat[2],-quat[3],quat[0],quat[1]), dqt[1].getCol(i))*normFac;
440 1525464 : if (i!=0) {
441 1144098 : addDerivativeOnMatrixArgument( stored[4+i], 2, 4+i, index1, ind2, tempDot+(-bond[i]*normFac*normFac*yf), myvals );
442 : } else {
443 381366 : addDerivativeOnMatrixArgument( stored[4+i], 2, 4+i, index1, ind2, 0.0, myvals );
444 : }
445 :
446 :
447 : }
448 :
449 : //k component
450 : pref=1;
451 : pref2=1;
452 1906830 : for (unsigned i=0; i<4; i++) {
453 1525464 : if(i==2) {
454 : pref=-1;
455 : } else {
456 : pref=1;
457 : }
458 1525464 : if(i==1) {
459 : pref2=-1;
460 : } else {
461 : pref2=1;
462 : }
463 :
464 1525464 : myvals.addValue( getConstPntrToComponent(3)->getPositionInStream(), normFac*pref*quatTemp[i]*quat[(3-i)]);
465 1525464 : zf+=normFac*pref*quatTemp[i]*quat[(3-i)];
466 1525464 : if( doNotCalculateDerivatives() ) {
467 0 : continue ;
468 : }
469 1525464 : tempDot=(dotProduct(Vector4d(quat[3],quat[2],-quat[1],quat[0]), dqt[0].getCol(i)) + pref2*quatTemp[(3-i)])*normFac;
470 1525464 : addDerivativeOnVectorArgument( stored[i], 3, i, index1, tempDot, myvals);
471 : }
472 : //addDerivativeOnMatrixArgument( false, 3, 4, index1, ind2, 0.0 , myvals );
473 :
474 1906830 : for(unsigned i=0; i<4; ++i) {
475 1525464 : tempDot=dotProduct(Vector4d(quat[3],quat[2],-quat[1],quat[0]), dqt[1].getCol(i))*normFac;
476 1525464 : if (i!=0) {
477 1144098 : addDerivativeOnMatrixArgument( stored[4+i], 3, 4+i, index1, ind2, tempDot+(-bond[i]*normFac*normFac*zf), myvals );
478 : } else {
479 381366 : addDerivativeOnMatrixArgument( stored[4+i], 3, 4+i, index1, ind2, 0.0, myvals );
480 : }
481 :
482 :
483 : }
484 381366 : if( doNotCalculateDerivatives() ) {
485 : return ;
486 : }
487 :
488 1906830 : for(unsigned outcomp=0; outcomp<4; ++outcomp) {
489 1525464 : unsigned ostrn = getConstPntrToComponent(outcomp)->getPositionInStream();
490 7627320 : for(unsigned i=4; i<8; ++i) {
491 : bool found=false;
492 6101856 : for(unsigned j=4; j<i; ++j) {
493 4576392 : if( arg_deriv_starts[i]==arg_deriv_starts[j] ) {
494 : found=true;
495 : break;
496 : }
497 : }
498 6101856 : if( found || !stored[i] ) {
499 4576392 : continue;
500 : }
501 :
502 : unsigned istrn = getPntrToArgument(i)->getPositionInStream();
503 24407424 : for(unsigned k=0; k<myvals.getNumberActive(istrn); ++k) {
504 22881960 : unsigned kind=myvals.getActiveIndex(istrn,k);
505 22881960 : myvals.updateIndex( ostrn, kind );
506 : }
507 : }
508 : }
509 : }
510 :
511 1614 : void QuaternionBondProductMatrix::runEndOfRowJobs( const unsigned& ival, const std::vector<unsigned> & indices, MultiValue& myvals ) const {
512 1614 : if( doNotCalculateDerivatives() || !matrixChainContinues() ) {
513 : return ;
514 : }
515 :
516 8070 : for(unsigned j=0; j<getNumberOfComponents(); ++j) {
517 6456 : unsigned nmat = getConstPntrToComponent(j)->getPositionInMatrixStash(), nmat_ind = myvals.getNumberOfMatrixRowDerivatives( nmat );
518 : std::vector<unsigned>& matrix_indices( myvals.getMatrixRowDerivativeIndices( nmat ) );
519 : unsigned ntwo_atoms = myvals.getSplitIndex();
520 : // Quaternion
521 32280 : for(unsigned k=0; k<4; ++k) {
522 25824 : matrix_indices[nmat_ind] = arg_deriv_starts[k] + ival;
523 25824 : nmat_ind++;
524 : }
525 : // Loop over row of matrix
526 32280 : for(unsigned n=4; n<8; ++n) {
527 : bool found=false;
528 25824 : for(unsigned k=4; k<n; ++k) {
529 19368 : if( arg_deriv_starts[k]==arg_deriv_starts[n] ) {
530 : found=true;
531 : break;
532 : }
533 : }
534 25824 : if( found ) {
535 : continue;
536 : }
537 : unsigned istrn = getPntrToArgument(n)->getPositionInMatrixStash();
538 : std::vector<unsigned>& imat_indices( myvals.getMatrixRowDerivativeIndices( istrn ) );
539 4660320 : for(unsigned k=0; k<myvals.getNumberOfMatrixRowDerivatives( istrn ); ++k) {
540 4653864 : matrix_indices[nmat_ind + k] = arg_deriv_starts[n] + imat_indices[k];
541 : }
542 6456 : nmat_ind += myvals.getNumberOfMatrixRowDerivatives( getPntrToArgument(4)->getPositionInMatrixStash() );
543 : }
544 : myvals.setNumberOfMatrixRowDerivatives( nmat, nmat_ind );
545 : }
546 : }
547 :
548 : }
549 : }
|
