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 static ffx.potential.parameters.MultipoleType.t000;
41  import static ffx.potential.parameters.MultipoleType.t001;
42  import static ffx.potential.parameters.MultipoleType.t002;
43  import static ffx.potential.parameters.MultipoleType.t003;
44  import static ffx.potential.parameters.MultipoleType.t010;
45  import static ffx.potential.parameters.MultipoleType.t011;
46  import static ffx.potential.parameters.MultipoleType.t012;
47  import static ffx.potential.parameters.MultipoleType.t020;
48  import static ffx.potential.parameters.MultipoleType.t021;
49  import static ffx.potential.parameters.MultipoleType.t030;
50  import static ffx.potential.parameters.MultipoleType.t100;
51  import static ffx.potential.parameters.MultipoleType.t101;
52  import static ffx.potential.parameters.MultipoleType.t102;
53  import static ffx.potential.parameters.MultipoleType.t110;
54  import static ffx.potential.parameters.MultipoleType.t111;
55  import static ffx.potential.parameters.MultipoleType.t120;
56  import static ffx.potential.parameters.MultipoleType.t200;
57  import static ffx.potential.parameters.MultipoleType.t201;
58  import static ffx.potential.parameters.MultipoleType.t210;
59  import static ffx.potential.parameters.MultipoleType.t300;
60  import static java.lang.Math.PI;
61  import static java.lang.String.format;
62  import static org.apache.commons.math3.util.FastMath.sqrt;
63  
64  import edu.rit.pj.IntegerForLoop;
65  import edu.rit.pj.IntegerSchedule;
66  import edu.rit.pj.ParallelRegion;
67  import edu.rit.pj.ParallelTeam;
68  import edu.rit.pj.reduction.SharedDouble;
69  import ffx.crystal.Crystal;
70  import ffx.numerics.atomic.AtomicDoubleArray3D;
71  import ffx.numerics.multipole.MultipoleUtilities;
72  import ffx.potential.bonded.Atom;
73  import ffx.potential.extended.ExtendedSystem;
74  import ffx.potential.nonbonded.ReciprocalSpace;
75  import java.util.logging.Level;
76  import java.util.logging.Logger;
77  
78  /**
79   * Parallel evaluation of the PME reciprocal space energy and gradient.
80   *
81   * @author Michael J. Schnieders
82   * @since 1.0
83   */
84  public class ReciprocalEnergyRegion extends ParallelRegion {
85  
86    private static final Logger logger = Logger.getLogger(ReciprocalEnergyRegion.class.getName());
87  
88    private static final double SQRT_PI = sqrt(PI);
89    private static final int tensorCount = MultipoleUtilities.tensorCount(3);
90    private final double electric;
91    private final double aewald;
92    private final double aewald1;
93    private final double aewald2;
94    private final double aewald3;
95    private final double aewald4;
96    private final double oneThird;
97    private final double twoThirds;
98    private final int maxThreads;
99    private final SharedDouble inducedDipoleSelfEnergy;
100   private final SharedDouble inducedDipoleRecipEnergy;
101   private final PermanentReciprocalEnergyLoop[] permanentReciprocalEnergyLoop;
102   private final InducedDipoleReciprocalEnergyLoop[] inducedDipoleReciprocalEnergyLoop;
103   /** Dimensions of [nsymm][nAtoms][3] */
104   public double[][][] inducedDipole;
105   public double[][][] inducedDipoleCR;
106   /** Convenience reference to inducedDipole[0] */
107   private double[][] ind;
108   /** Convenience reference to inducedDipoleCR[0] */
109   private double[][] indCR;
110   private double nfftX, nfftY, nfftZ;
111   private double[][] multipole;
112   private double[][] fracMultipole;
113   private double[][] fracMultipolePhi;
114   private double[][] fracInducedDipolePhi;
115   private double[][] fracInducedDipolePhiCR;
116   private double chargeCorrectionEnergy;
117   private double permanentSelfEnergy;
118   private double permanentReciprocalEnergy;
119   /** An ordered array of atoms in the system. */
120   private Atom[] atoms;
121   /** Unit cell and spacegroup information. */
122   private Crystal crystal;
123   /**
124    * When computing the polarization energy at Lambda there are 3 pieces.
125    *
126    * <p>1.) Upol(1) = The polarization energy computed normally (ie. system with ligand).
127    *
128    * <p>2.) Uenv = The polarization energy of the system without the ligand.
129    *
130    * <p>3.) Uligand = The polarization energy of the ligand by itself.
131    *
132    * <p>Upol(L) = L*Upol(1) + (1-L)*(Uenv + Uligand)
133    *
134    * <p>Set the "use" array to true for all atoms for part 1. Set the "use" array to true for all
135    * atoms except the ligand for part 2. Set the "use" array to true only for the ligand atoms for
136    * part 3.
137    *
138    * <p>The "use" array can also be employed to turn off atoms for computing the electrostatic
139    * energy of sub-structures.
140    */
141   private boolean[] use;
142   /** Dimensions of [nsymm][nAtoms][10] */
143   private double[][][] globalMultipole;
144   private double[][][] globalFracMultipole;
145   private double[][][] titrationMultipole;
146   private double[][][] tautomerMultipole;
147   private ExtendedSystem extendedSystem = null;
148   private boolean esvTerm = false;
149   private double[][] cartMultipolePhi;
150   private double[][] cartesianDipolePhi;
151   private double[][] cartesianDipolePhiCR;
152   /** Reciprocal space instance. */
153   private ReciprocalSpace reciprocalSpace;
154   private Polarization polarization;
155   /** Atomic Gradient array. */
156   private AtomicDoubleArray3D grad;
157   /** Atomic Torque array. */
158   private AtomicDoubleArray3D torque;
159   /** Partial derivative of the gradient with respect to Lambda. */
160   private AtomicDoubleArray3D lambdaGrad;
161   /** Partial derivative of the torque with respect to Lambda. */
162   private AtomicDoubleArray3D lambdaTorque;
163   /** If true, compute coordinate gradient. */
164   private boolean gradient;
165   /**
166    * If lambdaTerm is true, some ligand atom interactions with the environment are being turned
167    * on/off.
168    */
169   private boolean lambdaTerm;
170   /** Partial derivative with respect to Lambda. */
171   private SharedDouble shareddEdLambda;
172   /** Second partial derivative with respect to Lambda. */
173   private SharedDouble sharedd2EdLambda2;
174   // Alchemical control
175   private double permanentScale;
176   private double dlPowPerm;
177   private double d2lPowPerm;
178   private double polarizationScale;
179   private double dlPowPol;
180   private double d2lPowPol;
181   private double dEdLSign;
182 
183   public ReciprocalEnergyRegion(int nt, double aewald, double electric) {
184     permanentReciprocalEnergyLoop = new PermanentReciprocalEnergyLoop[nt];
185     inducedDipoleReciprocalEnergyLoop = new InducedDipoleReciprocalEnergyLoop[nt];
186     inducedDipoleSelfEnergy = new SharedDouble();
187     inducedDipoleRecipEnergy = new SharedDouble();
188     maxThreads = nt;
189     this.electric = electric;
190     this.aewald = aewald;
191     aewald1 = -electric * aewald / SQRT_PI;
192     aewald2 = 2.0 * aewald * aewald;
193     aewald3 = -2.0 / 3.0 * electric * aewald * aewald * aewald / SQRT_PI;
194     aewald4 = -2.0 * aewald3;
195     oneThird = 1.0 / 3.0;
196     twoThirds = 2.0 / 3.0;
197   }
198 
199   /**
200    * Execute the ReciprocalEnergyRegion with the passed ParallelTeam.
201    *
202    * @param parallelTeam The ParallelTeam instance to execute with.
203    */
204   public void executeWith(ParallelTeam parallelTeam) {
205     try {
206       parallelTeam.execute(this);
207     } catch (Exception e) {
208       String message = " Exception computing the electrostatic energy.\n";
209       logger.log(Level.SEVERE, message, e);
210     }
211   }
212 
213   @Override
214   public void finish() {
215     /*
216      The permanent multipole self energy contributions are large
217      enough that rounding differences that result from threads
218      finishing in different orders removes deterministic behavior.
219     */
220     permanentSelfEnergy = 0.0;
221     permanentReciprocalEnergy = 0.0;
222     double totalCharge = 0.0;
223     for (int i = 0; i < maxThreads; i++) {
224       totalCharge += permanentReciprocalEnergyLoop[i].totalCharge;
225       permanentSelfEnergy += permanentReciprocalEnergyLoop[i].eSelf;
226       permanentReciprocalEnergy += permanentReciprocalEnergyLoop[i].eRecip;
227     }
228 
229     if (totalCharge == 0.0) {
230       chargeCorrectionEnergy = 0.0;
231     } else {
232       Crystal unitCell = crystal.getUnitCell();
233       int nSymm = unitCell.getNumSymOps();
234       // Compute the total charge for the unit cell from the asymmetric unit total.
235       totalCharge = totalCharge * nSymm;
236       double denom = unitCell.volume * aewald * aewald;
237       chargeCorrectionEnergy = -0.5 * electric * PI * (totalCharge * totalCharge) / denom;
238       // Normalize the unit cell correction to the asymmetric unit.
239       chargeCorrectionEnergy = chargeCorrectionEnergy / nSymm;
240       if (lambdaTerm) {
241         shareddEdLambda.addAndGet(chargeCorrectionEnergy * dlPowPerm * dEdLSign);
242         sharedd2EdLambda2.addAndGet(chargeCorrectionEnergy * d2lPowPerm * dEdLSign);
243       }
244       chargeCorrectionEnergy = chargeCorrectionEnergy * permanentScale;
245     }
246   }
247 
248   public double getInducedDipoleReciprocalEnergy() {
249     return inducedDipoleRecipEnergy.get();
250   }
251 
252   public double getInducedDipoleSelfEnergy() {
253     return inducedDipoleSelfEnergy.get();
254   }
255 
256   public double getPermanentReciprocalEnergy() {
257     return permanentReciprocalEnergy;
258   }
259 
260   public double getPermanentSelfEnergy() {
261     return permanentSelfEnergy;
262   }
263 
264   public double getPermanentChargeCorrectionEnergy() {
265     return chargeCorrectionEnergy;
266   }
267 
268   public void init(
269       Atom[] atoms,
270       Crystal crystal,
271       boolean gradient,
272       boolean lambdaTerm,
273       boolean esvTerm,
274       boolean[] use,
275       double[][][] globalMultipole,
276       double[][][] globalFracMultipole,
277       double[][][] titrationMultipole,
278       double[][][] tautomerMultipole,
279       double[][] cartMultipolePhi,
280       double[][] fracMultipolePhi,
281       Polarization polarization,
282       double[][][] inducedDipole,
283       double[][][] inducedDipoleCR,
284       double[][] cartesianDipolePhi,
285       double[][] cartesianDipolePhiCR,
286       double[][] fracInducedDipolePhi,
287       double[][] fracInducedDipolePhiCR,
288       ReciprocalSpace reciprocalSpace,
289       AlchemicalParameters alchemicalParameters,
290       ExtendedSystem extendedSystem,
291       AtomicDoubleArray3D grad,
292       AtomicDoubleArray3D torque,
293       AtomicDoubleArray3D lambdaGrad,
294       AtomicDoubleArray3D lambdaTorque,
295       SharedDouble shareddEdLambda,
296       SharedDouble sharedd2EdLambda2
297       ) {
298     this.atoms = atoms;
299     this.crystal = crystal;
300     this.gradient = gradient;
301     this.lambdaTerm = lambdaTerm;
302     this.esvTerm = esvTerm;
303     this.use = use;
304     // Permanent
305     this.globalMultipole = globalMultipole;
306     this.globalFracMultipole = globalFracMultipole;
307     this.titrationMultipole = titrationMultipole;
308     this.tautomerMultipole = tautomerMultipole;
309     this.cartMultipolePhi = cartMultipolePhi;
310     this.fracMultipolePhi = fracMultipolePhi;
311     // Polarization
312     this.polarization = polarization;
313     this.inducedDipole = inducedDipole;
314     this.inducedDipoleCR = inducedDipoleCR;
315     this.cartesianDipolePhi = cartesianDipolePhi;
316     this.cartesianDipolePhiCR = cartesianDipolePhiCR;
317     this.fracInducedDipolePhi = fracInducedDipolePhi;
318     this.fracInducedDipolePhiCR = fracInducedDipolePhiCR;
319     this.reciprocalSpace = reciprocalSpace;
320     this.extendedSystem = extendedSystem;
321     // Alchemical control
322     this.permanentScale = alchemicalParameters.permanentScale;
323     this.dlPowPerm = alchemicalParameters.dlPowPerm;
324     this.d2lPowPerm = alchemicalParameters.d2lPowPerm;
325     this.polarizationScale = alchemicalParameters.polarizationScale;
326     this.dlPowPol = alchemicalParameters.dlPowPol;
327     this.d2lPowPol = alchemicalParameters.d2lPowPol;
328     this.dEdLSign = alchemicalParameters.dEdLSign;
329     // Output
330     this.grad = grad;
331     this.torque = torque;
332     this.lambdaGrad = lambdaGrad;
333     this.lambdaTorque = lambdaTorque;
334     this.shareddEdLambda = shareddEdLambda;
335     this.sharedd2EdLambda2 = sharedd2EdLambda2;
336   }
337 
338   @Override
339   public void run() throws Exception {
340     int threadIndex = getThreadIndex();
341     if (permanentReciprocalEnergyLoop[threadIndex] == null) {
342       permanentReciprocalEnergyLoop[threadIndex] = new PermanentReciprocalEnergyLoop();
343       inducedDipoleReciprocalEnergyLoop[threadIndex] = new InducedDipoleReciprocalEnergyLoop();
344     }
345     try {
346       int nAtoms = atoms.length;
347       execute(0, nAtoms - 1, permanentReciprocalEnergyLoop[threadIndex]);
348       if (polarization != Polarization.NONE) {
349         execute(0, nAtoms - 1, inducedDipoleReciprocalEnergyLoop[threadIndex]);
350       }
351     } catch (Exception e) {
352       String message =
353           "Fatal exception computing the real space field in thread " + threadIndex + "\n";
354       logger.log(Level.SEVERE, message, e);
355     }
356   }
357 
358   @Override
359   public void start() {
360     multipole = globalMultipole[0];
361     fracMultipole = globalFracMultipole[0];
362     ind = inducedDipole[0];
363     indCR = inducedDipoleCR[0];
364     inducedDipoleSelfEnergy.set(0.0);
365     inducedDipoleRecipEnergy.set(0.0);
366     nfftX = reciprocalSpace.getXDim();
367     nfftY = reciprocalSpace.getYDim();
368     nfftZ = reciprocalSpace.getZDim();
369   }
370 
371   private class PermanentReciprocalEnergyLoop extends IntegerForLoop {
372 
373     int threadID;
374     double totalCharge;
375     double eSelf;
376     double eRecip;
377 
378     @Override
379     public void finish() {
380       eSelf *= permanentScale;
381       eRecip *= permanentScale * 0.5 * electric;
382     }
383 
384     @Override
385     public void run(int lb, int ub) {
386 
387       // Permanent multipole self energy and gradient.
388       for (int i = lb; i <= ub; i++) {
389         if (use[i]) {
390           double[] in = globalMultipole[0][i];
391           totalCharge += in[t000];
392           double cii = in[t000] * in[t000];
393           double dii = in[t100] * in[t100] + in[t010] * in[t010] + in[t001] * in[t001];
394           double qii = in[t200] * in[t200] + in[t020] * in[t020] + in[t002] * in[t002]
395               + 2.0 * (in[t110] * in[t110] + in[t101] * in[t101] + in[t011] * in[t011]);
396           eSelf += aewald1 * (cii + aewald2 * (dii / 3.0 + 2.0 * aewald2 * qii / 45.0));
397 
398           if (esvTerm && extendedSystem.isTitrating(i)) {
399             double[] indot = titrationMultipole[0][i];
400             double ciidot = indot[t000] * in[t000];
401             double diidot = indot[t100] * in[t100] + indot[t010] * in[t010] + indot[t001] * in[t001];
402             double qiidot = indot[t200] * in[t200] + indot[t020] * in[t020] + indot[t002] * in[t002]
403                 + 2.0 * (indot[t110] * in[t110] + indot[t101] * in[t101] + indot[t011] * in[t011]);
404             double eSelfTitrDot =
405                 2.0 * aewald1 * (ciidot + aewald2 * (diidot / 3.0 + 2.0 * aewald2 * qiidot / 45.0));
406             double eSelfTautDot = 0.0;
407             if (extendedSystem.isTautomerizing(i)) {
408               indot = tautomerMultipole[0][i];
409               ciidot = indot[t000] * in[t000];
410               diidot = indot[t100] * in[t100] + indot[t010] * in[t010] + indot[t001] * in[t001];
411               qiidot = indot[t200] * in[t200] + indot[t020] * in[t020] + indot[t002] * in[t002]
412                   + 2.0 * (indot[t110] * in[t110] + indot[t101] * in[t101] + indot[t011] * in[t011]);
413               eSelfTautDot = 2.0 * aewald1 * (ciidot + aewald2 * (diidot / 3.0
414                   + 2.0 * aewald2 * qiidot / 45.0));
415             }
416             double factor = permanentScale;
417             extendedSystem.addPermElecDeriv(i, eSelfTitrDot * factor, eSelfTautDot * factor);
418           }
419         }
420       }
421       if (lambdaTerm) {
422         shareddEdLambda.addAndGet(eSelf * dlPowPerm * dEdLSign);
423         sharedd2EdLambda2.addAndGet(eSelf * d2lPowPerm * dEdLSign);
424       }
425 
426       // Permanent multipole reciprocal space energy and gradient.
427       final double[][] recip = crystal.getUnitCell().A;
428 
429       double dUdL = 0.0;
430       double d2UdL2 = 0.0;
431       for (int i = lb; i <= ub; i++) {
432         if (use[i]) {
433           final double[] phi = cartMultipolePhi[i];
434           final double[] mpole = multipole[i];
435           final double[] fmpole = fracMultipole[i];
436           double e = mpole[t000] * phi[t000]
437               + mpole[t100] * phi[t100] + mpole[t010] * phi[t010] + mpole[t001] * phi[t001]
438               + oneThird * (mpole[t200] * phi[t200] + mpole[t020] * phi[t020]
439               + mpole[t002] * phi[t002]
440               + 2.0 * (mpole[t110] * phi[t110] + mpole[t101] * phi[t101] + mpole[t011] * phi[t011]));
441           eRecip += e;
442           if (esvTerm && extendedSystem.isTitrating(i)) {
443             double[] mpoleDot = titrationMultipole[0][i];
444             double edotTitr = mpoleDot[t000] * phi[t000]
445                 + mpoleDot[t100] * phi[t100] + mpoleDot[t010] * phi[t010]
446                 + mpoleDot[t001] * phi[t001]
447                 + oneThird * (mpoleDot[t200] * phi[t200] + mpoleDot[t020] * phi[t020]
448                 + mpoleDot[t002] * phi[t002]
449                 + 2.0 * (mpoleDot[t110] * phi[t110] + mpoleDot[t101] * phi[t101]
450                 + mpoleDot[t011] * phi[t011]));
451             double edotTaut = 0.0;
452             if (extendedSystem.isTautomerizing(i)) {
453               mpoleDot = tautomerMultipole[0][i];
454               edotTaut = mpoleDot[t000] * phi[t000]
455                   + mpoleDot[t100] * phi[t100] + mpoleDot[t010] * phi[t010]
456                   + mpoleDot[t001] * phi[t001]
457                   + oneThird * (mpoleDot[t200] * phi[t200] + mpoleDot[t020] * phi[t020]
458                   + mpoleDot[t002] * phi[t002]
459                   + 2.0 * (mpoleDot[t110] * phi[t110] + mpoleDot[t101] * phi[t101]
460                   + mpoleDot[t011] * phi[t011]));
461             }
462             double factor = permanentScale * electric;
463             extendedSystem.addPermElecDeriv(i, edotTitr * factor, edotTaut * factor);
464           }
465           if (gradient || lambdaTerm) {
466             final double[] fPhi = fracMultipolePhi[i];
467             double gx = fmpole[t000] * fPhi[t100]
468                 + fmpole[t100] * fPhi[t200] + fmpole[t010] * fPhi[t110] + fmpole[t001] * fPhi[t101]
469                 + fmpole[t200] * fPhi[t300] + fmpole[t020] * fPhi[t120] + fmpole[t002] * fPhi[t102]
470                 + fmpole[t110] * fPhi[t210] + fmpole[t101] * fPhi[t201] + fmpole[t011] * fPhi[t111];
471             double gy = fmpole[t000] * fPhi[t010]
472                 + fmpole[t100] * fPhi[t110] + fmpole[t010] * fPhi[t020] + fmpole[t001] * fPhi[t011]
473                 + fmpole[t200] * fPhi[t210] + fmpole[t020] * fPhi[t030] + fmpole[t002] * fPhi[t012]
474                 + fmpole[t110] * fPhi[t120] + fmpole[t101] * fPhi[t111] + fmpole[t011] * fPhi[t021];
475             double gz = fmpole[t000] * fPhi[t001]
476                 + fmpole[t100] * fPhi[t101] + fmpole[t010] * fPhi[t011] + fmpole[t001] * fPhi[t002]
477                 + fmpole[t200] * fPhi[t201] + fmpole[t020] * fPhi[t021] + fmpole[t002] * fPhi[t003]
478                 + fmpole[t110] * fPhi[t111] + fmpole[t101] * fPhi[t102] + fmpole[t011] * fPhi[t012];
479             gx *= nfftX;
480             gy *= nfftY;
481             gz *= nfftZ;
482             final double dfx = recip[0][0] * gx + recip[0][1] * gy + recip[0][2] * gz;
483             final double dfy = recip[1][0] * gx + recip[1][1] * gy + recip[1][2] * gz;
484             final double dfz = recip[2][0] * gx + recip[2][1] * gy + recip[2][2] * gz;
485             // Compute dipole torques
486             double tqx = -mpole[t010] * phi[t001] + mpole[t001] * phi[t010];
487             double tqy = -mpole[t001] * phi[t100] + mpole[t100] * phi[t001];
488             double tqz = -mpole[t100] * phi[t010] + mpole[t010] * phi[t100];
489             // Compute quadrupole torques
490             tqx -= twoThirds * (
491                 mpole[t110] * phi[t101] + mpole[t020] * phi[t011] + mpole[t011] * phi[t002]
492                     - mpole[t101] * phi[t110] - mpole[t011] * phi[t020] - mpole[t002] * phi[t011]);
493             tqy -= twoThirds * (
494                 mpole[t101] * phi[t200] + mpole[t011] * phi[t110] + mpole[t002] * phi[t101]
495                     - mpole[t200] * phi[t101] - mpole[t110] * phi[t011] - mpole[t101] * phi[t002]);
496             tqz -= twoThirds * (
497                 mpole[t200] * phi[t110] + mpole[t110] * phi[t020] + mpole[t101] * phi[t011]
498                     - mpole[t110] * phi[t200] - mpole[t020] * phi[t110] - mpole[t011] * phi[t101]);
499             if (gradient) {
500               double factor = permanentScale * electric;
501               grad.add(threadID, i, factor * dfx, factor * dfy, factor * dfz);
502               torque.add(threadID, i, factor * tqx, factor * tqy, factor * tqz);
503             }
504             if (lambdaTerm) {
505               dUdL += dEdLSign * dlPowPerm * e;
506               d2UdL2 += dEdLSign * d2lPowPerm * e;
507               double factor = dEdLSign * dlPowPerm * electric;
508               lambdaGrad.add(threadID, i, factor * dfx, factor * dfy, factor * dfz);
509               lambdaTorque.add(threadID, i, factor * tqx, factor * tqy, factor * tqz);
510             }
511           }
512         }
513       }
514 
515       if (lambdaTerm) {
516         shareddEdLambda.addAndGet(0.5 * dUdL * electric);
517         sharedd2EdLambda2.addAndGet(0.5 * d2UdL2 * electric);
518       }
519     }
520 
521     @Override
522     public IntegerSchedule schedule() {
523       return IntegerSchedule.fixed();
524     }
525 
526     @Override
527     public void start() {
528       totalCharge = 0.0;
529       eSelf = 0.0;
530       eRecip = 0.0;
531       threadID = getThreadIndex();
532     }
533   }
534 
535   private class InducedDipoleReciprocalEnergyLoop extends IntegerForLoop {
536 
537     private final double[] sfPhi = new double[tensorCount];
538     private final double[] sPhi = new double[tensorCount];
539     private double eSelf;
540     private double eRecip;
541     private int threadID;
542 
543     @Override
544     public void finish() {
545       inducedDipoleSelfEnergy.addAndGet(polarizationScale * eSelf);
546       inducedDipoleRecipEnergy.addAndGet(polarizationScale * eRecip);
547     }
548 
549     @Override
550     public void run(int lb, int ub) throws Exception {
551       // Induced dipole self energy and gradient.
552       for (int i = lb; i <= ub; i++) {
553         if (use[i]) {
554           final double[] indi = ind[i];
555           final double[] multipolei = multipole[i];
556           final double dix = multipolei[t100];
557           final double diy = multipolei[t010];
558           final double diz = multipolei[t001];
559           final double dii = indi[0] * dix + indi[1] * diy + indi[2] * diz;
560           eSelf += aewald3 * dii;
561         }
562       }
563       if (lambdaTerm) {
564         shareddEdLambda.addAndGet(dEdLSign * dlPowPol * eSelf);
565         sharedd2EdLambda2.addAndGet(dEdLSign * d2lPowPol * eSelf);
566       }
567       if (gradient || esvTerm) {
568         for (int i = lb; i <= ub; i++) {
569           if (use[i]) {
570             final double[] indi = ind[i];
571             final double[] indpi = indCR[i];
572             final double[] multipolei = multipole[i];
573             final double dix = multipolei[t100];
574             final double diy = multipolei[t010];
575             final double diz = multipolei[t001];
576             final double uix = 0.5 * (indi[0] + indpi[0]);
577             final double uiy = 0.5 * (indi[1] + indpi[1]);
578             final double uiz = 0.5 * (indi[2] + indpi[2]);
579             if (gradient) {
580               final double tix = aewald4 * (diy * uiz - diz * uiy);
581               final double tiy = aewald4 * (diz * uix - dix * uiz);
582               final double tiz = aewald4 * (dix * uiy - diy * uix);
583               torque.add(threadID, i, polarizationScale * tix, polarizationScale * tiy, polarizationScale * tiz);
584               if (lambdaTerm) {
585                 double factor = dEdLSign * dlPowPol;
586                 lambdaTorque.add(threadID, i, factor * tix, factor * tiy, factor * tiz);
587               }
588             }
589             if (esvTerm && extendedSystem.isTitrating(i)) {
590               double[] titrDot = titrationMultipole[0][i];
591               double indiDotMdot = uix * titrDot[t100] + uiy * titrDot[t010] + uiz * titrDot[t001];
592               double dIndSelfTitrdLi = polarizationScale * 2.0 * aewald3 * indiDotMdot;
593               double dIndSelfTautdLi = 0.0;
594               if (extendedSystem.isTautomerizing(i)) {
595                 double[] tautDot = tautomerMultipole[0][i];
596                 indiDotMdot = uix * tautDot[t100] + uiy * tautDot[t010] + uiz * tautDot[t001];
597                 dIndSelfTautdLi = polarizationScale * 2.0 * aewald3 * indiDotMdot;
598               }
599               extendedSystem.addIndElecDeriv(i, dIndSelfTitrdLi, dIndSelfTautdLi);
600             }
601           }
602         }
603       }
604 
605       // Induced dipole reciprocal space energy and gradient.
606       for (int i = lb; i <= ub; i++) {
607         double[] fracInd = new double[3];
608         double[] fracIndCR = new double[3];
609         if (use[i]) {
610           final double[] fPhi = fracMultipolePhi[i];
611           reciprocalSpace.cartToFracInducedDipole(inducedDipole[0][i], inducedDipoleCR[0][i], fracInd, fracIndCR);
612           final double indx = fracInd[0];
613           final double indy = fracInd[1];
614           final double indz = fracInd[2];
615           eRecip += indx * fPhi[t100] + indy * fPhi[t010] + indz * fPhi[t001];
616           if (gradient || esvTerm) {
617             final double[] iPhi = cartesianDipolePhi[i];
618             final double[] iCRPhi = cartesianDipolePhiCR[i];
619             final double[] fiPhi = fracInducedDipolePhi[i];
620             final double[] fiCRPhi = fracInducedDipolePhiCR[i];
621             final double[] mpolei = multipole[i];
622             final double[] fmpolei = fracMultipole[i];
623             final double inpx = fracIndCR[0];
624             final double inpy = fracIndCR[1];
625             final double inpz = fracIndCR[2];
626             final double insx = indx + inpx;
627             final double insy = indy + inpy;
628             final double insz = indz + inpz;
629             for (int t = 0; t < tensorCount; t++) {
630               sPhi[t] = 0.5 * (iPhi[t] + iCRPhi[t]);
631               sfPhi[t] = fiPhi[t] + fiCRPhi[t];
632             }
633             double gx = insx * fPhi[t200] + insy * fPhi[t110] + insz * fPhi[t101];
634             double gy = insx * fPhi[t110] + insy * fPhi[t020] + insz * fPhi[t011];
635             double gz = insx * fPhi[t101] + insy * fPhi[t011] + insz * fPhi[t002];
636             if (polarization == Polarization.MUTUAL) {
637               gx += indx * fiCRPhi[t200] + inpx * fiPhi[t200] + indy * fiCRPhi[t110]
638                   + inpy * fiPhi[t110] + indz * fiCRPhi[t101] + inpz * fiPhi[t101];
639               gy += indx * fiCRPhi[t110] + inpx * fiPhi[t110] + indy * fiCRPhi[t020]
640                   + inpy * fiPhi[t020] + indz * fiCRPhi[t011] + inpz * fiPhi[t011];
641               gz += indx * fiCRPhi[t101] + inpx * fiPhi[t101] + indy * fiCRPhi[t011]
642                   + inpy * fiPhi[t011] + indz * fiCRPhi[t002] + inpz * fiPhi[t002];
643             }
644             gx += fmpolei[t000] * sfPhi[t100] + fmpolei[t100] * sfPhi[t200]
645                 + fmpolei[t010] * sfPhi[t110] + fmpolei[t001] * sfPhi[t101]
646                 + fmpolei[t200] * sfPhi[t300] + fmpolei[t020] * sfPhi[t120]
647                 + fmpolei[t002] * sfPhi[t102] + fmpolei[t110] * sfPhi[t210]
648                 + fmpolei[t101] * sfPhi[t201] + fmpolei[t011] * sfPhi[t111];
649             gy += fmpolei[t000] * sfPhi[t010] + fmpolei[t100] * sfPhi[t110]
650                 + fmpolei[t010] * sfPhi[t020] + fmpolei[t001] * sfPhi[t011]
651                 + fmpolei[t200] * sfPhi[t210] + fmpolei[t020] * sfPhi[t030]
652                 + fmpolei[t002] * sfPhi[t012] + fmpolei[t110] * sfPhi[t120]
653                 + fmpolei[t101] * sfPhi[t111] + fmpolei[t011] * sfPhi[t021];
654             gz += fmpolei[t000] * sfPhi[t001] + fmpolei[t100] * sfPhi[t101]
655                 + fmpolei[t010] * sfPhi[t011] + fmpolei[t001] * sfPhi[t002]
656                 + fmpolei[t200] * sfPhi[t201] + fmpolei[t020] * sfPhi[t021]
657                 + fmpolei[t002] * sfPhi[t003] + fmpolei[t110] * sfPhi[t111]
658                 + fmpolei[t101] * sfPhi[t102] + fmpolei[t011] * sfPhi[t012];
659             gx *= nfftX;
660             gy *= nfftY;
661             gz *= nfftZ;
662             double[][] recip = crystal.getUnitCell().A;
663             double dfx = recip[0][0] * gx + recip[0][1] * gy + recip[0][2] * gz;
664             double dfy = recip[1][0] * gx + recip[1][1] * gy + recip[1][2] * gz;
665             double dfz = recip[2][0] * gx + recip[2][1] * gy + recip[2][2] * gz;
666             dfx *= 0.5 * electric;
667             dfy *= 0.5 * electric;
668             dfz *= 0.5 * electric;
669             // Compute dipole torques
670             double tqx = -mpolei[t010] * sPhi[t001] + mpolei[t001] * sPhi[t010];
671             double tqy = -mpolei[t001] * sPhi[t100] + mpolei[t100] * sPhi[t001];
672             double tqz = -mpolei[t100] * sPhi[t010] + mpolei[t010] * sPhi[t100];
673             // Compute quadrupole torques
674             tqx -= twoThirds * (
675                 mpolei[t110] * sPhi[t101] + mpolei[t020] * sPhi[t011] + mpolei[t011] * sPhi[t002]
676                     - mpolei[t101] * sPhi[t110] - mpolei[t011] * sPhi[t020]
677                     - mpolei[t002] * sPhi[t011]);
678             tqy -= twoThirds * (
679                 mpolei[t101] * sPhi[t200] + mpolei[t011] * sPhi[t110] + mpolei[t002] * sPhi[t101]
680                     - mpolei[t200] * sPhi[t101] - mpolei[t110] * sPhi[t011]
681                     - mpolei[t101] * sPhi[t002]);
682             tqz -= twoThirds * (
683                 mpolei[t200] * sPhi[t110] + mpolei[t110] * sPhi[t020] + mpolei[t101] * sPhi[t011]
684                     - mpolei[t110] * sPhi[t200] - mpolei[t020] * sPhi[t110]
685                     - mpolei[t011] * sPhi[t101]);
686             tqx *= electric;
687             tqy *= electric;
688             tqz *= electric;
689             grad.add(threadID, i, polarizationScale * dfx, polarizationScale * dfy,
690                 polarizationScale * dfz);
691             torque.add(threadID, i, polarizationScale * tqx, polarizationScale * tqy,
692                 polarizationScale * tqz);
693             if (lambdaTerm) {
694               double factor = dEdLSign * dlPowPol;
695               lambdaGrad.add(threadID, i, factor * dfx, factor * dfy, factor * dfz);
696               lambdaTorque.add(threadID, i, factor * tqx, factor * tqy, factor * tqz);
697             }
698             if (esvTerm && extendedSystem.isTitrating(i)) {
699               final double[] mTitrDot = titrationMultipole[0][i];
700               double eTitrDot = mTitrDot[t000] * sPhi[t000]
701                   + mTitrDot[t100] * sPhi[t100] + mTitrDot[t010] * sPhi[t010] + mTitrDot[t001] * sPhi[t001]
702                   + oneThird * (mTitrDot[t200] * sPhi[t200] + mTitrDot[t020] * sPhi[t020] + mTitrDot[t002] * sPhi[t002]
703                   + 2.0 * (mTitrDot[t110] * sPhi[t110] + mTitrDot[t101] * sPhi[t101] + mTitrDot[t011] * sPhi[t011]));
704               double eTautDot = 0.0;
705               if (extendedSystem.isTautomerizing(i)) {
706                 final double[] mTautDot = tautomerMultipole[0][i];
707                 eTautDot = mTautDot[t000] * sPhi[t000]
708                     + mTautDot[t100] * sPhi[t100] + mTautDot[t010] * sPhi[t010] + mTautDot[t001] * sPhi[t001]
709                     + oneThird * (mTautDot[t200] * sPhi[t200] + mTautDot[t020] * sPhi[t020] + mTautDot[t002] * sPhi[t002]
710                     + 2.0 * (mTautDot[t110] * sPhi[t110] + mTautDot[t101] * sPhi[t101] + mTautDot[t011] * sPhi[t011]));
711               }
712               double factor = polarizationScale * electric;
713               extendedSystem.addIndElecDeriv(i, eTitrDot * factor, eTautDot * factor);
714             }
715           }
716         }
717       }
718       eRecip *= 0.5 * electric;
719       if (lambdaTerm) {
720         shareddEdLambda.addAndGet(dEdLSign * dlPowPol * eRecip);
721         sharedd2EdLambda2.addAndGet(dEdLSign * d2lPowPol * eRecip);
722       }
723     }
724 
725     @Override
726     public IntegerSchedule schedule() {
727       return IntegerSchedule.fixed();
728     }
729 
730     @Override
731     public void start() {
732       eSelf = 0.0;
733       eRecip = 0.0;
734       threadID = getThreadIndex();
735     }
736   }
737 }