1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38 package ffx.numerics.multipole;
39
40 import org.junit.Test;
41 import org.junit.runner.RunWith;
42 import org.junit.runners.Parameterized;
43 import org.junit.runners.Parameterized.Parameters;
44
45 import java.util.Arrays;
46 import java.util.Collection;
47 import java.util.logging.Logger;
48
49 import static ffx.numerics.multipole.MultipoleTensorTest.*;
50 import static java.lang.String.format;
51 import static org.apache.commons.math3.util.FastMath.pow;
52 import static org.junit.Assert.assertEquals;
53
54 @RunWith(Parameterized.class)
55 public class AMOEBAPlusMultipoleTensorGlobalTest {
56
57
58
59 private static final Logger logger = Logger.getLogger(AMOEBAPlusMultipoleTensorGlobalTest.class.getName());
60 private final double tolerance = 1.0e-13;
61 private final double fdTolerance = 1.0e-6;
62 private final double[] r = new double[3];
63 private final double[] r2 = new double[3];
64 private final int order;
65 private final String info;
66 private final Operator operator;
67
68 public AMOEBAPlusMultipoleTensorGlobalTest(
69 String info,
70 int order,
71 Operator operator) {
72 this.info = info;
73 this.order = order;
74 r[0] = QkXYZ[0] - QiXYZ[0];
75 r[1] = QkXYZ[1] - QiXYZ[1];
76 r[2] = QkXYZ[2] - QiXYZ[2];
77 r2[0] = QjXYZ[0] - QkXYZ[0];
78 r2[1] = QjXYZ[1] - QkXYZ[1];
79 r2[2] = QjXYZ[2] - QkXYZ[2];
80 this.operator = operator;
81 }
82
83 @Parameters
84 public static Collection<Object[]> data() {
85 return Arrays.asList(
86 new Object[][]{
87 {"Order 5 Aplus Direct", 5, Operator.AMOEBA_PLUS_OVERLAP_FIELD}
88 });
89 }
90
91 @Test
92 public void permanentMultipoleEnergyAndGradTest() {
93 double delta = 1.0e-8;
94 double delta2 = 2.0 * delta;
95 double[] Fi = new double[3];
96 double[] Fk = new double[3];
97 double[] Ti = new double[3];
98 double[] Tk = new double[3];
99
100
101 PolarizableMultipole mI = new PolarizableMultipole(QiAmoebaP, QiInduced, QiInduced, Zi);
102 PolarizableMultipole mK = new PolarizableMultipole(QkAmoebaP, QkInduced, QkInduced, Zk);
103 MultipoleTensor multipoleTensor = new AmoebaPlusOverlapTensorGlobal(order, chargePenAlphaOx, chargePenAlphaOx);
104 multipoleTensor.generateTensor(r);
105 double eOverlap = multipoleTensor.multipoleEnergyAndGradient(mI, mK, Fi, Fk, Ti, Tk);
106
107
108 double aX = Fk[0];
109 double aY = Fk[1];
110 double aZ = Fk[2];
111
112 r[0] += delta;
113 multipoleTensor.generateTensor(r);
114 double posX = multipoleTensor.multipoleEnergy(mI, mK);
115 r[0] -= delta2;
116 multipoleTensor.generateTensor(r);
117 double negX = multipoleTensor.multipoleEnergy(mI, mK);
118 r[0] += delta;
119 r[1] += delta;
120 multipoleTensor.generateTensor(r);
121 double posY = multipoleTensor.multipoleEnergy(mI, mK);
122 r[1] -= delta2;
123 multipoleTensor.generateTensor(r);
124 double negY = multipoleTensor.multipoleEnergy(mI, mK);
125 r[1] += delta;
126 r[2] += delta;
127 multipoleTensor.generateTensor(r);
128 double posZ = multipoleTensor.multipoleEnergy(mI, mK);
129 r[2] -= delta2;
130 multipoleTensor.generateTensor(r);
131 double negZ = multipoleTensor.multipoleEnergy(mI, mK);
132 r[2] += delta;
133
134 double expect = aX;
135 double actual = (posX - negX) / delta2;
136 assertEquals(info + " OO Overlap Force X", expect, actual, fdTolerance);
137 expect = aY;
138 actual = (posY - negY) / delta2;
139 assertEquals(info + " OO Overlap Force Y", expect, actual, fdTolerance);
140 expect = aZ;
141 actual = (posZ - negZ) / delta2;
142 assertEquals(info + " OO Overlap Force Z", expect, actual, fdTolerance);
143 Arrays.fill(Fi, 0);
144 Arrays.fill(Fk, 0);
145
146
147 AmoebaPlusDampTensorGlobal multipoleTensorDamp = new AmoebaPlusDampTensorGlobal(3, chargePenAlphaOx, chargePenAlphaOx);
148 multipoleTensorDamp.generateTensor(r);
149 double eProton = multipoleTensorDamp.coreInteractionAndGradient(mI, mK, Fi, Fk);
150 double e = eOverlap + eProton;
151 assertEquals(" Oxygen-Oxygen Perm Multipole Coulomb Energy", amoebaPlusMPoleEnergyOO, e, tolerance);
152
153
154 aX = Fk[0];
155 aY = Fk[1];
156 aZ = Fk[2];
157
158 r[0] += delta;
159 multipoleTensorDamp.generateTensor(r);
160 posX = multipoleTensorDamp.coreInteraction(mI, mK);
161 r[0] -= delta2;
162 multipoleTensorDamp.generateTensor(r);
163 negX = multipoleTensorDamp.coreInteraction(mI, mK);
164 r[0] += delta;
165 r[1] += delta;
166 multipoleTensorDamp.generateTensor(r);
167 posY = multipoleTensorDamp.coreInteraction(mI, mK);
168 r[1] -= delta2;
169 multipoleTensorDamp.generateTensor(r);
170 negY = multipoleTensorDamp.coreInteraction(mI, mK);
171 r[1] += delta;
172 r[2] += delta;
173 multipoleTensorDamp.generateTensor(r);
174 posZ = multipoleTensorDamp.coreInteraction(mI, mK);
175 r[2] -= delta2;
176 multipoleTensorDamp.generateTensor(r);
177 negZ = multipoleTensorDamp.coreInteraction(mI, mK);
178 r[2] += delta;
179
180 expect = aX;
181 actual = (posX - negX) / delta2;
182 assertEquals(info + " OO Core Force X", expect, actual, fdTolerance);
183 expect = aY;
184 actual = (posY - negY) / delta2;
185 assertEquals(info + " OO Core Force Y", expect, actual, fdTolerance);
186 expect = aZ;
187 actual = (posZ - negZ) / delta2;
188 assertEquals(info + " OO Core Force Z", expect, actual, fdTolerance);
189 Arrays.fill(Fi, 0);
190 Arrays.fill(Fk, 0);
191
192
193 mI = new PolarizableMultipole(QkAmoebaP, QkInduced, QkInduced, Zk);
194 mK = new PolarizableMultipole(QjAmoebaP, QjInduced, QjInduced, Zj);
195 multipoleTensor = new AmoebaPlusOverlapTensorGlobal(order, chargePenAlphaOx, chargePenAlphaHyd);
196 multipoleTensor.generateTensor(r2);
197 eOverlap = multipoleTensor.multipoleEnergyAndGradient(mI, mK, Fi, Fk, Ti, Tk);
198
199
200 aX = Fk[0];
201 aY = Fk[1];
202 aZ = Fk[2];
203
204 r2[0] += delta;
205 multipoleTensor.generateTensor(r2);
206 posX = multipoleTensor.multipoleEnergy(mI, mK);
207 r2[0] -= delta2;
208 multipoleTensor.generateTensor(r2);
209 negX = multipoleTensor.multipoleEnergy(mI, mK);
210 r2[0] += delta;
211 r2[1] += delta;
212 multipoleTensor.generateTensor(r2);
213 posY = multipoleTensor.multipoleEnergy(mI, mK);
214 r2[1] -= delta2;
215 multipoleTensor.generateTensor(r2);
216 negY = multipoleTensor.multipoleEnergy(mI, mK);
217 r2[1] += delta;
218 r2[2] += delta;
219 multipoleTensor.generateTensor(r2);
220 posZ = multipoleTensor.multipoleEnergy(mI, mK);
221 r2[2] -= delta2;
222 multipoleTensor.generateTensor(r2);
223 negZ = multipoleTensor.multipoleEnergy(mI, mK);
224 r2[2] += delta;
225
226 expect = aX;
227 actual = (posX - negX) / delta2;
228 assertEquals(info + " OH Overlap Force X", expect, actual, fdTolerance);
229 expect = aY;
230 actual = (posY - negY) / delta2;
231 assertEquals(info + " OH Overlap Force Y", expect, actual, fdTolerance);
232 expect = aZ;
233 actual = (posZ - negZ) / delta2;
234 assertEquals(info + " OH Overlap Force Z", expect, actual, fdTolerance);
235 Arrays.fill(Fi, 0);
236 Arrays.fill(Fk, 0);
237
238
239
240 multipoleTensorDamp = new AmoebaPlusDampTensorGlobal(3, chargePenAlphaOx, chargePenAlphaHyd);
241 multipoleTensorDamp.generateTensor(r2);
242
243 eProton = multipoleTensorDamp.coreInteractionAndGradient(mI, mK, Fi, Fk);
244 e = eOverlap + eProton;
245 assertEquals(" Hydrogen-Oxygen Perm Multipole Damped Coulomb Energy", amoebaPlusMPoleEnergyOH, e, tolerance);
246
247
248 aX = Fk[0];
249 aY = Fk[1];
250 aZ = Fk[2];
251
252 r2[0] += delta;
253 multipoleTensorDamp.generateTensor(r2);
254 posX = multipoleTensorDamp.coreInteraction(mI, mK);
255 r2[0] -= delta2;
256 multipoleTensorDamp.generateTensor(r2);
257 negX = multipoleTensorDamp.coreInteraction(mI, mK);
258 r2[0] += delta;
259 r2[1] += delta;
260 multipoleTensorDamp.generateTensor(r2);
261 posY = multipoleTensorDamp.coreInteraction(mI, mK);
262 r2[1] -= delta2;
263 multipoleTensorDamp.generateTensor(r2);
264 negY = multipoleTensorDamp.coreInteraction(mI, mK);
265 r2[1] += delta;
266 r2[2] += delta;
267 multipoleTensorDamp.generateTensor(r2);
268 posZ = multipoleTensorDamp.coreInteraction(mI, mK);
269 r2[2] -= delta2;
270 multipoleTensorDamp.generateTensor(r2);
271 negZ = multipoleTensorDamp.coreInteraction(mI, mK);
272 r2[2] += delta;
273
274 expect = aX;
275 actual = (posX - negX) / delta2;
276 assertEquals(info + " OH Core Force X", expect, actual, fdTolerance);
277 expect = aY;
278 actual = (posY - negY) / delta2;
279 assertEquals(info + " OH Core Force Y", expect, actual, fdTolerance);
280 expect = aZ;
281 actual = (posZ - negZ) / delta2;
282 assertEquals(info + " OH Core Force Z", expect, actual, fdTolerance);
283 }
284
285 @Test
286 public void inducedDipoleEnergyAndGradient() {
287 double[] Fi = new double[3];
288 double[] Ti = new double[3];
289 double[] Tk = new double[3];
290 double thole = .7;
291 double aiak = 1.0 / (0.99595940317287035 * 0.99595940317287035);
292
293
294 PolarizableMultipole mI = new PolarizableMultipole(QiAmoebaP, QiInduced, QiInduced, Zi);
295 PolarizableMultipole mK = new PolarizableMultipole(QkAmoebaP, QkInduced, QkInduced, Zk);
296 MultipoleTensor multipoleTensor = new TholeTensorGlobal(4, thole, aiak, true);
297 multipoleTensor.generateTensor(r);
298
299 double e = multipoleTensor.polarizationEnergyAndGradient(mI, mK, scaleMutual,
300 scaleMutual, 0.0, Fi, Ti, Tk);
301 assertEquals(" Oxygen-Oxygen Thole-Damped Polarization Energy", amoebaPlusIndDipoleEnergyOO, e, tolerance);
302 double[] aplusIndGrad = new double[]
303 {-7.8921873374791535E-004,
304 -6.3709119990021249E-004,
305 -2.4873367831826638E-021};
306 assertEquals(info + " Polarization GradX", aplusIndGrad[0], Fi[0], tolerance);
307 assertEquals(info + " Polarization GradY", aplusIndGrad[1], Fi[1], tolerance);
308 assertEquals(info + " Polarization GradZ", aplusIndGrad[2], Fi[2], tolerance);
309 }
310 }