1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38 package ffx.potential.nonbonded.implicit;
39
40 import edu.rit.pj.IntegerForLoop;
41 import edu.rit.pj.ParallelRegion;
42 import edu.rit.pj.ParallelTeam;
43 import edu.rit.pj.reduction.SharedDouble;
44 import edu.rit.pj.reduction.SharedInteger;
45 import ffx.crystal.Crystal;
46 import ffx.crystal.SymOp;
47 import ffx.numerics.atomic.AtomicDoubleArray;
48 import ffx.numerics.atomic.AtomicDoubleArray.AtomicDoubleArrayImpl;
49 import ffx.numerics.atomic.AtomicDoubleArray3D;
50 import ffx.numerics.multipole.GKEnergyQI;
51 import ffx.numerics.multipole.PolarizableMultipole;
52 import ffx.numerics.multipole.QIFrame;
53 import ffx.potential.bonded.Atom;
54 import ffx.potential.nonbonded.GeneralizedKirkwood.NonPolarModel;
55 import ffx.potential.nonbonded.pme.Polarization;
56
57 import java.util.List;
58 import java.util.logging.Level;
59 import java.util.logging.Logger;
60
61 import static ffx.numerics.multipole.GKSource.cn;
62 import static ffx.potential.parameters.MultipoleType.t000;
63 import static ffx.potential.parameters.MultipoleType.t001;
64 import static ffx.potential.parameters.MultipoleType.t002;
65 import static ffx.potential.parameters.MultipoleType.t010;
66 import static ffx.potential.parameters.MultipoleType.t011;
67 import static ffx.potential.parameters.MultipoleType.t020;
68 import static ffx.potential.parameters.MultipoleType.t100;
69 import static ffx.potential.parameters.MultipoleType.t101;
70 import static ffx.potential.parameters.MultipoleType.t110;
71 import static ffx.potential.parameters.MultipoleType.t200;
72 import static java.lang.String.format;
73 import static org.apache.commons.math3.util.FastMath.exp;
74 import static org.apache.commons.math3.util.FastMath.sqrt;
75
76
77
78
79
80
81
82 public class GKEnergyRegion extends ParallelRegion {
83
84 private static final Logger logger = Logger.getLogger(GKEnergyRegion.class.getName());
85
86 protected static final double oneThird = 1.0 / 3.0;
87
88 protected final double electric;
89
90 protected final Polarization polarization;
91 protected final NonPolarModel nonPolarModel;
92
93 private final double soluteDielectric;
94 private final double solventDielectric;
95
96
97
98
99
100
101
102 protected final double dOffset = 0.09;
103
104 protected final double surfaceTension;
105
106 protected final double gkc;
107
108 private final double fc;
109
110 private final double fd;
111
112 private final double fq;
113
114 private final double probe;
115 private final SharedDouble sharedPermanentGKEnergy;
116 private final SharedDouble sharedPolarizationGKEnergy;
117 private final SharedDouble sharedGKEnergy;
118 private final SharedInteger sharedInteractions;
119 private final IntegerForLoop[] gkEnergyLoop;
120
121 private Atom[] atoms;
122
123 private double[][][] inducedDipole;
124
125 private double[][][] inducedDipoleCR;
126
127 private double[][][] globalMultipole;
128
129 private Crystal crystal;
130
131 private double[][][] sXYZ;
132
133 private int[][][] neighborLists;
134
135 private boolean[] use = null;
136
137 private double cut2;
138
139 private double[] baseRadius;
140
141 private double[] born;
142
143 private boolean gradient = false;
144
145 private AtomicDoubleArray3D grad;
146
147 private AtomicDoubleArray3D torque;
148
149 private AtomicDoubleArray sharedBornGrad;
150
151 private AtomicDoubleArray selfEnergy;
152
153 private AtomicDoubleArray crossEnergy;
154
155
156 public GKEnergyRegion(
157 int nt, Polarization polarization, NonPolarModel nonPolarModel,
158 double surfaceTension, double probe, double electric,
159 double soluteDieletric, double solventDieletric, double gkc, boolean gkQI) {
160
161
162 this.soluteDielectric = soluteDieletric;
163 this.solventDielectric = solventDieletric;
164 fc = cn(0, soluteDieletric, solventDieletric);
165 fd = cn(1, soluteDieletric, solventDieletric);
166 fq = cn(2, soluteDieletric, solventDieletric);
167 this.gkc = gkc;
168
169 this.polarization = polarization;
170 this.nonPolarModel = nonPolarModel;
171 this.surfaceTension = surfaceTension;
172 this.probe = probe;
173
174 this.electric = electric;
175
176 gkEnergyLoop = new IntegerForLoop[nt];
177 for (int i = 0; i < nt; i++) {
178 if (gkQI) {
179 gkEnergyLoop[i] = new GKEnergyLoopQI();
180 } else {
181 gkEnergyLoop[i] = new GKEnergyLoop();
182 }
183 }
184 sharedPermanentGKEnergy = new SharedDouble();
185 sharedPolarizationGKEnergy = new SharedDouble();
186 sharedGKEnergy = new SharedDouble();
187 sharedInteractions = new SharedInteger();
188 }
189
190 public double getEnergy() {
191 return sharedGKEnergy.get();
192 }
193
194 public double getPermanentEnergy() {
195 return sharedPermanentGKEnergy.get();
196 }
197
198 public double getPolarizationEnergy() {
199 return sharedPolarizationGKEnergy.get();
200 }
201
202 public AtomicDoubleArray getSelfEnergy() {
203 int nAtoms = atoms.length;
204 selfEnergy.reduce(0, nAtoms - 1);
205 return selfEnergy;
206 }
207
208 public int getInteractions() {
209 return sharedInteractions.get();
210 }
211
212 public void init(
213 Atom[] atoms,
214 double[][][] globalMultipole,
215 double[][][] inducedDipole,
216 double[][][] inducedDipoleCR,
217 Crystal crystal,
218 double[][][] sXYZ,
219 int[][][] neighborLists,
220 boolean[] use,
221 double cut2,
222 double[] baseRadius,
223 double[] born,
224 boolean gradient,
225 ParallelTeam parallelTeam,
226 AtomicDoubleArray3D grad,
227 AtomicDoubleArray3D torque,
228 AtomicDoubleArray sharedBornGrad) {
229
230 this.atoms = atoms;
231 this.globalMultipole = globalMultipole;
232 this.inducedDipole = inducedDipole;
233 this.inducedDipoleCR = inducedDipoleCR;
234 this.crystal = crystal;
235 this.sXYZ = sXYZ;
236 this.neighborLists = neighborLists;
237 this.use = use;
238 this.cut2 = cut2;
239 this.baseRadius = baseRadius;
240 this.born = born;
241 this.gradient = gradient;
242
243 this.grad = grad;
244 this.torque = torque;
245 this.sharedBornGrad = sharedBornGrad;
246
247 int nAtoms = atoms.length;
248 int nThreads = gkEnergyLoop.length;
249 if (selfEnergy == null || selfEnergy.size() != atoms.length) {
250 selfEnergy = AtomicDoubleArray.atomicDoubleArrayFactory(
251 AtomicDoubleArrayImpl.MULTI, nThreads, nAtoms);
252 crossEnergy = AtomicDoubleArray.atomicDoubleArrayFactory(
253 AtomicDoubleArrayImpl.MULTI, nThreads, nAtoms);
254 } else {
255 selfEnergy.reset(parallelTeam, 0, nAtoms - 1);
256 crossEnergy.reset(parallelTeam, 0, nAtoms - 1);
257 }
258 }
259
260 @Override
261 public void run() {
262 try {
263 int nAtoms = atoms.length;
264 int threadIndex = getThreadIndex();
265 execute(0, nAtoms - 1, gkEnergyLoop[threadIndex]);
266 } catch (Exception e) {
267 String message = "Fatal exception computing GK Energy in thread " + getThreadIndex() + "\n";
268 logger.log(Level.SEVERE, message, e);
269 }
270 }
271
272 @Override
273 public void start() {
274 sharedPermanentGKEnergy.set(0.0);
275 sharedPolarizationGKEnergy.set(0.0);
276 sharedGKEnergy.set(0.0);
277 sharedInteractions.set(0);
278 }
279
280 @Override
281 public void finish() {
282 if (logger.isLoggable(Level.FINER)) {
283 int nAtoms = atoms.length;
284 selfEnergy.reduce(0, nAtoms - 1);
285 crossEnergy.reduce(0, nAtoms - 1);
286 logger.finer(" Generalized Kirkwood Self-Energies and Cross-Energies\n");
287 double selfSum = 0.0;
288 double crossSum = 0.0;
289 for (int i = 0; i < nAtoms; i++) {
290 double self = selfEnergy.get(i);
291 double cross = crossEnergy.get(i);
292 logger.finer(format("GKSELF %5d %16.8f %16.8f", i, self, cross));
293 selfSum += self;
294 crossSum += cross;
295 }
296 logger.finer(format(" %16.8f %16.8f %16.8f\n",
297 selfSum, crossSum, selfSum + crossSum));
298 }
299 }
300
301
302
303
304
305
306 private class GKEnergyLoop extends IntegerForLoop {
307
308 private final double[][] a;
309 private final double[][] b;
310 private final double[] gc;
311 private final double[] gux;
312 private final double[] guy;
313 private final double[] guz;
314 private final double[] gqxx;
315 private final double[] gqyy;
316 private final double[] gqzz;
317 private final double[] gqxy;
318 private final double[] gqxz;
319 private final double[] gqyz;
320 private final double[] dx_local;
321 private double ci, uxi, uyi, uzi, qxxi, qxyi, qxzi, qyyi, qyzi, qzzi;
322 private double ck, uxk, uyk, uzk, qxxk, qxyk, qxzk, qyyk, qyzk, qzzk;
323 private double dxi, dyi, dzi, pxi, pyi, pzi, sxi, syi, szi;
324 private double dxk, dyk, dzk, pxk, pyk, pzk, sxk, syk, szk;
325 private double xr, yr, zr, xr2, yr2, zr2, rbi, rbk;
326 private double xi, yi, zi;
327 private double dedxi, dedyi, dedzi;
328 private double dborni;
329 private double trqxi, trqyi, trqzi;
330 private int count;
331 private int iSymm;
332 private int threadID;
333 private final double[][] transOp;
334 private double gkPermanentEnergy;
335 private double gkPolarizationEnergy;
336 private double gkEnergy;
337
338 GKEnergyLoop() {
339 a = new double[6][4];
340 b = new double[5][3];
341 gc = new double[31];
342 gux = new double[31];
343 guy = new double[31];
344 guz = new double[31];
345 gqxx = new double[31];
346 gqyy = new double[31];
347 gqzz = new double[31];
348 gqxy = new double[31];
349 gqxz = new double[31];
350 gqyz = new double[31];
351 dx_local = new double[3];
352 transOp = new double[3][3];
353 }
354
355 @Override
356 public void finish() {
357 sharedInteractions.addAndGet(count);
358 sharedGKEnergy.addAndGet(gkEnergy);
359 sharedPermanentGKEnergy.addAndGet(gkPermanentEnergy);
360 sharedPolarizationGKEnergy.addAndGet(gkPolarizationEnergy);
361 }
362
363 @Override
364 public void run(int lb, int ub) {
365
366 double[] x = sXYZ[0][0];
367 double[] y = sXYZ[0][1];
368 double[] z = sXYZ[0][2];
369
370 int nSymm = crystal.spaceGroup.symOps.size();
371 List<SymOp> symOps = crystal.spaceGroup.symOps;
372
373 for (iSymm = 0; iSymm < nSymm; iSymm++) {
374 SymOp symOp = symOps.get(iSymm);
375 crystal.getTransformationOperator(symOp, transOp);
376 for (int i = lb; i <= ub; i++) {
377 if (!use[i]) {
378 continue;
379 }
380
381
382 dedxi = 0.0;
383 dedyi = 0.0;
384 dedzi = 0.0;
385 dborni = 0.0;
386 trqxi = 0.0;
387 trqyi = 0.0;
388 trqzi = 0.0;
389
390 xi = x[i];
391 yi = y[i];
392 zi = z[i];
393 final double[] multipolei = globalMultipole[0][i];
394 ci = multipolei[t000];
395 uxi = multipolei[t100];
396 uyi = multipolei[t010];
397 uzi = multipolei[t001];
398 qxxi = multipolei[t200] * oneThird;
399 qyyi = multipolei[t020] * oneThird;
400 qzzi = multipolei[t002] * oneThird;
401 qxyi = multipolei[t110] * oneThird;
402 qxzi = multipolei[t101] * oneThird;
403 qyzi = multipolei[t011] * oneThird;
404 dxi = inducedDipole[0][i][0];
405 dyi = inducedDipole[0][i][1];
406 dzi = inducedDipole[0][i][2];
407 pxi = inducedDipoleCR[0][i][0];
408 pyi = inducedDipoleCR[0][i][1];
409 pzi = inducedDipoleCR[0][i][2];
410 sxi = dxi + pxi;
411 syi = dyi + pyi;
412 szi = dzi + pzi;
413 rbi = born[i];
414 int[] list = neighborLists[iSymm][i];
415 for (int k : list) {
416 if (!use[k]) {
417 continue;
418 }
419 interaction(i, k);
420 }
421 if (iSymm == 0) {
422
423 interaction(i, i);
424
425
426
427
428
429
430 switch (nonPolarModel) {
431 case BORN_SOLV:
432 case BORN_CAV_DISP:
433 double r = baseRadius[i] + dOffset + probe;
434 double ratio = (baseRadius[i] + dOffset) / born[i];
435 ratio *= ratio;
436 ratio *= (ratio * ratio);
437 double saTerm = surfaceTension * r * r * ratio / 6.0;
438 gkEnergy += saTerm;
439 sharedBornGrad.sub(threadID, i, 6.0 * saTerm / born[i]);
440 break;
441 default:
442 break;
443 }
444 }
445 if (gradient) {
446 grad.add(threadID, i, dedxi, dedyi, dedzi);
447 torque.add(threadID, i, trqxi, trqyi, trqzi);
448 sharedBornGrad.add(threadID, i, dborni);
449 }
450 }
451 }
452 }
453
454 @Override
455 public void start() {
456 gkEnergy = 0.0;
457 gkPermanentEnergy = 0.0;
458 gkPolarizationEnergy = 0.0;
459 count = 0;
460 threadID = getThreadIndex();
461 }
462
463 private void interaction(int i, int k) {
464 dx_local[0] = sXYZ[iSymm][0][k] - xi;
465 dx_local[1] = sXYZ[iSymm][1][k] - yi;
466 dx_local[2] = sXYZ[iSymm][2][k] - zi;
467 double r2 = crystal.image(dx_local);
468 if (r2 > cut2) {
469 return;
470 }
471 xr = dx_local[0];
472 yr = dx_local[1];
473 zr = dx_local[2];
474 xr2 = xr * xr;
475 yr2 = yr * yr;
476 zr2 = zr * zr;
477 rbk = born[k];
478
479 final double[] multipolek = globalMultipole[iSymm][k];
480 ck = multipolek[t000];
481 uxk = multipolek[t100];
482 uyk = multipolek[t010];
483 uzk = multipolek[t001];
484 qxxk = multipolek[t200] * oneThird;
485 qyyk = multipolek[t020] * oneThird;
486 qzzk = multipolek[t002] * oneThird;
487 qxyk = multipolek[t110] * oneThird;
488 qxzk = multipolek[t101] * oneThird;
489 qyzk = multipolek[t011] * oneThird;
490 dxk = inducedDipole[iSymm][k][0];
491 dyk = inducedDipole[iSymm][k][1];
492 dzk = inducedDipole[iSymm][k][2];
493 pxk = inducedDipoleCR[iSymm][k][0];
494 pyk = inducedDipoleCR[iSymm][k][1];
495 pzk = inducedDipoleCR[iSymm][k][2];
496 sxk = dxk + pxk;
497 syk = dyk + pyk;
498 szk = dzk + pzk;
499 final double rb2 = rbi * rbk;
500 final double expterm = exp(-r2 / (gkc * rb2));
501 final double expc = expterm / gkc;
502 final double expc1 = 1.0 - expc;
503 final double expcr = r2 * expterm / (gkc * gkc * rb2 * rb2);
504 final double dexpc = -2.0 / (gkc * rb2);
505 double expcdexpc = -expc * dexpc;
506 final double dexpcr = 2.0 / (gkc * rb2 * rb2);
507 final double dgfdr = 0.5 * expterm * (1.0 + r2 / (rb2 * gkc));
508 final double gf2 = 1.0 / (r2 + rb2 * expterm);
509 final double gf = sqrt(gf2);
510 final double gf3 = gf2 * gf;
511 final double gf5 = gf3 * gf2;
512 final double gf7 = gf5 * gf2;
513 final double gf9 = gf7 * gf2;
514 final double gf11 = gf9 * gf2;
515
516
517 a[0][0] = gf;
518 a[1][0] = -gf3;
519 a[2][0] = 3.0 * gf5;
520 a[3][0] = -15.0 * gf7;
521 a[4][0] = 105.0 * gf9;
522 a[5][0] = -945.0 * gf11;
523
524
525 a[0][1] = expc1 * a[1][0];
526 a[1][1] = expc1 * a[2][0];
527 a[2][1] = expc1 * a[3][0];
528 a[3][1] = expc1 * a[4][0];
529 a[4][1] = expc1 * a[5][0];
530
531
532 a[0][2] = expc1 * a[1][1] + expcdexpc * a[1][0];
533 a[1][2] = expc1 * a[2][1] + expcdexpc * a[2][0];
534 a[2][2] = expc1 * a[3][1] + expcdexpc * a[3][0];
535 a[3][2] = expc1 * a[4][1] + expcdexpc * a[4][0];
536
537 if (gradient) {
538
539
540 expcdexpc = 2.0 * expcdexpc;
541 a[0][3] = expc1 * a[1][2] + expcdexpc * a[1][1];
542 a[1][3] = expc1 * a[2][2] + expcdexpc * a[2][1];
543 a[2][3] = expc1 * a[3][2] + expcdexpc * a[3][1];
544 expcdexpc = -expc * dexpc * dexpc;
545 a[0][3] = a[0][3] + expcdexpc * a[1][0];
546 a[1][3] = a[1][3] + expcdexpc * a[2][0];
547 a[2][3] = a[2][3] + expcdexpc * a[3][0];
548
549
550 b[0][0] = dgfdr * a[1][0];
551 b[1][0] = dgfdr * a[2][0];
552 b[2][0] = dgfdr * a[3][0];
553 b[3][0] = dgfdr * a[4][0];
554 b[4][0] = dgfdr * a[5][0];
555
556
557 b[0][1] = b[1][0] - expcr * a[1][0] - expc * b[1][0];
558 b[1][1] = b[2][0] - expcr * a[2][0] - expc * b[2][0];
559 b[2][1] = b[3][0] - expcr * a[3][0] - expc * b[3][0];
560 b[3][1] = b[4][0] - expcr * a[4][0] - expc * b[4][0];
561
562
563 b[0][2] = b[1][1] - (expcr * (a[1][1] + dexpc * a[1][0])
564 + expc * (b[1][1] + dexpcr * a[1][0] + dexpc * b[1][0]));
565 b[1][2] = b[2][1] - (expcr * (a[2][1] + dexpc * a[2][0])
566 + expc * (b[2][1] + dexpcr * a[2][0] + dexpc * b[2][0]));
567 b[2][2] = b[3][1] - (expcr * (a[3][1] + dexpc * a[3][0])
568 + expc * (b[3][1] + dexpcr * a[3][0] + dexpc * b[3][0]));
569
570
571 b[0][0] = electric * fc * b[0][0];
572 b[0][1] = electric * fc * b[0][1];
573 b[0][2] = electric * fc * b[0][2];
574 b[1][0] = electric * fd * b[1][0];
575 b[1][1] = electric * fd * b[1][1];
576 b[1][2] = electric * fd * b[1][2];
577 b[2][0] = electric * fq * b[2][0];
578 b[2][1] = electric * fq * b[2][1];
579 b[2][2] = electric * fq * b[2][2];
580 }
581
582
583 a[0][0] = electric * fc * a[0][0];
584 a[0][1] = electric * fc * a[0][1];
585 a[0][2] = electric * fc * a[0][2];
586 a[0][3] = electric * fc * a[0][3];
587 a[1][0] = electric * fd * a[1][0];
588 a[1][1] = electric * fd * a[1][1];
589 a[1][2] = electric * fd * a[1][2];
590 a[1][3] = electric * fd * a[1][3];
591 a[2][0] = electric * fq * a[2][0];
592 a[2][1] = electric * fq * a[2][1];
593 a[2][2] = electric * fq * a[2][2];
594 a[2][3] = electric * fq * a[2][3];
595
596
597 energyTensors();
598
599
600 double eik = energy(i, k);
601 gkEnergy += eik;
602 count++;
603
604
605
606
607
608
609
610
611
612
613 if (gradient) {
614
615 gradientTensors();
616
617
618 permanentEnergyGradient(i, k);
619
620 if (polarization != Polarization.NONE) {
621
622 polarizationEnergyGradient(i, k);
623 }
624 }
625 }
626
627 private void energyTensors() {
628
629 gc[1] = a[0][0];
630 gux[1] = xr * a[1][0];
631 guy[1] = yr * a[1][0];
632 guz[1] = zr * a[1][0];
633 gqxx[1] = xr2 * a[2][0];
634 gqyy[1] = yr2 * a[2][0];
635 gqzz[1] = zr2 * a[2][0];
636 gqxy[1] = xr * yr * a[2][0];
637 gqxz[1] = xr * zr * a[2][0];
638 gqyz[1] = yr * zr * a[2][0];
639
640
641 gc[2] = xr * a[0][1];
642 gc[3] = yr * a[0][1];
643 gc[4] = zr * a[0][1];
644 gux[2] = a[1][0] + xr2 * a[1][1];
645 gux[3] = xr * yr * a[1][1];
646 gux[4] = xr * zr * a[1][1];
647 guy[2] = gux[3];
648 guy[3] = a[1][0] + yr2 * a[1][1];
649 guy[4] = yr * zr * a[1][1];
650 guz[2] = gux[4];
651 guz[3] = guy[4];
652 guz[4] = a[1][0] + zr2 * a[1][1];
653 gqxx[2] = xr * (2.0 * a[2][0] + xr2 * a[2][1]);
654 gqxx[3] = yr * xr2 * a[2][1];
655 gqxx[4] = zr * xr2 * a[2][1];
656 gqyy[2] = xr * yr2 * a[2][1];
657 gqyy[3] = yr * (2.0 * a[2][0] + yr2 * a[2][1]);
658 gqyy[4] = zr * yr2 * a[2][1];
659 gqzz[2] = xr * zr2 * a[2][1];
660 gqzz[3] = yr * zr2 * a[2][1];
661 gqzz[4] = zr * (2.0 * a[2][0] + zr2 * a[2][1]);
662 gqxy[2] = yr * (a[2][0] + xr2 * a[2][1]);
663 gqxy[3] = xr * (a[2][0] + yr2 * a[2][1]);
664 gqxy[4] = zr * xr * yr * a[2][1];
665 gqxz[2] = zr * (a[2][0] + xr2 * a[2][1]);
666 gqxz[3] = gqxy[4];
667 gqxz[4] = xr * (a[2][0] + zr2 * a[2][1]);
668 gqyz[2] = gqxy[4];
669 gqyz[3] = zr * (a[2][0] + yr2 * a[2][1]);
670 gqyz[4] = yr * (a[2][0] + zr2 * a[2][1]);
671
672
673 gc[5] = a[0][1] + xr2 * a[0][2];
674 gc[6] = xr * yr * a[0][2];
675 gc[7] = xr * zr * a[0][2];
676 gc[8] = a[0][1] + yr2 * a[0][2];
677 gc[9] = yr * zr * a[0][2];
678 gc[10] = a[0][1] + zr2 * a[0][2];
679 gux[5] = xr * (3.0 * a[1][1] + xr2 * a[1][2]);
680 gux[6] = yr * (a[1][1] + xr2 * a[1][2]);
681 gux[7] = zr * (a[1][1] + xr2 * a[1][2]);
682 gux[8] = xr * (a[1][1] + yr2 * a[1][2]);
683 gux[9] = zr * xr * yr * a[1][2];
684 gux[10] = xr * (a[1][1] + zr2 * a[1][2]);
685 guy[5] = yr * (a[1][1] + xr2 * a[1][2]);
686 guy[6] = xr * (a[1][1] + yr2 * a[1][2]);
687 guy[7] = gux[9];
688 guy[8] = yr * (3.0 * a[1][1] + yr2 * a[1][2]);
689 guy[9] = zr * (a[1][1] + yr2 * a[1][2]);
690 guy[10] = yr * (a[1][1] + zr2 * a[1][2]);
691 guz[5] = zr * (a[1][1] + xr2 * a[1][2]);
692 guz[6] = gux[9];
693 guz[7] = xr * (a[1][1] + zr2 * a[1][2]);
694 guz[8] = zr * (a[1][1] + yr2 * a[1][2]);
695 guz[9] = yr * (a[1][1] + zr2 * a[1][2]);
696 guz[10] = zr * (3.0 * a[1][1] + zr2 * a[1][2]);
697 gqxx[5] = 2.0 * a[2][0] + xr2 * (5.0 * a[2][1] + xr2 * a[2][2]);
698 gqxx[6] = yr * xr * (2.0 * a[2][1] + xr2 * a[2][2]);
699 gqxx[7] = zr * xr * (2.0 * a[2][1] + xr2 * a[2][2]);
700 gqxx[8] = xr2 * (a[2][1] + yr2 * a[2][2]);
701 gqxx[9] = zr * yr * xr2 * a[2][2];
702 gqxx[10] = xr2 * (a[2][1] + zr2 * a[2][2]);
703 gqyy[5] = yr2 * (a[2][1] + xr2 * a[2][2]);
704 gqyy[6] = xr * yr * (2.0 * a[2][1] + yr2 * a[2][2]);
705 gqyy[7] = xr * zr * yr2 * a[2][2];
706 gqyy[8] = 2.0 * a[2][0] + yr2 * (5.0 * a[2][1] + yr2 * a[2][2]);
707 gqyy[9] = yr * zr * (2.0 * a[2][1] + yr2 * a[2][2]);
708 gqyy[10] = yr2 * (a[2][1] + zr2 * a[2][2]);
709 gqzz[5] = zr2 * (a[2][1] + xr2 * a[2][2]);
710 gqzz[6] = xr * yr * zr2 * a[2][2];
711 gqzz[7] = xr * zr * (2.0 * a[2][1] + zr2 * a[2][2]);
712 gqzz[8] = zr2 * (a[2][1] + yr2 * a[2][2]);
713 gqzz[9] = yr * zr * (2.0 * a[2][1] + zr2 * a[2][2]);
714 gqzz[10] = 2.0 * a[2][0] + zr2 * (5.0 * a[2][1] + zr2 * a[2][2]);
715 gqxy[5] = xr * yr * (3.0 * a[2][1] + xr2 * a[2][2]);
716 gqxy[6] = a[2][0] + (xr2 + yr2) * a[2][1] + xr2 * yr2 * a[2][2];
717 gqxy[7] = zr * yr * (a[2][1] + xr2 * a[2][2]);
718 gqxy[8] = xr * yr * (3.0 * a[2][1] + yr2 * a[2][2]);
719 gqxy[9] = zr * xr * (a[2][1] + yr2 * a[2][2]);
720 gqxy[10] = xr * yr * (a[2][1] + zr2 * a[2][2]);
721 gqxz[5] = xr * zr * (3.0 * a[2][1] + xr2 * a[2][2]);
722 gqxz[6] = yr * zr * (a[2][1] + xr2 * a[2][2]);
723 gqxz[7] = a[2][0] + (xr2 + zr2) * a[2][1] + xr2 * zr2 * a[2][2];
724 gqxz[8] = xr * zr * (a[2][1] + yr2 * a[2][2]);
725 gqxz[9] = xr * yr * (a[2][1] + zr2 * a[2][2]);
726 gqxz[10] = xr * zr * (3.0 * a[2][1] + zr2 * a[2][2]);
727 gqyz[5] = zr * yr * (a[2][1] + xr2 * a[2][2]);
728 gqyz[6] = xr * zr * (a[2][1] + yr2 * a[2][2]);
729 gqyz[7] = xr * yr * (a[2][1] + zr2 * a[2][2]);
730 gqyz[8] = yr * zr * (3.0 * a[2][1] + yr2 * a[2][2]);
731 gqyz[9] = a[2][0] + (yr2 + zr2) * a[2][1] + yr2 * zr2 * a[2][2];
732 gqyz[10] = yr * zr * (3.0 * a[2][1] + zr2 * a[2][2]);
733 }
734
735 private double energy(int i, int k) {
736
737
738 double esym =
739 ci * ck * gc[1]
740 - (uxi * (uxk * gux[2] + uyk * guy[2] + uzk * guz[2])
741 + uyi * (uxk * gux[3] + uyk * guy[3] + uzk * guz[3])
742 + uzi * (uxk * gux[4] + uyk * guy[4] + uzk * guz[4]));
743 double ewi =
744 ci * (uxk * gc[2] + uyk * gc[3] + uzk * gc[4])
745 - ck * (uxi * gux[1] + uyi * guy[1] + uzi * guz[1])
746 + ci
747 * (qxxk * gc[5]
748 + qyyk * gc[8]
749 + qzzk * gc[10]
750 + 2.0 * (qxyk * gc[6] + qxzk * gc[7] + qyzk * gc[9]))
751 + ck
752 * (qxxi * gqxx[1]
753 + qyyi * gqyy[1]
754 + qzzi * gqzz[1]
755 + 2.0 * (qxyi * gqxy[1] + qxzi * gqxz[1] + qyzi * gqyz[1]))
756 - uxi
757 * (qxxk * gux[5]
758 + qyyk * gux[8]
759 + qzzk * gux[10]
760 + 2.0 * (qxyk * gux[6] + qxzk * gux[7] + qyzk * gux[9]))
761 - uyi
762 * (qxxk * guy[5]
763 + qyyk * guy[8]
764 + qzzk * guy[10]
765 + 2.0 * (qxyk * guy[6] + qxzk * guy[7] + qyzk * guy[9]))
766 - uzi
767 * (qxxk * guz[5]
768 + qyyk * guz[8]
769 + qzzk * guz[10]
770 + 2.0 * (qxyk * guz[6] + qxzk * guz[7] + qyzk * guz[9]))
771 + uxk
772 * (qxxi * gqxx[2]
773 + qyyi * gqyy[2]
774 + qzzi * gqzz[2]
775 + 2.0 * (qxyi * gqxy[2] + qxzi * gqxz[2] + qyzi * gqyz[2]))
776 + uyk
777 * (qxxi * gqxx[3]
778 + qyyi * gqyy[3]
779 + qzzi * gqzz[3]
780 + 2.0 * (qxyi * gqxy[3] + qxzi * gqxz[3] + qyzi * gqyz[3]))
781 + uzk
782 * (qxxi * gqxx[4]
783 + qyyi * gqyy[4]
784 + qzzi * gqzz[4]
785 + 2.0 * (qxyi * gqxy[4] + qxzi * gqxz[4] + qyzi * gqyz[4]))
786 + qxxi
787 * (qxxk * gqxx[5]
788 + qyyk * gqxx[8]
789 + qzzk * gqxx[10]
790 + 2.0 * (qxyk * gqxx[6] + qxzk * gqxx[7] + qyzk * gqxx[9]))
791 + qyyi
792 * (qxxk * gqyy[5]
793 + qyyk * gqyy[8]
794 + qzzk * gqyy[10]
795 + 2.0 * (qxyk * gqyy[6] + qxzk * gqyy[7] + qyzk * gqyy[9]))
796 + qzzi
797 * (qxxk * gqzz[5]
798 + qyyk * gqzz[8]
799 + qzzk * gqzz[10]
800 + 2.0 * (qxyk * gqzz[6] + qxzk * gqzz[7] + qyzk * gqzz[9]))
801 + 2.0
802 * (qxyi
803 * (qxxk * gqxy[5]
804 + qyyk * gqxy[8]
805 + qzzk * gqxy[10]
806 + 2.0 * (qxyk * gqxy[6] + qxzk * gqxy[7] + qyzk * gqxy[9]))
807 + qxzi
808 * (qxxk * gqxz[5]
809 + qyyk * gqxz[8]
810 + qzzk * gqxz[10]
811 + 2.0 * (qxyk * gqxz[6] + qxzk * gqxz[7] + qyzk * gqxz[9]))
812 + qyzi
813 * (qxxk * gqyz[5]
814 + qyyk * gqyz[8]
815 + qzzk * gqyz[10]
816 + 2.0 * (qxyk * gqyz[6] + qxzk * gqyz[7] + qyzk * gqyz[9])));
817 double ewk =
818 ci * (uxk * gux[1] + uyk * guy[1] + uzk * guz[1])
819 - ck * (uxi * gc[2] + uyi * gc[3] + uzi * gc[4])
820 + ci
821 * (qxxk * gqxx[1]
822 + qyyk * gqyy[1]
823 + qzzk * gqzz[1]
824 + 2.0 * (qxyk * gqxy[1] + qxzk * gqxz[1] + qyzk * gqyz[1]))
825 + ck
826 * (qxxi * gc[5]
827 + qyyi * gc[8]
828 + qzzi * gc[10]
829 + 2.0 * (qxyi * gc[6] + qxzi * gc[7] + qyzi * gc[9]))
830 - uxi
831 * (qxxk * gqxx[2]
832 + qyyk * gqyy[2]
833 + qzzk * gqzz[2]
834 + 2.0 * (qxyk * gqxy[2] + qxzk * gqxz[2] + qyzk * gqyz[2]))
835 - uyi
836 * (qxxk * gqxx[3]
837 + qyyk * gqyy[3]
838 + qzzk * gqzz[3]
839 + 2.0 * (qxyk * gqxy[3] + qxzk * gqxz[3] + qyzk * gqyz[3]))
840 - uzi
841 * (qxxk * gqxx[4]
842 + qyyk * gqyy[4]
843 + qzzk * gqzz[4]
844 + 2.0 * (qxyk * gqxy[4] + qxzk * gqxz[4] + qyzk * gqyz[4]))
845 + uxk
846 * (qxxi * gux[5]
847 + qyyi * gux[8]
848 + qzzi * gux[10]
849 + 2.0 * (qxyi * gux[6] + qxzi * gux[7] + qyzi * gux[9]))
850 + uyk
851 * (qxxi * guy[5]
852 + qyyi * guy[8]
853 + qzzi * guy[10]
854 + 2.0 * (qxyi * guy[6] + qxzi * guy[7] + qyzi * guy[9]))
855 + uzk
856 * (qxxi * guz[5]
857 + qyyi * guz[8]
858 + qzzi * guz[10]
859 + 2.0 * (qxyi * guz[6] + qxzi * guz[7] + qyzi * guz[9]))
860 + qxxi
861 * (qxxk * gqxx[5]
862 + qyyk * gqyy[5]
863 + qzzk * gqzz[5]
864 + 2.0 * (qxyk * gqxy[5] + qxzk * gqxz[5] + qyzk * gqyz[5]))
865 + qyyi
866 * (qxxk * gqxx[8]
867 + qyyk * gqyy[8]
868 + qzzk * gqzz[8]
869 + 2.0 * (qxyk * gqxy[8] + qxzk * gqxz[8] + qyzk * gqyz[8]))
870 + qzzi
871 * (qxxk * gqxx[10]
872 + qyyk * gqyy[10]
873 + qzzk * gqzz[10]
874 + 2.0 * (qxyk * gqxy[10] + qxzk * gqxz[10] + qyzk * gqyz[10]))
875 + 2.0
876 * (qxyi
877 * (qxxk * gqxx[6]
878 + qyyk * gqyy[6]
879 + qzzk * gqzz[6]
880 + 2.0 * (qxyk * gqxy[6] + qxzk * gqxz[6] + qyzk * gqyz[6]))
881 + qxzi
882 * (qxxk * gqxx[7]
883 + qyyk * gqyy[7]
884 + qzzk * gqzz[7]
885 + 2.0 * (qxyk * gqxy[7] + qxzk * gqxz[7] + qyzk * gqyz[7]))
886 + qyzi
887 * (qxxk * gqxx[9]
888 + qyyk * gqyy[9]
889 + qzzk * gqzz[9]
890 + 2.0 * (qxyk * gqxy[9] + qxzk * gqxz[9] + qyzk * gqyz[9])));
891 double e = esym + 0.5 * (ewi + ewk);
892 double ei = 0.0;
893
894
895
896 if (polarization != Polarization.NONE) {
897 double esymi =
898 -uxi * (dxk * gux[2] + dyk * guy[2] + dzk * guz[2])
899 - uyi * (dxk * gux[3] + dyk * guy[3] + dzk * guz[3])
900 - uzi * (dxk * gux[4] + dyk * guy[4] + dzk * guz[4])
901 - uxk * (dxi * gux[2] + dyi * guy[2] + dzi * guz[2])
902 - uyk * (dxi * gux[3] + dyi * guy[3] + dzi * guz[3])
903 - uzk * (dxi * gux[4] + dyi * guy[4] + dzi * guz[4]);
904 double ewii =
905 ci * (dxk * gc[2] + dyk * gc[3] + dzk * gc[4])
906 - ck * (dxi * gux[1] + dyi * guy[1] + dzi * guz[1])
907 - dxi
908 * (qxxk * gux[5]
909 + qyyk * gux[8]
910 + qzzk * gux[10]
911 + 2.0 * (qxyk * gux[6] + qxzk * gux[7] + qyzk * gux[9]))
912 - dyi
913 * (qxxk * guy[5]
914 + qyyk * guy[8]
915 + qzzk * guy[10]
916 + 2.0 * (qxyk * guy[6] + qxzk * guy[7] + qyzk * guy[9]))
917 - dzi
918 * (qxxk * guz[5]
919 + qyyk * guz[8]
920 + qzzk * guz[10]
921 + 2.0 * (qxyk * guz[6] + qxzk * guz[7] + qyzk * guz[9]))
922 + dxk
923 * (qxxi * gqxx[2]
924 + qyyi * gqyy[2]
925 + qzzi * gqzz[2]
926 + 2.0 * (qxyi * gqxy[2] + qxzi * gqxz[2] + qyzi * gqyz[2]))
927 + dyk
928 * (qxxi * gqxx[3]
929 + qyyi * gqyy[3]
930 + qzzi * gqzz[3]
931 + 2.0 * (qxyi * gqxy[3] + qxzi * gqxz[3] + qyzi * gqyz[3]))
932 + dzk
933 * (qxxi * gqxx[4]
934 + qyyi * gqyy[4]
935 + qzzi * gqzz[4]
936 + 2.0 * (qxyi * gqxy[4] + qxzi * gqxz[4] + qyzi * gqyz[4]));
937 double ewki =
938 ci * (dxk * gux[1] + dyk * guy[1] + dzk * guz[1])
939 - ck * (dxi * gc[2] + dyi * gc[3] + dzi * gc[4])
940 - dxi
941 * (qxxk * gqxx[2]
942 + qyyk * gqyy[2]
943 + qzzk * gqzz[2]
944 + 2.0 * (qxyk * gqxy[2] + qxzk * gqxz[2] + qyzk * gqyz[2]))
945 - dyi
946 * (qxxk * gqxx[3]
947 + qyyk * gqyy[3]
948 + qzzk * gqzz[3]
949 + 2.0 * (qxyk * gqxy[3] + qxzk * gqxz[3] + qyzk * gqyz[3]))
950 - dzi
951 * (qxxk * gqxx[4]
952 + qyyk * gqyy[4]
953 + qzzk * gqzz[4]
954 + 2.0 * (qxyk * gqxy[4] + qxzk * gqxz[4] + qyzk * gqyz[4]))
955 + dxk
956 * (qxxi * gux[5]
957 + qyyi * gux[8]
958 + qzzi * gux[10]
959 + 2.0 * (qxyi * gux[6] + qxzi * gux[7] + qyzi * gux[9]))
960 + dyk
961 * (qxxi * guy[5]
962 + qyyi * guy[8]
963 + qzzi * guy[10]
964 + 2.0 * (qxyi * guy[6] + qxzi * guy[7] + qyzi * guy[9]))
965 + dzk
966 * (qxxi * guz[5]
967 + qyyi * guz[8]
968 + qzzi * guz[10]
969 + 2.0 * (qxyi * guz[6] + qxzi * guz[7] + qyzi * guz[9]));
970 ei = 0.5 * (esymi + 0.5 * (ewii + ewki));
971 }
972
973 if (i == k) {
974 e *= 0.5;
975 ei *= 0.5;
976 selfEnergy.add(threadID, i, e + ei);
977 } else {
978 double half = 0.5 * (e + ei);
979 crossEnergy.add(threadID, i, half);
980 crossEnergy.add(threadID, k, half);
981 }
982
983 gkPermanentEnergy += e;
984 gkPolarizationEnergy += ei;
985
986 return e + ei;
987 }
988
989 private void gradientTensors() {
990
991 gc[21] = b[0][0];
992 gux[21] = xr * b[1][0];
993 guy[21] = yr * b[1][0];
994 guz[21] = zr * b[1][0];
995 gqxx[21] = xr2 * b[2][0];
996 gqyy[21] = yr2 * b[2][0];
997 gqzz[21] = zr2 * b[2][0];
998 gqxy[21] = xr * yr * b[2][0];
999 gqxz[21] = xr * zr * b[2][0];
1000 gqyz[21] = yr * zr * b[2][0];
1001
1002
1003 gc[22] = xr * b[0][1];
1004 gc[23] = yr * b[0][1];
1005 gc[24] = zr * b[0][1];
1006 gux[22] = b[1][0] + xr2 * b[1][1];
1007 gux[23] = xr * yr * b[1][1];
1008 gux[24] = xr * zr * b[1][1];
1009 guy[22] = gux[23];
1010 guy[23] = b[1][0] + yr2 * b[1][1];
1011 guy[24] = yr * zr * b[1][1];
1012 guz[22] = gux[24];
1013 guz[23] = guy[24];
1014 guz[24] = b[1][0] + zr2 * b[1][1];
1015 gqxx[22] = xr * (2.0 * b[2][0] + xr2 * b[2][1]);
1016 gqxx[23] = yr * xr2 * b[2][1];
1017 gqxx[24] = zr * xr2 * b[2][1];
1018 gqyy[22] = xr * yr2 * b[2][1];
1019 gqyy[23] = yr * (2.0 * b[2][0] + yr2 * b[2][1]);
1020 gqyy[24] = zr * yr2 * b[2][1];
1021 gqzz[22] = xr * zr2 * b[2][1];
1022 gqzz[23] = yr * zr2 * b[2][1];
1023 gqzz[24] = zr * (2.0 * b[2][0] + zr2 * b[2][1]);
1024 gqxy[22] = yr * (b[2][0] + xr2 * b[2][1]);
1025 gqxy[23] = xr * (b[2][0] + yr2 * b[2][1]);
1026 gqxy[24] = zr * xr * yr * b[2][1];
1027 gqxz[22] = zr * (b[2][0] + xr2 * b[2][1]);
1028 gqxz[23] = gqxy[24];
1029 gqxz[24] = xr * (b[2][0] + zr2 * b[2][1]);
1030 gqyz[22] = gqxy[24];
1031 gqyz[23] = zr * (b[2][0] + yr2 * b[2][1]);
1032 gqyz[24] = yr * (b[2][0] + zr2 * b[2][1]);
1033
1034
1035 gc[25] = b[0][1] + xr2 * b[0][2];
1036 gc[26] = xr * yr * b[0][2];
1037 gc[27] = xr * zr * b[0][2];
1038 gc[28] = b[0][1] + yr2 * b[0][2];
1039 gc[29] = yr * zr * b[0][2];
1040 gc[30] = b[0][1] + zr2 * b[0][2];
1041 gux[25] = xr * (3.0 * b[1][1] + xr2 * b[1][2]);
1042 gux[26] = yr * (b[1][1] + xr2 * b[1][2]);
1043 gux[27] = zr * (b[1][1] + xr2 * b[1][2]);
1044 gux[28] = xr * (b[1][1] + yr2 * b[1][2]);
1045 gux[29] = zr * xr * yr * b[1][2];
1046 gux[30] = xr * (b[1][1] + zr2 * b[1][2]);
1047 guy[25] = yr * (b[1][1] + xr2 * b[1][2]);
1048 guy[26] = xr * (b[1][1] + yr2 * b[1][2]);
1049 guy[27] = gux[29];
1050 guy[28] = yr * (3.0 * b[1][1] + yr2 * b[1][2]);
1051 guy[29] = zr * (b[1][1] + yr2 * b[1][2]);
1052 guy[30] = yr * (b[1][1] + zr2 * b[1][2]);
1053 guz[25] = zr * (b[1][1] + xr2 * b[1][2]);
1054 guz[26] = gux[29];
1055 guz[27] = xr * (b[1][1] + zr2 * b[1][2]);
1056 guz[28] = zr * (b[1][1] + yr2 * b[1][2]);
1057 guz[29] = yr * (b[1][1] + zr2 * b[1][2]);
1058 guz[30] = zr * (3.0 * b[1][1] + zr2 * b[1][2]);
1059 gqxx[25] = 2.0 * b[2][0] + xr2 * (5.0 * b[2][1] + xr2 * b[2][2]);
1060 gqxx[26] = yr * xr * (2.0 * b[2][1] + xr2 * b[2][2]);
1061 gqxx[27] = zr * xr * (2.0 * b[2][1] + xr2 * b[2][2]);
1062 gqxx[28] = xr2 * (b[2][1] + yr2 * b[2][2]);
1063 gqxx[29] = zr * yr * xr2 * b[2][2];
1064 gqxx[30] = xr2 * (b[2][1] + zr2 * b[2][2]);
1065 gqyy[25] = yr2 * (b[2][1] + xr2 * b[2][2]);
1066 gqyy[26] = xr * yr * (2.0 * b[2][1] + yr2 * b[2][2]);
1067 gqyy[27] = xr * zr * yr2 * b[2][2];
1068 gqyy[28] = 2.0 * b[2][0] + yr2 * (5.0 * b[2][1] + yr2 * b[2][2]);
1069 gqyy[29] = yr * zr * (2.0 * b[2][1] + yr2 * b[2][2]);
1070 gqyy[30] = yr2 * (b[2][1] + zr2 * b[2][2]);
1071 gqzz[25] = zr2 * (b[2][1] + xr2 * b[2][2]);
1072 gqzz[26] = xr * yr * zr2 * b[2][2];
1073 gqzz[27] = xr * zr * (2.0 * b[2][1] + zr2 * b[2][2]);
1074 gqzz[28] = zr2 * (b[2][1] + yr2 * b[2][2]);
1075 gqzz[29] = yr * zr * (2.0 * b[2][1] + zr2 * b[2][2]);
1076 gqzz[30] = 2.0 * b[2][0] + zr2 * (5.0 * b[2][1] + zr2 * b[2][2]);
1077 gqxy[25] = xr * yr * (3.0 * b[2][1] + xr2 * b[2][2]);
1078 gqxy[26] = b[2][0] + (xr2 + yr2) * b[2][1] + xr2 * yr2 * b[2][2];
1079 gqxy[27] = zr * yr * (b[2][1] + xr2 * b[2][2]);
1080 gqxy[28] = xr * yr * (3.0 * b[2][1] + yr2 * b[2][2]);
1081 gqxy[29] = zr * xr * (b[2][1] + yr2 * b[2][2]);
1082 gqxy[30] = xr * yr * (b[2][1] + zr2 * b[2][2]);
1083 gqxz[25] = xr * zr * (3.0 * b[2][1] + xr2 * b[2][2]);
1084 gqxz[26] = yr * zr * (b[2][1] + xr2 * b[2][2]);
1085 gqxz[27] = b[2][0] + (xr2 + zr2) * b[2][1] + xr2 * zr2 * b[2][2];
1086 gqxz[28] = xr * zr * (b[2][1] + yr2 * b[2][2]);
1087 gqxz[29] = xr * yr * (b[2][1] + zr2 * b[2][2]);
1088 gqxz[30] = xr * zr * (3.0 * b[2][1] + zr2 * b[2][2]);
1089 gqyz[25] = zr * yr * (b[2][1] + xr2 * b[2][2]);
1090 gqyz[26] = xr * zr * (b[2][1] + yr2 * b[2][2]);
1091 gqyz[27] = xr * yr * (b[2][1] + zr2 * b[2][2]);
1092 gqyz[28] = yr * zr * (3.0 * b[2][1] + yr2 * b[2][2]);
1093 gqyz[29] = b[2][0] + (yr2 + zr2) * b[2][1] + yr2 * zr2 * b[2][2];
1094 gqyz[30] = yr * zr * (3.0 * b[2][1] + zr2 * b[2][2]);
1095
1096
1097 gc[11] = xr * (3.0 * a[0][2] + xr2 * a[0][3]);
1098 gc[12] = yr * (a[0][2] + xr2 * a[0][3]);
1099 gc[13] = zr * (a[0][2] + xr2 * a[0][3]);
1100 gc[14] = xr * (a[0][2] + yr2 * a[0][3]);
1101 gc[15] = xr * yr * zr * a[0][3];
1102 gc[16] = xr * (a[0][2] + zr2 * a[0][3]);
1103 gc[17] = yr * (3.0 * a[0][2] + yr2 * a[0][3]);
1104 gc[18] = zr * (a[0][2] + yr2 * a[0][3]);
1105 gc[19] = yr * (a[0][2] + zr2 * a[0][3]);
1106 gc[20] = zr * (3.0 * a[0][2] + zr2 * a[0][3]);
1107
1108 gux[11] = 3.0 * a[1][1] + xr2 * (6.0 * a[1][2] + xr2 * a[1][3]);
1109 gux[12] = xr * yr * (3.0 * a[1][2] + xr2 * a[1][3]);
1110 gux[13] = xr * zr * (3.0 * a[1][2] + xr2 * a[1][3]);
1111 gux[14] = a[1][1] + (xr2 + yr2) * a[1][2] + xr2 * yr2 * a[1][3];
1112 gux[15] = yr * zr * (a[1][2] + xr2 * a[1][3]);
1113 gux[16] = a[1][1] + (xr2 + zr2) * a[1][2] + xr2 * zr2 * a[1][3];
1114 gux[17] = xr * yr * (3.0 * a[1][2] + yr2 * a[1][3]);
1115 gux[18] = xr * zr * (a[1][2] + yr2 * a[1][3]);
1116 gux[19] = xr * yr * (a[1][2] + zr2 * a[1][3]);
1117 gux[20] = xr * zr * (3.0 * a[1][2] + zr2 * a[1][3]);
1118
1119 guy[11] = gux[12];
1120 guy[12] = gux[14];
1121 guy[13] = gux[15];
1122 guy[14] = gux[17];
1123 guy[15] = gux[18];
1124 guy[16] = gux[19];
1125 guy[17] = 3.0 * a[1][1] + yr2 * (6.0 * a[1][2] + yr2 * a[1][3]);
1126 guy[18] = yr * zr * (3.0 * a[1][2] + yr2 * a[1][3]);
1127 guy[19] = a[1][1] + (yr2 + zr2) * a[1][2] + yr2 * zr2 * a[1][3];
1128 guy[20] = yr * zr * (3.0 * a[1][2] + zr2 * a[1][3]);
1129
1130 guz[11] = gux[13];
1131 guz[12] = gux[15];
1132 guz[13] = gux[16];
1133 guz[14] = gux[18];
1134 guz[15] = gux[19];
1135 guz[16] = gux[20];
1136 guz[17] = guy[18];
1137 guz[18] = guy[19];
1138 guz[19] = guy[20];
1139 guz[20] = 3.0 * a[1][1] + zr2 * (6.0 * a[1][2] + zr2 * a[1][3]);
1140
1141 gqxx[11] = xr * (12.0 * a[2][1] + xr2 * (9.0 * a[2][2] + xr2 * a[2][3]));
1142 gqxx[12] = yr * (2.0 * a[2][1] + xr2 * (5.0 * a[2][2] + xr2 * a[2][3]));
1143 gqxx[13] = zr * (2.0 * a[2][1] + xr2 * (5.0 * a[2][2] + xr2 * a[2][3]));
1144 gqxx[14] = xr * (2.0 * a[2][1] + yr2 * 2.0 * a[2][2] + xr2 * (a[2][2] + yr2 * a[2][3]));
1145 gqxx[15] = xr * yr * zr * (2.0 * a[2][2] + xr2 * a[2][3]);
1146 gqxx[16] = xr * (2.0 * a[2][1] + zr2 * 2.0 * a[2][2] + xr2 * (a[2][2] + zr2 * a[2][3]));
1147 gqxx[17] = yr * xr2 * (3.0 * a[2][2] + yr2 * a[2][3]);
1148 gqxx[18] = zr * xr2 * (a[2][2] + yr2 * a[2][3]);
1149 gqxx[19] = yr * xr2 * (a[2][2] + zr2 * a[2][3]);
1150 gqxx[20] = zr * xr2 * (3.0 * a[2][2] + zr2 * a[2][3]);
1151
1152 gqxy[11] = yr * (3.0 * a[2][1] + xr2 * (6.0 * a[2][2] + xr2 * a[2][3]));
1153 gqxy[12] = xr * (3.0 * (a[2][1] + yr2 * a[2][2]) + xr2 * (a[2][2] + yr2 * a[2][3]));
1154 gqxy[13] = xr * yr * zr * (3.0 * a[2][2] + xr2 * a[2][3]);
1155 gqxy[14] = yr * (3.0 * (a[2][1] + xr2 * a[2][2]) + yr2 * (a[2][2] + xr2 * a[2][3]));
1156 gqxy[15] = zr * (a[2][1] + (yr2 + xr2) * a[2][2] + yr2 * xr2 * a[2][3]);
1157 gqxy[16] = yr * (a[2][1] + (xr2 + zr2) * a[2][2] + xr2 * zr2 * a[2][3]);
1158 gqxy[17] = xr * (3.0 * (a[2][1] + yr2 * a[2][2]) + yr2 * (3.0 * a[2][2] + yr2 * a[2][3]));
1159 gqxy[18] = xr * yr * zr * (3.0 * a[2][2] + yr2 * a[2][3]);
1160 gqxy[19] = xr * (a[2][1] + (yr2 + zr2) * a[2][2] + yr2 * zr2 * a[2][3]);
1161 gqxy[20] = xr * yr * zr * (3.0 * a[2][2] + zr2 * a[2][3]);
1162
1163 gqxz[11] = zr * (3.0 * a[2][1] + xr2 * (6.0 * a[2][2] + xr2 * a[2][3]));
1164 gqxz[12] = xr * yr * zr * (3.0 * a[2][2] + xr2 * a[2][3]);
1165 gqxz[13] = xr * (3.0 * (a[2][1] + zr2 * a[2][2]) + xr2 * (a[2][2] + zr2 * a[2][3]));
1166 gqxz[14] = zr * (a[2][1] + (xr2 + yr2) * a[2][2] + xr2 * yr2 * a[2][3]);
1167 gqxz[15] = yr * (a[2][1] + (xr2 + zr2) * a[2][2] + zr2 * xr2 * a[2][3]);
1168 gqxz[16] = zr * (3.0 * (a[2][1] + xr2 * a[2][2]) + zr2 * (a[2][2] + xr2 * a[2][3]));
1169 gqxz[17] = xr * yr * zr * (3.0 * a[2][2] + yr2 * a[2][3]);
1170 gqxz[18] = xr * (a[2][1] + (zr2 + yr2) * a[2][2] + zr2 * yr2 * a[2][3]);
1171 gqxz[19] = xr * yr * zr * (3.0 * a[2][2] + zr2 * a[2][3]);
1172 gqxz[20] = xr * (3.0 * a[2][1] + zr2 * (6.0 * a[2][2] + zr2 * a[2][3]));
1173
1174 gqyy[11] = xr * yr2 * (3.0 * a[2][2] + xr2 * a[2][3]);
1175 gqyy[12] = yr * (2.0 * a[2][1] + xr2 * 2.0 * a[2][2] + yr2 * (a[2][2] + xr2 * a[2][3]));
1176 gqyy[13] = zr * yr2 * (a[2][2] + xr2 * a[2][3]);
1177 gqyy[14] = xr * (2.0 * a[2][1] + yr2 * (5.0 * a[2][2] + yr2 * a[2][3]));
1178 gqyy[15] = xr * yr * zr * (2.0 * a[2][2] + yr2 * a[2][3]);
1179 gqyy[16] = xr * yr2 * (a[2][2] + zr2 * a[2][3]);
1180 gqyy[17] = yr * (12.0 * a[2][1] + yr2 * (9.0 * a[2][2] + yr2 * a[2][3]));
1181 gqyy[18] = zr * (2.0 * a[2][1] + yr2 * (5.0 * a[2][2] + yr2 * a[2][3]));
1182 gqyy[19] = yr * (2.0 * a[2][1] + zr2 * 2.0 * a[2][2] + yr2 * (a[2][2] + zr2 * a[2][3]));
1183 gqyy[20] = zr * yr2 * (3.0 * a[2][2] + zr2 * a[2][3]);
1184
1185 gqyz[11] = xr * yr * zr * (3.0 * a[2][2] + xr2 * a[2][3]);
1186 gqyz[12] = zr * (a[2][1] + (xr2 + yr2) * a[2][2] + xr2 * yr2 * a[2][3]);
1187 gqyz[13] = yr * (a[2][1] + (xr2 + zr2) * a[2][2] + xr2 * zr2 * a[2][3]);
1188 gqyz[14] = xr * yr * zr * (3.0 * a[2][2] + yr2 * a[2][3]);
1189 gqyz[15] = xr * (a[2][1] + (yr2 + zr2) * a[2][2] + yr2 * zr2 * a[2][3]);
1190 gqyz[16] = xr * yr * zr * (3.0 * a[2][2] + zr2 * a[2][3]);
1191 gqyz[17] = zr * (3.0 * a[2][1] + yr2 * (6.0 * a[2][2] + yr2 * a[2][3]));
1192 gqyz[18] = yr * (3.0 * (a[2][1] + zr2 * a[2][2]) + yr2 * (a[2][2] + zr2 * a[2][3]));
1193 gqyz[19] = zr * (3.0 * (a[2][1] + yr2 * a[2][2]) + zr2 * (a[2][2] + yr2 * a[2][3]));
1194 gqyz[20] = yr * (3.0 * a[2][1] + zr2 * (6.0 * a[2][2] + zr2 * a[2][3]));
1195
1196 gqzz[11] = xr * zr2 * (3.0 * a[2][2] + xr2 * a[2][3]);
1197 gqzz[12] = yr * (zr2 * a[2][2] + xr2 * (zr2 * a[2][3]));
1198 gqzz[13] = zr * (2.0 * a[2][1] + xr2 * 2.0 * a[2][2] + zr2 * (a[2][2] + xr2 * a[2][3]));
1199 gqzz[14] = xr * zr2 * (a[2][2] + yr2 * a[2][3]);
1200 gqzz[15] = xr * yr * zr * (2.0 * a[2][2] + zr2 * a[2][3]);
1201 gqzz[16] = xr * (2.0 * a[2][1] + zr2 * (5.0 * a[2][2] + zr2 * a[2][3]));
1202 gqzz[17] = yr * zr2 * (3.0 * a[2][2] + yr2 * a[2][3]);
1203 gqzz[18] = zr * (2.0 * a[2][1] + yr2 * 2.0 * a[2][2] + zr2 * (a[2][2] + yr2 * a[2][3]));
1204 gqzz[19] = yr * (2.0 * a[2][1] + zr2 * (5.0 * a[2][2] + zr2 * a[2][3]));
1205 gqzz[20] = zr * (12.0 * a[2][1] + zr2 * (9.0 * a[2][2] + zr2 * a[2][3]));
1206 }
1207
1208 private void permanentEnergyGradient(int i, int k) {
1209 final double desymdr =
1210 ci * ck * gc[21]
1211 - (uxi * (uxk * gux[22] + uyk * guy[22] + uzk * guz[22])
1212 + uyi * (uxk * gux[23] + uyk * guy[23] + uzk * guz[23])
1213 + uzi * (uxk * gux[24] + uyk * guy[24] + uzk * guz[24]));
1214 final double dewidr =
1215 ci * (uxk * gc[22] + uyk * gc[23] + uzk * gc[24])
1216 - ck * (uxi * gux[21] + uyi * guy[21] + uzi * guz[21])
1217 + ci
1218 * (qxxk * gc[25]
1219 + qyyk * gc[28]
1220 + qzzk * gc[30]
1221 + 2.0 * (qxyk * gc[26] + qxzk * gc[27] + qyzk * gc[29]))
1222 + ck
1223 * (qxxi * gqxx[21]
1224 + qyyi * gqyy[21]
1225 + qzzi * gqzz[21]
1226 + 2.0 * (qxyi * gqxy[21] + qxzi * gqxz[21] + qyzi * gqyz[21]))
1227 - uxi
1228 * (qxxk * gux[25]
1229 + qyyk * gux[28]
1230 + qzzk * gux[30]
1231 + 2.0 * (qxyk * gux[26] + qxzk * gux[27] + qyzk * gux[29]))
1232 - uyi
1233 * (qxxk * guy[25]
1234 + qyyk * guy[28]
1235 + qzzk * guy[30]
1236 + 2.0 * (qxyk * guy[26] + qxzk * guy[27] + qyzk * guy[29]))
1237 - uzi
1238 * (qxxk * guz[25]
1239 + qyyk * guz[28]
1240 + qzzk * guz[30]
1241 + 2.0 * (qxyk * guz[26] + qxzk * guz[27] + qyzk * guz[29]))
1242 + uxk
1243 * (qxxi * gqxx[22]
1244 + qyyi * gqyy[22]
1245 + qzzi * gqzz[22]
1246 + 2.0 * (qxyi * gqxy[22] + qxzi * gqxz[22] + qyzi * gqyz[22]))
1247 + uyk
1248 * (qxxi * gqxx[23]
1249 + qyyi * gqyy[23]
1250 + qzzi * gqzz[23]
1251 + 2.0 * (qxyi * gqxy[23] + qxzi * gqxz[23] + qyzi * gqyz[23]))
1252 + uzk
1253 * (qxxi * gqxx[24]
1254 + qyyi * gqyy[24]
1255 + qzzi * gqzz[24]
1256 + 2.0 * (qxyi * gqxy[24] + qxzi * gqxz[24] + qyzi * gqyz[24]))
1257 + qxxi
1258 * (qxxk * gqxx[25]
1259 + qyyk * gqxx[28]
1260 + qzzk * gqxx[30]
1261 + 2.0 * (qxyk * gqxx[26] + qxzk * gqxx[27] + qyzk * gqxx[29]))
1262 + qyyi
1263 * (qxxk * gqyy[25]
1264 + qyyk * gqyy[28]
1265 + qzzk * gqyy[30]
1266 + 2.0 * (qxyk * gqyy[26] + qxzk * gqyy[27] + qyzk * gqyy[29]))
1267 + qzzi
1268 * (qxxk * gqzz[25]
1269 + qyyk * gqzz[28]
1270 + qzzk * gqzz[30]
1271 + 2.0 * (qxyk * gqzz[26] + qxzk * gqzz[27] + qyzk * gqzz[29]))
1272 + 2.0
1273 * (qxyi
1274 * (qxxk * gqxy[25]
1275 + qyyk * gqxy[28]
1276 + qzzk * gqxy[30]
1277 + 2.0 * (qxyk * gqxy[26] + qxzk * gqxy[27] + qyzk * gqxy[29]))
1278 + qxzi
1279 * (qxxk * gqxz[25]
1280 + qyyk * gqxz[28]
1281 + qzzk * gqxz[30]
1282 + 2.0 * (qxyk * gqxz[26] + qxzk * gqxz[27] + qyzk * gqxz[29]))
1283 + qyzi
1284 * (qxxk * gqyz[25]
1285 + qyyk * gqyz[28]
1286 + qzzk * gqyz[30]
1287 + 2.0 * (qxyk * gqyz[26] + qxzk * gqyz[27] + qyzk * gqyz[29])));
1288 final double dewkdr =
1289 ci * (uxk * gux[21] + uyk * guy[21] + uzk * guz[21])
1290 - ck * (uxi * gc[22] + uyi * gc[23] + uzi * gc[24])
1291 + ci
1292 * (qxxk * gqxx[21]
1293 + qyyk * gqyy[21]
1294 + qzzk * gqzz[21]
1295 + 2.0 * (qxyk * gqxy[21] + qxzk * gqxz[21] + qyzk * gqyz[21]))
1296 + ck
1297 * (qxxi * gc[25]
1298 + qyyi * gc[28]
1299 + qzzi * gc[30]
1300 + 2.0 * (qxyi * gc[26] + qxzi * gc[27] + qyzi * gc[29]))
1301 - uxi
1302 * (qxxk * gqxx[22]
1303 + qyyk * gqyy[22]
1304 + qzzk * gqzz[22]
1305 + 2.0 * (qxyk * gqxy[22] + qxzk * gqxz[22] + qyzk * gqyz[22]))
1306 - uyi
1307 * (qxxk * gqxx[23]
1308 + qyyk * gqyy[23]
1309 + qzzk * gqzz[23]
1310 + 2.0 * (qxyk * gqxy[23] + qxzk * gqxz[23] + qyzk * gqyz[23]))
1311 - uzi
1312 * (qxxk * gqxx[24]
1313 + qyyk * gqyy[24]
1314 + qzzk * gqzz[24]
1315 + 2.0 * (qxyk * gqxy[24] + qxzk * gqxz[24] + qyzk * gqyz[24]))
1316 + uxk
1317 * (qxxi * gux[25]
1318 + qyyi * gux[28]
1319 + qzzi * gux[30]
1320 + 2.0 * (qxyi * gux[26] + qxzi * gux[27] + qyzi * gux[29]))
1321 + uyk
1322 * (qxxi * guy[25]
1323 + qyyi * guy[28]
1324 + qzzi * guy[30]
1325 + 2.0 * (qxyi * guy[26] + qxzi * guy[27] + qyzi * guy[29]))
1326 + uzk
1327 * (qxxi * guz[25]
1328 + qyyi * guz[28]
1329 + qzzi * guz[30]
1330 + 2.0 * (qxyi * guz[26] + qxzi * guz[27] + qyzi * guz[29]))
1331 + qxxi
1332 * (qxxk * gqxx[25]
1333 + qyyk * gqyy[25]
1334 + qzzk * gqzz[25]
1335 + 2.0 * (qxyk * gqxy[25] + qxzk * gqxz[25] + qyzk * gqyz[25]))
1336 + qyyi
1337 * (qxxk * gqxx[28]
1338 + qyyk * gqyy[28]
1339 + qzzk * gqzz[28]
1340 + 2.0 * (qxyk * gqxy[28] + qxzk * gqxz[28] + qyzk * gqyz[28]))
1341 + qzzi
1342 * (qxxk * gqxx[30]
1343 + qyyk * gqyy[30]
1344 + qzzk * gqzz[30]
1345 + 2.0 * (qxyk * gqxy[30] + qxzk * gqxz[30] + qyzk * gqyz[30]))
1346 + 2.0
1347 * (qxyi
1348 * (qxxk * gqxx[26]
1349 + qyyk * gqyy[26]
1350 + qzzk * gqzz[26]
1351 + 2.0 * (qxyk * gqxy[26] + qxzk * gqxz[26] + qyzk * gqyz[26]))
1352 + qxzi
1353 * (qxxk * gqxx[27]
1354 + qyyk * gqyy[27]
1355 + qzzk * gqzz[27]
1356 + 2.0 * (qxyk * gqxy[27] + qxzk * gqxz[27] + qyzk * gqyz[27]))
1357 + qyzi
1358 * (qxxk * gqxx[29]
1359 + qyyk * gqyy[29]
1360 + qzzk * gqzz[29]
1361 + 2.0 * (qxyk * gqxy[29] + qxzk * gqxz[29] + qyzk * gqyz[29])));
1362 final double dsumdr = desymdr + 0.5 * (dewidr + dewkdr);
1363 final double drbi = rbk * dsumdr;
1364 final double drbk = rbi * dsumdr;
1365
1366 double selfScale = 1.0;
1367 if (i == k) {
1368 if (iSymm == 0) {
1369 sharedBornGrad.add(threadID, i, drbi);
1370 return;
1371 } else {
1372 selfScale = 0.5;
1373 }
1374 }
1375
1376
1377 final double dedx = selfScale * dEdX();
1378 final double dedy = selfScale * dEdY();
1379 final double dedz = selfScale * dEdZ();
1380 dedxi -= dedx;
1381 dedyi -= dedy;
1382 dedzi -= dedz;
1383 dborni += selfScale * drbi;
1384
1385 final double dedxk = dedx * transOp[0][0] + dedy * transOp[1][0] + dedz * transOp[2][0];
1386 final double dedyk = dedx * transOp[0][1] + dedy * transOp[1][1] + dedz * transOp[2][1];
1387 final double dedzk = dedx * transOp[0][2] + dedy * transOp[1][2] + dedz * transOp[2][2];
1388 grad.add(threadID, k, dedxk, dedyk, dedzk);
1389 sharedBornGrad.add(threadID, k, selfScale * drbk);
1390 permanentEnergyTorque(i, k);
1391 }
1392
1393 private double dEdZ() {
1394 final double desymdz =
1395 ci * ck * gc[4]
1396 - (uxi * (uxk * gux[7] + uyk * guy[7] + uzk * guz[7])
1397 + uyi * (uxk * gux[9] + uyk * guy[9] + uzk * guz[9])
1398 + uzi * (uxk * gux[10] + uyk * guy[10] + uzk * guz[10]));
1399 final double dewidz =
1400 ci * (uxk * gc[7] + uyk * gc[9] + uzk * gc[10])
1401 - ck * (uxi * gux[4] + uyi * guy[4] + uzi * guz[4])
1402 + ci
1403 * (qxxk * gc[13]
1404 + qyyk * gc[18]
1405 + qzzk * gc[20]
1406 + 2.0 * (qxyk * gc[15] + qxzk * gc[16] + qyzk * gc[19]))
1407 + ck
1408 * (qxxi * gqxx[4]
1409 + qyyi * gqyy[4]
1410 + qzzi * gqzz[4]
1411 + 2.0 * (qxyi * gqxy[4] + qxzi * gqxz[4] + qyzi * gqyz[4]))
1412 - uxi
1413 * (qxxk * gux[13]
1414 + qyyk * gux[18]
1415 + qzzk * gux[20]
1416 + 2.0 * (qxyk * gux[15] + qxzk * gux[16] + qyzk * gux[19]))
1417 - uyi
1418 * (qxxk * guy[13]
1419 + qyyk * guy[18]
1420 + qzzk * guy[20]
1421 + 2.0 * (qxyk * guy[15] + qxzk * guy[16] + qyzk * guy[19]))
1422 - uzi
1423 * (qxxk * guz[13]
1424 + qyyk * guz[18]
1425 + qzzk * guz[20]
1426 + 2.0 * (qxyk * guz[15] + qxzk * guz[16] + qyzk * guz[19]))
1427 + uxk
1428 * (qxxi * gqxx[7]
1429 + qyyi * gqyy[7]
1430 + qzzi * gqzz[7]
1431 + 2.0 * (qxyi * gqxy[7] + qxzi * gqxz[7] + qyzi * gqyz[7]))
1432 + uyk
1433 * (qxxi * gqxx[9]
1434 + qyyi * gqyy[9]
1435 + qzzi * gqzz[9]
1436 + 2.0 * (qxyi * gqxy[9] + qxzi * gqxz[9] + qyzi * gqyz[9]))
1437 + uzk
1438 * (qxxi * gqxx[10]
1439 + qyyi * gqyy[10]
1440 + qzzi * gqzz[10]
1441 + 2.0 * (qxyi * gqxy[10] + qxzi * gqxz[10] + qyzi * gqyz[10]))
1442 + qxxi
1443 * (qxxk * gqxx[13]
1444 + qyyk * gqxx[18]
1445 + qzzk * gqxx[20]
1446 + 2.0 * (qxyk * gqxx[15] + qxzk * gqxx[16] + qyzk * gqxx[19]))
1447 + qyyi
1448 * (qxxk * gqyy[13]
1449 + qyyk * gqyy[18]
1450 + qzzk * gqyy[20]
1451 + 2.0 * (qxyk * gqyy[15] + qxzk * gqyy[16] + qyzk * gqyy[19]))
1452 + qzzi
1453 * (qxxk * gqzz[13]
1454 + qyyk * gqzz[18]
1455 + qzzk * gqzz[20]
1456 + 2.0 * (qxyk * gqzz[15] + qxzk * gqzz[16] + qyzk * gqzz[19]))
1457 + 2.0
1458 * (qxyi
1459 * (qxxk * gqxy[13]
1460 + qyyk * gqxy[18]
1461 + qzzk * gqxy[20]
1462 + 2.0 * (qxyk * gqxy[15] + qxzk * gqxy[16] + qyzk * gqxy[19]))
1463 + qxzi
1464 * (qxxk * gqxz[13]
1465 + qyyk * gqxz[18]
1466 + qzzk * gqxz[20]
1467 + 2.0 * (qxyk * gqxz[15] + qxzk * gqxz[16] + qyzk * gqxz[19]))
1468 + qyzi
1469 * (qxxk * gqyz[13]
1470 + qyyk * gqyz[18]
1471 + qzzk * gqyz[20]
1472 + 2.0 * (qxyk * gqyz[15] + qxzk * gqyz[16] + qyzk * gqyz[19])));
1473 final double dewkdz =
1474 ci * (uxk * gux[4] + uyk * guy[4] + uzk * guz[4])
1475 - ck * (uxi * gc[7] + uyi * gc[9] + uzi * gc[10])
1476 + ci
1477 * (qxxk * gqxx[4]
1478 + qyyk * gqyy[4]
1479 + qzzk * gqzz[4]
1480 + 2.0 * (qxyk * gqxy[4] + qxzk * gqxz[4] + qyzk * gqyz[4]))
1481 + ck
1482 * (qxxi * gc[13]
1483 + qyyi * gc[18]
1484 + qzzi * gc[20]
1485 + 2.0 * (qxyi * gc[15] + qxzi * gc[16] + qyzi * gc[19]))
1486 - uxi
1487 * (qxxk * gqxx[7]
1488 + qyyk * gqyy[7]
1489 + qzzk * gqzz[7]
1490 + 2.0 * (qxyk * gqxy[7] + qxzk * gqxz[7] + qyzk * gqyz[7]))
1491 - uyi
1492 * (qxxk * gqxx[9]
1493 + qyyk * gqyy[9]
1494 + qzzk * gqzz[9]
1495 + 2.0 * (qxyk * gqxy[9] + qxzk * gqxz[9] + qyzk * gqyz[9]))
1496 - uzi
1497 * (qxxk * gqxx[10]
1498 + qyyk * gqyy[10]
1499 + qzzk * gqzz[10]
1500 + 2.0 * (qxyk * gqxy[10] + qxzk * gqxz[10] + qyzk * gqyz[10]))
1501 + uxk
1502 * (qxxi * gux[13]
1503 + qyyi * gux[18]
1504 + qzzi * gux[20]
1505 + 2.0 * (qxyi * gux[15] + qxzi * gux[16] + qyzi * gux[19]))
1506 + uyk
1507 * (qxxi * guy[13]
1508 + qyyi * guy[18]
1509 + qzzi * guy[20]
1510 + 2.0 * (qxyi * guy[15] + qxzi * guy[16] + qyzi * guy[19]))
1511 + uzk
1512 * (qxxi * guz[13]
1513 + qyyi * guz[18]
1514 + qzzi * guz[20]
1515 + 2.0 * (qxyi * guz[15] + qxzi * guz[16] + qyzi * guz[19]))
1516 + qxxi
1517 * (qxxk * gqxx[13]
1518 + qyyk * gqyy[13]
1519 + qzzk * gqzz[13]
1520 + 2.0 * (qxyk * gqxy[13] + qxzk * gqxz[13] + qyzk * gqyz[13]))
1521 + qyyi
1522 * (qxxk * gqxx[18]
1523 + qyyk * gqyy[18]
1524 + qzzk * gqzz[18]
1525 + 2.0 * (qxyk * gqxy[18] + qxzk * gqxz[18] + qyzk * gqyz[18]))
1526 + qzzi
1527 * (qxxk * gqxx[20]
1528 + qyyk * gqyy[20]
1529 + qzzk * gqzz[20]
1530 + 2.0 * (qxyk * gqxy[20] + qxzk * gqxz[20] + qyzk * gqyz[20]))
1531 + 2.0
1532 * (qxyi
1533 * (qxxk * gqxx[15]
1534 + qyyk * gqyy[15]
1535 + qzzk * gqzz[15]
1536 + 2.0 * (qxyk * gqxy[15] + qxzk * gqxz[15] + qyzk * gqyz[15]))
1537 + qxzi
1538 * (qxxk * gqxx[16]
1539 + qyyk * gqyy[16]
1540 + qzzk * gqzz[16]
1541 + 2.0 * (qxyk * gqxy[16] + qxzk * gqxz[16] + qyzk * gqyz[16]))
1542 + qyzi
1543 * (qxxk * gqxx[19]
1544 + qyyk * gqyy[19]
1545 + qzzk * gqzz[19]
1546 + 2.0 * (qxyk * gqxy[19] + qxzk * gqxz[19] + qyzk * gqyz[19])));
1547 return desymdz + 0.5 * (dewidz + dewkdz);
1548 }
1549
1550 private double dEdY() {
1551 final double desymdy =
1552 ci * ck * gc[3]
1553 - (uxi * (uxk * gux[6] + uyk * guy[6] + uzk * guz[6])
1554 + uyi * (uxk * gux[8] + uyk * guy[8] + uzk * guz[8])
1555 + uzi * (uxk * gux[9] + uyk * guy[9] + uzk * guz[9]));
1556 final double dewidy =
1557 ci * (uxk * gc[6] + uyk * gc[8] + uzk * gc[9])
1558 - ck * (uxi * gux[3] + uyi * guy[3] + uzi * guz[3])
1559 + ci
1560 * (qxxk * gc[12]
1561 + qyyk * gc[17]
1562 + qzzk * gc[19]
1563 + 2.0 * (qxyk * gc[14] + qxzk * gc[15] + qyzk * gc[18]))
1564 + ck
1565 * (qxxi * gqxx[3]
1566 + qyyi * gqyy[3]
1567 + qzzi * gqzz[3]
1568 + 2.0 * (qxyi * gqxy[3] + qxzi * gqxz[3] + qyzi * gqyz[3]))
1569 - uxi
1570 * (qxxk * gux[12]
1571 + qyyk * gux[17]
1572 + qzzk * gux[19]
1573 + 2.0 * (qxyk * gux[14] + qxzk * gux[15] + qyzk * gux[18]))
1574 - uyi
1575 * (qxxk * guy[12]
1576 + qyyk * guy[17]
1577 + qzzk * guy[19]
1578 + 2.0 * (qxyk * guy[14] + qxzk * guy[15] + qyzk * guy[18]))
1579 - uzi
1580 * (qxxk * guz[12]
1581 + qyyk * guz[17]
1582 + qzzk * guz[19]
1583 + 2.0 * (qxyk * guz[14] + qxzk * guz[15] + qyzk * guz[18]))
1584 + uxk
1585 * (qxxi * gqxx[6]
1586 + qyyi * gqyy[6]
1587 + qzzi * gqzz[6]
1588 + 2.0 * (qxyi * gqxy[6] + qxzi * gqxz[6] + qyzi * gqyz[6]))
1589 + uyk
1590 * (qxxi * gqxx[8]
1591 + qyyi * gqyy[8]
1592 + qzzi * gqzz[8]
1593 + 2.0 * (qxyi * gqxy[8] + qxzi * gqxz[8] + qyzi * gqyz[8]))
1594 + uzk
1595 * (qxxi * gqxx[9]
1596 + qyyi * gqyy[9]
1597 + qzzi * gqzz[9]
1598 + 2.0 * (qxyi * gqxy[9] + qxzi * gqxz[9] + qyzi * gqyz[9]))
1599 + qxxi
1600 * (qxxk * gqxx[12]
1601 + qyyk * gqxx[17]
1602 + qzzk * gqxx[19]
1603 + 2.0 * (qxyk * gqxx[14] + qxzk * gqxx[15] + qyzk * gqxx[18]))
1604 + qyyi
1605 * (qxxk * gqyy[12]
1606 + qyyk * gqyy[17]
1607 + qzzk * gqyy[19]
1608 + 2.0 * (qxyk * gqyy[14] + qxzk * gqyy[15] + qyzk * gqyy[18]))
1609 + qzzi
1610 * (qxxk * gqzz[12]
1611 + qyyk * gqzz[17]
1612 + qzzk * gqzz[19]
1613 + 2.0 * (qxyk * gqzz[14] + qxzk * gqzz[15] + qyzk * gqzz[18]))
1614 + 2.0
1615 * (qxyi
1616 * (qxxk * gqxy[12]
1617 + qyyk * gqxy[17]
1618 + qzzk * gqxy[19]
1619 + 2.0 * (qxyk * gqxy[14] + qxzk * gqxy[15] + qyzk * gqxy[18]))
1620 + qxzi
1621 * (qxxk * gqxz[12]
1622 + qyyk * gqxz[17]
1623 + qzzk * gqxz[19]
1624 + 2.0 * (qxyk * gqxz[14] + qxzk * gqxz[15] + qyzk * gqxz[18]))
1625 + qyzi
1626 * (qxxk * gqyz[12]
1627 + qyyk * gqyz[17]
1628 + qzzk * gqyz[19]
1629 + 2.0 * (qxyk * gqyz[14] + qxzk * gqyz[15] + qyzk * gqyz[18])));
1630 final double dewkdy =
1631 ci * (uxk * gux[3] + uyk * guy[3] + uzk * guz[3])
1632 - ck * (uxi * gc[6] + uyi * gc[8] + uzi * gc[9])
1633 + ci
1634 * (qxxk * gqxx[3]
1635 + qyyk * gqyy[3]
1636 + qzzk * gqzz[3]
1637 + 2.0 * (qxyk * gqxy[3] + qxzk * gqxz[3] + qyzk * gqyz[3]))
1638 + ck
1639 * (qxxi * gc[12]
1640 + qyyi * gc[17]
1641 + qzzi * gc[19]
1642 + 2.0 * (qxyi * gc[14] + qxzi * gc[15] + qyzi * gc[18]))
1643 - uxi
1644 * (qxxk * gqxx[6]
1645 + qyyk * gqyy[6]
1646 + qzzk * gqzz[6]
1647 + 2.0 * (qxyk * gqxy[6] + qxzk * gqxz[6] + qyzk * gqyz[6]))
1648 - uyi
1649 * (qxxk * gqxx[8]
1650 + qyyk * gqyy[8]
1651 + qzzk * gqzz[8]
1652 + 2.0 * (qxyk * gqxy[8] + qxzk * gqxz[8] + qyzk * gqyz[8]))
1653 - uzi
1654 * (qxxk * gqxx[9]
1655 + qyyk * gqyy[9]
1656 + qzzk * gqzz[9]
1657 + 2.0 * (qxyk * gqxy[9] + qxzk * gqxz[9] + qyzk * gqyz[9]))
1658 + uxk
1659 * (qxxi * gux[12]
1660 + qyyi * gux[17]
1661 + qzzi * gux[19]
1662 + 2.0 * (qxyi * gux[14] + qxzi * gux[15] + qyzi * gux[18]))
1663 + uyk
1664 * (qxxi * guy[12]
1665 + qyyi * guy[17]
1666 + qzzi * guy[19]
1667 + 2.0 * (qxyi * guy[14] + qxzi * guy[15] + qyzi * guy[18]))
1668 + uzk
1669 * (qxxi * guz[12]
1670 + qyyi * guz[17]
1671 + qzzi * guz[19]
1672 + 2.0 * (qxyi * guz[14] + qxzi * guz[15] + qyzi * guz[18]))
1673 + qxxi
1674 * (qxxk * gqxx[12]
1675 + qyyk * gqyy[12]
1676 + qzzk * gqzz[12]
1677 + 2.0 * (qxyk * gqxy[12] + qxzk * gqxz[12] + qyzk * gqyz[12]))
1678 + qyyi
1679 * (qxxk * gqxx[17]
1680 + qyyk * gqyy[17]
1681 + qzzk * gqzz[17]
1682 + 2.0 * (qxyk * gqxy[17] + qxzk * gqxz[17] + qyzk * gqyz[17]))
1683 + qzzi
1684 * (qxxk * gqxx[19]
1685 + qyyk * gqyy[19]
1686 + qzzk * gqzz[19]
1687 + 2.0 * (qxyk * gqxy[19] + qxzk * gqxz[19] + qyzk * gqyz[19]))
1688 + 2.0
1689 * (qxyi
1690 * (qxxk * gqxx[14]
1691 + qyyk * gqyy[14]
1692 + qzzk * gqzz[14]
1693 + 2.0 * (qxyk * gqxy[14] + qxzk * gqxz[14] + qyzk * gqyz[14]))
1694 + qxzi
1695 * (qxxk * gqxx[15]
1696 + qyyk * gqyy[15]
1697 + qzzk * gqzz[15]
1698 + 2.0 * (qxyk * gqxy[15] + qxzk * gqxz[15] + qyzk * gqyz[15]))
1699 + qyzi
1700 * (qxxk * gqxx[18]
1701 + qyyk * gqyy[18]
1702 + qzzk * gqzz[18]
1703 + 2.0 * (qxyk * gqxy[18] + qxzk * gqxz[18] + qyzk * gqyz[18])));
1704
1705 return desymdy + 0.5 * (dewidy + dewkdy);
1706 }
1707
1708 private double dEdX() {
1709 final double desymdx =
1710 ci * ck * gc[2]
1711 - (uxi * (uxk * gux[5] + uyk * guy[5] + uzk * guz[5])
1712 + uyi * (uxk * gux[6] + uyk * guy[6] + uzk * guz[6])
1713 + uzi * (uxk * gux[7] + uyk * guy[7] + uzk * guz[7]));
1714 final double dewidx =
1715 ci * (uxk * gc[5] + uyk * gc[6] + uzk * gc[7])
1716 - ck * (uxi * gux[2] + uyi * guy[2] + uzi * guz[2])
1717 + ci
1718 * (qxxk * gc[11]
1719 + qyyk * gc[14]
1720 + qzzk * gc[16]
1721 + 2.0 * (qxyk * gc[12] + qxzk * gc[13] + qyzk * gc[15]))
1722 + ck
1723 * (qxxi * gqxx[2]
1724 + qyyi * gqyy[2]
1725 + qzzi * gqzz[2]
1726 + 2.0 * (qxyi * gqxy[2] + qxzi * gqxz[2] + qyzi * gqyz[2]))
1727 - uxi
1728 * (qxxk * gux[11]
1729 + qyyk * gux[14]
1730 + qzzk * gux[16]
1731 + 2.0 * (qxyk * gux[12] + qxzk * gux[13] + qyzk * gux[15]))
1732 - uyi
1733 * (qxxk * guy[11]
1734 + qyyk * guy[14]
1735 + qzzk * guy[16]
1736 + 2.0 * (qxyk * guy[12] + qxzk * guy[13] + qyzk * guy[15]))
1737 - uzi
1738 * (qxxk * guz[11]
1739 + qyyk * guz[14]
1740 + qzzk * guz[16]
1741 + 2.0 * (qxyk * guz[12] + qxzk * guz[13] + qyzk * guz[15]))
1742 + uxk
1743 * (qxxi * gqxx[5]
1744 + qyyi * gqyy[5]
1745 + qzzi * gqzz[5]
1746 + 2.0 * (qxyi * gqxy[5] + qxzi * gqxz[5] + qyzi * gqyz[5]))
1747 + uyk
1748 * (qxxi * gqxx[6]
1749 + qyyi * gqyy[6]
1750 + qzzi * gqzz[6]
1751 + 2.0 * (qxyi * gqxy[6] + qxzi * gqxz[6] + qyzi * gqyz[6]))
1752 + uzk
1753 * (qxxi * gqxx[7]
1754 + qyyi * gqyy[7]
1755 + qzzi * gqzz[7]
1756 + 2.0 * (qxyi * gqxy[7] + qxzi * gqxz[7] + qyzi * gqyz[7]))
1757 + qxxi
1758 * (qxxk * gqxx[11]
1759 + qyyk * gqxx[14]
1760 + qzzk * gqxx[16]
1761 + 2.0 * (qxyk * gqxx[12] + qxzk * gqxx[13] + qyzk * gqxx[15]))
1762 + qyyi
1763 * (qxxk * gqyy[11]
1764 + qyyk * gqyy[14]
1765 + qzzk * gqyy[16]
1766 + 2.0 * (qxyk * gqyy[12] + qxzk * gqyy[13] + qyzk * gqyy[15]))
1767 + qzzi
1768 * (qxxk * gqzz[11]
1769 + qyyk * gqzz[14]
1770 + qzzk * gqzz[16]
1771 + 2.0 * (qxyk * gqzz[12] + qxzk * gqzz[13] + qyzk * gqzz[15]))
1772 + 2.0
1773 * (qxyi
1774 * (qxxk * gqxy[11]
1775 + qyyk * gqxy[14]
1776 + qzzk * gqxy[16]
1777 + 2.0 * (qxyk * gqxy[12] + qxzk * gqxy[13] + qyzk * gqxy[15]))
1778 + qxzi
1779 * (qxxk * gqxz[11]
1780 + qyyk * gqxz[14]
1781 + qzzk * gqxz[16]
1782 + 2.0 * (qxyk * gqxz[12] + qxzk * gqxz[13] + qyzk * gqxz[15]))
1783 + qyzi
1784 * (qxxk * gqyz[11]
1785 + qyyk * gqyz[14]
1786 + qzzk * gqyz[16]
1787 + 2.0 * (qxyk * gqyz[12] + qxzk * gqyz[13] + qyzk * gqyz[15])));
1788 final double dewkdx =
1789 ci * (uxk * gux[2] + uyk * guy[2] + uzk * guz[2])
1790 - ck * (uxi * gc[5] + uyi * gc[6] + uzi * gc[7])
1791 + ci
1792 * (qxxk * gqxx[2]
1793 + qyyk * gqyy[2]
1794 + qzzk * gqzz[2]
1795 + 2.0 * (qxyk * gqxy[2] + qxzk * gqxz[2] + qyzk * gqyz[2]))
1796 + ck
1797 * (qxxi * gc[11]
1798 + qyyi * gc[14]
1799 + qzzi * gc[16]
1800 + 2.0 * (qxyi * gc[12] + qxzi * gc[13] + qyzi * gc[15]))
1801 - uxi
1802 * (qxxk * gqxx[5]
1803 + qyyk * gqyy[5]
1804 + qzzk * gqzz[5]
1805 + 2.0 * (qxyk * gqxy[5] + qxzk * gqxz[5] + qyzk * gqyz[5]))
1806 - uyi
1807 * (qxxk * gqxx[6]
1808 + qyyk * gqyy[6]
1809 + qzzk * gqzz[6]
1810 + 2.0 * (qxyk * gqxy[6] + qxzk * gqxz[6] + qyzk * gqyz[6]))
1811 - uzi
1812 * (qxxk * gqxx[7]
1813 + qyyk * gqyy[7]
1814 + qzzk * gqzz[7]
1815 + 2.0 * (qxyk * gqxy[7] + qxzk * gqxz[7] + qyzk * gqyz[7]))
1816 + uxk
1817 * (qxxi * gux[11]
1818 + qyyi * gux[14]
1819 + qzzi * gux[16]
1820 + 2.0 * (qxyi * gux[12] + qxzi * gux[13] + qyzi * gux[15]))
1821 + uyk
1822 * (qxxi * guy[11]
1823 + qyyi * guy[14]
1824 + qzzi * guy[16]
1825 + 2.0 * (qxyi * guy[12] + qxzi * guy[13] + qyzi * guy[15]))
1826 + uzk
1827 * (qxxi * guz[11]
1828 + qyyi * guz[14]
1829 + qzzi * guz[16]
1830 + 2.0 * (qxyi * guz[12] + qxzi * guz[13] + qyzi * guz[15]))
1831 + qxxi
1832 * (qxxk * gqxx[11]
1833 + qyyk * gqyy[11]
1834 + qzzk * gqzz[11]
1835 + 2.0 * (qxyk * gqxy[11] + qxzk * gqxz[11] + qyzk * gqyz[11]))
1836 + qyyi
1837 * (qxxk * gqxx[14]
1838 + qyyk * gqyy[14]
1839 + qzzk * gqzz[14]
1840 + 2.0 * (qxyk * gqxy[14] + qxzk * gqxz[14] + qyzk * gqyz[14]))
1841 + qzzi
1842 * (qxxk * gqxx[16]
1843 + qyyk * gqyy[16]
1844 + qzzk * gqzz[16]
1845 + 2.0 * (qxyk * gqxy[16] + qxzk * gqxz[16] + qyzk * gqyz[16]))
1846 + 2.0
1847 * (qxyi
1848 * (qxxk * gqxx[12]
1849 + qyyk * gqyy[12]
1850 + qzzk * gqzz[12]
1851 + 2.0 * (qxyk * gqxy[12] + qxzk * gqxz[12] + qyzk * gqyz[12]))
1852 + qxzi
1853 * (qxxk * gqxx[13]
1854 + qyyk * gqyy[13]
1855 + qzzk * gqzz[13]
1856 + 2.0 * (qxyk * gqxy[13] + qxzk * gqxz[13] + qyzk * gqyz[13]))
1857 + qyzi
1858 * (qxxk * gqxx[15]
1859 + qyyk * gqyy[15]
1860 + qzzk * gqzz[15]
1861 + 2.0 * (qxyk * gqxy[15] + qxzk * gqxz[15] + qyzk * gqyz[15])));
1862
1863 return desymdx + 0.5 * (dewidx + dewkdx);
1864 }
1865
1866 private void permanentEnergyTorque(int i, int k) {
1867
1868 final double ix =
1869 uxk * gux[2]
1870 + uyk * gux[3]
1871 + uzk * gux[4]
1872 + 0.5
1873 * (ck * gux[1]
1874 + qxxk * gux[5]
1875 + qyyk * gux[8]
1876 + qzzk * gux[10]
1877 + 2.0 * (qxyk * gux[6] + qxzk * gux[7] + qyzk * gux[9])
1878 + ck * gc[2]
1879 + qxxk * gqxx[2]
1880 + qyyk * gqyy[2]
1881 + qzzk * gqzz[2]
1882 + 2.0 * (qxyk * gqxy[2] + qxzk * gqxz[2] + qyzk * gqyz[2]));
1883 final double iy =
1884 uxk * guy[2]
1885 + uyk * guy[3]
1886 + uzk * guy[4]
1887 + 0.5
1888 * (ck * guy[1]
1889 + qxxk * guy[5]
1890 + qyyk * guy[8]
1891 + qzzk * guy[10]
1892 + 2.0 * (qxyk * guy[6] + qxzk * guy[7] + qyzk * guy[9])
1893 + ck * gc[3]
1894 + qxxk * gqxx[3]
1895 + qyyk * gqyy[3]
1896 + qzzk * gqzz[3]
1897 + 2.0 * (qxyk * gqxy[3] + qxzk * gqxz[3] + qyzk * gqyz[3]));
1898 final double iz =
1899 uxk * guz[2]
1900 + uyk * guz[3]
1901 + uzk * guz[4]
1902 + 0.5
1903 * (ck * guz[1]
1904 + qxxk * guz[5]
1905 + qyyk * guz[8]
1906 + qzzk * guz[10]
1907 + 2.0 * (qxyk * guz[6] + qxzk * guz[7] + qyzk * guz[9])
1908 + ck * gc[4]
1909 + qxxk * gqxx[4]
1910 + qyyk * gqyy[4]
1911 + qzzk * gqzz[4]
1912 + 2.0 * (qxyk * gqxy[4] + qxzk * gqxz[4] + qyzk * gqyz[4]));
1913 final double kx =
1914 uxi * gux[2]
1915 + uyi * gux[3]
1916 + uzi * gux[4]
1917 - 0.5
1918 * (ci * gux[1]
1919 + qxxi * gux[5]
1920 + qyyi * gux[8]
1921 + qzzi * gux[10]
1922 + 2.0 * (qxyi * gux[6] + qxzi * gux[7] + qyzi * gux[9])
1923 + ci * gc[2]
1924 + qxxi * gqxx[2]
1925 + qyyi * gqyy[2]
1926 + qzzi * gqzz[2]
1927 + 2.0 * (qxyi * gqxy[2] + qxzi * gqxz[2] + qyzi * gqyz[2]));
1928 final double ky =
1929 uxi * guy[2]
1930 + uyi * guy[3]
1931 + uzi * guy[4]
1932 - 0.5
1933 * (ci * guy[1]
1934 + qxxi * guy[5]
1935 + qyyi * guy[8]
1936 + qzzi * guy[10]
1937 + 2.0 * (qxyi * guy[6] + qxzi * guy[7] + qyzi * guy[9])
1938 + ci * gc[3]
1939 + qxxi * gqxx[3]
1940 + qyyi * gqyy[3]
1941 + qzzi * gqzz[3]
1942 + 2.0 * (qxyi * gqxy[3] + qxzi * gqxz[3] + qyzi * gqyz[3]));
1943 final double kz =
1944 uxi * guz[2]
1945 + uyi * guz[3]
1946 + uzi * guz[4]
1947 - 0.5
1948 * (ci * guz[1]
1949 + qxxi * guz[5]
1950 + qyyi * guz[8]
1951 + qzzi * guz[10]
1952 + 2.0 * (qxyi * guz[6] + qxzi * guz[7] + qyzi * guz[9])
1953 + ci * gc[4]
1954 + qxxi * gqxx[4]
1955 + qyyi * gqyy[4]
1956 + qzzi * gqzz[4]
1957 + 2.0 * (qxyi * gqxy[4] + qxzi * gqxz[4] + qyzi * gqyz[4]));
1958 double tix = uyi * iz - uzi * iy;
1959 double tiy = uzi * ix - uxi * iz;
1960 double tiz = uxi * iy - uyi * ix;
1961 double tkx = uyk * kz - uzk * ky;
1962 double tky = uzk * kx - uxk * kz;
1963 double tkz = uxk * ky - uyk * kx;
1964
1965
1966 final double ixx =
1967 -0.5
1968 * (ck * gqxx[1]
1969 + uxk * gqxx[2]
1970 + uyk * gqxx[3]
1971 + uzk * gqxx[4]
1972 + qxxk * gqxx[5]
1973 + qyyk * gqxx[8]
1974 + qzzk * gqxx[10]
1975 + 2.0 * (qxyk * gqxx[6] + qxzk * gqxx[7] + qyzk * gqxx[9])
1976 + ck * gc[5]
1977 + uxk * gux[5]
1978 + uyk * guy[5]
1979 + uzk * guz[5]
1980 + qxxk * gqxx[5]
1981 + qyyk * gqyy[5]
1982 + qzzk * gqzz[5]
1983 + 2.0 * (qxyk * gqxy[5] + qxzk * gqxz[5] + qyzk * gqyz[5]));
1984 final double ixy =
1985 -0.5
1986 * (ck * gqxy[1]
1987 + uxk * gqxy[2]
1988 + uyk * gqxy[3]
1989 + uzk * gqxy[4]
1990 + qxxk * gqxy[5]
1991 + qyyk * gqxy[8]
1992 + qzzk * gqxy[10]
1993 + 2.0 * (qxyk * gqxy[6] + qxzk * gqxy[7] + qyzk * gqxy[9])
1994 + ck * gc[6]
1995 + uxk * gux[6]
1996 + uyk * guy[6]
1997 + uzk * guz[6]
1998 + qxxk * gqxx[6]
1999 + qyyk * gqyy[6]
2000 + qzzk * gqzz[6]
2001 + 2.0 * (qxyk * gqxy[6] + qxzk * gqxz[6] + qyzk * gqyz[6]));
2002 final double ixz =
2003 -0.5
2004 * (ck * gqxz[1]
2005 + uxk * gqxz[2]
2006 + uyk * gqxz[3]
2007 + uzk * gqxz[4]
2008 + qxxk * gqxz[5]
2009 + qyyk * gqxz[8]
2010 + qzzk * gqxz[10]
2011 + 2.0 * (qxyk * gqxz[6] + qxzk * gqxz[7] + qyzk * gqxz[9])
2012 + ck * gc[7]
2013 + uxk * gux[7]
2014 + uyk * guy[7]
2015 + uzk * guz[7]
2016 + qxxk * gqxx[7]
2017 + qyyk * gqyy[7]
2018 + qzzk * gqzz[7]
2019 + 2.0 * (qxyk * gqxy[7] + qxzk * gqxz[7] + qyzk * gqyz[7]));
2020 final double iyy =
2021 -0.5
2022 * (ck * gqyy[1]
2023 + uxk * gqyy[2]
2024 + uyk * gqyy[3]
2025 + uzk * gqyy[4]
2026 + qxxk * gqyy[5]
2027 + qyyk * gqyy[8]
2028 + qzzk * gqyy[10]
2029 + 2.0 * (qxyk * gqyy[6] + qxzk * gqyy[7] + qyzk * gqyy[9])
2030 + ck * gc[8]
2031 + uxk * gux[8]
2032 + uyk * guy[8]
2033 + uzk * guz[8]
2034 + qxxk * gqxx[8]
2035 + qyyk * gqyy[8]
2036 + qzzk * gqzz[8]
2037 + 2.0 * (qxyk * gqxy[8] + qxzk * gqxz[8] + qyzk * gqyz[8]));
2038 final double iyz =
2039 -0.5
2040 * (ck * gqyz[1]
2041 + uxk * gqyz[2]
2042 + uyk * gqyz[3]
2043 + uzk * gqyz[4]
2044 + qxxk * gqyz[5]
2045 + qyyk * gqyz[8]
2046 + qzzk * gqyz[10]
2047 + 2.0 * (qxyk * gqyz[6] + qxzk * gqyz[7] + qyzk * gqyz[9])
2048 + ck * gc[9]
2049 + uxk * gux[9]
2050 + uyk * guy[9]
2051 + uzk * guz[9]
2052 + qxxk * gqxx[9]
2053 + qyyk * gqyy[9]
2054 + qzzk * gqzz[9]
2055 + 2.0 * (qxyk * gqxy[9] + qxzk * gqxz[9] + qyzk * gqyz[9]));
2056 final double izz =
2057 -0.5
2058 * (ck * gqzz[1]
2059 + uxk * gqzz[2]
2060 + uyk * gqzz[3]
2061 + uzk * gqzz[4]
2062 + qxxk * gqzz[5]
2063 + qyyk * gqzz[8]
2064 + qzzk * gqzz[10]
2065 + 2.0 * (qxyk * gqzz[6] + qxzk * gqzz[7] + qyzk * gqzz[9])
2066 + ck * gc[10]
2067 + uxk * gux[10]
2068 + uyk * guy[10]
2069 + uzk * guz[10]
2070 + qxxk * gqxx[10]
2071 + qyyk * gqyy[10]
2072 + qzzk * gqzz[10]
2073 + 2.0 * (qxyk * gqxy[10] + qxzk * gqxz[10] + qyzk * gqyz[10]));
2074 final double iyx = ixy;
2075 final double izx = ixz;
2076 final double izy = iyz;
2077 final double kxx =
2078 -0.5
2079 * (ci * gqxx[1]
2080 - uxi * gqxx[2]
2081 - uyi * gqxx[3]
2082 - uzi * gqxx[4]
2083 + qxxi * gqxx[5]
2084 + qyyi * gqxx[8]
2085 + qzzi * gqxx[10]
2086 + 2.0 * (qxyi * gqxx[6] + qxzi * gqxx[7] + qyzi * gqxx[9])
2087 + ci * gc[5]
2088 - uxi * gux[5]
2089 - uyi * guy[5]
2090 - uzi * guz[5]
2091 + qxxi * gqxx[5]
2092 + qyyi * gqyy[5]
2093 + qzzi * gqzz[5]
2094 + 2.0 * (qxyi * gqxy[5] + qxzi * gqxz[5] + qyzi * gqyz[5]));
2095 double kxy =
2096 -0.5
2097 * (ci * gqxy[1]
2098 - uxi * gqxy[2]
2099 - uyi * gqxy[3]
2100 - uzi * gqxy[4]
2101 + qxxi * gqxy[5]
2102 + qyyi * gqxy[8]
2103 + qzzi * gqxy[10]
2104 + 2.0 * (qxyi * gqxy[6] + qxzi * gqxy[7] + qyzi * gqxy[9])
2105 + ci * gc[6]
2106 - uxi * gux[6]
2107 - uyi * guy[6]
2108 - uzi * guz[6]
2109 + qxxi * gqxx[6]
2110 + qyyi * gqyy[6]
2111 + qzzi * gqzz[6]
2112 + 2.0 * (qxyi * gqxy[6] + qxzi * gqxz[6] + qyzi * gqyz[6]));
2113 double kxz =
2114 -0.5
2115 * (ci * gqxz[1]
2116 - uxi * gqxz[2]
2117 - uyi * gqxz[3]
2118 - uzi * gqxz[4]
2119 + qxxi * gqxz[5]
2120 + qyyi * gqxz[8]
2121 + qzzi * gqxz[10]
2122 + 2.0 * (qxyi * gqxz[6] + qxzi * gqxz[7] + qyzi * gqxz[9])
2123 + ci * gc[7]
2124 - uxi * gux[7]
2125 - uyi * guy[7]
2126 - uzi * guz[7]
2127 + qxxi * gqxx[7]
2128 + qyyi * gqyy[7]
2129 + qzzi * gqzz[7]
2130 + 2.0 * (qxyi * gqxy[7] + qxzi * gqxz[7] + qyzi * gqyz[7]));
2131 double kyy =
2132 -0.5
2133 * (ci * gqyy[1]
2134 - uxi * gqyy[2]
2135 - uyi * gqyy[3]
2136 - uzi * gqyy[4]
2137 + qxxi * gqyy[5]
2138 + qyyi * gqyy[8]
2139 + qzzi * gqyy[10]
2140 + 2.0 * (qxyi * gqyy[6] + qxzi * gqyy[7] + qyzi * gqyy[9])
2141 + ci * gc[8]
2142 - uxi * gux[8]
2143 - uyi * guy[8]
2144 - uzi * guz[8]
2145 + qxxi * gqxx[8]
2146 + qyyi * gqyy[8]
2147 + qzzi * gqzz[8]
2148 + 2.0 * (qxyi * gqxy[8] + qxzi * gqxz[8] + qyzi * gqyz[8]));
2149 double kyz =
2150 -0.5
2151 * (ci * gqyz[1]
2152 - uxi * gqyz[2]
2153 - uyi * gqyz[3]
2154 - uzi * gqyz[4]
2155 + qxxi * gqyz[5]
2156 + qyyi * gqyz[8]
2157 + qzzi * gqyz[10]
2158 + 2.0 * (qxyi * gqyz[6] + qxzi * gqyz[7] + qyzi * gqyz[9])
2159 + ci * gc[9]
2160 - uxi * gux[9]
2161 - uyi * guy[9]
2162 - uzi * guz[9]
2163 + qxxi * gqxx[9]
2164 + qyyi * gqyy[9]
2165 + qzzi * gqzz[9]
2166 + 2.0 * (qxyi * gqxy[9] + qxzi * gqxz[9] + qyzi * gqyz[9]));
2167 double kzz =
2168 -0.5
2169 * (ci * gqzz[1]
2170 - uxi * gqzz[2]
2171 - uyi * gqzz[3]
2172 - uzi * gqzz[4]
2173 + qxxi * gqzz[5]
2174 + qyyi * gqzz[8]
2175 + qzzi * gqzz[10]
2176 + 2.0 * (qxyi * gqzz[6] + qxzi * gqzz[7] + qyzi * gqzz[9])
2177 + ci * gc[10]
2178 - uxi * gux[10]
2179 - uyi * guy[10]
2180 - uzi * guz[10]
2181 + qxxi * gqxx[10]
2182 + qyyi * gqyy[10]
2183 + qzzi * gqzz[10]
2184 + 2.0 * (qxyi * gqxy[10] + qxzi * gqxz[10] + qyzi * gqyz[10]));
2185 double kyx = kxy;
2186 double kzx = kxz;
2187 double kzy = kyz;
2188 tix += 2.0 * (qxyi * ixz + qyyi * iyz + qyzi * izz - qxzi * ixy - qyzi * iyy - qzzi * izy);
2189 tiy += 2.0 * (qxzi * ixx + qyzi * iyx + qzzi * izx - qxxi * ixz - qxyi * iyz - qxzi * izz);
2190 tiz += 2.0 * (qxxi * ixy + qxyi * iyy + qxzi * izy - qxyi * ixx - qyyi * iyx - qyzi * izx);
2191 tkx += 2.0 * (qxyk * kxz + qyyk * kyz + qyzk * kzz - qxzk * kxy - qyzk * kyy - qzzk * kzy);
2192 tky += 2.0 * (qxzk * kxx + qyzk * kyx + qzzk * kzx - qxxk * kxz - qxyk * kyz - qxzk * kzz);
2193 tkz += 2.0 * (qxxk * kxy + qxyk * kyy + qxzk * kzy - qxyk * kxx - qyyk * kyx - qyzk * kzx);
2194 if (i == k) {
2195 double selfScale = 0.5;
2196 tix *= selfScale;
2197 tiy *= selfScale;
2198 tiz *= selfScale;
2199 tkx *= selfScale;
2200 tky *= selfScale;
2201 tkz *= selfScale;
2202 }
2203 trqxi += tix;
2204 trqyi += tiy;
2205 trqzi += tiz;
2206
2207 final double rtkx = tkx * transOp[0][0] + tky * transOp[1][0] + tkz * transOp[2][0];
2208 final double rtky = tkx * transOp[0][1] + tky * transOp[1][1] + tkz * transOp[2][1];
2209 final double rtkz = tkx * transOp[0][2] + tky * transOp[1][2] + tkz * transOp[2][2];
2210 torque.add(threadID, k, rtkx, rtky, rtkz);
2211 }
2212
2213 private void polarizationEnergyGradient(int i, int k) {
2214
2215
2216 final double dpsymdx =
2217 -uxi * (sxk * gux[5] + syk * guy[5] + szk * guz[5])
2218 - uyi * (sxk * gux[6] + syk * guy[6] + szk * guz[6])
2219 - uzi * (sxk * gux[7] + syk * guy[7] + szk * guz[7])
2220 - uxk * (sxi * gux[5] + syi * guy[5] + szi * guz[5])
2221 - uyk * (sxi * gux[6] + syi * guy[6] + szi * guz[6])
2222 - uzk * (sxi * gux[7] + syi * guy[7] + szi * guz[7]);
2223 final double dpwidx =
2224 ci * (sxk * gc[5] + syk * gc[6] + szk * gc[7])
2225 - ck * (sxi * gux[2] + syi * guy[2] + szi * guz[2])
2226 - sxi
2227 * (qxxk * gux[11]
2228 + qyyk * gux[14]
2229 + qzzk * gux[16]
2230 + 2.0 * (qxyk * gux[12] + qxzk * gux[13] + qyzk * gux[15]))
2231 - syi
2232 * (qxxk * guy[11]
2233 + qyyk * guy[14]
2234 + qzzk * guy[16]
2235 + 2.0 * (qxyk * guy[12] + qxzk * guy[13] + qyzk * guy[15]))
2236 - szi
2237 * (qxxk * guz[11]
2238 + qyyk * guz[14]
2239 + qzzk * guz[16]
2240 + 2.0 * (qxyk * guz[12] + qxzk * guz[13] + qyzk * guz[15]))
2241 + sxk
2242 * (qxxi * gqxx[5]
2243 + qyyi * gqyy[5]
2244 + qzzi * gqzz[5]
2245 + 2.0 * (qxyi * gqxy[5] + qxzi * gqxz[5] + qyzi * gqyz[5]))
2246 + syk
2247 * (qxxi * gqxx[6]
2248 + qyyi * gqyy[6]
2249 + qzzi * gqzz[6]
2250 + 2.0 * (qxyi * gqxy[6] + qxzi * gqxz[6] + qyzi * gqyz[6]))
2251 + szk
2252 * (qxxi * gqxx[7]
2253 + qyyi * gqyy[7]
2254 + qzzi * gqzz[7]
2255 + 2.0 * (qxyi * gqxy[7] + qxzi * gqxz[7] + qyzi * gqyz[7]));
2256 final double dpwkdx =
2257 ci * (sxk * gux[2] + syk * guy[2] + szk * guz[2])
2258 - ck * (sxi * gc[5] + syi * gc[6] + szi * gc[7])
2259 - sxi
2260 * (qxxk * gqxx[5]
2261 + qyyk * gqyy[5]
2262 + qzzk * gqzz[5]
2263 + 2.0 * (qxyk * gqxy[5] + qxzk * gqxz[5] + qyzk * gqyz[5]))
2264 - syi
2265 * (qxxk * gqxx[6]
2266 + qyyk * gqyy[6]
2267 + qzzk * gqzz[6]
2268 + 2.0 * (qxyk * gqxy[6] + qxzk * gqxz[6] + qyzk * gqyz[6]))
2269 - szi
2270 * (qxxk * gqxx[7]
2271 + qyyk * gqyy[7]
2272 + qzzk * gqzz[7]
2273 + 2.0 * (qxyk * gqxy[7] + qxzk * gqxz[7] + qyzk * gqyz[7]))
2274 + sxk
2275 * (qxxi * gux[11]
2276 + qyyi * gux[14]
2277 + qzzi * gux[16]
2278 + 2.0 * (qxyi * gux[12] + qxzi * gux[13] + qyzi * gux[15]))
2279 + syk
2280 * (qxxi * guy[11]
2281 + qyyi * guy[14]
2282 + qzzi * guy[16]
2283 + 2.0 * (qxyi * guy[12] + qxzi * guy[13] + qyzi * guy[15]))
2284 + szk
2285 * (qxxi * guz[11]
2286 + qyyi * guz[14]
2287 + qzzi * guz[16]
2288 + 2.0 * (qxyi * guz[12] + qxzi * guz[13] + qyzi * guz[15]));
2289 final double dpsymdy =
2290 -uxi * (sxk * gux[6] + syk * guy[6] + szk * guz[6])
2291 - uyi * (sxk * gux[8] + syk * guy[8] + szk * guz[8])
2292 - uzi * (sxk * gux[9] + syk * guy[9] + szk * guz[9])
2293 - uxk * (sxi * gux[6] + syi * guy[6] + szi * guz[6])
2294 - uyk * (sxi * gux[8] + syi * guy[8] + szi * guz[8])
2295 - uzk * (sxi * gux[9] + syi * guy[9] + szi * guz[9]);
2296 final double dpwidy =
2297 ci * (sxk * gc[6] + syk * gc[8] + szk * gc[9])
2298 - ck * (sxi * gux[3] + syi * guy[3] + szi * guz[3])
2299 - sxi
2300 * (qxxk * gux[12]
2301 + qyyk * gux[17]
2302 + qzzk * gux[19]
2303 + 2.0 * (qxyk * gux[14] + qxzk * gux[15] + qyzk * gux[18]))
2304 - syi
2305 * (qxxk * guy[12]
2306 + qyyk * guy[17]
2307 + qzzk * guy[19]
2308 + 2.0 * (qxyk * guy[14] + qxzk * guy[15] + qyzk * guy[18]))
2309 - szi
2310 * (qxxk * guz[12]
2311 + qyyk * guz[17]
2312 + qzzk * guz[19]
2313 + 2.0 * (qxyk * guz[14] + qxzk * guz[15] + qyzk * guz[18]))
2314 + sxk
2315 * (qxxi * gqxx[6]
2316 + qyyi * gqyy[6]
2317 + qzzi * gqzz[6]
2318 + 2.0 * (qxyi * gqxy[6] + qxzi * gqxz[6] + qyzi * gqyz[6]))
2319 + syk
2320 * (qxxi * gqxx[8]
2321 + qyyi * gqyy[8]
2322 + qzzi * gqzz[8]
2323 + 2.0 * (qxyi * gqxy[8] + qxzi * gqxz[8] + qyzi * gqyz[8]))
2324 + szk
2325 * (qxxi * gqxx[9]
2326 + qyyi * gqyy[9]
2327 + qzzi * gqzz[9]
2328 + 2.0 * (qxyi * gqxy[9] + qxzi * gqxz[9] + qyzi * gqyz[9]));
2329 final double dpwkdy =
2330 ci * (sxk * gux[3] + syk * guy[3] + szk * guz[3])
2331 - ck * (sxi * gc[6] + syi * gc[8] + szi * gc[9])
2332 - sxi
2333 * (qxxk * gqxx[6]
2334 + qyyk * gqyy[6]
2335 + qzzk * gqzz[6]
2336 + 2.0 * (qxyk * gqxy[6] + qxzk * gqxz[6] + qyzk * gqyz[6]))
2337 - syi
2338 * (qxxk * gqxx[8]
2339 + qyyk * gqyy[8]
2340 + qzzk * gqzz[8]
2341 + 2.0 * (qxyk * gqxy[8] + qxzk * gqxz[8] + qyzk * gqyz[8]))
2342 - szi
2343 * (qxxk * gqxx[9]
2344 + qyyk * gqyy[9]
2345 + qzzk * gqzz[9]
2346 + 2.0 * (qxyk * gqxy[9] + qxzk * gqxz[9] + qyzk * gqyz[9]))
2347 + sxk
2348 * (qxxi * gux[12]
2349 + qyyi * gux[17]
2350 + qzzi * gux[19]
2351 + 2.0 * (qxyi * gux[14] + qxzi * gux[15] + qyzi * gux[18]))
2352 + syk
2353 * (qxxi * guy[12]
2354 + qyyi * guy[17]
2355 + qzzi * guy[19]
2356 + 2.0 * (qxyi * guy[14] + qxzi * guy[15] + qyzi * guy[18]))
2357 + szk
2358 * (qxxi * guz[12]
2359 + qyyi * guz[17]
2360 + qzzi * guz[19]
2361 + 2.0 * (qxyi * guz[14] + qxzi * guz[15] + qyzi * guz[18]));
2362 final double dpsymdz =
2363 -uxi * (sxk * gux[7] + syk * guy[7] + szk * guz[7])
2364 - uyi * (sxk * gux[9] + syk * guy[9] + szk * guz[9])
2365 - uzi * (sxk * gux[10] + syk * guy[10] + szk * guz[10])
2366 - uxk * (sxi * gux[7] + syi * guy[7] + szi * guz[7])
2367 - uyk * (sxi * gux[9] + syi * guy[9] + szi * guz[9])
2368 - uzk * (sxi * gux[10] + syi * guy[10] + szi * guz[10]);
2369 final double dpwidz =
2370 ci * (sxk * gc[7] + syk * gc[9] + szk * gc[10])
2371 - ck * (sxi * gux[4] + syi * guy[4] + szi * guz[4])
2372 - sxi
2373 * (qxxk * gux[13]
2374 + qyyk * gux[18]
2375 + qzzk * gux[20]
2376 + 2.0 * (qxyk * gux[15] + qxzk * gux[16] + qyzk * gux[19]))
2377 - syi
2378 * (qxxk * guy[13]
2379 + qyyk * guy[18]
2380 + qzzk * guy[20]
2381 + 2.0 * (qxyk * guy[15] + qxzk * guy[16] + qyzk * guy[19]))
2382 - szi
2383 * (qxxk * guz[13]
2384 + qyyk * guz[18]
2385 + qzzk * guz[20]
2386 + 2.0 * (qxyk * guz[15] + qxzk * guz[16] + qyzk * guz[19]))
2387 + sxk
2388 * (qxxi * gqxx[7]
2389 + qyyi * gqyy[7]
2390 + qzzi * gqzz[7]
2391 + 2.0 * (qxyi * gqxy[7] + qxzi * gqxz[7] + qyzi * gqyz[7]))
2392 + syk
2393 * (qxxi * gqxx[9]
2394 + qyyi * gqyy[9]
2395 + qzzi * gqzz[9]
2396 + 2.0 * (qxyi * gqxy[9] + qxzi * gqxz[9] + qyzi * gqyz[9]))
2397 + szk
2398 * (qxxi * gqxx[10]
2399 + qyyi * gqyy[10]
2400 + qzzi * gqzz[10]
2401 + 2.0 * (qxyi * gqxy[10] + qxzi * gqxz[10] + qyzi * gqyz[10]));
2402 final double dpwkdz =
2403 ci * (sxk * gux[4] + syk * guy[4] + szk * guz[4])
2404 - ck * (sxi * gc[7] + syi * gc[9] + szi * gc[10])
2405 - sxi
2406 * (qxxk * gqxx[7]
2407 + qyyk * gqyy[7]
2408 + qzzk * gqzz[7]
2409 + 2.0 * (qxyk * gqxy[7] + qxzk * gqxz[7] + qyzk * gqyz[7]))
2410 - syi
2411 * (qxxk * gqxx[9]
2412 + qyyk * gqyy[9]
2413 + qzzk * gqzz[9]
2414 + 2.0 * (qxyk * gqxy[9] + qxzk * gqxz[9] + qyzk * gqyz[9]))
2415 - szi
2416 * (qxxk * gqxx[10]
2417 + qyyk * gqyy[10]
2418 + qzzk * gqzz[10]
2419 + 2.0 * (qxyk * gqxy[10] + qxzk * gqxz[10] + qyzk * gqyz[10]))
2420 + sxk
2421 * (qxxi * gux[13]
2422 + qyyi * gux[18]
2423 + qzzi * gux[20]
2424 + 2.0 * (qxyi * gux[15] + qxzi * gux[16] + qyzi * gux[19]))
2425 + syk
2426 * (qxxi * guy[13]
2427 + qyyi * guy[18]
2428 + qzzi * guy[20]
2429 + 2.0 * (qxyi * guy[15] + qxzi * guy[16] + qyzi * guy[19]))
2430 + szk
2431 * (qxxi * guz[13]
2432 + qyyi * guz[18]
2433 + qzzi * guz[20]
2434 + 2.0 * (qxyi * guz[15] + qxzi * guz[16] + qyzi * guz[19]));
2435
2436
2437
2438 final double dsymdr =
2439 -uxi * (sxk * gux[22] + syk * guy[22] + szk * guz[22])
2440 - uyi * (sxk * gux[23] + syk * guy[23] + szk * guz[23])
2441 - uzi * (sxk * gux[24] + syk * guy[24] + szk * guz[24])
2442 - uxk * (sxi * gux[22] + syi * guy[22] + szi * guz[22])
2443 - uyk * (sxi * gux[23] + syi * guy[23] + szi * guz[23])
2444 - uzk * (sxi * gux[24] + syi * guy[24] + szi * guz[24]);
2445 final double dwipdr =
2446 ci * (sxk * gc[22] + syk * gc[23] + szk * gc[24])
2447 - ck * (sxi * gux[21] + syi * guy[21] + szi * guz[21])
2448 - sxi
2449 * (qxxk * gux[25]
2450 + qyyk * gux[28]
2451 + qzzk * gux[30]
2452 + 2.0 * (qxyk * gux[26] + qxzk * gux[27] + qyzk * gux[29]))
2453 - syi
2454 * (qxxk * guy[25]
2455 + qyyk * guy[28]
2456 + qzzk * guy[30]
2457 + 2.0 * (qxyk * guy[26] + qxzk * guy[27] + qyzk * guy[29]))
2458 - szi
2459 * (qxxk * guz[25]
2460 + qyyk * guz[28]
2461 + qzzk * guz[30]
2462 + 2.0 * (qxyk * guz[26] + qxzk * guz[27] + qyzk * guz[29]))
2463 + sxk
2464 * (qxxi * gqxx[22]
2465 + qyyi * gqyy[22]
2466 + qzzi * gqzz[22]
2467 + 2.0 * (qxyi * gqxy[22] + qxzi * gqxz[22] + qyzi * gqyz[22]))
2468 + syk
2469 * (qxxi * gqxx[23]
2470 + qyyi * gqyy[23]
2471 + qzzi * gqzz[23]
2472 + 2.0 * (qxyi * gqxy[23] + qxzi * gqxz[23] + qyzi * gqyz[23]))
2473 + szk
2474 * (qxxi * gqxx[24]
2475 + qyyi * gqyy[24]
2476 + qzzi * gqzz[24]
2477 + 2.0 * (qxyi * gqxy[24] + qxzi * gqxz[24] + qyzi * gqyz[24]));
2478 final double dwkpdr =
2479 ci * (sxk * gux[21] + syk * guy[21] + szk * guz[21])
2480 - ck * (sxi * gc[22] + syi * gc[23] + szi * gc[24])
2481 - sxi
2482 * (qxxk * gqxx[22]
2483 + qyyk * gqyy[22]
2484 + qzzk * gqzz[22]
2485 + 2.0 * (qxyk * gqxy[22] + qxzk * gqxz[22] + qyzk * gqyz[22]))
2486 - syi
2487 * (qxxk * gqxx[23]
2488 + qyyk * gqyy[23]
2489 + qzzk * gqzz[23]
2490 + 2.0 * (qxyk * gqxy[23] + qxzk * gqxz[23] + qyzk * gqyz[23]))
2491 - szi
2492 * (qxxk * gqxx[24]
2493 + qyyk * gqyy[24]
2494 + qzzk * gqzz[24]
2495 + 2.0 * (qxyk * gqxy[24] + qxzk * gqxz[24] + qyzk * gqyz[24]))
2496 + sxk
2497 * (qxxi * gux[25]
2498 + qyyi * gux[28]
2499 + qzzi * gux[30]
2500 + 2.0 * (qxyi * gux[26] + qxzi * gux[27] + qyzi * gux[29]))
2501 + syk
2502 * (qxxi * guy[25]
2503 + qyyi * guy[28]
2504 + qzzi * guy[30]
2505 + 2.0 * (qxyi * guy[26] + qxzi * guy[27] + qyzi * guy[29]))
2506 + szk
2507 * (qxxi * guz[25]
2508 + qyyi * guz[28]
2509 + qzzi * guz[30]
2510 + 2.0 * (qxyi * guz[26] + qxzi * guz[27] + qyzi * guz[29]));
2511 double dpdx = 0.5 * (dpsymdx + 0.5 * (dpwidx + dpwkdx));
2512 double dpdy = 0.5 * (dpsymdy + 0.5 * (dpwidy + dpwkdy));
2513 double dpdz = 0.5 * (dpsymdz + 0.5 * (dpwidz + dpwkdz));
2514 double dsumdri = dsymdr + 0.5 * (dwipdr + dwkpdr);
2515 double dbi = 0.5 * rbk * dsumdri;
2516 double dbk = 0.5 * rbi * dsumdri;
2517 if (polarization == Polarization.MUTUAL) {
2518 dpdx -=
2519 0.5
2520 * (dxi * (pxk * gux[5] + pyk * gux[6] + pzk * gux[7])
2521 + dyi * (pxk * guy[5] + pyk * guy[6] + pzk * guy[7])
2522 + dzi * (pxk * guz[5] + pyk * guz[6] + pzk * guz[7])
2523 + dxk * (pxi * gux[5] + pyi * gux[6] + pzi * gux[7])
2524 + dyk * (pxi * guy[5] + pyi * guy[6] + pzi * guy[7])
2525 + dzk * (pxi * guz[5] + pyi * guz[6] + pzi * guz[7]));
2526 dpdy -=
2527 0.5
2528 * (dxi * (pxk * gux[6] + pyk * gux[8] + pzk * gux[9])
2529 + dyi * (pxk * guy[6] + pyk * guy[8] + pzk * guy[9])
2530 + dzi * (pxk * guz[6] + pyk * guz[8] + pzk * guz[9])
2531 + dxk * (pxi * gux[6] + pyi * gux[8] + pzi * gux[9])
2532 + dyk * (pxi * guy[6] + pyi * guy[8] + pzi * guy[9])
2533 + dzk * (pxi * guz[6] + pyi * guz[8] + pzi * guz[9]));
2534 dpdz -=
2535 0.5
2536 * (dxi * (pxk * gux[7] + pyk * gux[9] + pzk * gux[10])
2537 + dyi * (pxk * guy[7] + pyk * guy[9] + pzk * guy[10])
2538 + dzi * (pxk * guz[7] + pyk * guz[9] + pzk * guz[10])
2539 + dxk * (pxi * gux[7] + pyi * gux[9] + pzi * gux[10])
2540 + dyk * (pxi * guy[7] + pyi * guy[9] + pzi * guy[10])
2541 + dzk * (pxi * guz[7] + pyi * guz[9] + pzi * guz[10]));
2542 final double duvdr =
2543 dxi * (pxk * gux[22] + pyk * gux[23] + pzk * gux[24])
2544 + dyi * (pxk * guy[22] + pyk * guy[23] + pzk * guy[24])
2545 + dzi * (pxk * guz[22] + pyk * guz[23] + pzk * guz[24])
2546 + dxk * (pxi * gux[22] + pyi * gux[23] + pzi * gux[24])
2547 + dyk * (pxi * guy[22] + pyi * guy[23] + pzi * guy[24])
2548 + dzk * (pxi * guz[22] + pyi * guz[23] + pzi * guz[24]);
2549 dbi -= 0.5 * rbk * duvdr;
2550 dbk -= 0.5 * rbi * duvdr;
2551 }
2552
2553
2554 if (i == k && iSymm == 0) {
2555 dborni += dbi;
2556 } else {
2557 if (i == k) {
2558 dpdx *= 0.5;
2559 dpdy *= 0.5;
2560 dpdz *= 0.5;
2561 dbi *= 0.5;
2562 dbk *= 0.5;
2563 }
2564 dedxi -= dpdx;
2565 dedyi -= dpdy;
2566 dedzi -= dpdz;
2567 dborni += dbi;
2568 final double rdpdx = dpdx * transOp[0][0] + dpdy * transOp[1][0] + dpdz * transOp[2][0];
2569 final double rdpdy = dpdx * transOp[0][1] + dpdy * transOp[1][1] + dpdz * transOp[2][1];
2570 final double rdpdz = dpdx * transOp[0][2] + dpdy * transOp[1][2] + dpdz * transOp[2][2];
2571 grad.add(threadID, k, rdpdx, rdpdy, rdpdz);
2572 sharedBornGrad.add(threadID, k, dbk);
2573 }
2574 polarizationEnergyTorque(i, k);
2575 }
2576
2577 private void polarizationEnergyTorque(int i, int k) {
2578 double fix = 0.5 * (sxk * gux[2] + syk * guy[2] + szk * guz[2]);
2579 double fiy = 0.5 * (sxk * gux[3] + syk * guy[3] + szk * guz[3]);
2580 double fiz = 0.5 * (sxk * gux[4] + syk * guy[4] + szk * guz[4]);
2581 double fkx = 0.5 * (sxi * gux[2] + syi * guy[2] + szi * guz[2]);
2582 double fky = 0.5 * (sxi * gux[3] + syi * guy[3] + szi * guz[3]);
2583 double fkz = 0.5 * (sxi * gux[4] + syi * guy[4] + szi * guz[4]);
2584 double fixx =
2585 -0.25
2586 * ((sxk * gqxx[2] + syk * gqxx[3] + szk * gqxx[4])
2587 + (sxk * gux[5] + syk * guy[5] + szk * guz[5]));
2588 double fixy =
2589 -0.25
2590 * ((sxk * gqxy[2] + syk * gqxy[3] + szk * gqxy[4])
2591 + (sxk * gux[6] + syk * guy[6] + szk * guz[6]));
2592 double fixz =
2593 -0.25
2594 * ((sxk * gqxz[2] + syk * gqxz[3] + szk * gqxz[4])
2595 + (sxk * gux[7] + syk * guy[7] + szk * guz[7]));
2596 double fiyy =
2597 -0.25
2598 * ((sxk * gqyy[2] + syk * gqyy[3] + szk * gqyy[4])
2599 + (sxk * gux[8] + syk * guy[8] + szk * guz[8]));
2600 double fiyz =
2601 -0.25
2602 * ((sxk * gqyz[2] + syk * gqyz[3] + szk * gqyz[4])
2603 + (sxk * gux[9] + syk * guy[9] + szk * guz[9]));
2604 double fizz =
2605 -0.25
2606 * ((sxk * gqzz[2] + syk * gqzz[3] + szk * gqzz[4])
2607 + (sxk * gux[10] + syk * guy[10] + szk * guz[10]));
2608 double fiyx = fixy;
2609 double fizx = fixz;
2610 double fizy = fiyz;
2611 double fkxx =
2612 0.25
2613 * ((sxi * gqxx[2] + syi * gqxx[3] + szi * gqxx[4])
2614 + (sxi * gux[5] + syi * guy[5] + szi * guz[5]));
2615 double fkxy =
2616 0.25
2617 * ((sxi * gqxy[2] + syi * gqxy[3] + szi * gqxy[4])
2618 + (sxi * gux[6] + syi * guy[6] + szi * guz[6]));
2619 double fkxz =
2620 0.25
2621 * ((sxi * gqxz[2] + syi * gqxz[3] + szi * gqxz[4])
2622 + (sxi * gux[7] + syi * guy[7] + szi * guz[7]));
2623 double fkyy =
2624 0.25
2625 * ((sxi * gqyy[2] + syi * gqyy[3] + szi * gqyy[4])
2626 + (sxi * gux[8] + syi * guy[8] + szi * guz[8]));
2627 double fkyz =
2628 0.25
2629 * ((sxi * gqyz[2] + syi * gqyz[3] + szi * gqyz[4])
2630 + (sxi * gux[9] + syi * guy[9] + szi * guz[9]));
2631 double fkzz =
2632 0.25
2633 * ((sxi * gqzz[2] + syi * gqzz[3] + szi * gqzz[4])
2634 + (sxi * gux[10] + syi * guy[10] + szi * guz[10]));
2635 double fkyx = fkxy;
2636 double fkzx = fkxz;
2637 double fkzy = fkyz;
2638 if (i == k) {
2639 fix *= 0.5;
2640 fiy *= 0.5;
2641 fiz *= 0.5;
2642 fkx *= 0.5;
2643 fky *= 0.5;
2644 fkz *= 0.5;
2645 fixx *= 0.5;
2646 fixy *= 0.5;
2647 fixz *= 0.5;
2648 fiyx *= 0.5;
2649 fiyy *= 0.5;
2650 fiyz *= 0.5;
2651 fizx *= 0.5;
2652 fizy *= 0.5;
2653 fizz *= 0.5;
2654 fkxx *= 0.5;
2655 fkxy *= 0.5;
2656 fkxz *= 0.5;
2657 fkyx *= 0.5;
2658 fkyy *= 0.5;
2659 fkyz *= 0.5;
2660 fkzx *= 0.5;
2661 fkzy *= 0.5;
2662 fkzz *= 0.5;
2663 }
2664
2665
2666 double tix = uyi * fiz - uzi * fiy;
2667 double tiy = uzi * fix - uxi * fiz;
2668 double tiz = uxi * fiy - uyi * fix;
2669 double tkx = uyk * fkz - uzk * fky;
2670 double tky = uzk * fkx - uxk * fkz;
2671 double tkz = uxk * fky - uyk * fkx;
2672
2673
2674 tix +=
2675 2.0 * (qxyi * fixz + qyyi * fiyz + qyzi * fizz - qxzi * fixy - qyzi * fiyy - qzzi * fizy);
2676 tiy +=
2677 2.0 * (qxzi * fixx + qyzi * fiyx + qzzi * fizx - qxxi * fixz - qxyi * fiyz - qxzi * fizz);
2678 tiz +=
2679 2.0 * (qxxi * fixy + qxyi * fiyy + qxzi * fizy - qxyi * fixx - qyyi * fiyx - qyzi * fizx);
2680 tkx +=
2681 2.0 * (qxyk * fkxz + qyyk * fkyz + qyzk * fkzz - qxzk * fkxy - qyzk * fkyy - qzzk * fkzy);
2682 tky +=
2683 2.0 * (qxzk * fkxx + qyzk * fkyx + qzzk * fkzx - qxxk * fkxz - qxyk * fkyz - qxzk * fkzz);
2684 tkz +=
2685 2.0 * (qxxk * fkxy + qxyk * fkyy + qxzk * fkzy - qxyk * fkxx - qyyk * fkyx - qyzk * fkzx);
2686 trqxi += tix;
2687 trqyi += tiy;
2688 trqzi += tiz;
2689
2690 final double rx = tkx;
2691 final double ry = tky;
2692 final double rz = tkz;
2693 tkx = rx * transOp[0][0] + ry * transOp[1][0] + rz * transOp[2][0];
2694 tky = rx * transOp[0][1] + ry * transOp[1][1] + rz * transOp[2][1];
2695 tkz = rx * transOp[0][2] + ry * transOp[1][2] + rz * transOp[2][2];
2696 torque.add(threadID, k, tkx, tky, tkz);
2697 }
2698 }
2699
2700
2701
2702
2703
2704
2705 private class GKEnergyLoopQI extends IntegerForLoop {
2706
2707 private final double[] dx_local;
2708 private final double[] gradI = new double[3];
2709 private final double[] torqueI = new double[3];
2710 private final double[] torqueK = new double[3];
2711 private final double[] gI = new double[3];
2712 private final double[] tI = new double[3];
2713 private final double[] tK = new double[3];
2714
2715 private double rbi;
2716 PolarizableMultipole mI;
2717 PolarizableMultipole mK;
2718 QIFrame qiFrame;
2719 GKEnergyQI gkEnergyQI;
2720 private double xi, yi, zi;
2721 private double dedxi, dedyi, dedzi;
2722 private double dborni;
2723 private double trqxi, trqyi, trqzi;
2724 private int count;
2725 private int iSymm;
2726 private int threadID;
2727 private final double[][] transOp;
2728 private double gkEnergy;
2729 private double gkPermanentEnergy;
2730 private double gkPolarizationEnergy;
2731
2732 GKEnergyLoopQI() {
2733 mI = new PolarizableMultipole();
2734 mK = new PolarizableMultipole();
2735 qiFrame = new QIFrame();
2736 dx_local = new double[3];
2737 transOp = new double[3][3];
2738 }
2739
2740 @Override
2741 public void start() {
2742 gkEnergy = 0.0;
2743 gkPermanentEnergy = 0.0;
2744 gkPolarizationEnergy = 0.0;
2745 count = 0;
2746 threadID = getThreadIndex();
2747 gkEnergyQI = new GKEnergyQI(soluteDielectric, solventDielectric, gkc, gradient);
2748 }
2749
2750 @Override
2751 public void finish() {
2752 sharedInteractions.addAndGet(count);
2753 sharedGKEnergy.addAndGet(gkEnergy);
2754 sharedPermanentGKEnergy.addAndGet(gkPermanentEnergy);
2755 sharedPolarizationGKEnergy.addAndGet(gkPolarizationEnergy);
2756 }
2757
2758 @Override
2759 public void run(int lb, int ub) {
2760
2761 double[] x = sXYZ[0][0];
2762 double[] y = sXYZ[0][1];
2763 double[] z = sXYZ[0][2];
2764
2765 int nSymm = crystal.spaceGroup.symOps.size();
2766 List<SymOp> symOps = crystal.spaceGroup.symOps;
2767
2768 for (iSymm = 0; iSymm < nSymm; iSymm++) {
2769 SymOp symOp = symOps.get(iSymm);
2770 crystal.getTransformationOperator(symOp, transOp);
2771 for (int i = lb; i <= ub; i++) {
2772 if (!use[i]) {
2773 continue;
2774 }
2775
2776
2777 dedxi = 0.0;
2778 dedyi = 0.0;
2779 dedzi = 0.0;
2780 trqxi = 0.0;
2781 trqyi = 0.0;
2782 trqzi = 0.0;
2783 dborni = 0.0;
2784
2785
2786 xi = x[i];
2787 yi = y[i];
2788 zi = z[i];
2789 rbi = born[i];
2790 int[] list = neighborLists[iSymm][i];
2791 for (int k : list) {
2792 if (!use[k]) {
2793 continue;
2794 }
2795 interaction(i, k);
2796 }
2797 if (iSymm == 0) {
2798
2799 interaction(i, i);
2800
2801
2802
2803
2804
2805
2806 switch (nonPolarModel) {
2807 case BORN_SOLV:
2808 case BORN_CAV_DISP:
2809 double r = baseRadius[i] + dOffset + probe;
2810 double ratio = (baseRadius[i] + dOffset) / born[i];
2811 ratio *= ratio;
2812 ratio *= (ratio * ratio);
2813 double saTerm = surfaceTension * r * r * ratio / 6.0;
2814 gkEnergy += saTerm;
2815 sharedBornGrad.sub(threadID, i, 6.0 * saTerm / born[i]);
2816 break;
2817 default:
2818 break;
2819 }
2820 }
2821 if (gradient) {
2822 grad.add(threadID, i, dedxi, dedyi, dedzi);
2823 torque.add(threadID, i, trqxi, trqyi, trqzi);
2824 sharedBornGrad.add(threadID, i, dborni);
2825 }
2826 }
2827 }
2828 }
2829
2830 private void interaction(int i, int k) {
2831 dx_local[0] = sXYZ[iSymm][0][k] - xi;
2832 dx_local[1] = sXYZ[iSymm][1][k] - yi;
2833 dx_local[2] = sXYZ[iSymm][2][k] - zi;
2834 double r2 = crystal.image(dx_local);
2835 if (r2 > cut2) {
2836 return;
2837 }
2838
2839
2840 mI.setPermanentMultipole(globalMultipole[0][i]);
2841 mI.setInducedDipole(inducedDipole[0][i], inducedDipoleCR[0][i]);
2842
2843 mK.setPermanentMultipole(globalMultipole[iSymm][k]);
2844 mK.setInducedDipole(inducedDipole[iSymm][k], inducedDipoleCR[iSymm][k]);
2845
2846 double selfScale = 1.0;
2847 if (i == k) {
2848 selfScale = 0.5;
2849 } else {
2850 qiFrame.setAndRotate(dx_local, mI, mK);
2851 }
2852
2853
2854 double rbk = born[k];
2855 gkEnergyQI.initPotential(dx_local, r2, rbi, rbk);
2856
2857
2858 double eik;
2859 if (!gradient) {
2860
2861 eik = electric * selfScale * gkEnergyQI.multipoleEnergy(mI, mK);
2862 gkPermanentEnergy += eik;
2863 if (polarization != Polarization.NONE) {
2864
2865 double ep = electric * selfScale * gkEnergyQI.polarizationEnergy(mI, mK);
2866 gkPolarizationEnergy += ep;
2867 eik += ep;
2868 }
2869 } else {
2870 if (i == k) {
2871
2872 eik = electric * selfScale * gkEnergyQI.multipoleEnergy(mI, mK);
2873 gkPermanentEnergy += eik;
2874
2875 if (polarization != Polarization.NONE) {
2876
2877
2878
2879 double mutualMask = 1.0;
2880 if (polarization == Polarization.DIRECT) {
2881 mutualMask = 0.0;
2882 }
2883 double ep =
2884 electric * selfScale * gkEnergyQI.polarizationEnergyAndGradient(mI, mK, mutualMask,
2885 gradI, torqueI, torqueK);
2886 gkPolarizationEnergy += ep;
2887 eik += ep;
2888
2889 trqxi += electric * selfScale * torqueI[0];
2890 trqyi += electric * selfScale * torqueI[1];
2891 trqzi += electric * selfScale * torqueI[2];
2892 double tkx = electric * selfScale * torqueK[0];
2893 double tky = electric * selfScale * torqueK[1];
2894 double tkz = electric * selfScale * torqueK[2];
2895 final double rtkx = tkx * transOp[0][0] + tky * transOp[1][0] + tkz * transOp[2][0];
2896 final double rtky = tkx * transOp[0][1] + tky * transOp[1][1] + tkz * transOp[2][1];
2897 final double rtkz = tkx * transOp[0][2] + tky * transOp[1][2] + tkz * transOp[2][2];
2898 torque.add(threadID, k, rtkx, rtky, rtkz);
2899 }
2900
2901 gkEnergyQI.initBorn(dx_local, r2, rbi, rbk);
2902 double db = gkEnergyQI.multipoleEnergyBornGrad(mI, mK);
2903 if (polarization != Polarization.NONE) {
2904 db += gkEnergyQI.polarizationEnergyBornGrad(mI, mK, polarization == Polarization.MUTUAL);
2905 }
2906 dborni += electric * rbi * db;
2907 } else {
2908
2909 eik = electric * gkEnergyQI.multipoleEnergyAndGradient(mI, mK, gradI, torqueI, torqueK);
2910 gkPermanentEnergy += eik;
2911 if (polarization != Polarization.NONE) {
2912 double mutualMask = 1.0;
2913 if (polarization == Polarization.DIRECT) {
2914 mutualMask = 0.0;
2915 }
2916
2917 double ep =
2918 electric * gkEnergyQI.polarizationEnergyAndGradient(mI, mK, mutualMask, gI, tI, tK);
2919 eik += ep;
2920 gkPolarizationEnergy += ep;
2921 for (int j = 0; j < 3; j++) {
2922 gradI[j] += gI[j];
2923 torqueI[j] += tI[j];
2924 torqueK[j] += tK[j];
2925 }
2926 }
2927
2928
2929 for (int j = 0; j < 3; j++) {
2930 gradI[j] *= electric;
2931 torqueI[j] *= electric;
2932 torqueK[j] *= electric;
2933 }
2934
2935
2936 qiFrame.toGlobal(gradI);
2937 qiFrame.toGlobal(torqueI);
2938 qiFrame.toGlobal(torqueK);
2939
2940
2941 dedxi += gradI[0];
2942 dedyi += gradI[1];
2943 dedzi += gradI[2];
2944 trqxi += torqueI[0];
2945 trqyi += torqueI[1];
2946 trqzi += torqueI[2];
2947
2948
2949 double dedx = -gradI[0];
2950 double dedy = -gradI[1];
2951 double dedz = -gradI[2];
2952
2953 final double dedxk = dedx * transOp[0][0] + dedy * transOp[1][0] + dedz * transOp[2][0];
2954 final double dedyk = dedx * transOp[0][1] + dedy * transOp[1][1] + dedz * transOp[2][1];
2955 final double dedzk = dedx * transOp[0][2] + dedy * transOp[1][2] + dedz * transOp[2][2];
2956 double tkx = torqueK[0];
2957 double tky = torqueK[1];
2958 double tkz = torqueK[2];
2959 final double rtkx = tkx * transOp[0][0] + tky * transOp[1][0] + tkz * transOp[2][0];
2960 final double rtky = tkx * transOp[0][1] + tky * transOp[1][1] + tkz * transOp[2][1];
2961 final double rtkz = tkx * transOp[0][2] + tky * transOp[1][2] + tkz * transOp[2][2];
2962 grad.add(threadID, k, dedxk, dedyk, dedzk);
2963 torque.add(threadID, k, rtkx, rtky, rtkz);
2964
2965
2966 gkEnergyQI.initBorn(dx_local, r2, rbi, rbk);
2967 double db = gkEnergyQI.multipoleEnergyBornGrad(mI, mK);
2968 if (polarization != Polarization.NONE) {
2969 db += gkEnergyQI.polarizationEnergyBornGrad(mI, mK, polarization == Polarization.MUTUAL);
2970 }
2971 db *= electric;
2972 dborni += rbk * db;
2973 sharedBornGrad.add(threadID, k, rbi * db);
2974 }
2975 }
2976
2977 if (i == k) {
2978 selfEnergy.add(threadID, i, eik);
2979 } else {
2980 double half = 0.5 * eik;
2981 crossEnergy.add(threadID, i, half);
2982 crossEnergy.add(threadID, k, half);
2983 }
2984
2985 gkEnergy += eik;
2986 count++;
2987 }
2988 }
2989 }