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.pme;
39
40 import static ffx.numerics.special.Erf.erfc;
41 import static org.apache.commons.math3.util.FastMath.pow;
42 import static org.apache.commons.math3.util.FastMath.sqrt;
43
44 import ffx.utilities.FFXProperty;
45 import ffx.utilities.PropertyGroup;
46
47 /**
48 * Mutable Particle Mesh Ewald constants.
49 */
50 public class EwaldParameters {
51
52 /**
53 * Default cutoff values for PME under periodic boundary conditions.
54 */
55 public static final double DEFAULT_EWALD_CUTOFF = 7.0;
56 /**
57 * The sqrt of PI.
58 */
59 private static final double SQRT_PI = sqrt(Math.PI);
60
61 /**
62 * The default Ewald coefficient.
63 * <br>
64 * For charged systems, the converged Ewald electrostatic energy is a function of the
65 * Ewald coefficient. For this reason, we've chosen to use a default value of 0.545,
66 * for all real space Ewald cutoff values.
67 * <br>
68 * In this way, systemically more accurate values for the real space cutoff,
69 * b-spline order and reciprocal space grid will converge the total electrostatic energy.
70 */
71 public static final double DEFAULT_EWALD_COEFFICIENT = 0.545;
72
73 @FFXProperty(name = "ewald-alpha", propertyGroup = PropertyGroup.ParticleMeshEwald, defaultValue = "0.545",
74 description = """
75 Sets the value of the Ewald coefficient, which controls the width of the Gaussian screening charges during
76 particle mesh Ewald summation for multipole electrostatics. In the absence of the ewald-alpha keyword,
77 the default value is 0.545, which is appropriate for most applications.
78 """)
79 public double aewald;
80 public double aewald3;
81 public double an0;
82 public double an1;
83 public double an2;
84 public double an3;
85 public double an4;
86 public double an5;
87 public double off;
88 public double off2;
89
90 public EwaldParameters(double cutoff, double aewald) {
91 setEwaldParameters(cutoff, aewald);
92 }
93
94 /**
95 * Determine the real space Ewald parameters and permanent multipole self energy.
96 *
97 * @param off Real space cutoff.
98 * @param aewald Ewald convergence parameter (0.0 turns off reciprocal space).
99 */
100 public void setEwaldParameters(double off, double aewald) {
101 this.off = off;
102 this.aewald = aewald;
103 off2 = off * off;
104 double alsq2 = 2.0 * aewald * aewald;
105 double piEwald = Double.POSITIVE_INFINITY;
106 if (aewald > 0.0) {
107 piEwald = 1.0 / (SQRT_PI * aewald);
108 }
109 aewald3 = 4.0 / 3.0 * pow(aewald, 3.0) / SQRT_PI;
110 if (aewald > 0.0) {
111 an0 = alsq2 * piEwald;
112 an1 = alsq2 * an0;
113 an2 = alsq2 * an1;
114 an3 = alsq2 * an2;
115 an4 = alsq2 * an3;
116 an5 = alsq2 * an4;
117 } else {
118 an0 = 0.0;
119 an1 = 0.0;
120 an2 = 0.0;
121 an3 = 0.0;
122 an4 = 0.0;
123 an5 = 0.0;
124 }
125 }
126
127 /**
128 * A precision of 1.0e-8 results in an Ewald coefficient that ensures continuity in the real space
129 * gradient, but at the cost of increased amplitudes for high frequency reciprocal space structure
130 * factors.
131 */
132 private double ewaldCoefficient(double cutoff, double precision) {
133
134 double eps = 1.0e-8;
135 if (precision < 1.0e-1) {
136 eps = precision;
137 }
138
139 /*
140 * Get an approximate value from cutoff and tolerance.
141 */
142 double ratio = eps + 1.0;
143 double x = 0.5;
144 int i = 0;
145 // Larger values lead to a more "delta-function-like" Gaussian
146 while (ratio >= eps) {
147 i++;
148 x *= 2.0;
149 ratio = erfc(x * cutoff) / cutoff;
150 }
151 /*
152 * Use a binary search to refine the coefficient.
153 */
154 int k = i + 60;
155 double xlo = 0.0;
156 double xhi = x;
157 for (int j = 0; j < k; j++) {
158 x = (xlo + xhi) / 2.0;
159 ratio = erfc(x * cutoff) / cutoff;
160 if (ratio >= eps) {
161 xlo = x;
162 } else {
163 xhi = x;
164 }
165 }
166
167 return x;
168 }
169
170 /**
171 * Determine the Ewald real space cutoff given the Ewald coefficient and a target precision.
172 *
173 * @param coeff The Ewald coefficient in use.
174 * @param maxCutoff The maximum cutoff.
175 * @param eps The target precision.
176 * @return The determined real space Ewald cutoff.
177 */
178 public static double ewaldCutoff(double coeff, double maxCutoff, double eps) {
179 // Set the tolerance value; use of 1.0d-8 requires strict convergence of the real Space sum.
180 double ratio = erfc(coeff * maxCutoff) / maxCutoff;
181
182 if (ratio > eps) {
183 return maxCutoff;
184 }
185
186 // Use a binary search to refine the coefficient.
187 double xlo = 0.0;
188 double xhi = maxCutoff;
189 double cutoff = 0.0;
190 for (int j = 0; j < 100; j++) {
191 cutoff = (xlo + xhi) / 2.0;
192 ratio = erfc(coeff * cutoff) / cutoff;
193 if (ratio >= eps) {
194 xlo = cutoff;
195 } else {
196 xhi = cutoff;
197 }
198 }
199 return cutoff;
200 }
201 }