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