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.ParallelSection;
44  import edu.rit.pj.ParallelTeam;
45  import edu.rit.pj.reduction.SharedInteger;
46  import edu.rit.util.Range;
47  import ffx.crystal.Crystal;
48  import ffx.crystal.SymOp;
49  import ffx.numerics.atomic.AtomicDoubleArray3D;
50  import ffx.potential.bonded.Atom;
51  import ffx.potential.nonbonded.MaskingInterface;
52  import ffx.potential.nonbonded.ReciprocalSpace;
53  import ffx.potential.parameters.ForceField;
54  
55  import java.util.List;
56  import java.util.logging.Level;
57  import java.util.logging.Logger;
58  
59  import static ffx.numerics.special.Erf.erfc;
60  import static ffx.potential.parameters.MultipoleType.t000;
61  import static ffx.potential.parameters.MultipoleType.t001;
62  import static ffx.potential.parameters.MultipoleType.t002;
63  import static ffx.potential.parameters.MultipoleType.t010;
64  import static ffx.potential.parameters.MultipoleType.t011;
65  import static ffx.potential.parameters.MultipoleType.t020;
66  import static ffx.potential.parameters.MultipoleType.t100;
67  import static ffx.potential.parameters.MultipoleType.t101;
68  import static ffx.potential.parameters.MultipoleType.t110;
69  import static ffx.potential.parameters.MultipoleType.t200;
70  import static java.util.Arrays.copyOf;
71  import static java.util.Arrays.fill;
72  import static org.apache.commons.math3.util.FastMath.exp;
73  import static org.apache.commons.math3.util.FastMath.min;
74  import static org.apache.commons.math3.util.FastMath.sqrt;
75  
76  /**
77   * Parallel computation of the permanent field.
78   *
79   * <p>This class can be executed by a ParallelTeam with exactly 2 threads.
80   *
81   * <p>The Real Space and Reciprocal Space Sections will be run concurrently, each with the number of
82   * threads defined by their respective ParallelTeam instances.
83   *
84   * @author Michael J. Schnieders
85   * @since 1.0
86   */
87  public class PermanentFieldRegion extends ParallelRegion implements MaskingInterface {
88  
89    private static final Logger logger = Logger.getLogger(PermanentFieldRegion.class.getName());
90    /**
91     * Constant applied to multipole interactions.
92     */
93    private static final double oneThird = 1.0 / 3.0;
94    /**
95     * Specify inter-molecular softcore.
96     */
97    private final boolean intermolecularSoftcore;
98    /**
99     * Specify intra-molecular softcore.
100    */
101   private final boolean intramolecularSoftcore;
102   /**
103    * Dimensions of [nsymm][nAtoms][3]
104    */
105   public double[][][] inducedDipole;
106   public double[][][] inducedDipoleCR;
107   /**
108    * Polarization groups.
109    */
110   protected int[][] ip11;
111   /**
112    * An ordered array of atoms in the system.
113    */
114   private Atom[] atoms;
115   /**
116    * Unit cell and spacegroup information.
117    */
118   private Crystal crystal;
119   /**
120    * Dimensions of [nsymm][xyz][nAtoms].
121    */
122   private double[][][] coordinates;
123   /**
124    * Dimensions of [nsymm][nAtoms][10]
125    */
126   private double[][][] globalMultipole;
127   /**
128    * Neighbor lists, including atoms beyond the real space cutoff. [nsymm][nAtoms][nAllNeighbors]
129    */
130   private int[][][] neighborLists;
131   /**
132    * Neighbor lists, without atoms beyond the preconditioner cutoff.
133    * [nSymm][nAtoms][nIncludedNeighbors]
134    */
135   private int[][][] preconditionerLists;
136   /**
137    * Number of neighboring atoms within the preconditioner cutoff. [nSymm][nAtoms]
138    */
139   private int[][] preconditionerCounts;
140   /**
141    * When computing the polarization energy at Lambda there are 3 pieces.
142    *
143    * <p>1.) Upol(1) = The polarization energy computed normally (ie. system with ligand).
144    *
145    * <p>2.) Uenv = The polarization energy of the system without the ligand.
146    *
147    * <p>3.) Uligand = The polarization energy of the ligand by itself.
148    *
149    * <p>Upol(L) = L*Upol(1) + (1-L)*(Uenv + Uligand)
150    *
151    * <p>Set the "use" array to true for all atoms for part 1. Set the "use" array to true for all
152    * atoms except the ligand for part 2. Set the "use" array to true only for the ligand atoms for
153    * part 3.
154    *
155    * <p>The "use" array can also be employed to turn off atoms for computing the electrostatic
156    * energy of sub-structures.
157    */
158   private boolean[] use;
159   /**
160    * Molecule number for each atom.
161    */
162   private int[] molecule;
163   private double[] ipdamp;
164   private double[] thole;
165   /**
166    * Masking of 1-2, 1-3 and 1-4 interactions.
167    */
168   private int[][] mask12;
169   private int[][] mask13;
170   private int[][] mask14;
171   /**
172    * The current LambdaMode of this PME instance (or OFF for no lambda dependence).
173    */
174   private LambdaMode lambdaMode = LambdaMode.OFF;
175   /**
176    * Reciprocal space instance.
177    */
178   private ReciprocalSpace reciprocalSpace;
179   private boolean reciprocalSpaceTerm;
180   private double off2;
181   private double preconditionerCutoff;
182   private double an0, an1, an2;
183   private double aewald;
184   /**
185    * Neighbor lists, without atoms beyond the real space cutoff. [nSymm][nAtoms][nIncludedNeighbors]
186    */
187   private int[][][] realSpaceLists;
188   /**
189    * Number of neighboring atoms within the real space cutoff. [nSymm][nAtoms]
190    */
191   private int[][] realSpaceCounts;
192   private Range[] realSpaceRanges;
193   private IntegerSchedule permanentSchedule;
194   /**
195    * Field array.
196    */
197   private AtomicDoubleArray3D field;
198   /**
199    * Chain rule field array.
200    */
201   private AtomicDoubleArray3D fieldCR;
202   private ScaleParameters scaleParameters;
203   private final PermanentRealSpaceFieldSection permanentRealSpaceFieldSection;
204   private final PermanentReciprocalSection permanentReciprocalSection;
205   /**
206    * Timing variables.
207    */
208   private PMETimings pmeTimings;
209 
210   public PermanentFieldRegion(ParallelTeam pt, ForceField forceField, boolean lambdaTerm) {
211     permanentRealSpaceFieldSection = new PermanentRealSpaceFieldSection(pt);
212     permanentReciprocalSection = new PermanentReciprocalSection();
213 
214     // Flag to indicate application of an intermolecular softcore potential.
215     if (lambdaTerm) {
216       intermolecularSoftcore = forceField.getBoolean("INTERMOLECULAR_SOFTCORE", false);
217       intramolecularSoftcore = forceField.getBoolean("INTRAMOLECULAR_SOFTCORE", false);
218     } else {
219       intermolecularSoftcore = false;
220       intramolecularSoftcore = false;
221     }
222   }
223 
224   /**
225    * Apply permanent field masking rules.
226    *
227    * @param i     The atom whose masking rules should be applied.
228    * @param is14  True if atom i and the current atom are 1-4 to each other.
229    * @param masks One or more masking arrays.
230    */
231   @Override
232   public void applyMask(int i, boolean[] is14, double[]... masks) {
233     if (ip11[i] != null) {
234       double[] energyMask = masks[0];
235       var m12 = mask12[i];
236       for (int value : m12) {
237         energyMask[value] = scaleParameters.p12scale;
238       }
239       var m13 = mask13[i];
240       for (int value : m13) {
241         energyMask[value] = scaleParameters.p13scale;
242       }
243       var m14 = mask14[i];
244       for (int value : m14) {
245         energyMask[value] = scaleParameters.p14scale;
246         for (int k : ip11[i]) {
247           if (k == value) {
248             energyMask[value] = scaleParameters.intra14Scale * scaleParameters.p14scale;
249             break;
250           }
251         }
252       }
253       // Apply group based polarization masking rule.
254       double[] inductionMask = masks[1];
255       for (int index : ip11[i]) {
256         inductionMask[index] = scaleParameters.d11scale;
257       }
258     }
259   }
260 
261   public void initTimings() {
262     permanentRealSpaceFieldSection.initTimings();
263   }
264 
265   public long getRealSpacePermTime() {
266     return permanentRealSpaceFieldSection.time;
267   }
268 
269   public long getInitTime(int threadId) {
270     return permanentRealSpaceFieldSection.permanentRealSpaceFieldRegion.initializationLoop[threadId].time;
271   }
272 
273   public long getPermTime(int threadId) {
274     return permanentRealSpaceFieldSection.permanentRealSpaceFieldRegion.permanentRealSpaceFieldLoop[threadId].time;
275   }
276 
277   public void init(
278       Atom[] atoms,
279       Crystal crystal,
280       double[][][] coordinates,
281       double[][][] globalMultipole,
282       double[][][] inducedDipole,
283       double[][][] inducedDipoleCR,
284       int[][][] neighborLists,
285       ScaleParameters scaleParameters,
286       boolean[] use,
287       int[] molecule,
288       double[] ipdamp,
289       double[] thole,
290       int[][] ip11,
291       int[][] mask12,
292       int[][] mask13,
293       int[][] mask14,
294       LambdaMode lambdaMode,
295       boolean reciprocalSpaceTerm,
296       ReciprocalSpace reciprocalSpace,
297       EwaldParameters ewaldParameters,
298       PCGSolver pcgSolver,
299       IntegerSchedule permanentSchedule,
300       RealSpaceNeighborParameters realSpaceNeighborParameters,
301       AtomicDoubleArray3D field,
302       AtomicDoubleArray3D fieldCR) {
303     this.atoms = atoms;
304     this.crystal = crystal;
305     this.coordinates = coordinates;
306     this.globalMultipole = globalMultipole;
307     this.inducedDipole = inducedDipole;
308     this.inducedDipoleCR = inducedDipoleCR;
309     this.neighborLists = neighborLists;
310     this.scaleParameters = scaleParameters;
311     this.use = use;
312     this.molecule = molecule;
313     this.ipdamp = ipdamp;
314     this.thole = thole;
315     this.ip11 = ip11;
316     this.mask12 = mask12;
317     this.mask13 = mask13;
318     this.mask14 = mask14;
319     this.lambdaMode = lambdaMode;
320     this.reciprocalSpaceTerm = reciprocalSpaceTerm;
321     this.reciprocalSpace = reciprocalSpace;
322     if (pcgSolver != null) {
323       this.preconditionerCutoff = pcgSolver.getPreconditionerCutoff();
324       this.preconditionerLists = pcgSolver.getPreconditionerLists();
325       this.preconditionerCounts = pcgSolver.getPreconditionerCounts();
326     }
327     this.aewald = ewaldParameters.aewald;
328     this.an0 = ewaldParameters.an0;
329     this.an1 = ewaldParameters.an1;
330     this.an2 = ewaldParameters.an2;
331     this.off2 = ewaldParameters.off2;
332     this.permanentSchedule = permanentSchedule;
333     this.realSpaceLists = realSpaceNeighborParameters.realSpaceLists;
334     this.realSpaceCounts = realSpaceNeighborParameters.realSpaceCounts;
335     this.realSpaceRanges = realSpaceNeighborParameters.realSpaceRanges;
336     this.field = field;
337     this.fieldCR = fieldCR;
338   }
339 
340   /**
341    * Remove permanent field masking rules.
342    *
343    * @param i     The atom whose masking rules should be removed.
344    * @param is14  True if atom i and the current atom are 1-4 to each other.
345    * @param masks One or more masking arrays.
346    */
347   @Override
348   public void removeMask(int i, boolean[] is14, double[]... masks) {
349     if (ip11[i] != null) {
350       double[] energyMask = masks[0];
351       var m12 = mask12[i];
352       for (int value : m12) {
353         energyMask[value] = 1.0;
354       }
355       var m13 = mask13[i];
356       for (int value : m13) {
357         energyMask[value] = 1.0;
358       }
359       var m14 = mask14[i];
360       for (int value : m14) {
361         energyMask[value] = 1.0;
362         for (int k : ip11[i]) {
363           if (k == value) {
364             energyMask[value] = 1.0;
365             break;
366           }
367         }
368       }
369       // Apply group based polarization masking rule.
370       double[] inductionMask = masks[1];
371       for (int index : ip11[i]) {
372         inductionMask[index] = 1.0;
373       }
374     }
375   }
376 
377   @Override
378   public void run() {
379     try {
380       execute(permanentRealSpaceFieldSection, permanentReciprocalSection);
381     } catch (RuntimeException e) {
382       String message = "Runtime exception computing the permanent multipole field.\n";
383       logger.log(Level.WARNING, message, e);
384       throw e;
385     } catch (Exception e) {
386       String message = "Fatal exception computing the permanent multipole field.\n";
387       logger.log(Level.SEVERE, message, e);
388     }
389   }
390 
391   /**
392    * Computes the Permanent Multipole Real Space Field.
393    */
394   private class PermanentRealSpaceFieldSection extends ParallelSection {
395 
396     private final PermanentRealSpaceFieldRegion permanentRealSpaceFieldRegion;
397     private final ParallelTeam parallelTeam;
398     protected long time;
399 
400     PermanentRealSpaceFieldSection(ParallelTeam pt) {
401       this.parallelTeam = pt;
402       int nt = pt.getThreadCount();
403       permanentRealSpaceFieldRegion = new PermanentRealSpaceFieldRegion(nt);
404     }
405 
406     public void initTimings() {
407       time = 0;
408       permanentRealSpaceFieldRegion.initTimings();
409     }
410 
411     @Override
412     public void run() {
413       try {
414         time -= System.nanoTime();
415         parallelTeam.execute(permanentRealSpaceFieldRegion);
416         time += System.nanoTime();
417       } catch (RuntimeException e) {
418         String message = "Fatal exception computing the real space field.\n";
419         logger.log(Level.WARNING, message, e);
420       } catch (Exception e) {
421         String message = "Fatal exception computing the real space field.\n";
422         logger.log(Level.SEVERE, message, e);
423       }
424     }
425   }
426 
427   /**
428    * Compute the permanent multipole reciprocal space contribution to the electric potential, field,
429    * etc. using the number of threads specified by the ParallelTeam used to construct the
430    * ReciprocalSpace instance.
431    */
432   private class PermanentReciprocalSection extends ParallelSection {
433 
434     @Override
435     public void run() {
436       if (reciprocalSpaceTerm && aewald > 0.0) {
437         reciprocalSpace.performConvolution();
438       }
439     }
440   }
441 
442   private class PermanentRealSpaceFieldRegion extends ParallelRegion {
443 
444     private final InitializationLoop[] initializationLoop;
445     private final PermanentRealSpaceFieldLoop[] permanentRealSpaceFieldLoop;
446     private final SharedInteger sharedCount;
447     private final int threadCount;
448 
449     PermanentRealSpaceFieldRegion(int nt) {
450       threadCount = nt;
451       initializationLoop = new InitializationLoop[threadCount];
452       permanentRealSpaceFieldLoop = new PermanentRealSpaceFieldLoop[threadCount];
453       sharedCount = new SharedInteger();
454       for (int i = 0; i < threadCount; i++) {
455         initializationLoop[i] = new InitializationLoop();
456         permanentRealSpaceFieldLoop[i] = new PermanentRealSpaceFieldLoop();
457       }
458     }
459 
460     public void initTimings() {
461       for (int i = 0; i < threadCount; i++) {
462         permanentRealSpaceFieldLoop[i].time = 0;
463         initializationLoop[i].time = 0;
464       }
465     }
466 
467     @Override
468     public void finish() {
469       if (realSpaceRanges == null) {
470         logger.severe(" RealSpaceRange array is null");
471       }
472 
473       int nAtoms = atoms.length;
474 
475       // Load balancing.
476       int id = 0;
477       int goal = sharedCount.get() / threadCount;
478       int num = 0;
479       int start = 0;
480 
481       for (int i = 0; i < nAtoms; i++) {
482         List<SymOp> symOps = crystal.spaceGroup.symOps;
483         int nSymm = symOps.size();
484         for (int iSymm = 0; iSymm < nSymm; iSymm++) {
485           num += realSpaceCounts[iSymm][i];
486         }
487         if (num >= goal) {
488           // Last thread gets the remaining atoms in its range.
489           if (id == threadCount - 1) {
490             realSpaceRanges[id] = new Range(start, nAtoms - 1);
491             break;
492           }
493 
494           realSpaceRanges[id] = new Range(start, i);
495 
496           // Reset the count.
497           num = 0;
498 
499           // Next thread.
500           id++;
501 
502           // Next range starts at i+1.
503           start = i + 1;
504 
505           // Out of atoms. Threads remaining get a null range.
506           if (start == nAtoms) {
507             for (int j = id; j < threadCount; j++) {
508               realSpaceRanges[j] = null;
509             }
510             break;
511           }
512         } else if (i == nAtoms - 1) {
513 
514           // Last atom without reaching goal for current thread.
515           realSpaceRanges[id] = new Range(start, nAtoms - 1);
516           for (int j = id + 1; j < threadCount; j++) {
517             realSpaceRanges[j] = null;
518           }
519         }
520       }
521     }
522 
523     @Override
524     public void run() {
525       int threadIndex = getThreadIndex();
526       try {
527         int nAtoms = atoms.length;
528         execute(0, nAtoms - 1, initializationLoop[threadIndex]);
529         execute(0, nAtoms - 1, permanentRealSpaceFieldLoop[threadIndex]);
530       } catch (RuntimeException e) {
531         String message = "Runtime exception computing the real space field.\n";
532         logger.log(Level.SEVERE, message, e);
533       } catch (Exception e) {
534         String message =
535             "Fatal exception computing the real space field in thread " + getThreadIndex() + "\n";
536         logger.log(Level.SEVERE, message, e);
537       }
538     }
539 
540     @Override
541     public void start() {
542       sharedCount.set(0);
543     }
544 
545     private class InitializationLoop extends IntegerForLoop {
546 
547       protected long time;
548 
549       @Override
550       public void start() {
551         time -= System.nanoTime();
552       }
553 
554       @Override
555       public void finish() {
556         time += System.nanoTime();
557       }
558 
559       @Override
560       public void run(int lb, int ub) {
561         // Initialize the induced dipole arrays.
562         List<SymOp> symOps = crystal.spaceGroup.symOps;
563         int nSymm = symOps.size();
564         for (int iSymm = 0; iSymm < nSymm; iSymm++) {
565           double[][] ind0 = inducedDipole[0];
566           double[][] indCR0 = inducedDipoleCR[0];
567           for (int i = lb; i <= ub; i++) {
568             double[] ind = ind0[i];
569             double[] indCR = indCR0[i];
570             ind[0] = 0.0;
571             ind[1] = 0.0;
572             ind[2] = 0.0;
573             indCR[0] = 0.0;
574             indCR[1] = 0.0;
575             indCR[2] = 0.0;
576           }
577         }
578       }
579 
580       @Override
581       public IntegerSchedule schedule() {
582         return IntegerSchedule.fixed();
583       }
584 
585     }
586 
587     private class PermanentRealSpaceFieldLoop extends IntegerForLoop {
588 
589       protected long time;
590       private final double[] dx_local;
591       private final double[][] transOp;
592       private int threadID;
593       private double[] inductionMaskLocal;
594       private double[] energyMaskLocal;
595       private int count;
596 
597       PermanentRealSpaceFieldLoop() {
598         super();
599         dx_local = new double[3];
600         transOp = new double[3][3];
601       }
602 
603       @Override
604       public void finish() {
605         sharedCount.addAndGet(count);
606         time += System.nanoTime();
607       }
608 
609       @Override
610       public void run(int lb, int ub) {
611         int[][] lists = neighborLists[0];
612         int[][] ewalds = realSpaceLists[0];
613         int[] counts = realSpaceCounts[0];
614         int[][] preLists = preconditionerLists[0];
615         int[] preCounts = preconditionerCounts[0];
616         final double[] x = coordinates[0][0];
617         final double[] y = coordinates[0][1];
618         final double[] z = coordinates[0][2];
619         final double[][] mpole = globalMultipole[0];
620         // Loop over atom chunk.
621         for (int i = lb; i <= ub; i++) {
622           if (!use[i]) {
623             continue;
624           }
625 
626           // Zero out accumulation variables for the field at atom i.
627           double fix = 0.0;
628           double fiy = 0.0;
629           double fiz = 0.0;
630           double fixCR = 0.0;
631           double fiyCR = 0.0;
632           double fizCR = 0.0;
633 
634           final int moleculei = molecule[i];
635           final double pdi = ipdamp[i];
636           final double pti = thole[i];
637           final double xi = x[i];
638           final double yi = y[i];
639           final double zi = z[i];
640           final double[] globalMultipolei = mpole[i];
641           final double ci = globalMultipolei[0];
642           final double dix = globalMultipolei[t100];
643           final double diy = globalMultipolei[t010];
644           final double diz = globalMultipolei[t001];
645           final double qixx = globalMultipolei[t200] * oneThird;
646           final double qiyy = globalMultipolei[t020] * oneThird;
647           final double qizz = globalMultipolei[t002] * oneThird;
648           final double qixy = globalMultipolei[t110] * oneThird;
649           final double qixz = globalMultipolei[t101] * oneThird;
650           final double qiyz = globalMultipolei[t011] * oneThird;
651 
652           // Apply field masking rules.
653           applyMask(i, null, energyMaskLocal, inductionMaskLocal);
654 
655           // Loop over the neighbor list.
656           final int[] list = lists[i];
657           counts[i] = 0;
658           preCounts[i] = 0;
659           int[] ewald = ewalds[i];
660           int[] preList = preLists[i];
661           for (int k : list) {
662             if (!use[k]) {
663               continue;
664             }
665             if (lambdaMode == LambdaMode.VAPOR) {
666               boolean sameMolecule = (moleculei == molecule[k]);
667               if ((intermolecularSoftcore && !sameMolecule)
668                   || (intramolecularSoftcore && sameMolecule)) {
669                 continue;
670               }
671             }
672             final double xk = x[k];
673             final double yk = y[k];
674             final double zk = z[k];
675             dx_local[0] = xk - xi;
676             dx_local[1] = yk - yi;
677             dx_local[2] = zk - zi;
678             final double r2 = crystal.image(dx_local);
679             if (r2 <= off2) {
680               count++;
681               // Store a short neighbor list for the SCF.
682               if (ewald.length <= counts[i]) {
683                 int len = ewald.length;
684                 ewalds[i] = copyOf(ewald, len + 10);
685                 ewald = ewalds[i];
686               }
687               ewald[counts[i]++] = k;
688               final double xr = dx_local[0];
689               final double yr = dx_local[1];
690               final double zr = dx_local[2];
691               final double pdk = ipdamp[k];
692               final double ptk = thole[k];
693               final double[] globalMultipolek = mpole[k];
694               final double ck = globalMultipolek[t000];
695               final double dkx = globalMultipolek[t100];
696               final double dky = globalMultipolek[t010];
697               final double dkz = globalMultipolek[t001];
698               final double qkxx = globalMultipolek[t200] * oneThird;
699               final double qkyy = globalMultipolek[t020] * oneThird;
700               final double qkzz = globalMultipolek[t002] * oneThird;
701               final double qkxy = globalMultipolek[t110] * oneThird;
702               final double qkxz = globalMultipolek[t101] * oneThird;
703               final double qkyz = globalMultipolek[t011] * oneThird;
704               double r = sqrt(r2);
705 
706               // Store a short neighbor list for the SCF pre-conditioner.
707               if (r < preconditionerCutoff) {
708                 if (preList.length <= preCounts[i]) {
709                   int len = preList.length;
710                   preLists[i] = copyOf(preList, len + 10);
711                   preList = preLists[i];
712                 }
713                 preList[preCounts[i]++] = k;
714               }
715 
716               // Calculate the error function damping terms.
717               final double ralpha = aewald * r;
718               final double exp2a = exp(-ralpha * ralpha);
719               final double rr1 = 1.0 / r;
720               final double rr2 = rr1 * rr1;
721               final double bn0 = erfc(ralpha) * rr1;
722               final double bn1 = (bn0 + an0 * exp2a) * rr2;
723               final double bn2 = (3.0 * bn1 + an1 * exp2a) * rr2;
724               final double bn3 = (5.0 * bn2 + an2 * exp2a) * rr2;
725 
726               // Compute the error function scaled and unscaled terms.
727               double scale3 = 1.0;
728               double scale5 = 1.0;
729               double scale7 = 1.0;
730               double damp = pdi * pdk;
731               final double pgamma = min(pti, ptk);
732               final double rdamp = r * damp;
733               damp = -pgamma * rdamp * rdamp * rdamp;
734               if (damp > -50.0) {
735                 double expdamp = exp(damp);
736                 scale3 = 1.0 - expdamp;
737                 scale5 = 1.0 - expdamp * (1.0 - damp);
738                 scale7 = 1.0 - expdamp * (1.0 - damp + 0.6 * damp * damp);
739               }
740               final double scale = inductionMaskLocal[k];
741               final double scalep = energyMaskLocal[k];
742               // Thole damping multiplied by the group-based mask.
743               final double dsc3 = scale3 * scale;
744               final double dsc5 = scale5 * scale;
745               final double dsc7 = scale7 * scale;
746               // Thole damping multiplied by the energy mask.
747               final double psc3 = scale3 * scalep;
748               final double psc5 = scale5 * scalep;
749               final double psc7 = scale7 * scalep;
750               final double rr3 = rr1 * rr2;
751               final double rr5 = 3.0 * rr3 * rr2;
752               final double rr7 = 5.0 * rr5 * rr2;
753               // 1.0 minus induction masks and Thole damping.
754               final double drr3 = (1.0 - dsc3) * rr3;
755               final double drr5 = (1.0 - dsc5) * rr5;
756               final double drr7 = (1.0 - dsc7) * rr7;
757               // 1.0 minus energy masks and Thole damping.
758               final double prr3 = (1.0 - psc3) * rr3;
759               final double prr5 = (1.0 - psc5) * rr5;
760               final double prr7 = (1.0 - psc7) * rr7;
761               final double dir = dix * xr + diy * yr + diz * zr;
762               final double qix = 2.0 * (qixx * xr + qixy * yr + qixz * zr);
763               final double qiy = 2.0 * (qixy * xr + qiyy * yr + qiyz * zr);
764               final double qiz = 2.0 * (qixz * xr + qiyz * yr + qizz * zr);
765               final double qir = (qix * xr + qiy * yr + qiz * zr) * 0.5;
766               // Ewald field for atom k (no masking).
767               final double bn123i = bn1 * ci + bn2 * dir + bn3 * qir;
768               final double fkmx = xr * bn123i - bn1 * dix - bn2 * qix;
769               final double fkmy = yr * bn123i - bn1 * diy - bn2 * qiy;
770               final double fkmz = zr * bn123i - bn1 * diz - bn2 * qiz;
771               // Correct Ewald field for over-counted induction interactions.
772               final double ddr357i = drr3 * ci + drr5 * dir + drr7 * qir;
773               final double fkdx = xr * ddr357i - drr3 * dix - drr5 * qix;
774               final double fkdy = yr * ddr357i - drr3 * diy - drr5 * qiy;
775               final double fkdz = zr * ddr357i - drr3 * diz - drr5 * qiz;
776               field.add(threadID, k, fkmx - fkdx, fkmy - fkdy, fkmz - fkdz);
777               // Correct Ewald field for over-counted energy interactions.
778               final double prr357i = prr3 * ci + prr5 * dir + prr7 * qir;
779               final double fkpx = xr * prr357i - prr3 * dix - prr5 * qix;
780               final double fkpy = yr * prr357i - prr3 * diy - prr5 * qiy;
781               final double fkpz = zr * prr357i - prr3 * diz - prr5 * qiz;
782               fieldCR.add(threadID, k, fkmx - fkpx, fkmy - fkpy, fkmz - fkpz);
783               final double dkr = dkx * xr + dky * yr + dkz * zr;
784               final double qkx = 2.0 * (qkxx * xr + qkxy * yr + qkxz * zr);
785               final double qky = 2.0 * (qkxy * xr + qkyy * yr + qkyz * zr);
786               final double qkz = 2.0 * (qkxz * xr + qkyz * yr + qkzz * zr);
787               final double qkr = (qkx * xr + qky * yr + qkz * zr) * 0.5;
788               final double bn123k = bn1 * ck - bn2 * dkr + bn3 * qkr;
789               // Ewald field for atom i (no masking).
790               final double fimx = -xr * bn123k - bn1 * dkx + bn2 * qkx;
791               final double fimy = -yr * bn123k - bn1 * dky + bn2 * qky;
792               final double fimz = -zr * bn123k - bn1 * dkz + bn2 * qkz;
793               // Correct Ewald field for over-counted induction interactions.
794               final double drr357k = drr3 * ck - drr5 * dkr + drr7 * qkr;
795               final double fidx = -xr * drr357k - drr3 * dkx + drr5 * qkx;
796               final double fidy = -yr * drr357k - drr3 * dky + drr5 * qky;
797               final double fidz = -zr * drr357k - drr3 * dkz + drr5 * qkz;
798               fix += fimx - fidx;
799               fiy += fimy - fidy;
800               fiz += fimz - fidz;
801               // Correct Ewald field for over-counted energy interactions.
802               final double prr357k = prr3 * ck - prr5 * dkr + prr7 * qkr;
803               final double fipx = -xr * prr357k - prr3 * dkx + prr5 * qkx;
804               final double fipy = -yr * prr357k - prr3 * dky + prr5 * qky;
805               final double fipz = -zr * prr357k - prr3 * dkz + prr5 * qkz;
806               fixCR += fimx - fipx;
807               fiyCR += fimy - fipy;
808               fizCR += fimz - fipz;
809             }
810           }
811           // Add in field contributions at Atom i.
812           field.add(threadID, i, fix, fiy, fiz);
813           fieldCR.add(threadID, i, fixCR, fiyCR, fizCR);
814           // Remove field masking rules.
815           removeMask(i, null, energyMaskLocal, inductionMaskLocal);
816         }
817 
818         // Loop over symmetry mates.
819         List<SymOp> symOps = crystal.spaceGroup.symOps;
820         int nSymm = symOps.size();
821         for (int iSymm = 1; iSymm < nSymm; iSymm++) {
822           SymOp symOp = crystal.spaceGroup.getSymOp(iSymm);
823           crystal.getTransformationOperator(symOp, transOp);
824           lists = neighborLists[iSymm];
825           ewalds = realSpaceLists[iSymm];
826           counts = realSpaceCounts[iSymm];
827           preLists = preconditionerLists[iSymm];
828           preCounts = preconditionerCounts[iSymm];
829           double[] xs = coordinates[iSymm][0];
830           double[] ys = coordinates[iSymm][1];
831           double[] zs = coordinates[iSymm][2];
832           double[][] mpoles = globalMultipole[iSymm];
833 
834           // Loop over atoms in a chunk of the asymmetric unit.
835           for (int i = lb; i <= ub; i++) {
836             if (!use[i]) {
837               continue;
838             }
839             // Zero out accumulation variables for the field at atom i.
840             double fix = 0.0;
841             double fiy = 0.0;
842             double fiz = 0.0;
843             final double pdi = ipdamp[i];
844             final double pti = thole[i];
845             final double[] multipolei = mpole[i];
846             final double ci = multipolei[t000];
847             final double dix = multipolei[t100];
848             final double diy = multipolei[t010];
849             final double diz = multipolei[t001];
850             final double qixx = multipolei[t200] * oneThird;
851             final double qiyy = multipolei[t020] * oneThird;
852             final double qizz = multipolei[t002] * oneThird;
853             final double qixy = multipolei[t110] * oneThird;
854             final double qixz = multipolei[t101] * oneThird;
855             final double qiyz = multipolei[t011] * oneThird;
856             final double xi = x[i];
857             final double yi = y[i];
858             final double zi = z[i];
859 
860             // Loop over the neighbor list.
861             final int[] list = lists[i];
862             counts[i] = 0;
863             preCounts[i] = 0;
864             int[] ewald = ewalds[i];
865             int[] preList = preLists[i];
866             for (int k : list) {
867               if (!use[k]) {
868                 continue;
869               }
870               final double xk = xs[k];
871               final double yk = ys[k];
872               final double zk = zs[k];
873               dx_local[0] = xk - xi;
874               dx_local[1] = yk - yi;
875               dx_local[2] = zk - zi;
876               final double r2 = crystal.image(dx_local);
877               if (r2 <= off2) {
878                 count++;
879                 // Store a short neighbor list for the SCF.
880                 if (ewald.length <= counts[i]) {
881                   int len = ewald.length;
882                   ewalds[i] = copyOf(ewald, len + 10);
883                   ewald = ewalds[i];
884                 }
885                 ewald[counts[i]++] = k;
886                 double selfScale = 1.0;
887                 if (i == k) {
888                   selfScale = 0.5;
889                 }
890                 final double xr = dx_local[0];
891                 final double yr = dx_local[1];
892                 final double zr = dx_local[2];
893                 final double pdk = ipdamp[k];
894                 final double ptk = thole[k];
895                 final double[] multipolek = mpoles[k];
896                 final double ck = multipolek[t000];
897                 final double dkx = multipolek[t100];
898                 final double dky = multipolek[t010];
899                 final double dkz = multipolek[t001];
900                 final double qkxx = multipolek[t200] * oneThird;
901                 final double qkyy = multipolek[t020] * oneThird;
902                 final double qkzz = multipolek[t002] * oneThird;
903                 final double qkxy = multipolek[t110] * oneThird;
904                 final double qkxz = multipolek[t101] * oneThird;
905                 final double qkyz = multipolek[t011] * oneThird;
906                 final double r = sqrt(r2);
907                 if (r < preconditionerCutoff) {
908                   if (preList.length <= preCounts[i]) {
909                     int len = preList.length;
910                     preLists[i] = copyOf(preList, len + 10);
911                     preList = preLists[i];
912                   }
913                   preList[preCounts[i]++] = k;
914                 }
915 
916                 // Calculate the error function damping terms.
917                 final double ralpha = aewald * r;
918                 final double exp2a = exp(-ralpha * ralpha);
919                 final double rr1 = 1.0 / r;
920                 final double rr2 = rr1 * rr1;
921                 final double bn0 = erfc(ralpha) * rr1;
922                 final double bn1 = (bn0 + an0 * exp2a) * rr2;
923                 final double bn2 = (3.0 * bn1 + an1 * exp2a) * rr2;
924                 final double bn3 = (5.0 * bn2 + an2 * exp2a) * rr2;
925 
926                 // Compute the error function scaled and unscaled terms.
927                 double scale3 = 1.0;
928                 double scale5 = 1.0;
929                 double scale7 = 1.0;
930                 double damp = pdi * pdk;
931                 final double pgamma = min(pti, ptk);
932                 final double rdamp = r * damp;
933                 damp = -pgamma * rdamp * rdamp * rdamp;
934                 if (damp > -50.0) {
935                   double expdamp = exp(damp);
936                   scale3 = 1.0 - expdamp;
937                   scale5 = 1.0 - expdamp * (1.0 - damp);
938                   scale7 = 1.0 - expdamp * (1.0 - damp + 0.6 * damp * damp);
939                 }
940 
941                 final double dsc3 = scale3;
942                 final double dsc5 = scale5;
943                 final double dsc7 = scale7;
944                 final double rr3 = rr1 * rr2;
945                 final double rr5 = 3.0 * rr3 * rr2;
946                 final double rr7 = 5.0 * rr5 * rr2;
947                 final double drr3 = (1.0 - dsc3) * rr3;
948                 final double drr5 = (1.0 - dsc5) * rr5;
949                 final double drr7 = (1.0 - dsc7) * rr7;
950 
951                 final double dkr = dkx * xr + dky * yr + dkz * zr;
952                 final double qkx = 2.0 * (qkxx * xr + qkxy * yr + qkxz * zr);
953                 final double qky = 2.0 * (qkxy * xr + qkyy * yr + qkyz * zr);
954                 final double qkz = 2.0 * (qkxz * xr + qkyz * yr + qkzz * zr);
955                 final double qkr = (qkx * xr + qky * yr + qkz * zr) * 0.5;
956                 final double bn123k = bn1 * ck - bn2 * dkr + bn3 * qkr;
957                 final double drr357k = drr3 * ck - drr5 * dkr + drr7 * qkr;
958                 final double fimx = -xr * bn123k - bn1 * dkx + bn2 * qkx;
959                 final double fimy = -yr * bn123k - bn1 * dky + bn2 * qky;
960                 final double fimz = -zr * bn123k - bn1 * dkz + bn2 * qkz;
961                 final double fidx = -xr * drr357k - drr3 * dkx + drr5 * qkx;
962                 final double fidy = -yr * drr357k - drr3 * dky + drr5 * qky;
963                 final double fidz = -zr * drr357k - drr3 * dkz + drr5 * qkz;
964 
965                 final double dir = dix * xr + diy * yr + diz * zr;
966                 final double qix = 2.0 * (qixx * xr + qixy * yr + qixz * zr);
967                 final double qiy = 2.0 * (qixy * xr + qiyy * yr + qiyz * zr);
968                 final double qiz = 2.0 * (qixz * xr + qiyz * yr + qizz * zr);
969                 final double qir = (qix * xr + qiy * yr + qiz * zr) * 0.5;
970                 final double bn123i = bn1 * ci + bn2 * dir + bn3 * qir;
971                 final double ddr357i = drr3 * ci + drr5 * dir + drr7 * qir;
972                 final double fkmx = xr * bn123i - bn1 * dix - bn2 * qix;
973                 final double fkmy = yr * bn123i - bn1 * diy - bn2 * qiy;
974                 final double fkmz = zr * bn123i - bn1 * diz - bn2 * qiz;
975                 final double fkdx = xr * ddr357i - drr3 * dix - drr5 * qix;
976                 final double fkdy = yr * ddr357i - drr3 * diy - drr5 * qiy;
977                 final double fkdz = zr * ddr357i - drr3 * diz - drr5 * qiz;
978                 fix += selfScale * (fimx - fidx);
979                 fiy += selfScale * (fimy - fidy);
980                 fiz += selfScale * (fimz - fidz);
981                 final double xc = selfScale * (fkmx - fkdx);
982                 final double yc = selfScale * (fkmy - fkdy);
983                 final double zc = selfScale * (fkmz - fkdz);
984                 final double fkx = xc * transOp[0][0] + yc * transOp[1][0] + zc * transOp[2][0];
985                 final double fky = xc * transOp[0][1] + yc * transOp[1][1] + zc * transOp[2][1];
986                 final double fkz = xc * transOp[0][2] + yc * transOp[1][2] + zc * transOp[2][2];
987                 field.add(threadID, k, fkx, fky, fkz);
988                 fieldCR.add(threadID, k, fkx, fky, fkz);
989               }
990             }
991             field.add(threadID, i, fix, fiy, fiz);
992             fieldCR.add(threadID, i, fix, fiy, fiz);
993           }
994         }
995       }
996 
997       @Override
998       public IntegerSchedule schedule() {
999         return permanentSchedule;
1000       }
1001 
1002       @Override
1003       public void start() {
1004         time = -System.nanoTime();
1005         threadID = getThreadIndex();
1006         count = 0;
1007         int nAtoms = atoms.length;
1008         if (inductionMaskLocal == null || inductionMaskLocal.length < nAtoms) {
1009           inductionMaskLocal = new double[nAtoms];
1010           energyMaskLocal = new double[nAtoms];
1011           fill(inductionMaskLocal, 1.0);
1012           fill(energyMaskLocal, 1.0);
1013         }
1014       }
1015     }
1016   }
1017 }