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.implicit;
39
40 import static ffx.potential.nonbonded.VanDerWaalsForm.getCombinedEps;
41 import static ffx.potential.nonbonded.VanDerWaalsForm.getCombinedRadius;
42 import static org.apache.commons.math3.util.FastMath.PI;
43 import static org.apache.commons.math3.util.FastMath.max;
44 import static org.apache.commons.math3.util.FastMath.min;
45 import static org.apache.commons.math3.util.FastMath.pow;
46 import static org.apache.commons.math3.util.FastMath.sqrt;
47
48 import edu.rit.pj.IntegerForLoop;
49 import edu.rit.pj.ParallelRegion;
50 import edu.rit.pj.reduction.SharedDouble;
51 import ffx.crystal.Crystal;
52 import ffx.numerics.atomic.AtomicDoubleArray3D;
53 import ffx.potential.bonded.Atom;
54 import ffx.potential.nonbonded.VanDerWaalsForm.EPSILON_RULE;
55 import ffx.potential.nonbonded.VanDerWaalsForm.RADIUS_RULE;
56 import ffx.potential.parameters.ForceField;
57 import ffx.potential.parameters.VDWType;
58 import java.util.logging.Level;
59 import java.util.logging.Logger;
60
61 /**
62 * Parallel calculation of continuum dispersion energy via pairwise descreening.
63 *
64 * @author Michael J. Schnieders
65 * @since 1.0
66 */
67 public class DispersionRegion extends ParallelRegion {
68
69 /** The dispersion integral begins for each atom at: Rmin_ij + DISPERSION_OFFSET */
70 public static final double DEFAULT_DISPERSION_OFFSET = 1.056;
71 /** Each solute atom blocks dispersion interactions with solvent: Rmin + SOLUTE_OFFSET */
72 public static final double DEFAULT_SOLUTE_OFFSET = 0.0;
73
74 private static final Logger logger = Logger.getLogger(DispersionRegion.class.getName());
75 /** The dispersion integral HCT overlap scale factor. */
76 private static final double DEFAULT_DISP_OVERLAP_FACTOR = 0.75;
77 /**
78 * Conversion between solute-solvent dispersion enthalpy and solute-solvent dispersion free
79 * energy.
80 */
81 private static final double SLEVY = 1.0;
82 /** Number density of water (water per A^3). */
83 private static final double AWATER = 0.033428;
84 /** AMOEBA '03 Water oxygen epsilon. */
85 private static final double EPSO = 0.1100;
86 /** AMOEBA '03 Water hydrogen epsilon. */
87 private static final double EPSH = 0.0135;
88 /** AMOEBA '03 Water oxygen Rmin (A). */
89 private static final double RMINO = 1.7025;
90 /** AMOEBA '03 Water hydrogen Rmin (A). */
91 private static final double RMINH = 1.3275;
92 /** AMOEBA epsilon rule. */
93 private static final EPSILON_RULE epsilonRule = EPSILON_RULE.HHG;
94 /** AMOEBA radius rule. */
95 private static final RADIUS_RULE radiusRule = RADIUS_RULE.CUBIC_MEAN;
96
97 private final DispersionLoop[] dispersionLoop;
98 private final SharedDouble sharedDispersion;
99 /** An ordered array of atoms in the system. */
100 protected Atom[] atoms;
101 /** Periodic boundary conditions and symmetry. */
102 private Crystal crystal;
103 /** Flag to indicate if an atom should be included. */
104 private boolean[] use = null;
105 /** Neighbor lists for each atom and symmetry operator. */
106 private int[][][] neighborLists;
107 /** Cartesian coordinates of each atom. */
108 private double[] x, y, z;
109 /** GK cut-off distance squared. */
110 private double cut2;
111 /** Boolean flag to indicate GK will be scaled by the lambda state variable. */
112 private boolean gradient = false;
113 /** Gradient array for each thread. */
114 private AtomicDoubleArray3D grad;
115 /** Where the dispersion integral begins for each atom (A): Rmin + dispersionOffset */
116 private double dispersionOffest;
117 /**
118 * Each solute atom blocks dispersion interactions with solvent with this radius (A): Rmin +
119 * soluteOffset
120 */
121 private double soluteOffset;
122 /** The dispersion integral HCT overlap scale factor (unitless). */
123 private double dispersionOverlapFactor;
124
125 /**
126 * DispersionRegion constructor.
127 *
128 * @param nt Number of threads.
129 * @param atoms Atom array.
130 * @param forceField ForceField in use.
131 */
132 public DispersionRegion(int nt, Atom[] atoms, ForceField forceField) {
133 dispersionLoop = new DispersionLoop[nt];
134 for (int i = 0; i < nt; i++) {
135 dispersionLoop[i] = new DispersionLoop();
136 }
137 sharedDispersion = new SharedDouble();
138
139 dispersionOffest = forceField.getDouble("DISPERSION_OFFSET", DEFAULT_DISPERSION_OFFSET);
140 soluteOffset = forceField.getDouble("SOLUTE_OFFSET", DEFAULT_SOLUTE_OFFSET);
141 dispersionOverlapFactor =
142 forceField.getDouble("DISP_OVERLAP_FACTOR", DEFAULT_DISP_OVERLAP_FACTOR);
143
144 allocate(atoms);
145 }
146
147 /**
148 * Compute a Buffered-14-7 tail correction.
149 *
150 * @param r0 The separation distance where the intregal begins.
151 * @param eps The mixed eps value.
152 * @param rmin The mixed rmin value.
153 * @return The tail correction.
154 */
155 private static double tailCorrection(double r0, double eps, double rmin) {
156 if (r0 < rmin) {
157 double r03 = r0 * r0 * r0;
158 double rmin3 = rmin * rmin * rmin;
159 return -eps * 4.0 * PI * (rmin3 - r03) / 3.0 - eps * PI * 18.0 * rmin3 / 11.0;
160 } else {
161 double ri2 = r0 * r0;
162 double ri4 = ri2 * ri2;
163 double ri7 = r0 * ri2 * ri4;
164 double ri11 = ri7 * ri4;
165 double rmin2 = rmin * rmin;
166 double rmin4 = rmin2 * rmin2;
167 double rmin7 = rmin * rmin2 * rmin4;
168 return 2.0 * PI * eps * rmin7 * (2.0 * rmin7 - 11.0 * ri7) / (11.0 * ri11);
169 }
170 }
171
172 /**
173 * Allocate storage given the Atom array.
174 *
175 * @param atoms Atom array in use.
176 */
177 public void allocate(Atom[] atoms) {
178 this.atoms = atoms;
179 }
180
181 /**
182 * The dispersion integral begins offset from the vdW radius.
183 *
184 * @param dispersionOffest The dispersion integral offset.
185 */
186 public void setDispersionOffest(double dispersionOffest) {
187 this.dispersionOffest = dispersionOffest;
188 }
189
190 /**
191 * The dispersion integral begins offset from the vdW radius.
192 *
193 * @return the dispersion integral offset.
194 */
195 public double getDispersionOffset() {
196 return dispersionOffest;
197 }
198
199 public double getDispersionOverlapFactor() {
200 return dispersionOverlapFactor;
201 }
202
203 /**
204 * Set the dispersion overlap HCT scale factor.
205 *
206 * @param dispersionOverlapFactor The dispersion integral HCT scale factor.
207 */
208 public void setDispersionOverlapFactor(double dispersionOverlapFactor) {
209 this.dispersionOverlapFactor = dispersionOverlapFactor;
210 }
211
212 public double getEnergy() {
213 return sharedDispersion.get();
214 }
215
216 public double getSoluteOffset() {
217 return soluteOffset;
218 }
219
220 public void setSoluteOffset(double soluteOffset) {
221 this.soluteOffset = soluteOffset;
222 }
223
224 /**
225 * Initialize the DispersionRegion for energy calculation.
226 *
227 * @param atoms Atom array.
228 * @param crystal Crystal for periodic boundary conditions.
229 * @param use Flag to indicate an atom is to be used.
230 * @param neighborLists Neighbor-list for each atom.
231 * @param x X-coordinate array.
232 * @param y Y-coordinate array.
233 * @param z Z-coordinate array.
234 * @param cut2 The cut-off distance squared.
235 * @param gradient If true, compute the gradient.
236 * @param grad Array to store the gradient.
237 */
238 public void init(
239 Atom[] atoms,
240 Crystal crystal,
241 boolean[] use,
242 int[][][] neighborLists,
243 double[] x,
244 double[] y,
245 double[] z,
246 double cut2,
247 boolean gradient,
248 AtomicDoubleArray3D grad) {
249 this.atoms = atoms;
250 this.crystal = crystal;
251 this.use = use;
252 this.neighborLists = neighborLists;
253 this.x = x;
254 this.y = y;
255 this.z = z;
256 this.cut2 = cut2;
257 this.gradient = gradient;
258 this.grad = grad;
259 }
260
261 /** {@inheritDoc} */
262 @Override
263 public void run() {
264 try {
265 int nAtoms = atoms.length;
266 execute(0, nAtoms - 1, dispersionLoop[getThreadIndex()]);
267 } catch (Exception e) {
268 String message =
269 "Fatal exception computing Dispersion energy in thread " + getThreadIndex() + "\n";
270 logger.log(Level.SEVERE, message, e);
271 }
272 }
273
274 /** {@inheritDoc} */
275 @Override
276 public void start() {
277 sharedDispersion.set(0.0);
278 }
279
280 /**
281 * Compute Dispersion energy for a range of atoms via pairwise descreening.
282 *
283 * @since 1.0
284 */
285 private class DispersionLoop extends IntegerForLoop {
286
287 private final double[] dx_local;
288 private double edisp;
289 private double r, r2, r3;
290 private double xr, yr, zr;
291 private int threadID;
292
293 DispersionLoop() {
294 dx_local = new double[3];
295 }
296
297 @Override
298 public void finish() {
299 sharedDispersion.addAndGet(edisp);
300 }
301
302 @Override
303 public void run(int lb, int ub) {
304 for (int i = lb; i <= ub; i++) {
305 if (!use[i]) {
306 continue;
307 }
308
309 // Compute the limit of atom alone in solvent.
310 VDWType type = atoms[i].getVDWType();
311 double epsi = type.wellDepth;
312 double rmini = type.radius / 2.0;
313 double cDisp = 0.0;
314 if (rmini > 0.0 && epsi > 0.0) {
315 double emixo = getCombinedEps(EPSO, epsi, RMINO, rmini, epsilonRule);
316 double rmixo = getCombinedRadius(RMINO, rmini, radiusRule);
317 // Start of integration of dispersion for atom i with water oxygen.
318 double riO = rmixo / 2.0 + dispersionOffest;
319 cDisp = tailCorrection(riO, emixo, rmixo);
320 double emixh = getCombinedEps(EPSH, epsi, RMINH, rmini, epsilonRule);
321 double rmixh = getCombinedRadius(RMINH, rmini, radiusRule);
322 // Start of integration of dispersion for atom i with water hydrogen.
323 double riH = rmixh / 2.0 + dispersionOffest;
324 cDisp += 2.0 * tailCorrection(riH, emixh, rmixh);
325 cDisp = SLEVY * AWATER * cDisp;
326 }
327
328 edisp += cDisp;
329
330 // Now descreen over neighbors.
331 double sum = 0.0;
332 final double xi = x[i];
333 final double yi = y[i];
334 final double zi = z[i];
335 int[] list = neighborLists[0][i];
336 for (int k : list) {
337 if (!use[k] || i == k) {
338 continue;
339 }
340 if (atoms[k].getVDWType().radius > 0.0) {
341 dx_local[0] = xi - x[k];
342 dx_local[1] = yi - y[k];
343 dx_local[2] = zi - z[k];
344 r2 = crystal.image(dx_local);
345 if (r2 > cut2) {
346 continue;
347 }
348 xr = dx_local[0];
349 yr = dx_local[1];
350 zr = dx_local[2];
351 r = sqrt(r2);
352 r3 = r * r2;
353 // Atom i descreened by atom k.
354 sum += removeSoluteDispersion(i, k);
355 // Flip the sign on {xr, yr, zr};
356 xr = -xr;
357 yr = -yr;
358 zr = -zr;
359 // Atom k descreened by atom i.
360 sum += removeSoluteDispersion(k, i);
361 }
362 }
363 // Subtract descreening.
364 edisp -= SLEVY * AWATER * sum;
365 }
366 }
367
368 @Override
369 public void start() {
370 threadID = getThreadIndex();
371 edisp = 0;
372 }
373
374 /**
375 * Remove solute dispersion between atom i and solvent due to blocking of solveny by atom k.
376 *
377 * @param i Atom i index.
378 * @param k Atom k index.
379 * @return Reduction of dispersion.
380 */
381 private double removeSoluteDispersion(int i, int k) {
382
383 // Van der Waals parameters for atom i.
384 VDWType type = atoms[i].getVDWType();
385 double epsi = type.wellDepth;
386 double rmini = type.radius / 2.0;
387
388 // Van der Waals radius for atom k.
389 double rmink = atoms[k].getVDWType().radius / 2.0;
390
391 // Parameters for atom i with water oxygen.
392 double emixo = getCombinedEps(EPSO, epsi, RMINO, rmini, epsilonRule);
393 double rmixo = getCombinedRadius(RMINO, rmini, radiusRule);
394 // Start of integration of dispersion for atom i with water oxygen.
395 double riO = rmixo / 2.0 + dispersionOffest;
396 double nO = 1.0;
397
398 // Parameters for atom i with water hydrogen.
399 double emixh = getCombinedEps(EPSH, epsi, RMINH, rmini, epsilonRule);
400 double rmixh = getCombinedRadius(RMINH, rmini, radiusRule);
401 // Start of integration of dispersion for atom i with water oxygen.
402 double riH = rmixh / 2.0 + dispersionOffest;
403 double nH = 2.0;
404
405 // Atom k blocks the interaction of atom i with solvent.
406 double sk = (rmink + soluteOffset) * dispersionOverlapFactor;
407 return interact(i, k, nO, riO, sk, rmixo, emixo) + interact(i, k, nH, riH, sk, rmixh, emixh);
408 }
409
410 /**
411 * Solvent dispersion blocking of atom i by atom k.
412 *
413 * @param i Index of atom i.
414 * @param k Index of atom k.
415 * @param factor Apply this factor to the interaction.
416 * @param ri Start of the dispersion integral.
417 * @param sk Size of the solvent blocking atom.
418 * @param rmix Combined Rmin for atom i and solvent atom.
419 * @param emix Combined Eps for atom i and solvent atom.
420 * @return Contribution of atom k to removing dispersion.
421 */
422 private double interact(
423 int i, int k, double factor, double ri, double sk, double rmix, double emix) {
424 double sum = 0.0;
425 // Nothing to do if the integral begins beyond r + sk (i.e. atom k does not exclude solvent)
426 if (ri < r + sk) {
427 // Zero out the derivative contribution of atom k.
428 double de = 0.0;
429 double sk2 = sk * sk;
430 // Compute the maximum of 1) the beginning of the integral and 2) closest edge of atom K.
431 double iStart = max(ri, r - sk);
432 // Use this as the lower limit for integrating the constant eps value below Rmin.
433 double lik = iStart;
434 // Interaction with water from lik to Rmin; nothing to do if the lower limit is greater than
435 // Rmin.
436 if (lik < rmix) {
437 double lik2 = lik * lik;
438 double lik3 = lik2 * lik;
439 double lik4 = lik3 * lik;
440 // Upper limit is the minimum of Rmin and the farthest edge of atom K.
441 double uik = min(r + sk, rmix);
442 double uik2 = uik * uik;
443 double uik3 = uik2 * uik;
444 double uik4 = uik3 * uik;
445 sum = integralBeforeRMin(emix, r, r2, sk2, lik2, lik3, lik4, uik2, uik3, uik4);
446 if (gradient) {
447 de =
448 integralBeforeRminDerivative(
449 ri, emix, rmix, r, r2, r3, sk, sk2, lik, lik2, lik3, uik, uik2, uik3);
450 }
451 }
452 // Upper limit the variable part of Uwca always the farthest edge of atom K.
453 double uik = r + sk;
454 // Interaction with water beyond Rmin, from lik to uik = r + sk.
455 if (uik > rmix) {
456 // Start the integral at the max of 1) iStart and 2) Rmin.
457 lik = max(iStart, rmix);
458 double lik2 = lik * lik;
459 double lik3 = lik2 * lik;
460 double lik4 = lik3 * lik;
461 double lik5 = lik4 * lik;
462 double lik6 = lik5 * lik;
463 double lik10 = lik5 * lik5;
464 double lik11 = lik10 * lik;
465 double lik12 = lik11 * lik;
466 double uik2 = uik * uik;
467 double uik3 = uik2 * uik;
468 double uik4 = uik3 * uik;
469 double uik5 = uik4 * uik;
470 double uik10 = uik5 * uik5;
471 double uik11 = uik10 * uik;
472 double uik12 = uik11 * uik;
473 double rmix7 = pow(rmix, 7);
474 sum +=
475 integratlAfterRmin(
476 emix, rmix7, r, r2, sk2, lik, lik2, lik3, lik4, lik5, lik10, lik11, lik12, uik,
477 uik2, uik3, uik4, uik5, uik10, uik11, uik12);
478 if (gradient) {
479 double lik13 = lik12 * lik;
480 double uik6 = uik5 * uik;
481 double uik13 = uik12 * uik;
482 de +=
483 integratlAfterRminDerivative(
484 ri, emix, rmix, rmix7, iStart, r, r2, r3, sk, sk2, lik, lik2, lik3, lik5, lik6,
485 lik12, lik13, uik, uik2, uik3, uik6, uik13);
486 }
487 }
488 // Increment the individual dispersion gradient components.
489 if (gradient) {
490 de = -de / r * factor * SLEVY * AWATER;
491 double dedx = de * xr;
492 double dedy = de * yr;
493 double dedz = de * zr;
494 grad.add(threadID, i, dedx, dedy, dedz);
495 grad.sub(threadID, k, dedx, dedy, dedz);
496 }
497 }
498 return factor * sum;
499 }
500
501 /**
502 * Integrate over the constant portion of the WCA disperion interaction Uwca(x) = eps; x < rmin.
503 *
504 * @param eps The well depth.
505 * @param rij The separation between the current atom and solvent blocking atom.
506 * @param rij2 The separation squared.
507 * @param rho2 The scaled size of the solvent blocking atom squared.
508 * @param lik2 The beginning of the integral squared.
509 * @param lik3 The beginning of the integral to the third.
510 * @param lik4 The beginning of the integral to the fourth.
511 * @param uik2 The end of the integral squared.
512 * @param uik3 The end of the integral to the third.
513 * @param uik4 The end of the integral to the fourth.
514 * @return The integral.
515 */
516 private double integralBeforeRMin(
517 double eps,
518 double rij,
519 double rij2,
520 double rho2,
521 double lik2,
522 double lik3,
523 double lik4,
524 double uik2,
525 double uik3,
526 double uik4) {
527 return -eps
528 * (4.0
529 * PI
530 * (3.0 * (lik4 - uik4) - 8.0 * rij * (lik3 - uik3) + 6.0 * (rij2 - rho2) * (lik2 - uik2))
531 / (48.0 * rij));
532 }
533
534 /**
535 * Derivative of the integral over the constant portion of the WCA disperion interaction Uwca(x)
536 * = eps; x < rmin.
537 *
538 * @param ri The beginning of the integral.
539 * @param eps The well depth.
540 * @param rmin The Rmin value.
541 * @param rij The separation between the current atom and solvent blocking atom.
542 * @param rij2 The separation squared.
543 * @param rij3 The separation to the third.
544 * @param rho The scaled size of the solvent blocking atom.
545 * @param rho2 The scaled size of the solvent blocking atom squared.
546 * @param lik The beginning of the integral.
547 * @param lik2 The beginning of the integral squared.
548 * @param lik3 The beginning of the integral to the third.
549 * @param uik The end of the integral.
550 * @param uik2 The end of the integral squared.
551 * @param uik3 The end of the integral to the third.
552 * @return The derivative of the integral.
553 */
554 private double integralBeforeRminDerivative(
555 double ri,
556 double eps,
557 double rmin,
558 double rij,
559 double rij2,
560 double rij3,
561 double rho,
562 double rho2,
563 double lik,
564 double lik2,
565 double lik3,
566 double uik,
567 double uik2,
568 double uik3) {
569 double dl;
570 if (ri > rij - rho) {
571 dl = (-lik2 + 2.0 * rij2 + 2.0 * rho2) * lik2;
572 } else {
573 dl =
574 (-lik3 + 4.0 * lik2 * rij - 6.0 * lik * rij2 + 2.0 * lik * rho2 + 4.0 * rij3
575 - 4.0 * rij * rho2)
576 * lik;
577 }
578 double du;
579 if (rij + rho > rmin) {
580 du = -(-uik2 + 2.0 * rij2 + 2.0 * rho2) * uik2;
581 } else {
582 du =
583 -(-uik3 + 4.0 * uik2 * rij - 6.0 * uik * rij2 + 2.0 * uik * rho2 + 4.0 * rij3
584 - 4.0 * rij * rho2)
585 * uik;
586 }
587 return -eps * PI * (dl + du) / (4.0 * rij2);
588 }
589
590 /**
591 * @param eps The well depth.
592 * @param rmin7 The rmin value to the seventh.
593 * @param rij The separation between the current atom and solvent blocking atom.
594 * @param rij2 The separation squared.
595 * @param rho2 The scaled size of the solvent blocking atom squared.
596 * @param lik The beginning of the integral.
597 * @param lik2 The beginning of the integral squared.
598 * @param lik3 The beginning of the integral to the third.
599 * @param lik4 The beginning of the integral to the fourth.
600 * @param lik5 The beginning of the integral to the fifth.
601 * @param lik10 The beginning of the integral to the tenth.
602 * @param lik11 The beginning of the integral to the eleventh.
603 * @param lik12 The beginning of the integral to the twelfth.
604 * @param uik The end of the integral.
605 * @param uik2 The end of the integral squared.
606 * @param uik3 The end of the integral to the third.
607 * @param uik4 The end of the integral to the fourth.
608 * @param uik5 The end of the integral to the fifth.
609 * @param uik10 The end of the integral to the tenth.
610 * @param uik11 The end of the integral to the eleventh.
611 * @param uik12 The end of the integral to the twelfth.
612 * @return The value of the integral.
613 */
614 private double integratlAfterRmin(
615 double eps,
616 double rmin7,
617 double rij,
618 double rij2,
619 double rho2,
620 double lik,
621 double lik2,
622 double lik3,
623 double lik4,
624 double lik5,
625 double lik10,
626 double lik11,
627 double lik12,
628 double uik,
629 double uik2,
630 double uik3,
631 double uik4,
632 double uik5,
633 double uik10,
634 double uik11,
635 double uik12) {
636 double er7 = eps * rmin7;
637 double term = (15.0 * uik * lik * rij * (uik4 - lik4)
638 - 10.0 * uik2 * lik2 * (uik3 - lik3)
639 + 6.0 * (rho2 - rij2) * (uik5 - lik5)) / (120.0 * rij * lik5 * uik5);
640 double term2 = (120.0 * uik * lik * rij * (uik11 - lik11)
641 - 66.0 * uik2 * lik2 * (uik10 - lik10)
642 + 55.0 * (rho2 - rij2) * (uik12 - lik12)) / (2640.0 * rij * lik12 * uik12);
643 double idisp = -2.0 * er7 * term;
644 double irep = er7 * rmin7 * term2;
645 return 4.0 * PI * (irep + idisp);
646 }
647
648 /**
649 * @param ri The beginning of the integral.
650 * @param eps The eps value.
651 * @param rmin The rmin value to the seventh.
652 * @param rmin7 The rmin value to the seventh.
653 * @param rij The separation between the current atom and solvent blocking atom.
654 * @param rij2 The separation squared.
655 * @param rij3 The separation cubed.
656 * @param rho The scaled size of the solvent blocking atom.
657 * @param rho2 The scaled size of the solvent blocking atom squared.
658 * @param lik The beginning of the integral.
659 * @param lik2 The beginning of the integral squared.
660 * @param lik3 The beginning of the integral to the third.
661 * @param lik5 The beginning of the integral to the fifth.
662 * @param lik6 The beginning of the integral to the sixth.
663 * @param lik12 The beginning of the integral to the twelfth.
664 * @param lik13 The beginning of the integral to the thirteenth.
665 * @param uik The end of the integral.
666 * @param uik2 The end of the integral squared.
667 * @param uik3 The end of the integral to the third.
668 * @param uik6 The end of the integral to the sixth.
669 * @param uik13 The end of the integral to the thirteenth.
670 * @return The value of the integral derivative.
671 */
672 private double integratlAfterRminDerivative(
673 double ri,
674 double eps,
675 double rmin,
676 double rmin7,
677 double rmax,
678 double rij,
679 double rij2,
680 double rij3,
681 double rho,
682 double rho2,
683 double lik,
684 double lik2,
685 double lik3,
686 double lik5,
687 double lik6,
688 double lik12,
689 double lik13,
690 double uik,
691 double uik2,
692 double uik3,
693 double uik6,
694 double uik13) {
695 double er7 = eps * rmin7;
696 double lowerTerm = lik2 * rij + rij3 - rij * rho2;
697 double upperTerm = uik2 * rij + rij3 - rij * rho2;
698
699 double dl;
700 if (ri > rij - rho || rmax < rmin) {
701 dl = -(-5.0 * lik2 + 3.0 * rij2 + 3.0 * rho2) / lik5;
702 } else {
703 dl = (5.0 * lik3 - 33.0 * lik * rij2 - 3.0 * lik * rho2 + 15.0 * lowerTerm) / lik6;
704 }
705 double du = -(5.0 * uik3 - 33.0 * uik * rij2 - 3.0 * uik * rho2 + 15.0 * upperTerm) / uik6;
706 double de = -2.0 * PI * er7 * (dl + du) / (15.0 * rij2);
707
708 if (ri > rij - rho || rmax < rmin) {
709 dl = -(-6.0 * lik2 + 5.0 * rij2 + 5.0 * rho2) / lik12;
710 } else {
711 dl = (6.0 * lik3 - 125.0 * lik * rij2 - 5.0 * lik * rho2 + 60.0 * lowerTerm) / lik13;
712 }
713 du = -(6.0 * uik3 - 125.0 * uik * rij2 - 5.0 * uik * rho2 + 60.0 * upperTerm) / uik13;
714 de += PI * er7 * rmin7 * (dl + du) / (60.0 * rij2);
715
716 return de;
717 }
718 }
719 }