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.numerics.multipole;
39  
40  import static org.apache.commons.math3.util.FastMath.exp;
41  
42  /**
43   * The TholeTensorQI class computes derivatives of Thole damping via recursion to order <= 4 for
44   * Cartesian multipoles in a quasi-internal frame.
45   *
46   * @author Michael J. Schnieders
47   * @see <a href="http://doi.org/10.1142/9789812830364_0002" target="_blank"> Matt Challacombe, Eric
48   *     Schwegler and Jan Almlof, Modern developments in Hartree-Fock theory: Fast methods for
49   *     computing the Coulomb matrix. Computational Chemistry: Review of Current Trends. pp. 53-107,
50   *     Ed. J. Leczszynski, World Scientifc, 1996. </a>
51   * @since 1.0
52   */
53  public class TholeTensorQI extends CoulombTensorQI {
54  
55    /** Constant <code>threeFifths=3.0 / 5.0</code> */
56    private static final double threeFifths = 3.0 / 5.0;
57  
58    /** Constant <code>oneThirtyFifth=1.0 / 35.0</code> */
59    private static final double oneThirtyFifth = 1.0 / 35.0;
60  
61    /**
62     * Thole damping parameter is set to min(pti,ptk)).
63     */
64    private double thole;
65  
66    /**
67     * AiAk parameter = 1/(alphaI^6*alphaK^6) where alpha is polarizability.
68     */
69    private double AiAk;
70  
71    /**
72     * Constructor for TholeTensorQI.
73     *
74     * @param order Tensor order.
75     * @param thole Thole damping parameter is set to min(pti,ptk)).
76     * @param AiAk parameter = 1/(alphaI^6*alphaK^6) where alpha is polarizability.
77     */
78    public TholeTensorQI(int order, double thole, double AiAk) {
79      super(order);
80      this.thole = thole;
81      this.AiAk = AiAk;
82      this.operator = OPERATOR.THOLE_FIELD;
83  
84      assert (order <= 4);
85    }
86  
87    /**
88     * Set Thole damping parameters
89     *
90     * @param thole a double.
91     * @param AiAk a double.
92     */
93    public void setThole(double thole, double AiAk) {
94      this.thole = thole;
95      this.AiAk = AiAk;
96    }
97  
98    /**
99     * Check if the Thole damping is exponential is greater than zero (or the interaction can be
100    * neglected).
101    *
102    * @param r The separation distance.
103    * @return True if -thole*u^3 is greater than -50.0.
104    */
105   public boolean checkThole(double r) {
106     double rAiAk = r * AiAk;
107     return (-thole * rAiAk * rAiAk * rAiAk > -50.0);
108   }
109 
110   /**
111    * Generate source terms for the Challacombe et al. recursion.
112    *
113    * @param T000 Location to store the source terms.
114    */
115   protected void source(double[] T000) {
116     // Compute the normal Coulomb auxiliary term.
117     double ir = 1.0 / R;
118     double ir2 = ir * ir;
119     for (int n = 0; n < o1; n++) {
120       T000[n] = coulombSource[n] * ir;
121       ir *= ir2;
122     }
123 
124     // Add the Thole damping terms: edamp = exp(-thole*u^3).
125     double u = R * AiAk;
126     double u3 = thole * u * u * u;
127     double u6 = u3 * u3;
128     double u9 = u6 * u3;
129     double expU3 = exp(-u3);
130 
131     // The zeroth order term is not calculated for Thole damping.
132     T000[0] = 0.0;
133     T000[1] *= expU3;
134     T000[2] *= (1.0 + u3) * expU3;
135     T000[3] *= (1.0 + u3 + threeFifths * u6) * expU3;
136     T000[4] *= (1.0 + u3 + (18.0 * u6 + 9.0 * u9) * oneThirtyFifth) * expU3;
137   }
138 
139 }