View Javadoc
1   // ******************************************************************************
2   //
3   // Title:       Force Field X.
4   // Description: Force Field X - Software for Molecular Biophysics.
5   // Copyright:   Copyright (c) Michael J. Schnieders 2001-2025.
6   //
7   // This file is part of Force Field X.
8   //
9   // Force Field X is free software; you can redistribute it and/or modify it
10  // under the terms of the GNU General Public License version 3 as published by
11  // the Free Software Foundation.
12  //
13  // Force Field X is distributed in the hope that it will be useful, but WITHOUT
14  // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
15  // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
16  // details.
17  //
18  // You should have received a copy of the GNU General Public License along with
19  // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
20  // Place, Suite 330, Boston, MA 02111-1307 USA
21  //
22  // Linking this library statically or dynamically with other modules is making a
23  // combined work based on this library. Thus, the terms and conditions of the
24  // GNU General Public License cover the whole combination.
25  //
26  // As a special exception, the copyright holders of this library give you
27  // permission to link this library with independent modules to produce an
28  // executable, regardless of the license terms of these independent modules, and
29  // to copy and distribute the resulting executable under terms of your choice,
30  // provided that you also meet, for each linked independent module, the terms
31  // and conditions of the license of that module. An independent module is a
32  // module which is not derived from or based on this library. If you modify this
33  // library, you may extend this exception to your version of the library, but
34  // you are not obligated to do so. If you do not wish to do so, delete this
35  // exception statement from your version.
36  //
37  // ******************************************************************************
38  package ffx.potential.nonbonded.pme;
39  
40  import edu.rit.pj.IntegerForLoop;
41  import edu.rit.pj.IntegerSchedule;
42  import edu.rit.pj.ParallelRegion;
43  import edu.rit.pj.ParallelTeam;
44  import edu.rit.pj.reduction.SharedDouble;
45  import edu.rit.pj.reduction.SharedInteger;
46  import ffx.crystal.Crystal;
47  import ffx.crystal.SymOp;
48  import ffx.numerics.atomic.AtomicDoubleArray3D;
49  import ffx.potential.bonded.Atom;
50  import ffx.potential.extended.ExtendedSystem;
51  import ffx.potential.nonbonded.MaskingInterface;
52  import ffx.potential.parameters.ForceField;
53  import ffx.potential.parameters.MultipoleType.MultipoleFrameDefinition;
54  import ffx.potential.utils.EnergyException;
55  
56  import java.util.List;
57  import java.util.logging.Level;
58  import java.util.logging.Logger;
59  
60  import static ffx.numerics.special.Erf.erfc;
61  import static ffx.potential.parameters.MultipoleType.t000;
62  import static ffx.potential.parameters.MultipoleType.t001;
63  import static ffx.potential.parameters.MultipoleType.t002;
64  import static ffx.potential.parameters.MultipoleType.t010;
65  import static ffx.potential.parameters.MultipoleType.t011;
66  import static ffx.potential.parameters.MultipoleType.t020;
67  import static ffx.potential.parameters.MultipoleType.t100;
68  import static ffx.potential.parameters.MultipoleType.t101;
69  import static ffx.potential.parameters.MultipoleType.t110;
70  import static ffx.potential.parameters.MultipoleType.t200;
71  import static java.lang.Double.isInfinite;
72  import static java.lang.Double.isNaN;
73  import static java.lang.String.format;
74  import static java.util.Arrays.fill;
75  import static org.apache.commons.math3.util.FastMath.exp;
76  import static org.apache.commons.math3.util.FastMath.min;
77  import static org.apache.commons.math3.util.FastMath.sqrt;
78  
79  /**
80   * Parallel evaluation of the PME real space energy and gradient.
81   *
82   * @author Michael J. Schnieders
83   * @since 1.0
84   */
85  public class RealSpaceEnergyRegion extends ParallelRegion implements MaskingInterface {
86    private static final Logger logger = Logger.getLogger(RealSpaceEnergyRegion.class.getName());
87    /**
88     * Constant applied to multipole interactions.
89     */
90    private static final double oneThird = 1.0 / 3.0;
91    private final int maxThreads;
92    private final double electric;
93    private final SharedInteger sharedInteractions;
94    private final ForceField.ELEC_FORM elecForm;
95    private final RealSpaceEnergyLoop[] realSpaceEnergyLoop;
96    private final FixedChargeEnergyLoop[] fixedChargeEnergyLoop;
97    /**
98     * Specify inter-molecular softcore.
99     */
100   private final boolean intermolecularSoftcore;
101   /**
102    * Specify intra-molecular softcore.
103    */
104   private final boolean intramolecularSoftcore;
105   /**
106    * Dimensions of [nsymm][nAtoms][3]
107    */
108   public double[][][] inducedDipole;
109   public double[][][] inducedDipoleCR;
110   /**
111    * Polarization groups.
112    */
113   protected int[][] ip11;
114   /**
115    * An ordered array of atoms in the system.
116    */
117   private Atom[] atoms;
118   /**
119    * Unit cell and spacegroup information.
120    */
121   private Crystal crystal;
122   /**
123    * Dimensions of [nsymm][xyz][nAtoms].
124    */
125   private double[][][] coordinates;
126   /**
127    * Multipole frame definition.
128    */
129   private MultipoleFrameDefinition[] frame;
130   /**
131    * Multipole frame defining atoms.
132    */
133   private int[][] axisAtom;
134   /**
135    * Dimensions of [nsymm][nAtoms][10]
136    */
137   private double[][][] globalMultipole;
138   private double[][][] titrationMultipole;
139   private double[][][] tautomerMultipole;
140 
141   private ExtendedSystem extendedSystem = null;
142   private boolean esvTerm = false;
143   /**
144    * When computing the polarization energy at Lambda there are 3 pieces.
145    *
146    * <p>1.) Upol(1) = The polarization energy computed normally (ie. system with ligand).
147    *
148    * <p>2.) Uenv = The polarization energy of the system without the ligand.
149    *
150    * <p>3.) Uligand = The polarization energy of the ligand by itself.
151    *
152    * <p>Upol(L) = L*Upol(1) + (1-L)*(Uenv + Uligand)
153    *
154    * <p>Set the "use" array to true for all atoms for part 1. Set the "use" array to true for all
155    * atoms except the ligand for part 2. Set the "use" array to true only for the ligand atoms for
156    * part 3.
157    *
158    * <p>The "use" array can also be employed to turn off atoms for computing the electrostatic
159    * energy of sub-structures.
160    */
161   private boolean[] use;
162   /**
163    * Molecule number for each atom.
164    */
165   private int[] molecule;
166   /**
167    * Flag for ligand atoms.
168    */
169   private boolean[] isSoft;
170   private double[] ipdamp;
171   private double[] thole;
172   /**
173    * Masking of 1-2, 1-3, 1-4 and 1-5 interactions.
174    */
175   private int[][] mask12;
176   private int[][] mask13;
177   private int[][] mask14;
178   private int[][] mask15;
179   /**
180    * Neighbor lists, without atoms beyond the real space cutoff. [nSymm][nAtoms][nIncludedNeighbors]
181    */
182   private int[][][] realSpaceLists;
183   /**
184    * Number of neighboring atoms within the real space cutoff. [nSymm][nAtoms]
185    */
186   private int[][] realSpaceCounts;
187   /**
188    * Pairwise schedule for load balancing.
189    */
190   private IntegerSchedule realSpaceSchedule;
191   private long[] realSpaceEnergyTime;
192   /**
193    * Atomic Gradient array.
194    */
195   private AtomicDoubleArray3D grad;
196   /**
197    * Atomic Torque array.
198    */
199   private AtomicDoubleArray3D torque;
200   /**
201    * Partial derivative of the gradient with respect to Lambda.
202    */
203   private AtomicDoubleArray3D lambdaGrad;
204   /**
205    * Partial derivative of the torque with respect to Lambda.
206    */
207   private AtomicDoubleArray3D lambdaTorque;
208   /**
209    * Partial derivative with respect to Lambda.
210    */
211   private SharedDouble shareddEdLambda;
212   /**
213    * Second partial derivative with respect to Lambda.
214    */
215   private SharedDouble sharedd2EdLambda2;
216   /**
217    * If true, compute coordinate gradient.
218    */
219   private boolean gradient;
220   /**
221    * If lambdaTerm is true, some ligand atom interactions with the environment are being turned
222    * on/off.
223    */
224   private boolean lambdaTerm;
225   /**
226    * If nnTerm is true, all intramolecular interactions are treated with a neural network.
227    */
228   private boolean nnTerm;
229 
230   /**
231    * The current LambdaMode of this PME instance (or OFF for no lambda dependence).
232    */
233   private LambdaMode lambdaMode = LambdaMode.OFF;
234   private Polarization polarization;
235   /**
236    * lAlpha = α*(1 - L)^2
237    */
238   private double lAlpha = 0.0;
239   private double dlAlpha = 0.0;
240   private double d2lAlpha = 0.0;
241   private double dEdLSign = 1.0;
242   /**
243    * lPowPerm = L^permanentLambdaExponent
244    */
245   private double dlPowPerm = 0.0;
246   private double d2lPowPerm = 0.0;
247   private boolean doPermanentRealSpace;
248   private double permanentScale = 1.0;
249   /**
250    * lPowPol = L^polarizationLambdaExponent
251    */
252   private double lPowPol = 1.0;
253   private double dlPowPol = 0.0;
254   private double d2lPowPol = 0.0;
255   private boolean doPolarization;
256   /**
257    * When computing the polarization energy at L there are 3 pieces.
258    *
259    * <p>1.) Upol(1) = The polarization energy computed normally (ie. system with ligand).
260    *
261    * <p>2.) Uenv = The polarization energy of the system without the ligand.
262    *
263    * <p>3.) Uligand = The polarization energy of the ligand by itself.
264    *
265    * <p>Upol(L) = L*Upol(1) + (1-L)*(Uenv + Uligand)
266    *
267    * <p>Set polarizationScale to L for part 1. Set polarizationScale to (1-L) for parts 2 & 3.
268    */
269   private double polarizationScale = 1.0;
270   // *************************************************************************
271   // Mutable Particle Mesh Ewald constants.
272   private double aewald;
273   private double an0;
274   private double an1;
275   private double an2;
276   private double an3;
277   private double an4;
278   private double an5;
279   private double permanentEnergy;
280   private double polarizationEnergy;
281   private ScaleParameters scaleParameters;
282 
283   public RealSpaceEnergyRegion(ForceField.ELEC_FORM elecForm, int nt, ForceField forceField, boolean lambdaTerm, double electric) {
284     maxThreads = nt;
285     this.electric = electric;
286     this.elecForm = elecForm;
287     sharedInteractions = new SharedInteger();
288     if (elecForm == ForceField.ELEC_FORM.FIXED_CHARGE) {
289       fixedChargeEnergyLoop = new FixedChargeEnergyLoop[nt];
290       realSpaceEnergyLoop = null;
291     } else {
292       fixedChargeEnergyLoop = null;
293       realSpaceEnergyLoop = new RealSpaceEnergyLoop[nt];
294     }
295 
296     // Flag to indicate application of an intermolecular softcore potential.
297     if (lambdaTerm) {
298       intermolecularSoftcore = forceField.getBoolean("INTERMOLECULAR_SOFTCORE", false);
299       intramolecularSoftcore = forceField.getBoolean("INTRAMOLECULAR_SOFTCORE", false);
300     } else {
301       intermolecularSoftcore = false;
302       intramolecularSoftcore = false;
303     }
304   }
305 
306   @Override
307   public void applyMask(int i, boolean[] is14, double[]... masks) {
308     if (ip11[i] != null) {
309       // AMOEBA
310       double[] permanentEnergyMask = masks[0];
311       double[] permanentFieldMask = masks[1];
312       for (int value : mask12[i]) {
313         permanentFieldMask[value] = scaleParameters.p12scale;
314         permanentEnergyMask[value] = scaleParameters.m12scale;
315       }
316       for (int value : mask13[i]) {
317         permanentFieldMask[value] = scaleParameters.p13scale;
318         permanentEnergyMask[value] = scaleParameters.m13scale;
319       }
320       for (int value : mask14[i]) {
321         permanentFieldMask[value] = scaleParameters.p14scale;
322         permanentEnergyMask[value] = scaleParameters.m14scale;
323         for (int k : ip11[i]) {
324           if (k == value) {
325             permanentFieldMask[value] = scaleParameters.intra14Scale * scaleParameters.p14scale;
326             break;
327           }
328         }
329       }
330       for (int value : mask15[i]) {
331         permanentEnergyMask[value] = scaleParameters.m15scale;
332       }
333       // Apply group based polarization masking rule.
334       double[] polarizationGroupMask = masks[2];
335       for (int j : ip11[i]) {
336         polarizationGroupMask[j] = scaleParameters.d11scale;
337       }
338     } else {
339       // Fixed Charge Force Fields.
340       double[] permanentEnergyMask = masks[0];
341       for (int value : mask12[i]) {
342         permanentEnergyMask[value] = scaleParameters.m12scale;
343       }
344       for (int value : mask13[i]) {
345         permanentEnergyMask[value] = scaleParameters.m13scale;
346       }
347       for (int value : mask14[i]) {
348         permanentEnergyMask[value] = scaleParameters.m14scale;
349       }
350     }
351   }
352 
353   /**
354    * Execute the RealSpaceEnergyRegion with the passed ParallelTeam.
355    *
356    * @param parallelTeam The ParallelTeam instance to execute with.
357    */
358   public void executeWith(ParallelTeam parallelTeam) {
359     try {
360       parallelTeam.execute(this);
361     } catch (Exception e) {
362       String message = " Exception computing the electrostatic energy.\n";
363       logger.log(Level.SEVERE, message, e);
364     }
365   }
366 
367   @Override
368   public void finish() {
369     permanentEnergy = 0.0;
370     polarizationEnergy = 0.0;
371     if (elecForm == ForceField.ELEC_FORM.FIXED_CHARGE) {
372       for (int i = 0; i < maxThreads; i++) {
373         permanentEnergy += fixedChargeEnergyLoop[i].permanentEnergy;
374       }
375     } else {
376       for (int i = 0; i < maxThreads; i++) {
377         permanentEnergy += realSpaceEnergyLoop[i].permanentEnergy;
378         polarizationEnergy += realSpaceEnergyLoop[i].inducedEnergy;
379       }
380     }
381     permanentEnergy *= electric;
382     polarizationEnergy *= electric;
383     if (isNaN(permanentEnergy)) {
384       throw new EnergyException(
385           format(" The permanent multipole energy is %16.8f", permanentEnergy));
386     }
387     if (isNaN(polarizationEnergy)) {
388       throw new EnergyException(
389           format(" The polarization energy is %16.8f", polarizationEnergy));
390     }
391   }
392 
393   public int getCount(int i) {
394     if (elecForm == ForceField.ELEC_FORM.FIXED_CHARGE) {
395       return fixedChargeEnergyLoop[i].getCount();
396     } else {
397       return realSpaceEnergyLoop[i].getCount();
398     }
399   }
400 
401   public int getInteractions() {
402     return sharedInteractions.get();
403   }
404 
405   public double getPermanentEnergy() {
406     return permanentEnergy;
407   }
408 
409   public double getPolarizationEnergy() {
410     return polarizationEnergy;
411   }
412 
413   public void init(
414       Atom[] atoms,
415       Crystal crystal,
416       ExtendedSystem extendedSystem,
417       boolean esvTerm,
418       double[][][] coordinates,
419       MultipoleFrameDefinition[] frame,
420       int[][] axisAtom,
421       double[][][] globalMultipole,
422       double[][][] titrationMultipole,
423       double[][][] tautomerMultipole,
424       double[][][] inducedDipole,
425       double[][][] inducedDipoleCR,
426       boolean[] use,
427       int[] molecule,
428       int[][] ip11,
429       int[][] mask12,
430       int[][] mask13,
431       int[][] mask14,
432       int[][] mask15,
433       boolean[] isSoft,
434       double[] ipdamp,
435       double[] thole,
436       RealSpaceNeighborParameters realSpaceNeighborParameters,
437       boolean gradient,
438       boolean lambdaTerm,
439       boolean nnTerm,
440       LambdaMode lambdaMode,
441       Polarization polarization,
442       EwaldParameters ewaldParameters,
443       ScaleParameters scaleParameters,
444       AlchemicalParameters alchemicalParameters,
445       long[] realSpaceEnergyTime,
446       AtomicDoubleArray3D grad,
447       AtomicDoubleArray3D torque,
448       AtomicDoubleArray3D lambdaGrad,
449       AtomicDoubleArray3D lambdaTorque,
450       SharedDouble shareddEdLambda,
451       SharedDouble sharedd2EdLambda2) {
452     // Input
453     this.atoms = atoms;
454     this.crystal = crystal;
455     this.extendedSystem = extendedSystem;
456     this.esvTerm = esvTerm;
457     this.coordinates = coordinates;
458     this.frame = frame;
459     this.axisAtom = axisAtom;
460     this.globalMultipole = globalMultipole;
461     this.titrationMultipole = titrationMultipole;
462     this.tautomerMultipole = tautomerMultipole;
463     this.inducedDipole = inducedDipole;
464     this.inducedDipoleCR = inducedDipoleCR;
465     this.use = use;
466     this.molecule = molecule;
467     this.ip11 = ip11;
468     this.mask12 = mask12;
469     this.mask13 = mask13;
470     this.mask14 = mask14;
471     this.mask15 = mask15;
472     this.isSoft = isSoft;
473     this.ipdamp = ipdamp;
474     this.thole = thole;
475     this.realSpaceLists = realSpaceNeighborParameters.realSpaceLists;
476     this.realSpaceCounts = realSpaceNeighborParameters.realSpaceCounts;
477     this.realSpaceSchedule = realSpaceNeighborParameters.realSpaceSchedule;
478     this.gradient = gradient;
479     this.lambdaTerm = lambdaTerm;
480     this.nnTerm = nnTerm;
481     this.lambdaMode = lambdaMode;
482     this.polarization = polarization;
483     this.lAlpha = alchemicalParameters.lAlpha;
484     this.dlAlpha = alchemicalParameters.dlAlpha;
485     this.d2lAlpha = alchemicalParameters.d2lAlpha;
486     this.dEdLSign = alchemicalParameters.dEdLSign;
487     this.dlPowPerm = alchemicalParameters.dlPowPerm;
488     this.d2lPowPerm = alchemicalParameters.d2lPowPerm;
489     this.doPermanentRealSpace = alchemicalParameters.doPermanentRealSpace;
490     this.permanentScale = alchemicalParameters.permanentScale;
491     this.lPowPol = alchemicalParameters.lPowPol;
492     this.dlPowPol = alchemicalParameters.dlPowPol;
493     this.d2lPowPol = alchemicalParameters.d2lPowPol;
494     this.doPolarization = alchemicalParameters.doPolarization;
495     this.polarizationScale = alchemicalParameters.polarizationScale;
496     this.aewald = ewaldParameters.aewald;
497     this.an0 = ewaldParameters.an0;
498     this.an1 = ewaldParameters.an1;
499     this.an2 = ewaldParameters.an2;
500     this.an3 = ewaldParameters.an3;
501     this.an4 = ewaldParameters.an4;
502     this.an5 = ewaldParameters.an5;
503     this.scaleParameters = scaleParameters;
504     // Output
505     this.realSpaceEnergyTime = realSpaceEnergyTime;
506     this.grad = grad;
507     this.torque = torque;
508     this.lambdaGrad = lambdaGrad;
509     this.lambdaTorque = lambdaTorque;
510     this.shareddEdLambda = shareddEdLambda;
511     this.sharedd2EdLambda2 = sharedd2EdLambda2;
512   }
513 
514   @Override
515   public void removeMask(int i, boolean[] is14, double[]... masks) {
516     if (ip11[i] != null) {
517       // AMOEBA
518       double[] permanentEnergyMask = masks[0];
519       double[] permanentFieldMask = masks[1];
520       for (int value : mask12[i]) {
521         permanentFieldMask[value] = 1.0;
522         permanentEnergyMask[value] = 1.0;
523       }
524       for (int value : mask13[i]) {
525         permanentFieldMask[value] = 1.0;
526         permanentEnergyMask[value] = 1.0;
527       }
528       for (int value : mask14[i]) {
529         permanentFieldMask[value] = 1.0;
530         permanentEnergyMask[value] = 1.0;
531         for (int k : ip11[i]) {
532           if (k == value) {
533             permanentFieldMask[value] = 1.0;
534             break;
535           }
536         }
537       }
538       for (int value : mask15[i]) {
539         permanentEnergyMask[value] = 1.0;
540       }
541       // Apply group based polarization masking rule.
542       double[] polarizationGroupMask = masks[2];
543       for (int j : ip11[i]) {
544         polarizationGroupMask[j] = 1.0;
545       }
546     } else {
547       // Fixed Charge Force Fields.
548       double[] permanentEnergyMask = masks[0];
549       for (int value : mask12[i]) {
550         permanentEnergyMask[value] = 1.0;
551       }
552       for (int value : mask13[i]) {
553         permanentEnergyMask[value] = 1.0;
554       }
555       for (int value : mask14[i]) {
556         permanentEnergyMask[value] = 1.0;
557       }
558     }
559   }
560 
561   @Override
562   public void run() {
563     int threadIndex = getThreadIndex();
564     if (elecForm == ForceField.ELEC_FORM.FIXED_CHARGE) {
565       if (fixedChargeEnergyLoop[threadIndex] == null) {
566         fixedChargeEnergyLoop[threadIndex] = new FixedChargeEnergyLoop();
567       }
568       try {
569         int nAtoms = atoms.length;
570         execute(0, nAtoms - 1, fixedChargeEnergyLoop[threadIndex]);
571       } catch (Exception e) {
572         String message =
573             "Fatal exception computing the real space energy in thread " + getThreadIndex() + "\n";
574         logger.log(Level.SEVERE, message, e);
575       }
576     } else {
577       if (realSpaceEnergyLoop[threadIndex] == null) {
578         realSpaceEnergyLoop[threadIndex] = new RealSpaceEnergyLoop();
579       }
580       try {
581         int nAtoms = atoms.length;
582         execute(0, nAtoms - 1, realSpaceEnergyLoop[threadIndex]);
583       } catch (Exception e) {
584         String message =
585             "Fatal exception computing the real space energy in thread " + getThreadIndex() + "\n";
586         logger.log(Level.SEVERE, message, e);
587       }
588     }
589   }
590 
591   @Override
592   public void start() {
593     sharedInteractions.set(0);
594   }
595 
596   /**
597    * Log the real space electrostatics interaction.
598    *
599    * @param i   Atom i.
600    * @param k   Atom j.
601    * @param r   The distance rij.
602    * @param eij The interaction energy.
603    * @since 1.0
604    */
605   private void log(int i, int k, double r, double eij, int count, double total) {
606     logger.info(format("%s  %6d  %6d  %12.10f  %12.10f %8d %16.10f",
607         "ELEC", atoms[i].getIndex(), atoms[k].getIndex(), r, eij, count, total));
608   }
609 
610   /**
611    * The Real Space Gradient Loop class contains methods and thread local variables to parallelize
612    * the evaluation of the real space permanent and polarization energies and gradients.
613    */
614   private class RealSpaceEnergyLoop extends IntegerForLoop {
615 
616     private final double[] dx_local;
617     private final double[][] rot_local;
618     private final Torque torques;
619     private double ci;
620     private double dix, diy, diz;
621     private double qixx, qiyy, qizz, qixy, qixz, qiyz;
622     private double ck;
623     private double dkx, dky, dkz;
624     private double qkxx, qkyy, qkzz, qkxy, qkxz, qkyz;
625     private double uix, uiy, uiz;
626     private double pix, piy, piz;
627     private double xr, yr, zr;
628     private double ukx, uky, ukz;
629     private double pkx, pky, pkz;
630     private double bn0, bn1, bn2, bn3, bn4, bn5, bn6;
631     private double rr1, rr3, rr5, rr7, rr9, rr11, rr13;
632     private double scale, scale3, scale5, scale7;
633     private double scalep, scaled;
634     private double ddsc3x, ddsc3y, ddsc3z;
635     private double ddsc5x, ddsc5y, ddsc5z;
636     private double ddsc7x, ddsc7y, ddsc7z;
637     private double l2;
638     private boolean soft;
639     private double selfScale;
640     private double permanentEnergy;
641     private double inducedEnergy;
642     private double dUdL, d2UdL2;
643     private int i, k, iSymm, count;
644     private double[] gxk_local, gyk_local, gzk_local;
645     private double[] txk_local, tyk_local, tzk_local;
646     private double[] lxk_local, lyk_local, lzk_local;
647     private double[] ltxk_local, ltyk_local, ltzk_local;
648     private double[] masking_local;
649     private double[] maskingp_local;
650     private double[] maskingd_local;
651     private int threadID;
652 
653     RealSpaceEnergyLoop() {
654       super();
655       dx_local = new double[3];
656       rot_local = new double[3][3];
657       torques = new Torque();
658     }
659 
660     @Override
661     public void finish() {
662       sharedInteractions.addAndGet(count);
663       if (lambdaTerm) {
664         shareddEdLambda.addAndGet(dUdL * electric);
665         sharedd2EdLambda2.addAndGet(d2UdL2 * electric);
666       }
667       realSpaceEnergyTime[getThreadIndex()] += System.nanoTime();
668     }
669 
670     public int getCount() {
671       return count;
672     }
673 
674     @Override
675     public void run(int lb, int ub) {
676       List<SymOp> symOps = crystal.spaceGroup.symOps;
677       int nSymm = symOps.size();
678       int nAtoms = atoms.length;
679       for (iSymm = 0; iSymm < nSymm; iSymm++) {
680         SymOp symOp = symOps.get(iSymm);
681         if (gradient) {
682           fill(gxk_local, 0.0);
683           fill(gyk_local, 0.0);
684           fill(gzk_local, 0.0);
685           fill(txk_local, 0.0);
686           fill(tyk_local, 0.0);
687           fill(tzk_local, 0.0);
688         }
689         if (lambdaTerm) {
690           fill(lxk_local, 0.0);
691           fill(lyk_local, 0.0);
692           fill(lzk_local, 0.0);
693           fill(ltxk_local, 0.0);
694           fill(ltyk_local, 0.0);
695           fill(ltzk_local, 0.0);
696         }
697         realSpaceChunk(lb, ub);
698         if (gradient) {
699           // Turn symmetry mate torques into gradients
700           int[] frameIndex = {-1, -1, -1, -1};
701           double[][] g = new double[4][3];
702           double[] trq = new double[3];
703           for (int i = 0; i < nAtoms; i++) {
704             fill(frameIndex, -1);
705             for (int j = 0; j < 4; j++) {
706               fill(g[j], 0.0);
707             }
708             trq[0] = txk_local[i];
709             trq[1] = tyk_local[i];
710             trq[2] = tzk_local[i];
711             torques.torque(i, iSymm, trq, frameIndex, g);
712             for (int j = 0; j < 4; j++) {
713               int index = frameIndex[j];
714               if (index >= 0) {
715                 double[] gj = g[j];
716                 gxk_local[index] += gj[0];
717                 gyk_local[index] += gj[1];
718                 gzk_local[index] += gj[2];
719               }
720             }
721           }
722           // Rotate symmetry mate gradients
723           if (iSymm != 0) {
724             crystal.applyTransSymRot(
725                 nAtoms, gxk_local, gyk_local, gzk_local, gxk_local, gyk_local, gzk_local, symOp,
726                 rot_local);
727           }
728           // Sum symmetry mate gradients into asymmetric unit gradients
729           for (int j = 0; j < nAtoms; j++) {
730             grad.add(threadID, j, gxk_local[j], gyk_local[j], gzk_local[j]);
731           }
732         }
733         if (lambdaTerm) {
734           // Turn symmetry mate torques into gradients
735           int[] frameIndex = {-1, -1, -1, -1};
736           double[][] g = new double[4][3];
737           double[] trq = new double[3];
738           for (int i = 0; i < nAtoms; i++) {
739             fill(frameIndex, -1);
740             for (int j = 0; j < 4; j++) {
741               fill(g[j], 0.0);
742             }
743             trq[0] = ltxk_local[i];
744             trq[1] = ltyk_local[i];
745             trq[2] = ltzk_local[i];
746             torques.torque(i, iSymm, trq, frameIndex, g);
747             for (int j = 0; j < 4; j++) {
748               int index = frameIndex[j];
749               if (index >= 0) {
750                 double[] gj = g[j];
751                 lxk_local[index] += gj[0];
752                 lyk_local[index] += gj[1];
753                 lzk_local[index] += gj[2];
754               }
755             }
756           }
757           // Rotate symmetry mate gradients
758           if (iSymm != 0) {
759             crystal.applyTransSymRot(
760                 nAtoms, lxk_local, lyk_local, lzk_local, lxk_local, lyk_local, lzk_local, symOp,
761                 rot_local);
762           }
763           // Sum symmetry mate gradients into asymmetric unit gradients
764           for (int j = 0; j < nAtoms; j++) {
765             lambdaGrad.add(threadID, j, lxk_local[j], lyk_local[j], lzk_local[j]);
766           }
767         }
768       }
769     }
770 
771     @Override
772     public IntegerSchedule schedule() {
773       return realSpaceSchedule;
774     }
775 
776     @Override
777     public void start() {
778       init();
779       threadID = getThreadIndex();
780       realSpaceEnergyTime[threadID] -= System.nanoTime();
781       permanentEnergy = 0.0;
782       inducedEnergy = 0.0;
783       count = 0;
784       if (lambdaTerm) {
785         dUdL = 0.0;
786         d2UdL2 = 0.0;
787       }
788       torques.init(axisAtom, frame, coordinates);
789     }
790 
791     private void init() {
792       int nAtoms = atoms.length;
793       if (masking_local == null || masking_local.length < nAtoms) {
794         txk_local = new double[nAtoms];
795         tyk_local = new double[nAtoms];
796         tzk_local = new double[nAtoms];
797         gxk_local = new double[nAtoms];
798         gyk_local = new double[nAtoms];
799         gzk_local = new double[nAtoms];
800         lxk_local = new double[nAtoms];
801         lyk_local = new double[nAtoms];
802         lzk_local = new double[nAtoms];
803         ltxk_local = new double[nAtoms];
804         ltyk_local = new double[nAtoms];
805         ltzk_local = new double[nAtoms];
806         masking_local = new double[nAtoms];
807         maskingp_local = new double[nAtoms];
808         maskingd_local = new double[nAtoms];
809         fill(masking_local, 1.0);
810         fill(maskingp_local, 1.0);
811         fill(maskingd_local, 1.0);
812       }
813     }
814 
815     /**
816      * Evaluate the real space permanent energy and polarization energy for a chunk of atoms.
817      *
818      * @param lb The lower bound of the chunk.
819      * @param ub The upper bound of the chunk.
820      */
821     private void realSpaceChunk(final int lb, final int ub) {
822       final double[] x = coordinates[0][0];
823       final double[] y = coordinates[0][1];
824       final double[] z = coordinates[0][2];
825       final double[][] mpole = globalMultipole[0];
826       final double[] zeropole = new double[10];
827       final double[][] ind = inducedDipole[0];
828       final double[][] indp = inducedDipoleCR[0];
829       final int[][] lists = realSpaceLists[iSymm];
830       final double[] neighborX = coordinates[iSymm][0];
831       final double[] neighborY = coordinates[iSymm][1];
832       final double[] neighborZ = coordinates[iSymm][2];
833       final double[][] neighborMultipole = globalMultipole[iSymm];
834       final double[][] neighborInducedDipole = inducedDipole[iSymm];
835       final double[][] neighborInducedDipolep = inducedDipoleCR[iSymm];
836       for (i = lb; i <= ub; i++) {
837         if (!use[i]) {
838           continue;
839         }
840         final int moleculei = molecule[i];
841         if (iSymm == 0) {
842           applyMask(i, null, masking_local, maskingp_local, maskingd_local);
843         }
844         final double xi = x[i];
845         final double yi = y[i];
846         final double zi = z[i];
847         final double[] globalMultipolei = mpole[i];
848         final double[] inducedDipolei = ind[i];
849         final double[] inducedDipolepi = indp[i];
850         setMultipoleI(globalMultipolei);
851         setInducedI(inducedDipolei);
852         setInducedpI(inducedDipolepi);
853         final boolean softi = isSoft[i];
854         final double pdi = ipdamp[i];
855         final double pti = thole[i];
856         final int[] list = lists[i];
857         final int npair = realSpaceCounts[iSymm][i];
858         for (int j = 0; j < npair; j++) {
859           k = list[j];
860           if (!use[k]) {
861             continue;
862           }
863           boolean sameMolecule = (moleculei == molecule[k]);
864           if (lambdaMode == LambdaMode.VAPOR) {
865             if ((intermolecularSoftcore && !sameMolecule)
866                 || (intramolecularSoftcore && sameMolecule)) {
867               continue;
868             }
869           }
870           selfScale = 1.0;
871           if (i == k) {
872             selfScale = 0.5;
873           }
874           double beta = 0.0;
875           l2 = 1.0;
876           soft = (softi || isSoft[k]);
877           if (soft && doPermanentRealSpace) {
878             beta = lAlpha;
879             l2 = permanentScale;
880           } else if (nnTerm) {
881             l2 = permanentScale;
882           }
883           final double xk = neighborX[k];
884           final double yk = neighborY[k];
885           final double zk = neighborZ[k];
886           dx_local[0] = xk - xi;
887           dx_local[1] = yk - yi;
888           dx_local[2] = zk - zi;
889           double r2 = crystal.image(dx_local);
890           xr = dx_local[0];
891           yr = dx_local[1];
892           zr = dx_local[2];
893           final double[] globalMultipolek = neighborMultipole[k];
894           final double[] inducedDipolek = neighborInducedDipole[k];
895           final double[] inducedDipolepk = neighborInducedDipolep[k];
896           setMultipoleK(globalMultipolek);
897           setInducedK(inducedDipolek);
898           setInducedpK(inducedDipolepk);
899           final double pdk = ipdamp[k];
900           final double ptk = thole[k];
901           scale = masking_local[k];
902           scalep = maskingp_local[k];
903           scaled = maskingd_local[k];
904           scale3 = 1.0;
905           scale5 = 1.0;
906           scale7 = 1.0;
907           double r = sqrt(r2 + beta);
908           double ralpha = aewald * r;
909           double exp2a = exp(-ralpha * ralpha);
910           rr1 = 1.0 / r;
911           double rr2 = rr1 * rr1;
912           bn0 = erfc(ralpha) * rr1;
913           bn1 = (bn0 + an0 * exp2a) * rr2;
914           bn2 = (3.0 * bn1 + an1 * exp2a) * rr2;
915           bn3 = (5.0 * bn2 + an2 * exp2a) * rr2;
916           bn4 = (7.0 * bn3 + an3 * exp2a) * rr2;
917           bn5 = (9.0 * bn4 + an4 * exp2a) * rr2;
918           bn6 = (11.0 * bn5 + an5 * exp2a) * rr2;
919           rr3 = rr1 * rr2;
920           rr5 = 3.0 * rr3 * rr2;
921           rr7 = 5.0 * rr5 * rr2;
922           rr9 = 7.0 * rr7 * rr2;
923           rr11 = 9.0 * rr9 * rr2;
924           rr13 = 11.0 * rr11 * rr2;
925           ddsc3x = 0.0;
926           ddsc3y = 0.0;
927           ddsc3z = 0.0;
928           ddsc5x = 0.0;
929           ddsc5y = 0.0;
930           ddsc5z = 0.0;
931           ddsc7x = 0.0;
932           ddsc7y = 0.0;
933           ddsc7z = 0.0;
934           double damp = pdi * pdk;
935           double thole = min(pti, ptk);
936           double rdamp = r * damp;
937           damp = -thole * rdamp * rdamp * rdamp;
938           if (damp > -50.0) {
939             final double expdamp = exp(damp);
940             scale3 = 1.0 - expdamp;
941             scale5 = 1.0 - expdamp * (1.0 - damp);
942             scale7 = 1.0 - expdamp * (1.0 - damp + 0.6 * damp * damp);
943             final double temp3 = -3.0 * damp * expdamp * rr2;
944             final double temp5 = -damp;
945             final double temp7 = -0.2 - 0.6 * damp;
946             ddsc3x = temp3 * xr;
947             ddsc3y = temp3 * yr;
948             ddsc3z = temp3 * zr;
949             ddsc5x = temp5 * ddsc3x;
950             ddsc5y = temp5 * ddsc3y;
951             ddsc5z = temp5 * ddsc3z;
952             ddsc7x = temp7 * ddsc5x;
953             ddsc7y = temp7 * ddsc5y;
954             ddsc7z = temp7 * ddsc5z;
955           }
956           if (doPermanentRealSpace) {
957             double ei = permanentPair(gradient, lambdaTerm);
958             if (isNaN(ei) || isInfinite(ei)) {
959               String message = format(" %s\n %s\n %s\n The permanent multipole energy between "
960                       + "atoms %d and %d (%d) is %16.8f at %16.8f A.",
961                   crystal.getUnitCell().toString(), atoms[i].toString(), atoms[k].toString(),
962                   i, k, iSymm, ei, r);
963               throw new EnergyException(message);
964             }
965             if (ei != 0.0) {
966               permanentEnergy += ei;
967               count++;
968               // log(i, k, r, ei * electric, count, permanentEnergy * electric);
969             }
970             if (esvTerm) {
971               if (extendedSystem.isTitrating(i)) {
972                 double titrdUdL = 0.0;
973                 double tautdUdL = 0.0;
974                 final double[][] titrMpole = titrationMultipole[0];
975                 final double[] titrMultipolei = titrMpole[i];
976                 setMultipoleI(titrMultipolei);
977                 titrdUdL = electric * permanentPair(false, false);
978                 if (extendedSystem.isTautomerizing(i)) {
979                   final double[][] tautMpole = tautomerMultipole[0];
980                   final double[] tautMultipolei = tautMpole[i];
981                   setMultipoleI(tautMultipolei);
982                   tautdUdL = electric * permanentPair(false, false);
983                 }
984                 extendedSystem.addPermElecDeriv(i, titrdUdL, tautdUdL);
985                 // Reset Multipoles between titration and tautomer ESVs
986                 setMultipoleI(globalMultipolei);
987               }
988               if (extendedSystem.isTitrating(k)) {
989                 double titrdUdL = 0.0;
990                 double tautdUdL = 0.0;
991                 final double[][] titrNeighborMpole = titrationMultipole[iSymm];
992                 final double[] titrMultipolek = titrNeighborMpole[k];
993                 setMultipoleK(titrMultipolek);
994                 titrdUdL = electric * permanentPair(false, false);
995                 if (extendedSystem.isTautomerizing(k)) {
996                   final double[][] tautNeighborMpole = tautomerMultipole[iSymm];
997                   final double[] tautMultipolek = tautNeighborMpole[k];
998                   setMultipoleK(tautMultipolek);
999                   tautdUdL = electric * permanentPair(false, false);
1000                 }
1001                 extendedSystem.addPermElecDeriv(k, titrdUdL, tautdUdL);
1002                 // Reset Multipoles between titration and tautomer ESVs
1003                 setMultipoleK(globalMultipolek);
1004               }
1005             }
1006           }
1007           if (polarization != Polarization.NONE && doPolarization) {
1008             // Polarization does not use the softcore tensors.
1009             if (soft && doPermanentRealSpace) {
1010               scale3 = 1.0;
1011               scale5 = 1.0;
1012               scale7 = 1.0;
1013               r = sqrt(r2);
1014               ralpha = aewald * r;
1015               exp2a = exp(-ralpha * ralpha);
1016               rr1 = 1.0 / r;
1017               rr2 = rr1 * rr1;
1018               bn0 = erfc(ralpha) * rr1;
1019               bn1 = (bn0 + an0 * exp2a) * rr2;
1020               bn2 = (3.0 * bn1 + an1 * exp2a) * rr2;
1021               bn3 = (5.0 * bn2 + an2 * exp2a) * rr2;
1022               bn4 = (7.0 * bn3 + an3 * exp2a) * rr2;
1023               bn5 = (9.0 * bn4 + an4 * exp2a) * rr2;
1024               bn6 = (11.0 * bn5 + an5 * exp2a) * rr2;
1025               rr3 = rr1 * rr2;
1026               rr5 = 3.0 * rr3 * rr2;
1027               rr7 = 5.0 * rr5 * rr2;
1028               rr9 = 7.0 * rr7 * rr2;
1029               rr11 = 9.0 * rr9 * rr2;
1030               ddsc3x = 0.0;
1031               ddsc3y = 0.0;
1032               ddsc3z = 0.0;
1033               ddsc5x = 0.0;
1034               ddsc5y = 0.0;
1035               ddsc5z = 0.0;
1036               ddsc7x = 0.0;
1037               ddsc7y = 0.0;
1038               ddsc7z = 0.0;
1039               damp = pdi * pdk;
1040               thole = min(pti, ptk);
1041               rdamp = r * damp;
1042               damp = -thole * rdamp * rdamp * rdamp;
1043               if (damp > -50.0) {
1044                 final double expdamp = exp(damp);
1045                 scale3 = 1.0 - expdamp;
1046                 scale5 = 1.0 - expdamp * (1.0 - damp);
1047                 scale7 = 1.0 - expdamp * (1.0 - damp + 0.6 * damp * damp);
1048                 final double temp3 = -3.0 * damp * expdamp * rr2;
1049                 final double temp5 = -damp;
1050                 final double temp7 = -0.2 - 0.6 * damp;
1051                 ddsc3x = temp3 * xr;
1052                 ddsc3y = temp3 * yr;
1053                 ddsc3z = temp3 * zr;
1054                 ddsc5x = temp5 * ddsc3x;
1055                 ddsc5y = temp5 * ddsc3y;
1056                 ddsc5z = temp5 * ddsc3z;
1057                 ddsc7x = temp7 * ddsc5x;
1058                 ddsc7y = temp7 * ddsc5y;
1059                 ddsc7z = temp7 * ddsc5z;
1060               }
1061             }
1062             double ei = polarizationPair(gradient, lambdaTerm);
1063             if (isNaN(ei) || isInfinite(ei)) {
1064               String message = format(" %s\n"
1065                       + " %s\n with induced dipole: %8.3f %8.3f %8.3f\n"
1066                       + " %s\n with induced dipole: %8.3f %8.3f %8.3f\n"
1067                       + " The polarization energy due to atoms "
1068                       + "%d and %d (%d) is %10.6f at %10.6f A.",
1069                   crystal.getUnitCell(), atoms[i], uix, uiy, uiz, atoms[k], ukx, uky, ukz, i + 1, k + 1, iSymm, ei, r);
1070               throw new EnergyException(message);
1071             }
1072             inducedEnergy += ei;
1073             if (esvTerm) {
1074               int esvI = extendedSystem.getTitrationESVIndex(i);
1075               int esvK = extendedSystem.getTitrationESVIndex(k);
1076               if (extendedSystem.isTitrating(i)) {
1077                 double titrdUdL = 0.0;
1078                 double tautdUdL = 0.0;
1079                 final double[][] titrMpole = titrationMultipole[0];
1080                 final double[] titrMultipolei = titrMpole[i];
1081                 setMultipoleI(titrMultipolei);
1082                 // Check if atom i and atom k are controlled by the same ESV.
1083                 if (extendedSystem.isTitrating(k) && esvI == esvK) {
1084                   final double[][] titrNeighborMpole = titrationMultipole[iSymm];
1085                   final double[] titrMultipolek = titrNeighborMpole[k];
1086                   setMultipoleK(titrMultipolek);
1087                 } else {
1088                   setMultipoleK(zeropole);
1089                 }
1090                 titrdUdL = polarizationPair(false, false);
1091                 // Swap induced dipoles and masking rules
1092                 setInducedI(inducedDipolepi);
1093                 setInducedK(inducedDipolepk);
1094                 scaled = maskingp_local[k];
1095                 scalep = maskingd_local[k];
1096                 titrdUdL += polarizationPair(false, false);
1097                 // Reset
1098                 setInducedI(inducedDipolei);
1099                 setInducedK(inducedDipolek);
1100                 scalep = maskingp_local[k];
1101                 scaled = maskingd_local[k];
1102 
1103                 if (extendedSystem.isTautomerizing(i)) {
1104                   final double[][] tautMpole = tautomerMultipole[0];
1105                   final double[] tautMultipolei = tautMpole[i];
1106                   setMultipoleI(tautMultipolei);
1107                   // Check if atom i and atom k are controlled by the same ESV.
1108                   if (extendedSystem.isTautomerizing(k) && esvI == esvK) {
1109                     final double[][] tautNeighborMpole = tautomerMultipole[iSymm];
1110                     final double[] tautMultipolek = tautNeighborMpole[k];
1111                     setMultipoleK(tautMultipolek);
1112                   } else {
1113                     setMultipoleK(zeropole);
1114                   }
1115                   tautdUdL = polarizationPair(false, false);
1116                   // Swap induced dipoles and masking
1117                   setInducedI(inducedDipolepi);
1118                   setInducedK(inducedDipolepk);
1119                   scaled = maskingp_local[k];
1120                   scalep = maskingd_local[k];
1121                   tautdUdL += polarizationPair(false, false);
1122                   // Reset induced dipoles and masking
1123                   setInducedI(inducedDipolei);
1124                   setInducedK(inducedDipolek);
1125                   scalep = maskingp_local[k];
1126                   scaled = maskingd_local[k];
1127                 }
1128                 extendedSystem.addIndElecDeriv(i, titrdUdL * electric, tautdUdL * electric);
1129                 // Reset permanent multipoles
1130                 setMultipoleI(globalMultipolei);
1131                 setMultipoleK(globalMultipolek);
1132               }
1133               // Collect further terms if atom i and atom k are controlled by different ESVs.
1134               if (extendedSystem.isTitrating(k) && esvI != esvK) {
1135                 double titrdUdL = 0.0;
1136                 double tautdUdL = 0.0;
1137                 final double[][] titrNeighborMpole = titrationMultipole[iSymm];
1138                 final double[] titrMultipolek = titrNeighborMpole[k];
1139                 setMultipoleK(titrMultipolek);
1140                 setMultipoleI(zeropole);
1141                 titrdUdL = polarizationPair(false, false);
1142                 // Swap induced dipoles and masking
1143                 setInducedI(inducedDipolepi);
1144                 setInducedK(inducedDipolepk);
1145                 scaled = maskingp_local[k];
1146                 scalep = maskingd_local[k];
1147                 titrdUdL += polarizationPair(false, false);
1148                 // Reset induced dipoles and masking
1149                 setInducedI(inducedDipolei);
1150                 setInducedK(inducedDipolek);
1151                 scalep = maskingp_local[k];
1152                 scaled = maskingd_local[k];
1153                 if (extendedSystem.isTautomerizing(k)) {
1154                   final double[][] tautNeighborMpole = tautomerMultipole[iSymm];
1155                   final double[] tautMultipolek = tautNeighborMpole[k];
1156                   setMultipoleK(tautMultipolek);
1157                   tautdUdL = polarizationPair(false, false);
1158                   // Swap induced dipoles and masking
1159                   setInducedI(inducedDipolepi);
1160                   setInducedK(inducedDipolepk);
1161                   scaled = maskingp_local[k];
1162                   scalep = maskingd_local[k];
1163                   tautdUdL += polarizationPair(false, false);
1164                   // Reset induced dipoles and masking
1165                   setInducedI(inducedDipolei);
1166                   setInducedK(inducedDipolek);
1167                   scalep = maskingp_local[k];
1168                   scaled = maskingd_local[k];
1169                 }
1170                 extendedSystem.addIndElecDeriv(k, titrdUdL * electric, tautdUdL * electric);
1171                 // Reset permanent multipoles
1172                 setMultipoleI(globalMultipolei);
1173                 setMultipoleK(globalMultipolek);
1174               }
1175             }
1176           }
1177         }
1178         if (iSymm == 0) {
1179           removeMask(i, null, masking_local, maskingp_local, maskingd_local);
1180         }
1181       }
1182     }
1183 
1184     /**
1185      * Evaluate the real space permanent energy for a pair of multipole sites.
1186      *
1187      * @return the permanent multipole energy.
1188      */
1189     private double permanentPair(boolean gradientLocal, boolean lambdaTermLocal) {
1190       final double dixdkx = diy * dkz - diz * dky;
1191       final double dixdky = diz * dkx - dix * dkz;
1192       final double dixdkz = dix * dky - diy * dkx;
1193       final double dixrx = diy * zr - diz * yr;
1194       final double dixry = diz * xr - dix * zr;
1195       final double dixrz = dix * yr - diy * xr;
1196       final double dkxrx = dky * zr - dkz * yr;
1197       final double dkxry = dkz * xr - dkx * zr;
1198       final double dkxrz = dkx * yr - dky * xr;
1199       final double qirx = qixx * xr + qixy * yr + qixz * zr;
1200       final double qiry = qixy * xr + qiyy * yr + qiyz * zr;
1201       final double qirz = qixz * xr + qiyz * yr + qizz * zr;
1202       final double qkrx = qkxx * xr + qkxy * yr + qkxz * zr;
1203       final double qkry = qkxy * xr + qkyy * yr + qkyz * zr;
1204       final double qkrz = qkxz * xr + qkyz * yr + qkzz * zr;
1205       final double qiqkrx = qixx * qkrx + qixy * qkry + qixz * qkrz;
1206       final double qiqkry = qixy * qkrx + qiyy * qkry + qiyz * qkrz;
1207       final double qiqkrz = qixz * qkrx + qiyz * qkry + qizz * qkrz;
1208       final double qkqirx = qkxx * qirx + qkxy * qiry + qkxz * qirz;
1209       final double qkqiry = qkxy * qirx + qkyy * qiry + qkyz * qirz;
1210       final double qkqirz = qkxz * qirx + qkyz * qiry + qkzz * qirz;
1211       final double qixqkx =
1212           qixy * qkxz + qiyy * qkyz + qiyz * qkzz - qixz * qkxy - qiyz * qkyy - qizz * qkyz;
1213       final double qixqky =
1214           qixz * qkxx + qiyz * qkxy + qizz * qkxz - qixx * qkxz - qixy * qkyz - qixz * qkzz;
1215       final double qixqkz =
1216           qixx * qkxy + qixy * qkyy + qixz * qkyz - qixy * qkxx - qiyy * qkxy - qiyz * qkxz;
1217       final double rxqirx = yr * qirz - zr * qiry;
1218       final double rxqiry = zr * qirx - xr * qirz;
1219       final double rxqirz = xr * qiry - yr * qirx;
1220       final double rxqkrx = yr * qkrz - zr * qkry;
1221       final double rxqkry = zr * qkrx - xr * qkrz;
1222       final double rxqkrz = xr * qkry - yr * qkrx;
1223       final double rxqikrx = yr * qiqkrz - zr * qiqkry;
1224       final double rxqikry = zr * qiqkrx - xr * qiqkrz;
1225       final double rxqikrz = xr * qiqkry - yr * qiqkrx;
1226       final double rxqkirx = yr * qkqirz - zr * qkqiry;
1227       final double rxqkiry = zr * qkqirx - xr * qkqirz;
1228       final double rxqkirz = xr * qkqiry - yr * qkqirx;
1229       final double qkrxqirx = qkry * qirz - qkrz * qiry;
1230       final double qkrxqiry = qkrz * qirx - qkrx * qirz;
1231       final double qkrxqirz = qkrx * qiry - qkry * qirx;
1232       final double qidkx = qixx * dkx + qixy * dky + qixz * dkz;
1233       final double qidky = qixy * dkx + qiyy * dky + qiyz * dkz;
1234       final double qidkz = qixz * dkx + qiyz * dky + qizz * dkz;
1235       final double qkdix = qkxx * dix + qkxy * diy + qkxz * diz;
1236       final double qkdiy = qkxy * dix + qkyy * diy + qkyz * diz;
1237       final double qkdiz = qkxz * dix + qkyz * diy + qkzz * diz;
1238       final double dixqkrx = diy * qkrz - diz * qkry;
1239       final double dixqkry = diz * qkrx - dix * qkrz;
1240       final double dixqkrz = dix * qkry - diy * qkrx;
1241       final double dkxqirx = dky * qirz - dkz * qiry;
1242       final double dkxqiry = dkz * qirx - dkx * qirz;
1243       final double dkxqirz = dkx * qiry - dky * qirx;
1244       final double rxqidkx = yr * qidkz - zr * qidky;
1245       final double rxqidky = zr * qidkx - xr * qidkz;
1246       final double rxqidkz = xr * qidky - yr * qidkx;
1247       final double rxqkdix = yr * qkdiz - zr * qkdiy;
1248       final double rxqkdiy = zr * qkdix - xr * qkdiz;
1249       final double rxqkdiz = xr * qkdiy - yr * qkdix;
1250 
1251       // Calculate the scalar products for permanent multipoles.
1252       final double sc2 = dix * dkx + diy * dky + diz * dkz;
1253       final double sc3 = dix * xr + diy * yr + diz * zr;
1254       final double sc4 = dkx * xr + dky * yr + dkz * zr;
1255       final double sc5 = qirx * xr + qiry * yr + qirz * zr;
1256       final double sc6 = qkrx * xr + qkry * yr + qkrz * zr;
1257       final double sc7 = qirx * dkx + qiry * dky + qirz * dkz;
1258       final double sc8 = qkrx * dix + qkry * diy + qkrz * diz;
1259       final double sc9 = qirx * qkrx + qiry * qkry + qirz * qkrz;
1260       final double sc10 =
1261           2.0 * (qixy * qkxy + qixz * qkxz + qiyz * qkyz) + qixx * qkxx + qiyy * qkyy + qizz * qkzz;
1262 
1263       // Calculate the gl functions for permanent multipoles.
1264       final double gl0 = ci * ck;
1265       final double gl1 = ck * sc3 - ci * sc4;
1266       final double gl2 = ci * sc6 + ck * sc5 - sc3 * sc4;
1267       final double gl3 = sc3 * sc6 - sc4 * sc5;
1268       final double gl4 = sc5 * sc6;
1269       final double gl5 = -4.0 * sc9;
1270       final double gl6 = sc2;
1271       final double gl7 = 2.0 * (sc7 - sc8);
1272       final double gl8 = 2.0 * sc10;
1273 
1274       // Compute the energy contributions for this interaction.
1275       final double scale1 = 1.0 - scale;
1276       final double ereal =
1277           gl0 * bn0 + (gl1 + gl6) * bn1 + (gl2 + gl7 + gl8) * bn2 + (gl3 + gl5) * bn3 + gl4 * bn4;
1278       final double efix =
1279           scale1
1280               * (gl0 * rr1
1281               + (gl1 + gl6) * rr3
1282               + (gl2 + gl7 + gl8) * rr5
1283               + (gl3 + gl5) * rr7
1284               + gl4 * rr9);
1285       final double e = selfScale * l2 * (ereal - efix);
1286       if (gradientLocal) {
1287         final double gf1 =
1288             bn1 * gl0 + bn2 * (gl1 + gl6) + bn3 * (gl2 + gl7 + gl8) + bn4 * (gl3 + gl5) + bn5 * gl4;
1289         final double gf2 = -ck * bn1 + sc4 * bn2 - sc6 * bn3;
1290         final double gf3 = ci * bn1 + sc3 * bn2 + sc5 * bn3;
1291         final double gf4 = 2.0 * bn2;
1292         final double gf5 = 2.0 * (-ck * bn2 + sc4 * bn3 - sc6 * bn4);
1293         final double gf6 = 2.0 * (-ci * bn2 - sc3 * bn3 - sc5 * bn4);
1294         final double gf7 = 4.0 * bn3;
1295 
1296         // Get the permanent force with screening.
1297         double ftm2x =
1298             gf1 * xr
1299                 + gf2 * dix
1300                 + gf3 * dkx
1301                 + gf4 * (qkdix - qidkx)
1302                 + gf5 * qirx
1303                 + gf6 * qkrx
1304                 + gf7 * (qiqkrx + qkqirx);
1305         double ftm2y =
1306             gf1 * yr
1307                 + gf2 * diy
1308                 + gf3 * dky
1309                 + gf4 * (qkdiy - qidky)
1310                 + gf5 * qiry
1311                 + gf6 * qkry
1312                 + gf7 * (qiqkry + qkqiry);
1313         double ftm2z =
1314             gf1 * zr
1315                 + gf2 * diz
1316                 + gf3 * dkz
1317                 + gf4 * (qkdiz - qidkz)
1318                 + gf5 * qirz
1319                 + gf6 * qkrz
1320                 + gf7 * (qiqkrz + qkqirz);
1321 
1322         // Get the permanent torque with screening.
1323         double ttm2x =
1324             -bn1 * dixdkx
1325                 + gf2 * dixrx
1326                 + gf4 * (dixqkrx + dkxqirx + rxqidkx - 2.0 * qixqkx)
1327                 - gf5 * rxqirx
1328                 - gf7 * (rxqikrx + qkrxqirx);
1329         double ttm2y =
1330             -bn1 * dixdky
1331                 + gf2 * dixry
1332                 + gf4 * (dixqkry + dkxqiry + rxqidky - 2.0 * qixqky)
1333                 - gf5 * rxqiry
1334                 - gf7 * (rxqikry + qkrxqiry);
1335         double ttm2z =
1336             -bn1 * dixdkz
1337                 + gf2 * dixrz
1338                 + gf4 * (dixqkrz + dkxqirz + rxqidkz - 2.0 * qixqkz)
1339                 - gf5 * rxqirz
1340                 - gf7 * (rxqikrz + qkrxqirz);
1341         double ttm3x =
1342             bn1 * dixdkx
1343                 + gf3 * dkxrx
1344                 - gf4 * (dixqkrx + dkxqirx + rxqkdix - 2.0 * qixqkx)
1345                 - gf6 * rxqkrx
1346                 - gf7 * (rxqkirx - qkrxqirx);
1347         double ttm3y =
1348             bn1 * dixdky
1349                 + gf3 * dkxry
1350                 - gf4 * (dixqkry + dkxqiry + rxqkdiy - 2.0 * qixqky)
1351                 - gf6 * rxqkry
1352                 - gf7 * (rxqkiry - qkrxqiry);
1353         double ttm3z =
1354             bn1 * dixdkz
1355                 + gf3 * dkxrz
1356                 - gf4 * (dixqkrz + dkxqirz + rxqkdiz - 2.0 * qixqkz)
1357                 - gf6 * rxqkrz
1358                 - gf7 * (rxqkirz - qkrxqirz);
1359 
1360         // Handle the case where scaling is used.
1361         if (scale1 != 0.0) {
1362           final double gfr1 =
1363               rr3 * gl0
1364                   + rr5 * (gl1 + gl6)
1365                   + rr7 * (gl2 + gl7 + gl8)
1366                   + rr9 * (gl3 + gl5)
1367                   + rr11 * gl4;
1368           final double gfr2 = -ck * rr3 + sc4 * rr5 - sc6 * rr7;
1369           final double gfr3 = ci * rr3 + sc3 * rr5 + sc5 * rr7;
1370           final double gfr4 = 2.0 * rr5;
1371           final double gfr5 = 2.0 * (-ck * rr5 + sc4 * rr7 - sc6 * rr9);
1372           final double gfr6 = 2.0 * (-ci * rr5 - sc3 * rr7 - sc5 * rr9);
1373           final double gfr7 = 4.0 * rr7;
1374 
1375           // Get the permanent force without screening.
1376           final double ftm2rx =
1377               gfr1 * xr
1378                   + gfr2 * dix
1379                   + gfr3 * dkx
1380                   + gfr4 * (qkdix - qidkx)
1381                   + gfr5 * qirx
1382                   + gfr6 * qkrx
1383                   + gfr7 * (qiqkrx + qkqirx);
1384           final double ftm2ry =
1385               gfr1 * yr
1386                   + gfr2 * diy
1387                   + gfr3 * dky
1388                   + gfr4 * (qkdiy - qidky)
1389                   + gfr5 * qiry
1390                   + gfr6 * qkry
1391                   + gfr7 * (qiqkry + qkqiry);
1392           final double ftm2rz =
1393               gfr1 * zr
1394                   + gfr2 * diz
1395                   + gfr3 * dkz
1396                   + gfr4 * (qkdiz - qidkz)
1397                   + gfr5 * qirz
1398                   + gfr6 * qkrz
1399                   + gfr7 * (qiqkrz + qkqirz);
1400 
1401           // Get the permanent torque without screening.
1402           final double ttm2rx =
1403               -rr3 * dixdkx
1404                   + gfr2 * dixrx
1405                   + gfr4 * (dixqkrx + dkxqirx + rxqidkx - 2.0 * qixqkx)
1406                   - gfr5 * rxqirx
1407                   - gfr7 * (rxqikrx + qkrxqirx);
1408           final double ttm2ry =
1409               -rr3 * dixdky
1410                   + gfr2 * dixry
1411                   + gfr4 * (dixqkry + dkxqiry + rxqidky - 2.0 * qixqky)
1412                   - gfr5 * rxqiry
1413                   - gfr7 * (rxqikry + qkrxqiry);
1414           final double ttm2rz =
1415               -rr3 * dixdkz
1416                   + gfr2 * dixrz
1417                   + gfr4 * (dixqkrz + dkxqirz + rxqidkz - 2.0 * qixqkz)
1418                   - gfr5 * rxqirz
1419                   - gfr7 * (rxqikrz + qkrxqirz);
1420           final double ttm3rx =
1421               rr3 * dixdkx
1422                   + gfr3 * dkxrx
1423                   - gfr4 * (dixqkrx + dkxqirx + rxqkdix - 2.0 * qixqkx)
1424                   - gfr6 * rxqkrx
1425                   - gfr7 * (rxqkirx - qkrxqirx);
1426           final double ttm3ry =
1427               rr3 * dixdky
1428                   + gfr3 * dkxry
1429                   - gfr4 * (dixqkry + dkxqiry + rxqkdiy - 2.0 * qixqky)
1430                   - gfr6 * rxqkry
1431                   - gfr7 * (rxqkiry - qkrxqiry);
1432           final double ttm3rz =
1433               rr3 * dixdkz
1434                   + gfr3 * dkxrz
1435                   - gfr4 * (dixqkrz + dkxqirz + rxqkdiz - 2.0 * qixqkz)
1436                   - gfr6 * rxqkrz
1437                   - gfr7 * (rxqkirz - qkrxqirz);
1438           ftm2x -= scale1 * ftm2rx;
1439           ftm2y -= scale1 * ftm2ry;
1440           ftm2z -= scale1 * ftm2rz;
1441           ttm2x -= scale1 * ttm2rx;
1442           ttm2y -= scale1 * ttm2ry;
1443           ttm2z -= scale1 * ttm2rz;
1444           ttm3x -= scale1 * ttm3rx;
1445           ttm3y -= scale1 * ttm3ry;
1446           ttm3z -= scale1 * ttm3rz;
1447         }
1448         double prefactor = electric * selfScale * l2;
1449         grad.add(threadID, i, prefactor * ftm2x, prefactor * ftm2y, prefactor * ftm2z);
1450         torque.add(threadID, i, prefactor * ttm2x, prefactor * ttm2y, prefactor * ttm2z);
1451 
1452 //        logger.info(format(" %d %d I Perm Force %17.15f %17.15f %17.15f", i, k, ftm2x, ftm2y, ftm2z));
1453 //        logger.info(format(" %d %d I Perm Torque %17.15f %17.15f %17.15f", i, k, ttm2x, ttm2y, ttm2z));
1454 //        logger.info(format(" %d %d K Perm Torque %17.15f %17.15f %17.15f", i, k, ttm3x, ttm3y, ttm3z));
1455 
1456         gxk_local[k] -= prefactor * ftm2x;
1457         gyk_local[k] -= prefactor * ftm2y;
1458         gzk_local[k] -= prefactor * ftm2z;
1459         txk_local[k] += prefactor * ttm3x;
1460         tyk_local[k] += prefactor * ttm3y;
1461         tzk_local[k] += prefactor * ttm3z;
1462 
1463         // This is dU/dL/dX for the first term of dU/dL: d[dlPow * ereal]/dx
1464         if (lambdaTerm && soft) {
1465           prefactor = electric * selfScale * dEdLSign * dlPowPerm;
1466           lambdaGrad.add(threadID, i, prefactor * ftm2x, prefactor * ftm2y, prefactor * ftm2z);
1467           lambdaTorque.add(threadID, i, prefactor * ttm2x, prefactor * ttm2y, prefactor * ttm2z);
1468           lxk_local[k] -= prefactor * ftm2x;
1469           lyk_local[k] -= prefactor * ftm2y;
1470           lzk_local[k] -= prefactor * ftm2z;
1471           ltxk_local[k] += prefactor * ttm3x;
1472           ltyk_local[k] += prefactor * ttm3y;
1473           ltzk_local[k] += prefactor * ttm3z;
1474         }
1475       }
1476       if (lambdaTermLocal && soft) {
1477         double dRealdL =
1478             gl0 * bn1 + (gl1 + gl6) * bn2 + (gl2 + gl7 + gl8) * bn3 + (gl3 + gl5) * bn4 + gl4 * bn5;
1479         double d2RealdL2 =
1480             gl0 * bn2 + (gl1 + gl6) * bn3 + (gl2 + gl7 + gl8) * bn4 + (gl3 + gl5) * bn5 + gl4 * bn6;
1481 
1482         dUdL += selfScale * (dEdLSign * dlPowPerm * ereal + l2 * dlAlpha * dRealdL);
1483         d2UdL2 +=
1484             selfScale
1485                 * (dEdLSign
1486                 * (d2lPowPerm * ereal
1487                 + dlPowPerm * dlAlpha * dRealdL
1488                 + dlPowPerm * dlAlpha * dRealdL)
1489                 + l2 * d2lAlpha * dRealdL
1490                 + l2 * dlAlpha * dlAlpha * d2RealdL2);
1491 
1492         double dFixdL =
1493             gl0 * rr3
1494                 + (gl1 + gl6) * rr5
1495                 + (gl2 + gl7 + gl8) * rr7
1496                 + (gl3 + gl5) * rr9
1497                 + gl4 * rr11;
1498         double d2FixdL2 =
1499             gl0 * rr5
1500                 + (gl1 + gl6) * rr7
1501                 + (gl2 + gl7 + gl8) * rr9
1502                 + (gl3 + gl5) * rr11
1503                 + gl4 * rr13;
1504         dFixdL *= scale1;
1505         d2FixdL2 *= scale1;
1506         dUdL -= selfScale * (dEdLSign * dlPowPerm * efix + l2 * dlAlpha * dFixdL);
1507         d2UdL2 -=
1508             selfScale
1509                 * (dEdLSign
1510                 * (d2lPowPerm * efix
1511                 + dlPowPerm * dlAlpha * dFixdL
1512                 + dlPowPerm * dlAlpha * dFixdL)
1513                 + l2 * d2lAlpha * dFixdL
1514                 + l2 * dlAlpha * dlAlpha * d2FixdL2);
1515 
1516         // Collect terms for dU/dL/dX for the second term of dU/dL: d[fL2*dfL1dL*dRealdL]/dX
1517         final double gf1 =
1518             bn2 * gl0 + bn3 * (gl1 + gl6) + bn4 * (gl2 + gl7 + gl8) + bn5 * (gl3 + gl5) + bn6 * gl4;
1519         final double gf2 = -ck * bn2 + sc4 * bn3 - sc6 * bn4;
1520         final double gf3 = ci * bn2 + sc3 * bn3 + sc5 * bn4;
1521         final double gf4 = 2.0 * bn3;
1522         final double gf5 = 2.0 * (-ck * bn3 + sc4 * bn4 - sc6 * bn5);
1523         final double gf6 = 2.0 * (-ci * bn3 - sc3 * bn4 - sc5 * bn5);
1524         final double gf7 = 4.0 * bn4;
1525 
1526         // Get the permanent force with screening.
1527         double ftm2x =
1528             gf1 * xr
1529                 + gf2 * dix
1530                 + gf3 * dkx
1531                 + gf4 * (qkdix - qidkx)
1532                 + gf5 * qirx
1533                 + gf6 * qkrx
1534                 + gf7 * (qiqkrx + qkqirx);
1535         double ftm2y =
1536             gf1 * yr
1537                 + gf2 * diy
1538                 + gf3 * dky
1539                 + gf4 * (qkdiy - qidky)
1540                 + gf5 * qiry
1541                 + gf6 * qkry
1542                 + gf7 * (qiqkry + qkqiry);
1543         double ftm2z =
1544             gf1 * zr
1545                 + gf2 * diz
1546                 + gf3 * dkz
1547                 + gf4 * (qkdiz - qidkz)
1548                 + gf5 * qirz
1549                 + gf6 * qkrz
1550                 + gf7 * (qiqkrz + qkqirz);
1551 
1552         // Get the permanent torque with screening.
1553         double ttm2x =
1554             -bn2 * dixdkx
1555                 + gf2 * dixrx
1556                 + gf4 * (dixqkrx + dkxqirx + rxqidkx - 2.0 * qixqkx)
1557                 - gf5 * rxqirx
1558                 - gf7 * (rxqikrx + qkrxqirx);
1559         double ttm2y =
1560             -bn2 * dixdky
1561                 + gf2 * dixry
1562                 + gf4 * (dixqkry + dkxqiry + rxqidky - 2.0 * qixqky)
1563                 - gf5 * rxqiry
1564                 - gf7 * (rxqikry + qkrxqiry);
1565         double ttm2z =
1566             -bn2 * dixdkz
1567                 + gf2 * dixrz
1568                 + gf4 * (dixqkrz + dkxqirz + rxqidkz - 2.0 * qixqkz)
1569                 - gf5 * rxqirz
1570                 - gf7 * (rxqikrz + qkrxqirz);
1571         double ttm3x =
1572             bn2 * dixdkx
1573                 + gf3 * dkxrx
1574                 - gf4 * (dixqkrx + dkxqirx + rxqkdix - 2.0 * qixqkx)
1575                 - gf6 * rxqkrx
1576                 - gf7 * (rxqkirx - qkrxqirx);
1577         double ttm3y =
1578             bn2 * dixdky
1579                 + gf3 * dkxry
1580                 - gf4 * (dixqkry + dkxqiry + rxqkdiy - 2.0 * qixqky)
1581                 - gf6 * rxqkry
1582                 - gf7 * (rxqkiry - qkrxqiry);
1583         double ttm3z =
1584             bn2 * dixdkz
1585                 + gf3 * dkxrz
1586                 - gf4 * (dixqkrz + dkxqirz + rxqkdiz - 2.0 * qixqkz)
1587                 - gf6 * rxqkrz
1588                 - gf7 * (rxqkirz - qkrxqirz);
1589 
1590         // Handle the case where scaling is used.
1591         if (scale1 != 0.0) {
1592           final double gfr1 =
1593               rr5 * gl0
1594                   + rr7 * (gl1 + gl6)
1595                   + rr9 * (gl2 + gl7 + gl8)
1596                   + rr11 * (gl3 + gl5)
1597                   + rr13 * gl4;
1598           final double gfr2 = -ck * rr5 + sc4 * rr7 - sc6 * rr9;
1599           final double gfr3 = ci * rr5 + sc3 * rr7 + sc5 * rr9;
1600           final double gfr4 = 2.0 * rr7;
1601           final double gfr5 = 2.0 * (-ck * rr7 + sc4 * rr9 - sc6 * rr11);
1602           final double gfr6 = 2.0 * (-ci * rr7 - sc3 * rr9 - sc5 * rr11);
1603           final double gfr7 = 4.0 * rr9;
1604 
1605           // Get the permanent force without screening.
1606           final double ftm2rx =
1607               gfr1 * xr
1608                   + gfr2 * dix
1609                   + gfr3 * dkx
1610                   + gfr4 * (qkdix - qidkx)
1611                   + gfr5 * qirx
1612                   + gfr6 * qkrx
1613                   + gfr7 * (qiqkrx + qkqirx);
1614           final double ftm2ry =
1615               gfr1 * yr
1616                   + gfr2 * diy
1617                   + gfr3 * dky
1618                   + gfr4 * (qkdiy - qidky)
1619                   + gfr5 * qiry
1620                   + gfr6 * qkry
1621                   + gfr7 * (qiqkry + qkqiry);
1622           final double ftm2rz =
1623               gfr1 * zr
1624                   + gfr2 * diz
1625                   + gfr3 * dkz
1626                   + gfr4 * (qkdiz - qidkz)
1627                   + gfr5 * qirz
1628                   + gfr6 * qkrz
1629                   + gfr7 * (qiqkrz + qkqirz);
1630 
1631           // Get the permanent torque without screening.
1632           final double ttm2rx =
1633               -rr5 * dixdkx
1634                   + gfr2 * dixrx
1635                   + gfr4 * (dixqkrx + dkxqirx + rxqidkx - 2.0 * qixqkx)
1636                   - gfr5 * rxqirx
1637                   - gfr7 * (rxqikrx + qkrxqirx);
1638           final double ttm2ry =
1639               -rr5 * dixdky
1640                   + gfr2 * dixry
1641                   + gfr4 * (dixqkry + dkxqiry + rxqidky - 2.0 * qixqky)
1642                   - gfr5 * rxqiry
1643                   - gfr7 * (rxqikry + qkrxqiry);
1644           final double ttm2rz =
1645               -rr5 * dixdkz
1646                   + gfr2 * dixrz
1647                   + gfr4 * (dixqkrz + dkxqirz + rxqidkz - 2.0 * qixqkz)
1648                   - gfr5 * rxqirz
1649                   - gfr7 * (rxqikrz + qkrxqirz);
1650           final double ttm3rx =
1651               rr5 * dixdkx
1652                   + gfr3 * dkxrx
1653                   - gfr4 * (dixqkrx + dkxqirx + rxqkdix - 2.0 * qixqkx)
1654                   - gfr6 * rxqkrx
1655                   - gfr7 * (rxqkirx - qkrxqirx);
1656           final double ttm3ry =
1657               rr5 * dixdky
1658                   + gfr3 * dkxry
1659                   - gfr4 * (dixqkry + dkxqiry + rxqkdiy - 2.0 * qixqky)
1660                   - gfr6 * rxqkry
1661                   - gfr7 * (rxqkiry - qkrxqiry);
1662           final double ttm3rz =
1663               rr5 * dixdkz
1664                   + gfr3 * dkxrz
1665                   - gfr4 * (dixqkrz + dkxqirz + rxqkdiz - 2.0 * qixqkz)
1666                   - gfr6 * rxqkrz
1667                   - gfr7 * (rxqkirz - qkrxqirz);
1668           ftm2x -= scale1 * ftm2rx;
1669           ftm2y -= scale1 * ftm2ry;
1670           ftm2z -= scale1 * ftm2rz;
1671           ttm2x -= scale1 * ttm2rx;
1672           ttm2y -= scale1 * ttm2ry;
1673           ttm2z -= scale1 * ttm2rz;
1674           ttm3x -= scale1 * ttm3rx;
1675           ttm3y -= scale1 * ttm3ry;
1676           ttm3z -= scale1 * ttm3rz;
1677         }
1678 
1679         // Add in dU/dL/dX for the second term of dU/dL: d[lPow*dlAlpha*dRealdL]/dX
1680         double prefactor = electric * selfScale * l2 * dlAlpha;
1681         lambdaGrad.add(threadID, i, prefactor * ftm2x, prefactor * ftm2y, prefactor * ftm2z);
1682         lambdaTorque.add(threadID, i, prefactor * ttm2x, prefactor * ttm2y, prefactor * ttm2z);
1683         lxk_local[k] -= prefactor * ftm2x;
1684         lyk_local[k] -= prefactor * ftm2y;
1685         lzk_local[k] -= prefactor * ftm2z;
1686         ltxk_local[k] += prefactor * ttm3x;
1687         ltyk_local[k] += prefactor * ttm3y;
1688         ltzk_local[k] += prefactor * ttm3z;
1689       }
1690 
1691       return e;
1692     }
1693 
1694     /**
1695      * Evaluate the polarization energy for a pair of polarizable multipole sites.
1696      *
1697      * @return the polarization energy.
1698      */
1699     private double polarizationPair(boolean gradientLocal, boolean lambdaTermLocal) {
1700       final double dsc3 = 1.0 - scale3 * scaled;
1701       final double dsc5 = 1.0 - scale5 * scaled;
1702       final double dsc7 = 1.0 - scale7 * scaled;
1703       final double psc3 = 1.0 - scale3 * scalep;
1704       final double psc5 = 1.0 - scale5 * scalep;
1705       final double psc7 = 1.0 - scale7 * scalep;
1706       final double usc3 = 1.0 - scale3;
1707       final double usc5 = 1.0 - scale5;
1708 
1709       final double dixukx = diy * ukz - diz * uky;
1710       final double dixuky = diz * ukx - dix * ukz;
1711       final double dixukz = dix * uky - diy * ukx;
1712       final double dkxuix = dky * uiz - dkz * uiy;
1713       final double dkxuiy = dkz * uix - dkx * uiz;
1714       final double dkxuiz = dkx * uiy - dky * uix;
1715       final double dixukpx = diy * pkz - diz * pky;
1716       final double dixukpy = diz * pkx - dix * pkz;
1717       final double dixukpz = dix * pky - diy * pkx;
1718       final double dkxuipx = dky * piz - dkz * piy;
1719       final double dkxuipy = dkz * pix - dkx * piz;
1720       final double dkxuipz = dkx * piy - dky * pix;
1721       final double dixrx = diy * zr - diz * yr;
1722       final double dixry = diz * xr - dix * zr;
1723       final double dixrz = dix * yr - diy * xr;
1724       final double dkxrx = dky * zr - dkz * yr;
1725       final double dkxry = dkz * xr - dkx * zr;
1726       final double dkxrz = dkx * yr - dky * xr;
1727       final double qirx = qixx * xr + qixy * yr + qixz * zr;
1728       final double qiry = qixy * xr + qiyy * yr + qiyz * zr;
1729       final double qirz = qixz * xr + qiyz * yr + qizz * zr;
1730       final double qkrx = qkxx * xr + qkxy * yr + qkxz * zr;
1731       final double qkry = qkxy * xr + qkyy * yr + qkyz * zr;
1732       final double qkrz = qkxz * xr + qkyz * yr + qkzz * zr;
1733       final double rxqirx = yr * qirz - zr * qiry;
1734       final double rxqiry = zr * qirx - xr * qirz;
1735       final double rxqirz = xr * qiry - yr * qirx;
1736       final double rxqkrx = yr * qkrz - zr * qkry;
1737       final double rxqkry = zr * qkrx - xr * qkrz;
1738       final double rxqkrz = xr * qkry - yr * qkrx;
1739       final double qiukx = qixx * ukx + qixy * uky + qixz * ukz;
1740       final double qiuky = qixy * ukx + qiyy * uky + qiyz * ukz;
1741       final double qiukz = qixz * ukx + qiyz * uky + qizz * ukz;
1742       final double qkuix = qkxx * uix + qkxy * uiy + qkxz * uiz;
1743       final double qkuiy = qkxy * uix + qkyy * uiy + qkyz * uiz;
1744       final double qkuiz = qkxz * uix + qkyz * uiy + qkzz * uiz;
1745       final double qiukpx = qixx * pkx + qixy * pky + qixz * pkz;
1746       final double qiukpy = qixy * pkx + qiyy * pky + qiyz * pkz;
1747       final double qiukpz = qixz * pkx + qiyz * pky + qizz * pkz;
1748       final double qkuipx = qkxx * pix + qkxy * piy + qkxz * piz;
1749       final double qkuipy = qkxy * pix + qkyy * piy + qkyz * piz;
1750       final double qkuipz = qkxz * pix + qkyz * piy + qkzz * piz;
1751       final double uixqkrx = uiy * qkrz - uiz * qkry;
1752       final double uixqkry = uiz * qkrx - uix * qkrz;
1753       final double uixqkrz = uix * qkry - uiy * qkrx;
1754       final double ukxqirx = uky * qirz - ukz * qiry;
1755       final double ukxqiry = ukz * qirx - ukx * qirz;
1756       final double ukxqirz = ukx * qiry - uky * qirx;
1757       final double uixqkrpx = piy * qkrz - piz * qkry;
1758       final double uixqkrpy = piz * qkrx - pix * qkrz;
1759       final double uixqkrpz = pix * qkry - piy * qkrx;
1760       final double ukxqirpx = pky * qirz - pkz * qiry;
1761       final double ukxqirpy = pkz * qirx - pkx * qirz;
1762       final double ukxqirpz = pkx * qiry - pky * qirx;
1763       final double rxqiukx = yr * qiukz - zr * qiuky;
1764       final double rxqiuky = zr * qiukx - xr * qiukz;
1765       final double rxqiukz = xr * qiuky - yr * qiukx;
1766       final double rxqkuix = yr * qkuiz - zr * qkuiy;
1767       final double rxqkuiy = zr * qkuix - xr * qkuiz;
1768       final double rxqkuiz = xr * qkuiy - yr * qkuix;
1769       final double rxqiukpx = yr * qiukpz - zr * qiukpy;
1770       final double rxqiukpy = zr * qiukpx - xr * qiukpz;
1771       final double rxqiukpz = xr * qiukpy - yr * qiukpx;
1772       final double rxqkuipx = yr * qkuipz - zr * qkuipy;
1773       final double rxqkuipy = zr * qkuipx - xr * qkuipz;
1774       final double rxqkuipz = xr * qkuipy - yr * qkuipx;
1775 
1776       // Calculate the scalar products for permanent multipoles.
1777       final double sc3 = dix * xr + diy * yr + diz * zr;
1778       final double sc4 = dkx * xr + dky * yr + dkz * zr;
1779       final double sc5 = qirx * xr + qiry * yr + qirz * zr;
1780       final double sc6 = qkrx * xr + qkry * yr + qkrz * zr;
1781 
1782       // Calculate the scalar products for polarization components.
1783       final double sci1 = uix * dkx + uiy * dky + uiz * dkz + dix * ukx + diy * uky + diz * ukz;
1784       final double sci3 = uix * xr + uiy * yr + uiz * zr;
1785       final double sci4 = ukx * xr + uky * yr + ukz * zr;
1786       final double sci7 = qirx * ukx + qiry * uky + qirz * ukz;
1787       final double sci8 = qkrx * uix + qkry * uiy + qkrz * uiz;
1788       final double scip1 = pix * dkx + piy * dky + piz * dkz + dix * pkx + diy * pky + diz * pkz;
1789       final double scip2 = uix * pkx + uiy * pky + uiz * pkz + pix * ukx + piy * uky + piz * ukz;
1790       final double scip3 = pix * xr + piy * yr + piz * zr;
1791       final double scip4 = pkx * xr + pky * yr + pkz * zr;
1792       final double scip7 = qirx * pkx + qiry * pky + qirz * pkz;
1793       final double scip8 = qkrx * pix + qkry * piy + qkrz * piz;
1794 
1795       // Calculate the gl functions for polarization components.
1796       final double gli1 = ck * sci3 - ci * sci4;
1797       final double gli2 = -sc3 * sci4 - sci3 * sc4;
1798       final double gli3 = sci3 * sc6 - sci4 * sc5;
1799       final double gli6 = sci1;
1800       final double gli7 = 2.0 * (sci7 - sci8);
1801       final double glip1 = ck * scip3 - ci * scip4;
1802       final double glip2 = -sc3 * scip4 - scip3 * sc4;
1803       final double glip3 = scip3 * sc6 - scip4 * sc5;
1804       final double glip6 = scip1;
1805       final double glip7 = 2.0 * (scip7 - scip8);
1806 
1807       // Compute the energy contributions for this interaction.
1808       final double ereal = (gli1 + gli6) * bn1 + (gli2 + gli7) * bn2 + gli3 * bn3;
1809       final double efix =
1810           (gli1 + gli6) * rr3 * psc3 + (gli2 + gli7) * rr5 * psc5 + gli3 * rr7 * psc7;
1811       final double e = selfScale * 0.5 * (ereal - efix);
1812 
1813       // logger.info(format(" %d %d Polarization Energy (real): %17.15f", i, k, 0.5 * ereal));
1814       // logger.info(format(" %d %d Polarization Energy (fix) : %17.15f", i, k, -0.5 * efix));
1815 
1816       if (!(gradientLocal || lambdaTermLocal)) {
1817         return polarizationScale * e;
1818       }
1819       boolean dorli = false;
1820       if (psc3 != 0.0 || dsc3 != 0.0 || usc3 != 0.0) {
1821         dorli = true;
1822       }
1823 
1824       // Get the induced force with screening.
1825       final double gfi1 =
1826           0.5 * bn2 * (gli1 + glip1 + gli6 + glip6)
1827               + 0.5 * bn2 * scip2
1828               + 0.5 * bn3 * (gli2 + glip2 + gli7 + glip7)
1829               - 0.5 * bn3 * (sci3 * scip4 + scip3 * sci4)
1830               + 0.5 * bn4 * (gli3 + glip3);
1831       final double gfi2 = -ck * bn1 + sc4 * bn2 - sc6 * bn3;
1832       final double gfi3 = ci * bn1 + sc3 * bn2 + sc5 * bn3;
1833       final double gfi4 = 2.0 * bn2;
1834       final double gfi5 = bn3 * (sci4 + scip4);
1835       final double gfi6 = -bn3 * (sci3 + scip3);
1836       double ftm2ix =
1837           gfi1 * xr
1838               + 0.5
1839               * (gfi2 * (uix + pix)
1840               + bn2 * (sci4 * pix + scip4 * uix)
1841               + gfi3 * (ukx + pkx)
1842               + bn2 * (sci3 * pkx + scip3 * ukx)
1843               + (sci4 + scip4) * bn2 * dix
1844               + (sci3 + scip3) * bn2 * dkx
1845               + gfi4 * (qkuix + qkuipx - qiukx - qiukpx))
1846               + gfi5 * qirx
1847               + gfi6 * qkrx;
1848       double ftm2iy =
1849           gfi1 * yr
1850               + 0.5
1851               * (gfi2 * (uiy + piy)
1852               + bn2 * (sci4 * piy + scip4 * uiy)
1853               + gfi3 * (uky + pky)
1854               + bn2 * (sci3 * pky + scip3 * uky)
1855               + (sci4 + scip4) * bn2 * diy
1856               + (sci3 + scip3) * bn2 * dky
1857               + gfi4 * (qkuiy + qkuipy - qiuky - qiukpy))
1858               + gfi5 * qiry
1859               + gfi6 * qkry;
1860       double ftm2iz =
1861           gfi1 * zr
1862               + 0.5
1863               * (gfi2 * (uiz + piz)
1864               + bn2 * (sci4 * piz + scip4 * uiz)
1865               + gfi3 * (ukz + pkz)
1866               + bn2 * (sci3 * pkz + scip3 * ukz)
1867               + (sci4 + scip4) * bn2 * diz
1868               + (sci3 + scip3) * bn2 * dkz
1869               + gfi4 * (qkuiz + qkuipz - qiukz - qiukpz))
1870               + gfi5 * qirz
1871               + gfi6 * qkrz;
1872 
1873       // Get the induced torque with screening.
1874       final double gti2 = 0.5 * bn2 * (sci4 + scip4);
1875       final double gti3 = 0.5 * bn2 * (sci3 + scip3);
1876       final double gti4 = gfi4;
1877       final double gti5 = gfi5;
1878       final double gti6 = gfi6;
1879       double ttm2ix =
1880           -0.5 * bn1 * (dixukx + dixukpx)
1881               + gti2 * dixrx
1882               - gti5 * rxqirx
1883               + 0.5 * gti4 * (ukxqirx + rxqiukx + ukxqirpx + rxqiukpx);
1884       double ttm2iy =
1885           -0.5 * bn1 * (dixuky + dixukpy)
1886               + gti2 * dixry
1887               - gti5 * rxqiry
1888               + 0.5 * gti4 * (ukxqiry + rxqiuky + ukxqirpy + rxqiukpy);
1889       double ttm2iz =
1890           -0.5 * bn1 * (dixukz + dixukpz)
1891               + gti2 * dixrz
1892               - gti5 * rxqirz
1893               + 0.5 * gti4 * (ukxqirz + rxqiukz + ukxqirpz + rxqiukpz);
1894       double ttm3ix =
1895           -0.5 * bn1 * (dkxuix + dkxuipx)
1896               + gti3 * dkxrx
1897               - gti6 * rxqkrx
1898               - 0.5 * gti4 * (uixqkrx + rxqkuix + uixqkrpx + rxqkuipx);
1899       double ttm3iy =
1900           -0.5 * bn1 * (dkxuiy + dkxuipy)
1901               + gti3 * dkxry
1902               - gti6 * rxqkry
1903               - 0.5 * gti4 * (uixqkry + rxqkuiy + uixqkrpy + rxqkuipy);
1904       double ttm3iz =
1905           -0.5 * bn1 * (dkxuiz + dkxuipz)
1906               + gti3 * dkxrz
1907               - gti6 * rxqkrz
1908               - 0.5 * gti4 * (uixqkrz + rxqkuiz + uixqkrpz + rxqkuipz);
1909       double ftm2rix = 0.0;
1910       double ftm2riy = 0.0;
1911       double ftm2riz = 0.0;
1912       double ttm2rix = 0.0;
1913       double ttm2riy = 0.0;
1914       double ttm2riz = 0.0;
1915       double ttm3rix = 0.0;
1916       double ttm3riy = 0.0;
1917       double ttm3riz = 0.0;
1918       if (dorli) {
1919         // Get the induced force without screening.
1920         final double gfri1 =
1921             0.5 * rr5 * ((gli1 + gli6) * psc3 + (glip1 + glip6) * dsc3 + scip2 * usc3)
1922                 + 0.5
1923                 * rr7
1924                 * ((gli7 + gli2) * psc5
1925                 + (glip7 + glip2) * dsc5
1926                 - (sci3 * scip4 + scip3 * sci4) * usc5)
1927                 + 0.5 * rr9 * (gli3 * psc7 + glip3 * dsc7);
1928         final double gfri4 = 2.0 * rr5;
1929         final double gfri5 = rr7 * (sci4 * psc7 + scip4 * dsc7);
1930         final double gfri6 = -rr7 * (sci3 * psc7 + scip3 * dsc7);
1931         ftm2rix =
1932             gfri1 * xr
1933                 + 0.5
1934                 * (-rr3 * ck * (uix * psc3 + pix * dsc3)
1935                 + rr5 * sc4 * (uix * psc5 + pix * dsc5)
1936                 - rr7 * sc6 * (uix * psc7 + pix * dsc7))
1937                 + (rr3 * ci * (ukx * psc3 + pkx * dsc3)
1938                 + rr5 * sc3 * (ukx * psc5 + pkx * dsc5)
1939                 + rr7 * sc5 * (ukx * psc7 + pkx * dsc7))
1940                 * 0.5
1941                 + rr5 * usc5 * (sci4 * pix + scip4 * uix + sci3 * pkx + scip3 * ukx) * 0.5
1942                 + 0.5 * (sci4 * psc5 + scip4 * dsc5) * rr5 * dix
1943                 + 0.5 * (sci3 * psc5 + scip3 * dsc5) * rr5 * dkx
1944                 + 0.5 * gfri4 * ((qkuix - qiukx) * psc5 + (qkuipx - qiukpx) * dsc5)
1945                 + gfri5 * qirx
1946                 + gfri6 * qkrx;
1947         ftm2riy =
1948             gfri1 * yr
1949                 + 0.5
1950                 * (-rr3 * ck * (uiy * psc3 + piy * dsc3)
1951                 + rr5 * sc4 * (uiy * psc5 + piy * dsc5)
1952                 - rr7 * sc6 * (uiy * psc7 + piy * dsc7))
1953                 + (rr3 * ci * (uky * psc3 + pky * dsc3)
1954                 + rr5 * sc3 * (uky * psc5 + pky * dsc5)
1955                 + rr7 * sc5 * (uky * psc7 + pky * dsc7))
1956                 * 0.5
1957                 + rr5 * usc5 * (sci4 * piy + scip4 * uiy + sci3 * pky + scip3 * uky) * 0.5
1958                 + 0.5 * (sci4 * psc5 + scip4 * dsc5) * rr5 * diy
1959                 + 0.5 * (sci3 * psc5 + scip3 * dsc5) * rr5 * dky
1960                 + 0.5 * gfri4 * ((qkuiy - qiuky) * psc5 + (qkuipy - qiukpy) * dsc5)
1961                 + gfri5 * qiry
1962                 + gfri6 * qkry;
1963         ftm2riz =
1964             gfri1 * zr
1965                 + 0.5
1966                 * (-rr3 * ck * (uiz * psc3 + piz * dsc3)
1967                 + rr5 * sc4 * (uiz * psc5 + piz * dsc5)
1968                 - rr7 * sc6 * (uiz * psc7 + piz * dsc7))
1969                 + (rr3 * ci * (ukz * psc3 + pkz * dsc3)
1970                 + rr5 * sc3 * (ukz * psc5 + pkz * dsc5)
1971                 + rr7 * sc5 * (ukz * psc7 + pkz * dsc7))
1972                 * 0.5
1973                 + rr5 * usc5 * (sci4 * piz + scip4 * uiz + sci3 * pkz + scip3 * ukz) * 0.5
1974                 + 0.5 * (sci4 * psc5 + scip4 * dsc5) * rr5 * diz
1975                 + 0.5 * (sci3 * psc5 + scip3 * dsc5) * rr5 * dkz
1976                 + 0.5 * gfri4 * ((qkuiz - qiukz) * psc5 + (qkuipz - qiukpz) * dsc5)
1977                 + gfri5 * qirz
1978                 + gfri6 * qkrz;
1979 
1980         // Get the induced torque without screening.
1981         final double gtri2 = 0.5 * rr5 * (sci4 * psc5 + scip4 * dsc5);
1982         final double gtri3 = 0.5 * rr5 * (sci3 * psc5 + scip3 * dsc5);
1983         final double gtri4 = gfri4;
1984         final double gtri5 = gfri5;
1985         final double gtri6 = gfri6;
1986         ttm2rix =
1987             -rr3 * (dixukx * psc3 + dixukpx * dsc3) * 0.5
1988                 + gtri2 * dixrx
1989                 - gtri5 * rxqirx
1990                 + gtri4 * ((ukxqirx + rxqiukx) * psc5 + (ukxqirpx + rxqiukpx) * dsc5) * 0.5;
1991         ttm2riy =
1992             -rr3 * (dixuky * psc3 + dixukpy * dsc3) * 0.5
1993                 + gtri2 * dixry
1994                 - gtri5 * rxqiry
1995                 + gtri4 * ((ukxqiry + rxqiuky) * psc5 + (ukxqirpy + rxqiukpy) * dsc5) * 0.5;
1996         ttm2riz =
1997             -rr3 * (dixukz * psc3 + dixukpz * dsc3) * 0.5
1998                 + gtri2 * dixrz
1999                 - gtri5 * rxqirz
2000                 + gtri4 * ((ukxqirz + rxqiukz) * psc5 + (ukxqirpz + rxqiukpz) * dsc5) * 0.5;
2001         ttm3rix =
2002             -rr3 * (dkxuix * psc3 + dkxuipx * dsc3) * 0.5
2003                 + gtri3 * dkxrx
2004                 - gtri6 * rxqkrx
2005                 - gtri4 * ((uixqkrx + rxqkuix) * psc5 + (uixqkrpx + rxqkuipx) * dsc5) * 0.5;
2006         ttm3riy =
2007             -rr3 * (dkxuiy * psc3 + dkxuipy * dsc3) * 0.5
2008                 + gtri3 * dkxry
2009                 - gtri6 * rxqkry
2010                 - gtri4 * ((uixqkry + rxqkuiy) * psc5 + (uixqkrpy + rxqkuipy) * dsc5) * 0.5;
2011         ttm3riz =
2012             -rr3 * (dkxuiz * psc3 + dkxuipz * dsc3) * 0.5
2013                 + gtri3 * dkxrz
2014                 - gtri6 * rxqkrz
2015                 - gtri4 * ((uixqkrz + rxqkuiz) * psc5 + (uixqkrpz + rxqkuipz) * dsc5) * 0.5;
2016       }
2017 
2018       // Account for partially excluded induced interactions.
2019       double temp3 = 0.5 * rr3 * ((gli1 + gli6) * scalep + (glip1 + glip6) * scaled);
2020       double temp5 = 0.5 * rr5 * ((gli2 + gli7) * scalep + (glip2 + glip7) * scaled);
2021       final double temp7 = 0.5 * rr7 * (gli3 * scalep + glip3 * scaled);
2022       final double fridmpx = temp3 * ddsc3x + temp5 * ddsc5x + temp7 * ddsc7x;
2023       final double fridmpy = temp3 * ddsc3y + temp5 * ddsc5y + temp7 * ddsc7y;
2024       final double fridmpz = temp3 * ddsc3z + temp5 * ddsc5z + temp7 * ddsc7z;
2025 
2026       // Find some scaling terms for induced-induced force.
2027       temp3 = 0.5 * rr3 * scip2;
2028       temp5 = -0.5 * rr5 * (sci3 * scip4 + scip3 * sci4);
2029       final double findmpx = temp3 * ddsc3x + temp5 * ddsc5x;
2030       final double findmpy = temp3 * ddsc3y + temp5 * ddsc5y;
2031       final double findmpz = temp3 * ddsc3z + temp5 * ddsc5z;
2032 
2033 //      logger.info(format(" %d %d Ewald Polarization grad (i):   %17.15f %17.15f %17.15f", i, k, ftm2ix, ftm2iy, ftm2iz));
2034 //      logger.info(format(" %d %d Ewald Polarization torque (i): %17.15f %17.15f %17.15f", i, k, ttm2ix, ttm2iy, ttm2iz));
2035 //      logger.info(format(" %d %d Ewald Polarization torque (k): %17.15f %17.15f %17.15f", i, k, ttm3ix, ttm3iy, ttm3iz));
2036 
2037       // Modify the forces for partially excluded interactions.
2038       ftm2ix = ftm2ix - fridmpx - findmpx;
2039       ftm2iy = ftm2iy - fridmpy - findmpy;
2040       ftm2iz = ftm2iz - fridmpz - findmpz;
2041 
2042       // Correction to convert mutual to direct polarization force.
2043       if (polarization == Polarization.DIRECT) {
2044         final double gfd = 0.5 * (bn2 * scip2 - bn3 * (scip3 * sci4 + sci3 * scip4));
2045         final double gfdr = 0.5 * (rr5 * scip2 * usc3 - rr7 * (scip3 * sci4 + sci3 * scip4) * usc5);
2046         ftm2ix = ftm2ix - gfd * xr - 0.5 * bn2 * (sci4 * pix + scip4 * uix + sci3 * pkx + scip3 * ukx);
2047         ftm2iy = ftm2iy - gfd * yr - 0.5 * bn2 * (sci4 * piy + scip4 * uiy + sci3 * pky + scip3 * uky);
2048         ftm2iz = ftm2iz - gfd * zr - 0.5 * bn2 * (sci4 * piz + scip4 * uiz + sci3 * pkz + scip3 * ukz);
2049         final double fdirx = gfdr * xr + 0.5 * usc5 * rr5 * (sci4 * pix + scip4 * uix + sci3 * pkx + scip3 * ukx);
2050         final double fdiry = gfdr * yr + 0.5 * usc5 * rr5 * (sci4 * piy + scip4 * uiy + sci3 * pky + scip3 * uky);
2051         final double fdirz = gfdr * zr + 0.5 * usc5 * rr5 * (sci4 * piz + scip4 * uiz + sci3 * pkz + scip3 * ukz);
2052         ftm2ix = ftm2ix + fdirx + findmpx;
2053         ftm2iy = ftm2iy + fdiry + findmpy;
2054         ftm2iz = ftm2iz + fdirz + findmpz;
2055       }
2056 
2057       // Correction for OPT induced dipoles.
2058       //            if (scfAlgorithm == ParticleMeshEwald.SCFAlgorithm.EPT) {
2059       //                double eptx = 0.0;
2060       //                double epty = 0.0;
2061       //                double eptz = 0.0;
2062       //                for (int jj = 0; jj < optOrder; jj++) {
2063       //                    double optix = optDipole[jj][i][0];
2064       //                    double optiy = optDipole[jj][i][1];
2065       //                    double optiz = optDipole[jj][i][2];
2066       //                    double optixp = optDipoleCR[jj][i][0];
2067       //                    double optiyp = optDipoleCR[jj][i][1];
2068       //                    double optizp = optDipoleCR[jj][i][2];
2069       //                    double uirm = optix * xr + optiy * yr + optiz * zr;
2070       //                    for (int mm = 0; mm < optOrder - jj; mm++) {
2071       //                        double optkx = optDipole[mm][k][0];
2072       //                        double optky = optDipole[mm][k][1];
2073       //                        double optkz = optDipole[mm][k][2];
2074       //                        double optkxp = optDipoleCR[mm][k][0];
2075       //                        double optkyp = optDipoleCR[mm][k][1];
2076       //                        double optkzp = optDipoleCR[mm][k][2];
2077       //                        double ukrm = optkx * xr + optky * yr + optkz * zr;
2078       //                        double term1 = bn2 - usc3 * rr5;
2079       //                        double term2 = bn3 - usc5 * rr7;
2080       //                        double term3 = usr5 + term1;
2081       //                        double term4 = rr3;
2082       //                        double term5 = -xr * term3 + ddsc3x * term4;
2083       //                        double term6 = -(bn2 - usc5 * rr5) + xr * xr * term2 - rr5 * xr *
2084       // ddsc5x;
2085       //                        double tixx = optix * term5 + uirm * term6;
2086       //                        double tkxx = optkx * term5 + ukrm * term6;
2087       //                        term5 = -yr * term3 + ddsc3y * term4;
2088       //                        term6 = -usr5 + yr * yr * term2 - rr5 * yr * ddsc5y;
2089       //                        double tiyy = optiy * term5 + uirm * term6;
2090       //                        double tkyy = optky * term5 + ukrm * term6;
2091       //                        term5 = -zr * term3 + ddsc3z * term4;
2092       //                        term6 = -usr5 + zr * zr * term2 - rr5 * zr * ddsc5z;
2093       //                        double tizz = optiz * term5 + uirm * term6;
2094       //                        double tkzz = optkz * term5 + ukrm * term6;
2095       //                        term4 = -usr5 * yr;
2096       //                        term5 = -xr * term1 + rr3 * ddsc3x;
2097       //                        term6 = xr * yr * term2 - rr5 * yr * ddsc5x;
2098       //                        double tixy = optix * term4 + optiy * term5 + uirm * term6;
2099       //                        double tkxy = optkx * term4 + optky * term5 + ukrm * term6;
2100       //                        term4 = -usr5 * zr;
2101       //                        term6 = xr * zr * term2 - rr5 * zr * ddsc5x;
2102       //                        double tixz = optix * term4 + optiz * term5 + uirm * term6;
2103       //                        double tkxz = optkx * term4 + optkz * term5 + ukrm * term6;
2104       //                        term5 = -yr * term1 + rr3 * ddsc3y;
2105       //                        term6 = yr * zr * term2 - rr5 * zr * ddsc5y;
2106       //                        double tiyz = optiy * term4 + optiz * term5 + uirm * term6;
2107       //                        double tkyz = optkz * term4 + optkz * term5 + ukrm * term6;
2108       //                        double depx = tixx * optkxp + tkxx * optixp
2109       //                                + tixy * optkyp + tkxy * optiyp
2110       //                                + tixz * optkzp + tkxz * optizp;
2111       //                        double depy = tixy * optkxp + tkxy * optixp
2112       //                                + tiyy * optkyp + tkyy * optiyp
2113       //                                + tiyz * optkzp + tkyz * optizp;
2114       //                        double depz = tixz * optkxp + tkxz * optixp
2115       //                                + tiyz * optkyp + tkyz * optiyp
2116       //                                + tizz * optkzp + tkzz * optizp;
2117       //                        double optCoefSum = optRegion.optCoefficientsSum[jj + mm + 1];
2118       //                        double fx = optCoefSum * depx;
2119       //                        double fy = optCoefSum * depy;
2120       //                        double fz = optCoefSum * depz;
2121       //                        eptx += fx;
2122       //                        epty += fy;
2123       //                        eptz += fz;
2124       //                    }
2125       //                }
2126       //                if (i == 0 && k == 1) {
2127       //                    logger.info(format(" %s %16.8f %16.8f %16.8f", "ept force", eptx, epty,
2128       // eptz));
2129       //                }
2130       //                ftm2ix += eptx;
2131       //                ftm2iy += epty;
2132       //                ftm2iz += eptz;
2133       //            }
2134 
2135       // Handle the case where scaling is used.
2136 
2137 //      logger.info(format(" %d %d Polarization Coulomb grad (i):   %17.15f %17.15f %17.15f", i, k,
2138 //          -ftm2rix - fridmpx - findmpx, -ftm2riy - fridmpy - findmpy, -ftm2riz - fridmpz - findmpz));
2139 //      logger.info(format(" %d %d Polarization Thole torque (i): %17.15f %17.15f %17.15f", i, k, -ttm2rix, -ttm2riy, -ttm2riz));
2140 //      logger.info(format(" %d %d Polarization Thole torque (k): %17.15f %17.15f %17.15f", i, k, -ttm3rix, -ttm3riy, -ttm3riz));
2141 
2142       ftm2ix = ftm2ix - ftm2rix;
2143       ftm2iy = ftm2iy - ftm2riy;
2144       ftm2iz = ftm2iz - ftm2riz;
2145       ttm2ix = ttm2ix - ttm2rix;
2146       ttm2iy = ttm2iy - ttm2riy;
2147       ttm2iz = ttm2iz - ttm2riz;
2148       ttm3ix = ttm3ix - ttm3rix;
2149       ttm3iy = ttm3iy - ttm3riy;
2150       ttm3iz = ttm3iz - ttm3riz;
2151 
2152       double scalar = electric * polarizationScale * selfScale;
2153 
2154       grad.add(threadID, i, scalar * ftm2ix, scalar * ftm2iy, scalar * ftm2iz);
2155       torque.add(threadID, i, scalar * ttm2ix, scalar * ttm2iy, scalar * ttm2iz);
2156       gxk_local[k] -= scalar * ftm2ix;
2157       gyk_local[k] -= scalar * ftm2iy;
2158       gzk_local[k] -= scalar * ftm2iz;
2159       txk_local[k] += scalar * ttm3ix;
2160       tyk_local[k] += scalar * ttm3iy;
2161       tzk_local[k] += scalar * ttm3iz;
2162       if (lambdaTermLocal) {
2163         dUdL += dEdLSign * dlPowPol * e;
2164         d2UdL2 += dEdLSign * d2lPowPol * e;
2165         scalar = electric * dEdLSign * dlPowPol * selfScale;
2166         lambdaGrad.add(threadID, i, scalar * ftm2ix, scalar * ftm2iy, scalar * ftm2iz);
2167         lambdaTorque.add(threadID, i, scalar * ttm2ix, scalar * ttm2iy, scalar * ttm2iz);
2168         lxk_local[k] -= scalar * ftm2ix;
2169         lyk_local[k] -= scalar * ftm2iy;
2170         lzk_local[k] -= scalar * ftm2iz;
2171         ltxk_local[k] += scalar * ttm3ix;
2172         ltyk_local[k] += scalar * ttm3iy;
2173         ltzk_local[k] += scalar * ttm3iz;
2174       }
2175       return polarizationScale * e;
2176     }
2177 
2178     private void setMultipoleI(double[] globalMultipolei) {
2179       ci = globalMultipolei[t000];
2180       dix = globalMultipolei[t100];
2181       diy = globalMultipolei[t010];
2182       diz = globalMultipolei[t001];
2183       qixx = globalMultipolei[t200] * oneThird;
2184       qiyy = globalMultipolei[t020] * oneThird;
2185       qizz = globalMultipolei[t002] * oneThird;
2186       qixy = globalMultipolei[t110] * oneThird;
2187       qixz = globalMultipolei[t101] * oneThird;
2188       qiyz = globalMultipolei[t011] * oneThird;
2189     }
2190 
2191     private void setMultipoleK(double[] globalMultipolek) {
2192       ck = globalMultipolek[t000];
2193       dkx = globalMultipolek[t100];
2194       dky = globalMultipolek[t010];
2195       dkz = globalMultipolek[t001];
2196       qkxx = globalMultipolek[t200] * oneThird;
2197       qkyy = globalMultipolek[t020] * oneThird;
2198       qkzz = globalMultipolek[t002] * oneThird;
2199       qkxy = globalMultipolek[t110] * oneThird;
2200       qkxz = globalMultipolek[t101] * oneThird;
2201       qkyz = globalMultipolek[t011] * oneThird;
2202     }
2203 
2204     private void setInducedI(double[] inducedDipolei) {
2205       uix = inducedDipolei[0];
2206       uiy = inducedDipolei[1];
2207       uiz = inducedDipolei[2];
2208     }
2209 
2210     private void setInducedK(double[] inducedDipolek) {
2211       ukx = inducedDipolek[0];
2212       uky = inducedDipolek[1];
2213       ukz = inducedDipolek[2];
2214     }
2215 
2216     private void setInducedpI(double[] inducedDipolepi) {
2217       pix = inducedDipolepi[0];
2218       piy = inducedDipolepi[1];
2219       piz = inducedDipolepi[2];
2220     }
2221 
2222     private void setInducedpK(double[] inducedDipolepk) {
2223       pkx = inducedDipolepk[0];
2224       pky = inducedDipolepk[1];
2225       pkz = inducedDipolepk[2];
2226     }
2227 
2228   }
2229 
2230   /**
2231    * The Real Space Gradient Loop class contains methods and thread local variables to parallelize
2232    * the evaluation of the real space permanent and polarization energies and gradients.
2233    */
2234   private class FixedChargeEnergyLoop extends IntegerForLoop {
2235 
2236     private final double[] dx_local;
2237     private final double[][] rot_local;
2238     private double ci;
2239     private double ck;
2240     private double xr, yr, zr;
2241     private double bn0, bn1, bn2;
2242     private double rr1, rr3, rr5;
2243     private double scale;
2244     private double l2;
2245     private boolean soft;
2246     private double selfScale;
2247     private double permanentEnergy;
2248     private double inducedEnergy;
2249     private double dUdL, d2UdL2;
2250     private int i, k, iSymm, count;
2251     private double[] gxk_local, gyk_local, gzk_local;
2252     private double[] lxk_local, lyk_local, lzk_local;
2253     private double[] masking_local;
2254     private int threadID;
2255 
2256     FixedChargeEnergyLoop() {
2257       super();
2258       dx_local = new double[3];
2259       rot_local = new double[3][3];
2260     }
2261 
2262     @Override
2263     public void finish() {
2264       sharedInteractions.addAndGet(count);
2265       if (lambdaTerm) {
2266         shareddEdLambda.addAndGet(dUdL * electric);
2267         sharedd2EdLambda2.addAndGet(d2UdL2 * electric);
2268       }
2269       realSpaceEnergyTime[getThreadIndex()] += System.nanoTime();
2270     }
2271 
2272     public int getCount() {
2273       return count;
2274     }
2275 
2276     @Override
2277     public void run(int lb, int ub) {
2278       List<SymOp> symOps = crystal.spaceGroup.symOps;
2279       int nSymm = symOps.size();
2280       int nAtoms = atoms.length;
2281       for (iSymm = 0; iSymm < nSymm; iSymm++) {
2282         SymOp symOp = symOps.get(iSymm);
2283         if (gradient) {
2284           fill(gxk_local, 0.0);
2285           fill(gyk_local, 0.0);
2286           fill(gzk_local, 0.0);
2287         }
2288         if (lambdaTerm) {
2289           fill(lxk_local, 0.0);
2290           fill(lyk_local, 0.0);
2291           fill(lzk_local, 0.0);
2292         }
2293         realSpaceChunk(lb, ub);
2294         if (gradient) {
2295           // Rotate symmetry mate gradients
2296           if (iSymm != 0) {
2297             crystal.applyTransSymRot(nAtoms, gxk_local, gyk_local, gzk_local, gxk_local, gyk_local, gzk_local, symOp, rot_local);
2298           }
2299           // Sum symmetry mate gradients into asymmetric unit gradients
2300           for (int j = 0; j < nAtoms; j++) {
2301             grad.add(threadID, j, gxk_local[j], gyk_local[j], gzk_local[j]);
2302           }
2303         }
2304         if (lambdaTerm) {
2305           // Rotate symmetry mate gradients
2306           if (iSymm != 0) {
2307             crystal.applyTransSymRot(
2308                 nAtoms, lxk_local, lyk_local, lzk_local, lxk_local, lyk_local, lzk_local, symOp,
2309                 rot_local);
2310           }
2311           // Sum symmetry mate gradients into asymmetric unit gradients
2312           for (int j = 0; j < nAtoms; j++) {
2313             lambdaGrad.add(threadID, j, lxk_local[j], lyk_local[j], lzk_local[j]);
2314           }
2315         }
2316       }
2317     }
2318 
2319     @Override
2320     public IntegerSchedule schedule() {
2321       return realSpaceSchedule;
2322     }
2323 
2324     @Override
2325     public void start() {
2326       init();
2327       threadID = getThreadIndex();
2328       realSpaceEnergyTime[threadID] -= System.nanoTime();
2329       permanentEnergy = 0.0;
2330       inducedEnergy = 0.0;
2331       count = 0;
2332       if (lambdaTerm) {
2333         dUdL = 0.0;
2334         d2UdL2 = 0.0;
2335       }
2336     }
2337 
2338     private void init() {
2339       int nAtoms = atoms.length;
2340       if (masking_local == null || masking_local.length < nAtoms) {
2341         gxk_local = new double[nAtoms];
2342         gyk_local = new double[nAtoms];
2343         gzk_local = new double[nAtoms];
2344         lxk_local = new double[nAtoms];
2345         lyk_local = new double[nAtoms];
2346         lzk_local = new double[nAtoms];
2347         masking_local = new double[nAtoms];
2348         fill(masking_local, 1.0);
2349       }
2350     }
2351 
2352     /**
2353      * Evaluate the real space permanent energy and polarization energy for a chunk of atoms.
2354      *
2355      * @param lb The lower bound of the chunk.
2356      * @param ub The upper bound of the chunk.
2357      */
2358     private void realSpaceChunk(final int lb, final int ub) {
2359       final double[] x = coordinates[0][0];
2360       final double[] y = coordinates[0][1];
2361       final double[] z = coordinates[0][2];
2362       final double[][] mpole = globalMultipole[0];
2363       final int[][] lists = realSpaceLists[iSymm];
2364       final double[] neighborX = coordinates[iSymm][0];
2365       final double[] neighborY = coordinates[iSymm][1];
2366       final double[] neighborZ = coordinates[iSymm][2];
2367       final double[][] neighborMultipole = globalMultipole[iSymm];
2368       for (i = lb; i <= ub; i++) {
2369         if (!use[i]) {
2370           continue;
2371         }
2372         final int moleculei = molecule[i];
2373         if (iSymm == 0) {
2374           applyMask(i, null, masking_local);
2375         }
2376         final double xi = x[i];
2377         final double yi = y[i];
2378         final double zi = z[i];
2379         final double[] globalMultipolei = mpole[i];
2380         setMultipoleI(globalMultipolei);
2381         final boolean softi = isSoft[i];
2382         final int[] list = lists[i];
2383         final int npair = realSpaceCounts[iSymm][i];
2384         for (int j = 0; j < npair; j++) {
2385           k = list[j];
2386           if (!use[k]) {
2387             continue;
2388           }
2389           boolean sameMolecule = (moleculei == molecule[k]);
2390           if (lambdaMode == LambdaMode.VAPOR) {
2391             if ((intermolecularSoftcore && !sameMolecule)
2392                 || (intramolecularSoftcore && sameMolecule)) {
2393               continue;
2394             }
2395           }
2396           selfScale = 1.0;
2397           if (i == k) {
2398             selfScale = 0.5;
2399           }
2400           double beta = 0.0;
2401           l2 = 1.0;
2402           soft = (softi || isSoft[k]);
2403           if (soft && doPermanentRealSpace) {
2404             beta = lAlpha;
2405             l2 = permanentScale;
2406           } else if (nnTerm) {
2407             l2 = permanentScale;
2408           }
2409           final double xk = neighborX[k];
2410           final double yk = neighborY[k];
2411           final double zk = neighborZ[k];
2412           dx_local[0] = xk - xi;
2413           dx_local[1] = yk - yi;
2414           dx_local[2] = zk - zi;
2415           double r2 = crystal.image(dx_local);
2416           xr = dx_local[0];
2417           yr = dx_local[1];
2418           zr = dx_local[2];
2419           final double[] globalMultipolek = neighborMultipole[k];
2420           setMultipoleK(globalMultipolek);
2421           scale = masking_local[k];
2422           double r = sqrt(r2 + beta);
2423           double ralpha = aewald * r;
2424           double exp2a = exp(-ralpha * ralpha);
2425           rr1 = 1.0 / r;
2426           double rr2 = rr1 * rr1;
2427           bn0 = erfc(ralpha) * rr1;
2428           bn1 = (bn0 + an0 * exp2a) * rr2;
2429           bn2 = (3.0 * bn1 + an1 * exp2a) * rr2;
2430           rr3 = rr1 * rr2;
2431           rr5 = 3.0 * rr3 * rr2;
2432           if (doPermanentRealSpace) {
2433             double ei = permanentPair(gradient, lambdaTerm);
2434             if (isNaN(ei) || isInfinite(ei)) {
2435               String message = format(" %s\n %s\n %s\n The permanent multipole energy between "
2436                       + "atoms %d and %d (%d) is %16.8f at %16.8f A.",
2437                   crystal.getUnitCell().toString(), atoms[i].toString(), atoms[k].toString(),
2438                   i, k, iSymm, ei, r);
2439               throw new EnergyException(message);
2440             }
2441             if (ei != 0.0) {
2442               permanentEnergy += ei;
2443               count++;
2444             }
2445             if (esvTerm) {
2446               if (extendedSystem.isTitrating(i)) {
2447                 double titrdUdL = 0.0;
2448                 double tautdUdL = 0.0;
2449                 final double[][] titrMpole = titrationMultipole[0];
2450                 final double[] titrMultipolei = titrMpole[i];
2451                 setMultipoleI(titrMultipolei);
2452                 titrdUdL = electric * permanentPair(false, false);
2453                 if (extendedSystem.isTautomerizing(i)) {
2454                   final double[][] tautMpole = tautomerMultipole[0];
2455                   final double[] tautMultipolei = tautMpole[i];
2456                   setMultipoleI(tautMultipolei);
2457                   tautdUdL = electric * permanentPair(false, false);
2458                 }
2459                 extendedSystem.addPermElecDeriv(i, titrdUdL, tautdUdL);
2460                 // Reset Multipoles between titration and tautomer ESVs
2461                 setMultipoleI(globalMultipolei);
2462               }
2463               if (extendedSystem.isTitrating(k)) {
2464                 double titrdUdL = 0.0;
2465                 double tautdUdL = 0.0;
2466                 final double[][] titrNeighborMpole = titrationMultipole[iSymm];
2467                 final double[] titrMultipolek = titrNeighborMpole[k];
2468                 setMultipoleK(titrMultipolek);
2469                 titrdUdL = electric * permanentPair(false, false);
2470                 if (extendedSystem.isTautomerizing(k)) {
2471                   final double[][] tautNeighborMpole = tautomerMultipole[iSymm];
2472                   final double[] tautMultipolek = tautNeighborMpole[k];
2473                   setMultipoleK(tautMultipolek);
2474                   tautdUdL = electric * permanentPair(false, false);
2475                 }
2476                 extendedSystem.addPermElecDeriv(k, titrdUdL, tautdUdL);
2477                 // Reset Multipoles between titration and tautomer ESVs
2478                 setMultipoleK(globalMultipolek);
2479               }
2480             }
2481           }
2482         }
2483         if (iSymm == 0) {
2484           removeMask(i, null, masking_local);
2485         }
2486       }
2487     }
2488 
2489     /**
2490      * Evaluate the real space permanent energy for a pair of multipole sites.
2491      *
2492      * @return the permanent multipole energy.
2493      */
2494     private double permanentPair(boolean gradientLocal, boolean lambdaTermLocal) {
2495 
2496       // Calculate the gl functions for permanent multipoles.
2497       final double gl0 = ci * ck;
2498 
2499       // Compute the energy contributions for this interaction.
2500       final double scale1 = 1.0 - scale;
2501       final double ereal = gl0 * bn0;
2502       final double efix = scale1 * (gl0 * rr1);
2503       final double e = selfScale * l2 * (ereal - efix);
2504       if (gradientLocal) {
2505         final double gf1 = bn1 * gl0;
2506         // Get the permanent force with screening.
2507         double ftm2x = gf1 * xr;
2508         double ftm2y = gf1 * yr;
2509         double ftm2z = gf1 * zr;
2510 
2511         // Handle the case where scaling is used.
2512         if (scale1 != 0.0) {
2513           final double gfr1 = rr3 * gl0;
2514           // Get the permanent force without screening.
2515           final double ftm2rx = gfr1 * xr;
2516           final double ftm2ry = gfr1 * yr;
2517           final double ftm2rz = gfr1 * zr;
2518           ftm2x -= scale1 * ftm2rx;
2519           ftm2y -= scale1 * ftm2ry;
2520           ftm2z -= scale1 * ftm2rz;
2521         }
2522         double prefactor = electric * selfScale * l2;
2523         grad.add(threadID, i, prefactor * ftm2x, prefactor * ftm2y, prefactor * ftm2z);
2524         gxk_local[k] -= prefactor * ftm2x;
2525         gyk_local[k] -= prefactor * ftm2y;
2526         gzk_local[k] -= prefactor * ftm2z;
2527 
2528         // This is dU/dL/dX for the first term of dU/dL: d[dlPow * ereal]/dx
2529         if (lambdaTerm && soft) {
2530           prefactor = electric * selfScale * dEdLSign * dlPowPerm;
2531           lambdaGrad.add(threadID, i, prefactor * ftm2x, prefactor * ftm2y, prefactor * ftm2z);
2532           lxk_local[k] -= prefactor * ftm2x;
2533           lyk_local[k] -= prefactor * ftm2y;
2534           lzk_local[k] -= prefactor * ftm2z;
2535         }
2536       }
2537 
2538       if (lambdaTermLocal && soft) {
2539         double dRealdL = gl0 * bn1;
2540         double d2RealdL2 = gl0 * bn2;
2541 
2542         dUdL += selfScale * (dEdLSign * dlPowPerm * ereal + l2 * dlAlpha * dRealdL);
2543         d2UdL2 += selfScale * (dEdLSign * (d2lPowPerm * ereal
2544             + dlPowPerm * dlAlpha * dRealdL
2545             + dlPowPerm * dlAlpha * dRealdL)
2546             + l2 * d2lAlpha * dRealdL
2547             + l2 * dlAlpha * dlAlpha * d2RealdL2);
2548 
2549         double dFixdL = gl0 * rr3;
2550         double d2FixdL2 = gl0 * rr5;
2551         dFixdL *= scale1;
2552         d2FixdL2 *= scale1;
2553         dUdL -= selfScale * (dEdLSign * dlPowPerm * efix + l2 * dlAlpha * dFixdL);
2554         d2UdL2 -= selfScale * (dEdLSign * (d2lPowPerm * efix
2555             + dlPowPerm * dlAlpha * dFixdL
2556             + dlPowPerm * dlAlpha * dFixdL)
2557             + l2 * d2lAlpha * dFixdL
2558             + l2 * dlAlpha * dlAlpha * d2FixdL2);
2559 
2560         // Collect terms for dU/dL/dX for the second term of dU/dL: d[fL2*dfL1dL*dRealdL]/dX
2561         final double gf1 = bn2 * gl0;
2562         // Get the permanent force with screening.
2563         double ftm2x = gf1 * xr;
2564         double ftm2y = gf1 * yr;
2565         double ftm2z = gf1 * zr;
2566 
2567         // Handle the case where scaling is used.
2568         if (scale1 != 0.0) {
2569           final double gfr1 = rr5 * gl0;
2570           // Get the permanent force without screening.
2571           final double ftm2rx = gfr1 * xr;
2572           final double ftm2ry = gfr1 * yr;
2573           final double ftm2rz = gfr1 * zr;
2574           ftm2x -= scale1 * ftm2rx;
2575           ftm2y -= scale1 * ftm2ry;
2576           ftm2z -= scale1 * ftm2rz;
2577         }
2578 
2579         // Add in dU/dL/dX for the second term of dU/dL: d[lPow*dlAlpha*dRealdL]/dX
2580         double prefactor = electric * selfScale * l2 * dlAlpha;
2581         lambdaGrad.add(threadID, i, prefactor * ftm2x, prefactor * ftm2y, prefactor * ftm2z);
2582         lxk_local[k] -= prefactor * ftm2x;
2583         lyk_local[k] -= prefactor * ftm2y;
2584         lzk_local[k] -= prefactor * ftm2z;
2585       }
2586       return e;
2587     }
2588 
2589     private void setMultipoleI(double[] globalMultipolei) {
2590       ci = globalMultipolei[t000];
2591     }
2592 
2593     private void setMultipoleK(double[] globalMultipolek) {
2594       ck = globalMultipolek[t000];
2595     }
2596 
2597   }
2598 
2599 }