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