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.switching;
39
40 import static java.lang.String.format;
41 import static org.apache.commons.math3.util.FastMath.max;
42 import static org.apache.commons.math3.util.FastMath.min;
43 import static org.apache.commons.math3.util.FastMath.pow;
44
45 /**
46 * The 6 coefficients of the multiplicative polynomial switch are unique given the distances "a" and
47 * "b". They are found by solving a system of 6 equations, which define the boundary conditions of
48 * the switch. <br> f(a) = 1 <br> f'(a) = f"(a) = 0 <br> f(b) = f'(b) = f"(b) = 0
49 *
50 * @author Michael J. Schnieders
51 */
52 public class MultiplicativeSwitch implements UnivariateSwitchingFunction {
53
54 private final double b;
55 private final double a;
56 private final double c0;
57 private final double c1;
58 private final double c2;
59 private final double c3;
60 private final double c4;
61 private final double c5;
62 private final double twoC2;
63 private final double threeC3;
64 private final double fourC4;
65 private final double fiveC5;
66
67 /**
68 * Constructs a MultiplicativeSwitch that starts at f(0)=1 and ends at f(1)=0. The switch smoothly
69 * interpolates from 1 to 0 across that range, with zero first and second derivatives at off and
70 * cut.
71 */
72 public MultiplicativeSwitch() {
73 this(0.0, 1.0);
74 }
75
76 /**
77 * Constructs a MultiplicativeSwitch that starts at f(a)=1 and ends at f(b)=0. The switch smoothly
78 * interpolates from 1 to 0 across that range, with zero first and second derivatives at off and
79 * cut.
80 *
81 * @param a f(a)=1
82 * @param b f(b)=0
83 */
84 public MultiplicativeSwitch(double a, double b) {
85
86 // f(a) = 1.0
87 this.a = a;
88 // f(b) = 0.0
89 this.b = b;
90
91 double a2 = a * a;
92 double b2 = b * b;
93
94 double denominator = pow(b - a, 5.0);
95 c0 = b * b2 * (b2 - 5.0 * a * b + 10.0 * a2) / denominator;
96 c1 = -30.0 * a2 * b2 / denominator;
97 c2 = 30.0 * b * a * (b + a) / denominator;
98 c3 = -10.0 * (a2 + 4.0 * a * b + b2) / denominator;
99 c4 = 15.0 * (a + b) / denominator;
100 c5 = -6.0 / denominator;
101 twoC2 = 2.0 * c2;
102 threeC3 = 3.0 * c3;
103 fourC4 = 4.0 * c4;
104 fiveC5 = 5.0 * c5;
105 }
106
107 /** {@inheritDoc} */
108 @Override
109 public boolean constantOutsideBounds() {
110 return false;
111 }
112
113 /**
114 * First derivative of the switching function at r.
115 *
116 * @param r r
117 * @param r2 r^2
118 * @param r3 r^3
119 * @param r4 r^4
120 * @return First derivative of switch at r
121 */
122 public double dtaper(double r, double r2, double r3, double r4) {
123 return fiveC5 * r4 + fourC4 * r3 + threeC3 * r2 + twoC2 * r + c1;
124 }
125
126 /**
127 * First derivative of the switching function at r.
128 *
129 * @param r r
130 * @return First derivative of switch at r
131 */
132 public double dtaper(double r) {
133 // Minimize number of multiply operations by storing r^2.
134 double r2 = r * r;
135 return dtaper(r, r2, r2 * r, r2 * r2);
136 }
137
138 /** {@inheritDoc} */
139 @Override
140 public double firstDerivative(double x) {
141 return dtaper(x);
142 }
143
144 /** {@inheritDoc} */
145 @Override
146 public int getHighestOrderZeroDerivative() {
147 return 2;
148 }
149
150 /** {@inheritDoc} */
151 @Override
152 public double getOneBound() {
153 return max(a, b);
154 }
155
156 /**
157 * Get the value where the switch starts.
158 *
159 * @return Switch start.
160 */
161 public double getSwitchEnd() {
162 return b;
163 }
164
165 /**
166 * Get the value where the switch starts.
167 *
168 * @return Switch start.
169 */
170 public double getSwitchStart() {
171 return a;
172 }
173
174 /** {@inheritDoc} */
175 @Override
176 public double getZeroBound() {
177 return min(b, a);
178 }
179
180 /** {@inheritDoc} */
181 @Override
182 public double nthDerivative(double x, int order) throws IllegalArgumentException {
183 if (order < 1) {
184 throw new IllegalArgumentException("Order must be >= 1");
185 }
186 switch (order) {
187 case 1:
188 return dtaper(x);
189 case 2:
190 return secondDerivative(x);
191 case 3:
192 double val = 60.0 * c5 * x * x;
193 val += 24.0 * c4 * x;
194 val += 6.0 * c3;
195 return val;
196 case 4:
197 val = 120.0 * c5 * x;
198 val += 24.0 * c4;
199 return val;
200 case 5:
201 return 120.0 * c5;
202 default:
203 return 0;
204 }
205 }
206
207 /** {@inheritDoc} */
208 @Override
209 public double secondDerivative(double x) {
210 double x2 = x * x;
211 double val = 20.0 * c5 * x2 * x;
212 val += 12.0 * c4 * x2;
213 val += 6.0 * c3 * x;
214 val += 2.0 * c2;
215 return val;
216 }
217
218 /** {@inheritDoc} */
219 @Override
220 public boolean symmetricToUnity() {
221 return true;
222 }
223
224 /**
225 * Value of the switching function at r.
226 *
227 * @param r r
228 * @param r2 r^2
229 * @param r3 r^3
230 * @param r4 r^4
231 * @param r5 r^5
232 * @return Value of switch at r
233 */
234 public double taper(double r, double r2, double r3, double r4, double r5) {
235 return c5 * r5 + c4 * r4 + c3 * r3 + c2 * r2 + c1 * r + c0;
236 }
237
238 /**
239 * Value of the switching function at r.
240 *
241 * @param r r
242 * @return Value of switch at r
243 */
244 public double taper(double r) {
245 // Minimize number of multiply operations by storing r^2, r^3.
246 double r2 = r * r;
247 double r3 = r2 * r;
248 return taper(r, r2, r3, r2 * r2, r3 * r2);
249 }
250
251 /** {@inheritDoc} */
252 @Override
253 public String toString() {
254 return format(
255 "Multiplicative switch of form f(x) = %8.4g*x^5 + "
256 + "%8.4g*x^4 + %8.4g*x^3 + %8.4g*x^2 + %8.4g*x + %8.4g",
257 c5, c4, c3, c2, c1, c0);
258 }
259
260 /** {@inheritDoc} */
261 @Override
262 public boolean validOutsideBounds() {
263 return false;
264 }
265
266 /** {@inheritDoc} */
267 @Override
268 public double valueAt(double x) throws IllegalArgumentException {
269 return taper(x);
270 }
271 }