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 ffx.crystal.Crystal;
46  import ffx.crystal.SymOp;
47  import ffx.numerics.atomic.AtomicDoubleArray3D;
48  import ffx.potential.bonded.Atom;
49  import ffx.potential.nonbonded.ReciprocalSpace;
50  import ffx.potential.parameters.ForceField;
51  
52  import java.util.List;
53  import java.util.logging.Level;
54  import java.util.logging.Logger;
55  
56  import static ffx.numerics.special.Erf.erfc;
57  import static org.apache.commons.math3.util.FastMath.exp;
58  import static org.apache.commons.math3.util.FastMath.min;
59  import static org.apache.commons.math3.util.FastMath.sqrt;
60  
61  /**
62   * Parallel calculation of the induced dipole field.
63   *
64   * <p>This region can be executed by a ParallelTeam with exactly 2 threads.
65   *
66   * <p>The Real Space and Reciprocal Space Sections will be run concurrently, each with the number of
67   * threads defined by their respective ParallelTeam instances.
68   *
69   * @author Michael J. Schnieders
70   * @since 1.0
71   */
72  public class InducedDipoleFieldRegion extends ParallelRegion {
73  
74    private static final Logger logger = Logger.getLogger(InducedDipoleFieldRegion.class.getName());
75    /**
76     * Specify inter-molecular softcore.
77     */
78    private final boolean intermolecularSoftcore;
79    /**
80     * Specify intra-molecular softcore.
81     */
82    private final boolean intramolecularSoftcore;
83    /**
84     * Dimensions of [nsymm][nAtoms][3]
85     */
86    public double[][][] inducedDipole;
87    public double[][][] inducedDipoleCR;
88    /**
89     * An ordered array of atoms in the system.
90     */
91    private Atom[] atoms;
92    /**
93     * Unit cell and spacegroup information.
94     */
95    private Crystal crystal;
96    /**
97     * When computing the polarization energy at Lambda there are 3 pieces.
98     *
99     * <p>1.) Upol(1) = The polarization energy computed normally (ie. system with ligand).
100    *
101    * <p>2.) Uenv = The polarization energy of the system without the ligand.
102    *
103    * <p>3.) Uligand = The polarization energy of the ligand by itself.
104    *
105    * <p>Upol(L) = L*Upol(1) + (1-L)*(Uenv + Uligand)
106    *
107    * <p>Set the "use" array to true for all atoms for part 1. Set the "use" array to true for all
108    * atoms except the ligand for part 2. Set the "use" array to true only for the ligand atoms for
109    * part 3.
110    *
111    * <p>The "use" array can also be employed to turn off atoms for computing the electrostatic
112    * energy of sub-structures.
113    */
114   private boolean[] use;
115   /**
116    * Molecule number for each atom.
117    */
118   private int[] molecule;
119   private double[] ipdamp;
120   private double[] thole;
121   /**
122    * Dimensions of [nsymm][xyz][nAtoms].
123    */
124   private double[][][] coordinates;
125   /**
126    * Neighbor lists, without atoms beyond the real space cutoff. [nSymm][nAtoms][nIncludedNeighbors]
127    */
128   private int[][][] realSpaceLists;
129   /**
130    * Number of neighboring atoms within the real space cutoff. [nSymm][nAtoms]
131    */
132   private int[][] realSpaceCounts;
133   /**
134    * Pairwise schedule for load balancing.
135    */
136   private IntegerSchedule realSpaceSchedule;
137   /**
138    * Reciprocal space instance.
139    */
140   private ReciprocalSpace reciprocalSpace;
141   private boolean reciprocalSpaceTerm;
142   /**
143    * The current LambdaMode of this PME instance (or OFF for no lambda dependence).
144    */
145   private LambdaMode lambdaMode = LambdaMode.OFF;
146   private double aewald;
147   private double an0, an1;
148   /**
149    * Field array.
150    */
151   private AtomicDoubleArray3D field;
152   /**
153    * Chain rule field array.
154    */
155   private AtomicDoubleArray3D fieldCR;
156 
157   private PMETimings pmeTimings;
158   private final InducedDipoleRealSpaceFieldSection inducedRealSpaceFieldSection;
159   private final InducedDipoleReciprocalFieldSection inducedReciprocalFieldSection;
160 
161 
162   public InducedDipoleFieldRegion(ParallelTeam parallelTeam,
163                                   ForceField forceField, boolean lambdaTerm) {
164     inducedRealSpaceFieldSection = new InducedDipoleRealSpaceFieldSection(parallelTeam);
165     inducedReciprocalFieldSection = new InducedDipoleReciprocalFieldSection();
166 
167     // Flag to indicate application of an intermolecular softcore potential.
168     if (lambdaTerm) {
169       intermolecularSoftcore = forceField.getBoolean("INTERMOLECULAR_SOFTCORE", false);
170       intramolecularSoftcore = forceField.getBoolean("INTRAMOLECULAR_SOFTCORE", false);
171     } else {
172       intermolecularSoftcore = false;
173       intramolecularSoftcore = false;
174     }
175   }
176 
177   /**
178    * Execute the InducedDipoleFieldRegion with the passed ParallelTeam.
179    *
180    * @param sectionTeam The ParallelTeam instance to execute with.
181    */
182   public void executeWith(ParallelTeam sectionTeam) {
183     try {
184       sectionTeam.execute(this);
185     } catch (RuntimeException e) {
186       String message = "Runtime exception computing the induced reciprocal space field.\n";
187       logger.log(Level.WARNING, message, e);
188       throw e;
189     } catch (Exception ex) {
190       String message = "Fatal exception computing the induced reciprocal space field.\n";
191       logger.log(Level.SEVERE, message, ex);
192     }
193   }
194 
195   public void init(
196       Atom[] atoms,
197       Crystal crystal,
198       boolean[] use,
199       int[] molecule,
200       double[] ipdamp,
201       double[] thole,
202       double[][][] coordinates,
203       RealSpaceNeighborParameters realSpaceNeighborParameters,
204       double[][][] inducedDipole,
205       double[][][] inducedDipoleCR,
206       boolean reciprocalSpaceTerm,
207       ReciprocalSpace reciprocalSpace,
208       LambdaMode lambdaMode,
209       EwaldParameters ewaldParameters,
210       AtomicDoubleArray3D field,
211       AtomicDoubleArray3D fieldCR,
212       PMETimings pmeTimings) {
213     // Input
214     this.atoms = atoms;
215     this.crystal = crystal;
216     this.use = use;
217     this.molecule = molecule;
218     this.ipdamp = ipdamp;
219     this.thole = thole;
220     this.coordinates = coordinates;
221     this.realSpaceLists = realSpaceNeighborParameters.realSpaceLists;
222     this.realSpaceCounts = realSpaceNeighborParameters.realSpaceCounts;
223     this.realSpaceSchedule = realSpaceNeighborParameters.realSpaceSchedule;
224     this.inducedDipole = inducedDipole;
225     this.inducedDipoleCR = inducedDipoleCR;
226     this.reciprocalSpaceTerm = reciprocalSpaceTerm;
227     this.reciprocalSpace = reciprocalSpace;
228     this.lambdaMode = lambdaMode;
229     this.aewald = ewaldParameters.aewald;
230     this.an0 = ewaldParameters.an0;
231     this.an1 = ewaldParameters.an1;
232     // Output
233     this.field = field;
234     this.fieldCR = fieldCR;
235     this.pmeTimings = pmeTimings;
236   }
237 
238   @Override
239   public void run() {
240     try {
241       if (reciprocalSpaceTerm && aewald > 0.0) {
242         execute(inducedRealSpaceFieldSection, inducedReciprocalFieldSection);
243       } else {
244         execute(inducedRealSpaceFieldSection);
245       }
246     } catch (Exception e) {
247       String message = "Fatal exception computing the induced dipole field.\n";
248       logger.log(Level.SEVERE, message, e);
249     }
250   }
251 
252   private class InducedDipoleRealSpaceFieldSection extends ParallelSection {
253 
254     private final InducedDipoleRealSpaceFieldRegion inducedDipoleRealSpaceFieldRegion;
255     private final ParallelTeam parallelTeam;
256 
257     InducedDipoleRealSpaceFieldSection(ParallelTeam parallelTeam) {
258       this.parallelTeam = parallelTeam;
259       int nThreads = parallelTeam.getThreadCount();
260       inducedDipoleRealSpaceFieldRegion = new InducedDipoleRealSpaceFieldRegion(nThreads);
261     }
262 
263     @Override
264     public void run() {
265       pmeTimings.realSpaceSCFTotalTime -= System.nanoTime();
266       try {
267         parallelTeam.execute(inducedDipoleRealSpaceFieldRegion);
268       } catch (Exception e) {
269         String message = "Fatal exception computing the real space field.\n";
270         logger.log(Level.SEVERE, message, e);
271       }
272       pmeTimings.realSpaceSCFTotalTime += System.nanoTime();
273     }
274   }
275 
276   private class InducedDipoleReciprocalFieldSection extends ParallelSection {
277 
278     @Override
279     public void run() {
280       reciprocalSpace.performConvolution();
281     }
282   }
283 
284   private class InducedDipoleRealSpaceFieldRegion extends ParallelRegion {
285 
286     private final InducedRealSpaceFieldLoop[] inducedRealSpaceFieldLoop;
287 
288     InducedDipoleRealSpaceFieldRegion(int threadCount) {
289       inducedRealSpaceFieldLoop = new InducedRealSpaceFieldLoop[threadCount];
290     }
291 
292     @Override
293     public void run() {
294       int threadIndex = getThreadIndex();
295       pmeTimings.realSpaceSCFTime[threadIndex] -= System.nanoTime();
296       if (inducedRealSpaceFieldLoop[threadIndex] == null) {
297         inducedRealSpaceFieldLoop[threadIndex] = new InducedRealSpaceFieldLoop();
298       }
299       try {
300         int nAtoms = atoms.length;
301         execute(0, nAtoms - 1, inducedRealSpaceFieldLoop[threadIndex]);
302       } catch (Exception e) {
303         String message = "Fatal exception computing the induced real space field in thread " + getThreadIndex() + "\n";
304         logger.log(Level.SEVERE, message, e);
305       }
306     }
307 
308     private class InducedRealSpaceFieldLoop extends IntegerForLoop {
309 
310       private int threadID;
311       private double[][] ind, indCR;
312       private double[] x, y, z;
313 
314       InducedRealSpaceFieldLoop() {
315       }
316 
317       @Override
318       public void run(int lb, int ub) {
319         final double[] dx = new double[3];
320         final double[][] transOp = new double[3][3];
321 
322         // Loop over a chunk of atoms.
323         int[][] lists = realSpaceLists[0];
324         int[] counts = realSpaceCounts[0];
325         for (int i = lb; i <= ub; i++) {
326           if (!use[i]) {
327             continue;
328           }
329           final int moleculei = molecule[i];
330           double fx = 0.0;
331           double fy = 0.0;
332           double fz = 0.0;
333           double px = 0.0;
334           double py = 0.0;
335           double pz = 0.0;
336           final double xi = x[i];
337           final double yi = y[i];
338           final double zi = z[i];
339           final double[] dipolei = ind[i];
340           final double uix = dipolei[0];
341           final double uiy = dipolei[1];
342           final double uiz = dipolei[2];
343           final double[] dipoleCRi = indCR[i];
344           final double pix = dipoleCRi[0];
345           final double piy = dipoleCRi[1];
346           final double piz = dipoleCRi[2];
347           final double pdi = ipdamp[i];
348           final double pti = thole[i];
349 
350           // Loop over the neighbor list.
351           final int[] list = lists[i];
352           final int npair = counts[i];
353           for (int j = 0; j < npair; j++) {
354             final int k = list[j];
355             if (!use[k]) {
356               continue;
357             }
358             boolean sameMolecule = (moleculei == molecule[k]);
359             if (lambdaMode == LambdaMode.VAPOR) {
360               if ((intermolecularSoftcore && !sameMolecule)
361                   || (intramolecularSoftcore && sameMolecule)) {
362                 continue;
363               }
364             }
365             final double pdk = ipdamp[k];
366             final double ptk = thole[k];
367             dx[0] = x[k] - xi;
368             dx[1] = y[k] - yi;
369             dx[2] = z[k] - zi;
370             final double r2 = crystal.image(dx);
371 
372             // Calculate the error function damping terms.
373             final double r = sqrt(r2);
374             final double rr1 = 1.0 / r;
375             final double rr2 = rr1 * rr1;
376             final double ralpha = aewald * r;
377             final double exp2a = exp(-ralpha * ralpha);
378             final double bn0 = erfc(ralpha) * rr1;
379             final double bn1 = (bn0 + an0 * exp2a) * rr2;
380             final double bn2 = (3.0 * bn1 + an1 * exp2a) * rr2;
381             double scale3 = 1.0;
382             double scale5 = 1.0;
383             double damp = pdi * pdk;
384             final double pgamma = min(pti, ptk);
385             final double rdamp = r * damp;
386             damp = -pgamma * rdamp * rdamp * rdamp;
387             if (damp > -50.0) {
388               final double expdamp = exp(damp);
389               scale3 = 1.0 - expdamp;
390               scale5 = 1.0 - expdamp * (1.0 - damp);
391             }
392             double rr3 = rr1 * rr2;
393             double rr5 = 3.0 * rr3 * rr2;
394             rr3 *= (1.0 - scale3);
395             rr5 *= (1.0 - scale5);
396             final double xr = dx[0];
397             final double yr = dx[1];
398             final double zr = dx[2];
399             final double[] dipolek = ind[k];
400             final double ukx = dipolek[0];
401             final double uky = dipolek[1];
402             final double ukz = dipolek[2];
403             final double ukr = ukx * xr + uky * yr + ukz * zr;
404             final double bn2ukr = bn2 * ukr;
405             final double fimx = -bn1 * ukx + bn2ukr * xr;
406             final double fimy = -bn1 * uky + bn2ukr * yr;
407             final double fimz = -bn1 * ukz + bn2ukr * zr;
408             final double rr5ukr = rr5 * ukr;
409             final double fidx = -rr3 * ukx + rr5ukr * xr;
410             final double fidy = -rr3 * uky + rr5ukr * yr;
411             final double fidz = -rr3 * ukz + rr5ukr * zr;
412             fx += (fimx - fidx);
413             fy += (fimy - fidy);
414             fz += (fimz - fidz);
415             final double[] dipolepk = indCR[k];
416             final double pkx = dipolepk[0];
417             final double pky = dipolepk[1];
418             final double pkz = dipolepk[2];
419             final double pkr = pkx * xr + pky * yr + pkz * zr;
420             final double bn2pkr = bn2 * pkr;
421             final double pimx = -bn1 * pkx + bn2pkr * xr;
422             final double pimy = -bn1 * pky + bn2pkr * yr;
423             final double pimz = -bn1 * pkz + bn2pkr * zr;
424             final double rr5pkr = rr5 * pkr;
425             final double pidx = -rr3 * pkx + rr5pkr * xr;
426             final double pidy = -rr3 * pky + rr5pkr * yr;
427             final double pidz = -rr3 * pkz + rr5pkr * zr;
428             px += (pimx - pidx);
429             py += (pimy - pidy);
430             pz += (pimz - pidz);
431             final double uir = uix * xr + uiy * yr + uiz * zr;
432             final double bn2uir = bn2 * uir;
433             final double fkmx = -bn1 * uix + bn2uir * xr;
434             final double fkmy = -bn1 * uiy + bn2uir * yr;
435             final double fkmz = -bn1 * uiz + bn2uir * zr;
436             final double rr5uir = rr5 * uir;
437             final double fkdx = -rr3 * uix + rr5uir * xr;
438             final double fkdy = -rr3 * uiy + rr5uir * yr;
439             final double fkdz = -rr3 * uiz + rr5uir * zr;
440             field.add(threadID, k, fkmx - fkdx, fkmy - fkdy, fkmz - fkdz);
441             final double pir = pix * xr + piy * yr + piz * zr;
442             final double bn2pir = bn2 * pir;
443             final double pkmx = -bn1 * pix + bn2pir * xr;
444             final double pkmy = -bn1 * piy + bn2pir * yr;
445             final double pkmz = -bn1 * piz + bn2pir * zr;
446             final double rr5pir = rr5 * pir;
447             final double pkdx = -rr3 * pix + rr5pir * xr;
448             final double pkdy = -rr3 * piy + rr5pir * yr;
449             final double pkdz = -rr3 * piz + rr5pir * zr;
450             fieldCR.add(threadID, k, pkmx - pkdx, pkmy - pkdy, pkmz - pkdz);
451           }
452           field.add(threadID, i, fx, fy, fz);
453           fieldCR.add(threadID, i, px, py, pz);
454         }
455 
456         // Loop over symmetry mates.
457         List<SymOp> symOps = crystal.spaceGroup.symOps;
458         int nSymm = symOps.size();
459         for (int iSymm = 1; iSymm < nSymm; iSymm++) {
460           SymOp symOp = crystal.spaceGroup.getSymOp(iSymm);
461           crystal.getTransformationOperator(symOp, transOp);
462           lists = realSpaceLists[iSymm];
463           counts = realSpaceCounts[iSymm];
464           final double[] xs = coordinates[iSymm][0];
465           final double[] ys = coordinates[iSymm][1];
466           final double[] zs = coordinates[iSymm][2];
467           final double[][] inds = inducedDipole[iSymm];
468           final double[][] indCRs = inducedDipoleCR[iSymm];
469 
470           // Loop over a chunk of atoms.
471           for (int i = lb; i <= ub; i++) {
472             if (!use[i]) {
473               continue;
474             }
475             double fx = 0.0;
476             double fy = 0.0;
477             double fz = 0.0;
478             double px = 0.0;
479             double py = 0.0;
480             double pz = 0.0;
481             final double xi = x[i];
482             final double yi = y[i];
483             final double zi = z[i];
484             final double[] dipolei = ind[i];
485             final double uix = dipolei[0];
486             final double uiy = dipolei[1];
487             final double uiz = dipolei[2];
488             final double[] dipoleCRi = indCR[i];
489             final double pix = dipoleCRi[0];
490             final double piy = dipoleCRi[1];
491             final double piz = dipoleCRi[2];
492             final double pdi = ipdamp[i];
493             final double pti = thole[i];
494 
495             // Loop over the neighbor list.
496             final int[] list = lists[i];
497             final int npair = counts[i];
498             for (int j = 0; j < npair; j++) {
499               final int k = list[j];
500               if (!use[k]) {
501                 continue;
502               }
503               double selfScale = 1.0;
504               if (i == k) {
505                 selfScale = 0.5;
506               }
507               final double pdk = ipdamp[k];
508               final double ptk = thole[k];
509               dx[0] = xs[k] - xi;
510               dx[1] = ys[k] - yi;
511               dx[2] = zs[k] - zi;
512               final double r2 = crystal.image(dx);
513 
514               // Calculate the error function damping terms.
515               final double r = sqrt(r2);
516               final double rr1 = 1.0 / r;
517               final double rr2 = rr1 * rr1;
518               final double ralpha = aewald * r;
519               final double exp2a = exp(-ralpha * ralpha);
520               final double bn0 = erfc(ralpha) * rr1;
521               final double bn1 = (bn0 + an0 * exp2a) * rr2;
522               final double bn2 = (3.0 * bn1 + an1 * exp2a) * rr2;
523               double scale3 = 1.0;
524               double scale5 = 1.0;
525               double damp = pdi * pdk;
526               final double pgamma = min(pti, ptk);
527               final double rdamp = r * damp;
528               damp = -pgamma * rdamp * rdamp * rdamp;
529               if (damp > -50.0) {
530                 final double expdamp = exp(damp);
531                 scale3 = 1.0 - expdamp;
532                 scale5 = 1.0 - expdamp * (1.0 - damp);
533               }
534               double rr3 = rr1 * rr2;
535               double rr5 = 3.0 * rr3 * rr2;
536               rr3 *= (1.0 - scale3);
537               rr5 *= (1.0 - scale5);
538               final double xr = dx[0];
539               final double yr = dx[1];
540               final double zr = dx[2];
541               final double[] dipolek = inds[k];
542               final double ukx = dipolek[0];
543               final double uky = dipolek[1];
544               final double ukz = dipolek[2];
545               final double[] dipolepk = indCRs[k];
546               final double pkx = dipolepk[0];
547               final double pky = dipolepk[1];
548               final double pkz = dipolepk[2];
549               final double ukr = ukx * xr + uky * yr + ukz * zr;
550               final double bn2ukr = bn2 * ukr;
551               final double fimx = -bn1 * ukx + bn2ukr * xr;
552               final double fimy = -bn1 * uky + bn2ukr * yr;
553               final double fimz = -bn1 * ukz + bn2ukr * zr;
554               final double rr5ukr = rr5 * ukr;
555               final double fidx = -rr3 * ukx + rr5ukr * xr;
556               final double fidy = -rr3 * uky + rr5ukr * yr;
557               final double fidz = -rr3 * ukz + rr5ukr * zr;
558               fx += selfScale * (fimx - fidx);
559               fy += selfScale * (fimy - fidy);
560               fz += selfScale * (fimz - fidz);
561               final double pkr = pkx * xr + pky * yr + pkz * zr;
562               final double bn2pkr = bn2 * pkr;
563               final double pimx = -bn1 * pkx + bn2pkr * xr;
564               final double pimy = -bn1 * pky + bn2pkr * yr;
565               final double pimz = -bn1 * pkz + bn2pkr * zr;
566               final double rr5pkr = rr5 * pkr;
567               final double pidx = -rr3 * pkx + rr5pkr * xr;
568               final double pidy = -rr3 * pky + rr5pkr * yr;
569               final double pidz = -rr3 * pkz + rr5pkr * zr;
570               px += selfScale * (pimx - pidx);
571               py += selfScale * (pimy - pidy);
572               pz += selfScale * (pimz - pidz);
573               final double uir = uix * xr + uiy * yr + uiz * zr;
574               final double bn2uir = bn2 * uir;
575               final double fkmx = -bn1 * uix + bn2uir * xr;
576               final double fkmy = -bn1 * uiy + bn2uir * yr;
577               final double fkmz = -bn1 * uiz + bn2uir * zr;
578               final double rr5uir = rr5 * uir;
579               final double fkdx = -rr3 * uix + rr5uir * xr;
580               final double fkdy = -rr3 * uiy + rr5uir * yr;
581               final double fkdz = -rr3 * uiz + rr5uir * zr;
582               double xc = selfScale * (fkmx - fkdx);
583               double yc = selfScale * (fkmy - fkdy);
584               double zc = selfScale * (fkmz - fkdz);
585               double kx = (xc * transOp[0][0] + yc * transOp[1][0] + zc * transOp[2][0]);
586               double ky = (xc * transOp[0][1] + yc * transOp[1][1] + zc * transOp[2][1]);
587               double kz = (xc * transOp[0][2] + yc * transOp[1][2] + zc * transOp[2][2]);
588               field.add(threadID, k, kx, ky, kz);
589               final double pir = pix * xr + piy * yr + piz * zr;
590               final double bn2pir = bn2 * pir;
591               final double pkmx = -bn1 * pix + bn2pir * xr;
592               final double pkmy = -bn1 * piy + bn2pir * yr;
593               final double pkmz = -bn1 * piz + bn2pir * zr;
594               final double rr5pir = rr5 * pir;
595               final double pkdx = -rr3 * pix + rr5pir * xr;
596               final double pkdy = -rr3 * piy + rr5pir * yr;
597               final double pkdz = -rr3 * piz + rr5pir * zr;
598               xc = selfScale * (pkmx - pkdx);
599               yc = selfScale * (pkmy - pkdy);
600               zc = selfScale * (pkmz - pkdz);
601               kx = (xc * transOp[0][0] + yc * transOp[1][0] + zc * transOp[2][0]);
602               ky = (xc * transOp[0][1] + yc * transOp[1][1] + zc * transOp[2][1]);
603               kz = (xc * transOp[0][2] + yc * transOp[1][2] + zc * transOp[2][2]);
604               fieldCR.add(threadID, k, kx, ky, kz);
605             }
606             field.add(threadID, i, fx, fy, fz);
607             fieldCR.add(threadID, i, px, py, pz);
608           }
609         }
610       }
611 
612       @Override
613       public IntegerSchedule schedule() {
614         return realSpaceSchedule;
615       }
616 
617       @Override
618       public void start() {
619         threadID = getThreadIndex();
620         x = coordinates[0][0];
621         y = coordinates[0][1];
622         z = coordinates[0][2];
623         ind = inducedDipole[0];
624         indCR = inducedDipoleCR[0];
625       }
626 
627       @Override
628       public void finish() {
629         pmeTimings.realSpaceSCFTime[threadID] += System.nanoTime();
630       }
631     }
632   }
633 }