View Javadoc
1   package ffx.numerics.multipole;
2   
3   import static java.lang.Math.abs;
4   import static org.apache.commons.math3.util.FastMath.exp;
5   
6   /**
7    * The AmoebaPlusDampTensorGlobal class computes derivatives of overlap via recursion to
8    * order <= 6 for Cartesian multipoles defined in AMOEBA+. These are the multipole overlap
9    * interactions. The zeroth order overlap term is 1-A*exp(-alpha1*r)-B*exp(-alpha2*r).
10   *
11   * @author Matthew J. Speranza
12   * @see <a href="http://doi.org/10.1142/9789812830364_0002" target="_blank"> Matt Challacombe, Eric
13   * Schwegler and Jan Almlof, Modern developments in Hartree-Fock theory: Fast methods for
14   * computing the Coulomb matrix. Computational Chemistry: Review of Current Trends. pp. 53-107,
15   * Ed. J. Leczszynski, World Scientifc, 1996. </a>
16   * @see <a href="https://doi.org/10.1039/C6CP06017J" target="_blank"> Rackers, Wang, Liu,
17   * Piquemal, Ren and Ponder, An optimized charge penetration model for use with the AMOEBA
18   * force field. Phys. Chem. Chem. Phys., 2017,19, 276-291. </a>
19   * @since 1.0
20   */
21  public class AmoebaPlusOverlapTensorGlobal extends CoulombTensorGlobal {
22  
23    private double aI;
24    private double aK;
25    private double A;
26    private double B;
27  
28  
29    /**
30     * Constructor for CoulombTensorGlobal.
31     *
32     * @param order  The order of the tensor.
33     * @param alphaI The Thole damping parameter for site I.
34     * @param alphaK The Thole damping parameter for site K.
35     */
36    public AmoebaPlusOverlapTensorGlobal(int order, double alphaI, double alphaK) {
37      super(order);
38      this.operator = Operator.AMOEBA_PLUS_OVERLAP_FIELD;
39      this.aI = alphaI;
40      this.aK = alphaK;
41      this.A = aK * aK / ((aK * aK) - (aI * aI));
42      this.B = aI * aI / ((aI * aI) - (aK * aK));
43      // Multipole-Multipole overlap only defined up to order 6 (OST - Quadrupole-Quadrupole)
44      assert (order <= 6);
45    }
46  
47    /**
48     * Generate source terms for the Challacombe et al. recursion.
49     *
50     * @param T000 Location to store the source terms.
51     */
52    @Override
53    protected void source(double[] T000) {
54      // Compute the normal Coulomb auxiliary term.
55      super.source(T000);
56  
57      // Add the Thole damping terms: edamp = exp(-thole*u^3).
58      overlapSource(aI, aK, A, B, R, T000);
59    }
60  
61    /**
62     * Generate source terms for the Challacombe et al. recursion.
63     *
64     * @param aI   The Thole damping parameter for site I.
65     * @param aK   The Thole damping parameter for site K.
66     * @param A    The combined damping parameter for site I.
67     * @param B    The combined damping parameter for site K.
68     * @param R    The separation distance.
69     * @param T000 Location to store the source terms.
70     */
71    protected static void overlapSource(double aI, double aK, double A, double B, double R, double[] T000) {
72      // Add the damping terms: edamp = 1-A*exp(-alphaI*r)-B*exp(-alphaK*r).
73      // 1/r^11 terms from tinker code & are not included in supp info of paper
74      if (abs(aI - aK) < 1e-7) { // ai = ak
75        double alpha = aI;
76        double alphaR = alpha * R;
77        double alphaR2 = alphaR * alphaR;
78        double alphaR4 = alphaR2 * alphaR2;
79        double expAR = exp(-alphaR);
80  
81        T000[0] *= 1 - (1 + alphaR / 2.0) * expAR; // 1/r
82        T000[1] *= 1 - (1 + alphaR + alphaR2 / 2.0) * expAR; // 1/r^3
83        T000[2] *= 1 - (1.0 + alphaR + alphaR2 / 2.0 + alphaR2 * alphaR / 6.0) * expAR; // 1/r^5
84        T000[3] *= 1 - (1.0 + alphaR + alphaR2 / 2.0 + alphaR2 * alphaR / 6.0 +
85            alphaR2 * alphaR2 / 30.0) * expAR; // 1/r^7
86        T000[4] *= 1 - (1.0 + alphaR + alphaR2 / 2.0 + alphaR2 * alphaR / 6.0 +
87            4.0 * alphaR4 / 105.0 + alphaR4 * alphaR / 210.0) * expAR; // 1/r^9
88        T000[5] *= 1 - (1.0 + alphaR + alphaR2 * .5 + alphaR2 * alphaR / 6.0 +
89            5.0 * alphaR4 / 126.0 + 2.0 * alphaR4 * alphaR / 315.0 +
90            alphaR4 * alphaR2 / 1890.0) * expAR; // 1/r^11
91  
92      } else { // ai != ak
93        assert (A != Double.POSITIVE_INFINITY && B != Double.POSITIVE_INFINITY);
94        double aIR = aI * R;
95        double aIR2 = aIR * aIR;
96        double aKR = aK * R;
97        double aKR2 = aKR * aKR;
98        double expAI = A * exp(-aIR);
99        double expAK = B * exp(-aKR);
100 
101       T000[0] *= 1 - expAI - expAK; // 1/r
102       T000[1] *= 1 - (1 + aIR) * expAI - (1 + aKR) * expAK; // 1/r^3
103       T000[2] *= 1 - (1.0 + aIR + aIR2 / 3.0) * expAI - (1.0 + aKR + aKR2 / 3.0) * expAK; // 1/r^5
104       T000[3] *= 1 - (1.0 + aIR + 2.0 * aIR2 / 5.0 + aIR * aIR2 / 15.0) * expAI - // 1/r^7
105           (1.0 + aKR + 2.0 * aKR2 / 5.0 + aKR * aKR2 / 15.0) * expAK;
106       T000[4] *= 1 - (1.0 + aIR + 3.0 * aIR2 / 7.0 + 2.0 * aIR * aIR2 / 21.0 + aIR2 * aIR2 / 105.0) * expAI - // 1/r^9
107           (1.0 + aKR + 3.0 * aKR2 / 7.0 + 2.0 * aKR * aKR2 / 21.0 + aKR2 * aKR2 / 105.0) * expAK;
108       T000[5] *= 1.0 - (1.0 + aIR + 4.0 * aIR2 / 9.0 + aIR2 * aI / 9.0 +
109           aIR2 * aIR2 / 63.0 + aIR2 * aIR2 * aIR / 945.0) * expAI
110           - (1.0 + aKR + 4.0 * aKR2 / 9.0 + aKR2 * aKR / 9.0 + aKR2 * aKR2 / 63.0 + aKR2 * aKR2 * aKR / 945.0) * expAK; // TODO: Get 6th order terms
111     }
112   }
113 }