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