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.implicit;
39  
40  import static ffx.potential.nonbonded.GeneralizedKirkwood.DEFAULT_GKC;
41  import static ffx.potential.parameters.MultipoleType.t000;
42  import static ffx.potential.parameters.MultipoleType.t001;
43  import static ffx.potential.parameters.MultipoleType.t002;
44  import static ffx.potential.parameters.MultipoleType.t010;
45  import static ffx.potential.parameters.MultipoleType.t011;
46  import static ffx.potential.parameters.MultipoleType.t020;
47  import static ffx.potential.parameters.MultipoleType.t100;
48  import static ffx.potential.parameters.MultipoleType.t101;
49  import static ffx.potential.parameters.MultipoleType.t110;
50  import static ffx.potential.parameters.MultipoleType.t200;
51  import static ffx.utilities.Constants.dWater;
52  import static org.apache.commons.math3.util.FastMath.exp;
53  import static org.apache.commons.math3.util.FastMath.sqrt;
54  
55  import edu.rit.pj.IntegerForLoop;
56  import edu.rit.pj.ParallelRegion;
57  import ffx.crystal.Crystal;
58  import ffx.crystal.SymOp;
59  import ffx.numerics.atomic.AtomicDoubleArray3D;
60  import ffx.potential.bonded.Atom;
61  import ffx.potential.parameters.ForceField;
62  import java.util.List;
63  import java.util.logging.Level;
64  import java.util.logging.Logger;
65  
66  /**
67   * Parallel computation of the Generalized Kirkwood permanent reaction field.
68   *
69   * @author Michael J. Schnieders
70   * @since 1.0
71   */
72  public class PermanentGKFieldRegion extends ParallelRegion {
73  
74    private static final Logger logger = Logger.getLogger(PermanentGKFieldRegion.class.getName());
75    /** Constant factor used with quadrupoles. */
76    private static final double oneThird = 1.0 / 3.0;
77    /** Empirical constant that controls the GK cross-term. */
78    private final double gkc;
79    /** Kirkwood monopole reaction field constant. */
80    private final double fc;
81    /** Kirkwood dipole reaction field constant. */
82    private final double fd;
83    /** Kirkwood quadrupole reaction field constant. */
84    private final double fq;
85  
86    private final PermanentGKFieldLoop[] permanentGKFieldLoop;
87    /** An ordered array of atoms in the system. */
88    protected Atom[] atoms;
89    /** Multipole moments for each symmetry operator. */
90    private double[][][] globalMultipole;
91    /** Periodic boundary conditions and symmetry. */
92    private Crystal crystal;
93    /** Atomic coordinates for each symmetry operator. */
94    private double[][][] sXYZ;
95    /** Neighbor lists for each atom and symmetry operator. */
96    private int[][][] neighborLists;
97    /** Flag to indicate if an atom should be included. */
98    private boolean[] use = null;
99    /** GK cut-off distance squared. */
100   private double cut2;
101   /** Born radius of each atom. */
102   private double[] born;
103   /** Atomic GK field array. */
104   private AtomicDoubleArray3D sharedGKField;
105 
106   public PermanentGKFieldRegion(int nt, ForceField forceField) {
107 
108     // Set the Kirkwood multipolar reaction field constants.
109     double epsilon = forceField.getDouble("GK_EPSILON", dWater);
110     fc = 1.0 * (1.0 - epsilon) / (0.0 + 1.0 * epsilon);
111     fd = 2.0 * (1.0 - epsilon) / (1.0 + 2.0 * epsilon);
112     fq = 3.0 * (1.0 - epsilon) / (2.0 + 3.0 * epsilon);
113 
114     gkc = forceField.getDouble("GKC", DEFAULT_GKC);
115 
116     permanentGKFieldLoop = new PermanentGKFieldLoop[nt];
117     for (int i = 0; i < nt; i++) {
118       permanentGKFieldLoop[i] = new PermanentGKFieldLoop();
119     }
120   }
121 
122   public void init(
123       Atom[] atoms,
124       double[][][] globalMultipole,
125       Crystal crystal,
126       double[][][] sXYZ,
127       int[][][] neighborLists,
128       boolean[] use,
129       double cut2,
130       double[] born,
131       AtomicDoubleArray3D sharedGKField) {
132     this.atoms = atoms;
133     this.globalMultipole = globalMultipole;
134     this.crystal = crystal;
135     this.sXYZ = sXYZ;
136     this.neighborLists = neighborLists;
137     this.use = use;
138     this.cut2 = cut2;
139     this.born = born;
140     this.sharedGKField = sharedGKField;
141   }
142 
143   @Override
144   public void run() {
145     try {
146       int threadIndex = getThreadIndex();
147       int nAtoms = atoms.length;
148       execute(0, nAtoms - 1, permanentGKFieldLoop[threadIndex]);
149     } catch (Exception e) {
150       String message = "Fatal exception computing GK Energy in thread " + getThreadIndex() + "\n";
151       logger.log(Level.SEVERE, message, e);
152     }
153   }
154 
155   /**
156    * Compute the Generalized Kirkwood permanent reaction field.
157    *
158    * @since 1.0
159    */
160   private class PermanentGKFieldLoop extends IntegerForLoop {
161 
162     private final double[][] a;
163     private final double[] gc;
164     private final double[] gux, guy, guz;
165     private final double[] gqxx, gqyy, gqzz;
166     private final double[] gqxy, gqxz, gqyz;
167     private final double[] dx_local;
168     private int threadID;
169     private double xi, yi, zi;
170     private double ci, uxi, uyi, uzi, qxxi, qxyi, qxzi, qyyi, qyzi, qzzi;
171     private double rbi;
172     private int iSymm;
173     private double[][] transOp;
174     // Extra padding to avert cache interference.
175     private long pad0, pad1, pad2, pad3, pad4, pad5, pad6, pad7;
176     private long pad8, pad9, pada, padb, padc, padd, pade, padf;
177 
178     PermanentGKFieldLoop() {
179       a = new double[4][3];
180       gc = new double[11];
181       gux = new double[11];
182       guy = new double[11];
183       guz = new double[11];
184       gqxx = new double[11];
185       gqyy = new double[11];
186       gqzz = new double[11];
187       gqxy = new double[11];
188       gqxz = new double[11];
189       gqyz = new double[11];
190       dx_local = new double[3];
191       transOp = new double[3][3];
192     }
193 
194     @Override
195     public void run(int lb, int ub) {
196 
197       threadID = getThreadIndex();
198 
199       double[] x = sXYZ[0][0];
200       double[] y = sXYZ[0][1];
201       double[] z = sXYZ[0][2];
202 
203       int nSymm = crystal.spaceGroup.symOps.size();
204       List<SymOp> symOps = crystal.spaceGroup.symOps;
205       for (iSymm = 0; iSymm < nSymm; iSymm++) {
206         SymOp symOp = symOps.get(iSymm);
207         crystal.getTransformationOperator(symOp, transOp);
208         for (int i = lb; i <= ub; i++) {
209           if (!use[i]) {
210             continue;
211           }
212           xi = x[i];
213           yi = y[i];
214           zi = z[i];
215           double[] multipolei = globalMultipole[0][i];
216           ci = multipolei[t000];
217           uxi = multipolei[t100];
218           uyi = multipolei[t010];
219           uzi = multipolei[t001];
220           qxxi = multipolei[t200] * oneThird;
221           qxyi = multipolei[t110] * oneThird;
222           qxzi = multipolei[t101] * oneThird;
223           qyyi = multipolei[t020] * oneThird;
224           qyzi = multipolei[t011] * oneThird;
225           qzzi = multipolei[t002] * oneThird;
226           rbi = born[i];
227           int[] list = neighborLists[iSymm][i];
228           for (int k : list) {
229             if (!use[k]) {
230               continue;
231             }
232             permanentGKField(i, k);
233           }
234 
235           // Include the self permanent reaction field, which is not in the neighbor list.
236           if (iSymm == 0) {
237             permanentGKField(i, i);
238           }
239         }
240       }
241     }
242 
243     private void permanentGKField(int i, int k) {
244       dx_local[0] = sXYZ[iSymm][0][k] - xi;
245       dx_local[1] = sXYZ[iSymm][1][k] - yi;
246       dx_local[2] = sXYZ[iSymm][2][k] - zi;
247       final double r2 = crystal.image(dx_local);
248       if (r2 > cut2) {
249         return;
250       }
251       double xr = dx_local[0];
252       double yr = dx_local[1];
253       double zr = dx_local[2];
254       double xr2 = xr * xr;
255       double yr2 = yr * yr;
256       double zr2 = zr * zr;
257       final double rbk = born[k];
258       final double[] multipolek = globalMultipole[iSymm][k];
259       final double ck = multipolek[t000];
260       final double uxk = multipolek[t100];
261       final double uyk = multipolek[t010];
262       final double uzk = multipolek[t001];
263       final double qxxk = multipolek[t200] * oneThird;
264       final double qxyk = multipolek[t110] * oneThird;
265       final double qxzk = multipolek[t101] * oneThird;
266       final double qyyk = multipolek[t020] * oneThird;
267       final double qyzk = multipolek[t011] * oneThird;
268       final double qzzk = multipolek[t002] * oneThird;
269       final double rb2 = rbi * rbk;
270       final double expterm = exp(-r2 / (gkc * rb2));
271       final double expc = expterm / gkc;
272       final double expc1 = 1.0 - expc;
273       final double dexpc = -2.0 / (gkc * rb2);
274       final double expcdexpc = -expc * dexpc;
275       final double gf2 = 1.0 / (r2 + rb2 * expterm);
276       final double gf = sqrt(gf2);
277       final double gf3 = gf2 * gf;
278       final double gf5 = gf3 * gf2;
279       final double gf7 = gf5 * gf2;
280 
281       // Reaction potential auxiliary terms.
282       a[0][0] = gf;
283       a[1][0] = -gf3;
284       a[2][0] = 3.0 * gf5;
285       a[3][0] = -15.0 * gf7;
286 
287       // Reaction potential gradient auxiliary terms.
288       a[0][1] = expc1 * a[1][0];
289       a[1][1] = expc1 * a[2][0];
290       a[2][1] = expc1 * a[3][0];
291       // 2nd reaction potential gradient auxiliary terms.
292       a[1][2] = expc1 * a[2][1] + expcdexpc * a[2][0];
293 
294       // Multiply the potential auxiliary terms by their dielectric functions.
295       a[0][1] = fc * a[0][1];
296       a[1][0] = fd * a[1][0];
297       a[1][1] = fd * a[1][1];
298       a[1][2] = fd * a[1][2];
299       a[2][0] = fq * a[2][0];
300       a[2][1] = fq * a[2][1];
301 
302       // Unweighted reaction potential tensor.
303       gux[1] = xr * a[1][0];
304       guy[1] = yr * a[1][0];
305       guz[1] = zr * a[1][0];
306 
307       // Unweighted reaction potential gradient tensor.
308       gc[2] = xr * a[0][1];
309       gc[3] = yr * a[0][1];
310       gc[4] = zr * a[0][1];
311       gux[2] = a[1][0] + xr2 * a[1][1];
312       gux[3] = xr * yr * a[1][1];
313       gux[4] = xr * zr * a[1][1];
314       guy[2] = gux[3];
315       guy[3] = a[1][0] + yr2 * a[1][1];
316       guy[4] = yr * zr * a[1][1];
317       guz[2] = gux[4];
318       guz[3] = guy[4];
319       guz[4] = a[1][0] + zr2 * a[1][1];
320       gqxx[2] = xr * (2.0 * a[2][0] + xr2 * a[2][1]);
321       gqxx[3] = yr * xr2 * a[2][1];
322       gqxx[4] = zr * xr2 * a[2][1];
323       gqyy[2] = xr * yr2 * a[2][1];
324       gqyy[3] = yr * (2.0 * a[2][0] + yr2 * a[2][1]);
325       gqyy[4] = zr * yr2 * a[2][1];
326       gqzz[2] = xr * zr2 * a[2][1];
327       gqzz[3] = yr * zr2 * a[2][1];
328       gqzz[4] = zr * (2.0 * a[2][0] + zr2 * a[2][1]);
329       gqxy[2] = yr * (a[2][0] + xr2 * a[2][1]);
330       gqxy[3] = xr * (a[2][0] + yr2 * a[2][1]);
331       gqxy[4] = zr * xr * yr * a[2][1];
332       gqxz[2] = zr * (a[2][0] + xr2 * a[2][1]);
333       gqxz[3] = gqxy[4];
334       gqxz[4] = xr * (a[2][0] + zr2 * a[2][1]);
335       gqyz[2] = gqxy[4];
336       gqyz[3] = zr * (a[2][0] + yr2 * a[2][1]);
337       gqyz[4] = yr * (a[2][0] + zr2 * a[2][1]);
338 
339       // Unweighted 2nd reaction potential gradient tensor.
340       gux[5] = xr * (3.0 * a[1][1] + xr2 * a[1][2]);
341       gux[6] = yr * (a[1][1] + xr2 * a[1][2]);
342       gux[7] = zr * (a[1][1] + xr2 * a[1][2]);
343       gux[8] = xr * (a[1][1] + yr2 * a[1][2]);
344       gux[9] = zr * xr * yr * a[1][2];
345       gux[10] = xr * (a[1][1] + zr2 * a[1][2]);
346       guy[5] = yr * (a[1][1] + xr2 * a[1][2]);
347       guy[6] = xr * (a[1][1] + yr2 * a[1][2]);
348       guy[7] = gux[9];
349       guy[8] = yr * (3.0 * a[1][1] + yr2 * a[1][2]);
350       guy[9] = zr * (a[1][1] + yr2 * a[1][2]);
351       guy[10] = yr * (a[1][1] + zr2 * a[1][2]);
352       guz[5] = zr * (a[1][1] + xr2 * a[1][2]);
353       guz[6] = gux[9];
354       guz[7] = xr * (a[1][1] + zr2 * a[1][2]);
355       guz[8] = zr * (a[1][1] + yr2 * a[1][2]);
356       guz[9] = yr * (a[1][1] + zr2 * a[1][2]);
357       guz[10] = zr * (3.0 * a[1][1] + zr2 * a[1][2]);
358 
359       // Generalized Kirkwood permanent reaction field.
360       double fix =
361           uxk * gux[2]
362               + uyk * gux[3]
363               + uzk * gux[4]
364               + 0.5
365                   * (ck * gux[1]
366                       + qxxk * gux[5]
367                       + qyyk * gux[8]
368                       + qzzk * gux[10]
369                       + 2.0 * (qxyk * gux[6] + qxzk * gux[7] + qyzk * gux[9]))
370               + 0.5
371                   * (ck * gc[2]
372                       + qxxk * gqxx[2]
373                       + qyyk * gqyy[2]
374                       + qzzk * gqzz[2]
375                       + 2.0 * (qxyk * gqxy[2] + qxzk * gqxz[2] + qyzk * gqyz[2]));
376       double fiy =
377           uxk * guy[2]
378               + uyk * guy[3]
379               + uzk * guy[4]
380               + 0.5
381                   * (ck * guy[1]
382                       + qxxk * guy[5]
383                       + qyyk * guy[8]
384                       + qzzk * guy[10]
385                       + 2.0 * (qxyk * guy[6] + qxzk * guy[7] + qyzk * guy[9]))
386               + 0.5
387                   * (ck * gc[3]
388                       + qxxk * gqxx[3]
389                       + qyyk * gqyy[3]
390                       + qzzk * gqzz[3]
391                       + 2.0 * (qxyk * gqxy[3] + qxzk * gqxz[3] + qyzk * gqyz[3]));
392       double fiz =
393           uxk * guz[2]
394               + uyk * guz[3]
395               + uzk * guz[4]
396               + 0.5
397                   * (ck * guz[1]
398                       + qxxk * guz[5]
399                       + qyyk * guz[8]
400                       + qzzk * guz[10]
401                       + 2.0 * (qxyk * guz[6] + qxzk * guz[7] + qyzk * guz[9]))
402               + 0.5
403                   * (ck * gc[4]
404                       + qxxk * gqxx[4]
405                       + qyyk * gqyy[4]
406                       + qzzk * gqzz[4]
407                       + 2.0 * (qxyk * gqxy[4] + qxzk * gqxz[4] + qyzk * gqyz[4]));
408       double fkx =
409           uxi * gux[2]
410               + uyi * gux[3]
411               + uzi * gux[4]
412               - 0.5
413                   * (ci * gux[1]
414                       + qxxi * gux[5]
415                       + qyyi * gux[8]
416                       + qzzi * gux[10]
417                       + 2.0 * (qxyi * gux[6] + qxzi * gux[7] + qyzi * gux[9]))
418               - 0.5
419                   * (ci * gc[2]
420                       + qxxi * gqxx[2]
421                       + qyyi * gqyy[2]
422                       + qzzi * gqzz[2]
423                       + 2.0 * (qxyi * gqxy[2] + qxzi * gqxz[2] + qyzi * gqyz[2]));
424       double fky =
425           uxi * guy[2]
426               + uyi * guy[3]
427               + uzi * guy[4]
428               - 0.5
429                   * (ci * guy[1]
430                       + qxxi * guy[5]
431                       + qyyi * guy[8]
432                       + qzzi * guy[10]
433                       + 2.0 * (qxyi * guy[6] + qxzi * guy[7] + qyzi * guy[9]))
434               - 0.5
435                   * (ci * gc[3]
436                       + qxxi * gqxx[3]
437                       + qyyi * gqyy[3]
438                       + qzzi * gqzz[3]
439                       + 2.0 * (qxyi * gqxy[3] + qxzi * gqxz[3] + qyzi * gqyz[3]));
440       double fkz =
441           uxi * guz[2]
442               + uyi * guz[3]
443               + uzi * guz[4]
444               - 0.5
445                   * (ci * guz[1]
446                       + qxxi * guz[5]
447                       + qyyi * guz[8]
448                       + qzzi * guz[10]
449                       + 2.0 * (qxyi * guz[6] + qxzi * guz[7] + qyzi * guz[9]))
450               - 0.5
451                   * (ci * gc[4]
452                       + qxxi * gqxx[4]
453                       + qyyi * gqyy[4]
454                       + qzzi * gqzz[4]
455                       + 2.0 * (qxyi * gqxy[4] + qxzi * gqxz[4] + qyzi * gqyz[4]));
456 
457       // Scale the self-field by half, such that it sums to one below.
458       if (i == k) {
459         fix *= 0.5;
460         fiy *= 0.5;
461         fiz *= 0.5;
462         fkx *= 0.5;
463         fky *= 0.5;
464         fkz *= 0.5;
465       }
466 
467       double xc = fkx;
468       double yc = fky;
469       double zc = fkz;
470       fkx = xc * transOp[0][0] + yc * transOp[1][0] + zc * transOp[2][0];
471       fky = xc * transOp[0][1] + yc * transOp[1][1] + zc * transOp[2][1];
472       fkz = xc * transOp[0][2] + yc * transOp[1][2] + zc * transOp[2][2];
473 
474       sharedGKField.add(threadID, i, fix, fiy, fiz);
475       sharedGKField.add(threadID, k, fkx, fky, fkz);
476     }
477   }
478 }