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 jdk.incubator.vector.DoubleVector;
41 import org.junit.Test;
42 import org.junit.runner.RunWith;
43 import org.junit.runners.Parameterized;
44 import org.junit.runners.Parameterized.Parameters;
45
46 import java.util.Arrays;
47 import java.util.Collection;
48
49 import static ffx.numerics.multipole.MultipoleTensorTest.Qi;
50 import static ffx.numerics.multipole.MultipoleTensorTest.Qk;
51 import static ffx.numerics.multipole.MultipoleTensorTest.Ui;
52 import static ffx.numerics.multipole.MultipoleTensorTest.UiEwald;
53 import static ffx.numerics.multipole.MultipoleTensorTest.Uk;
54 import static ffx.numerics.multipole.MultipoleTensorTest.UkEwald;
55 import static ffx.numerics.multipole.MultipoleTensorTest.permTorqueI;
56 import static ffx.numerics.multipole.MultipoleTensorTest.permTorqueK;
57 import static ffx.numerics.multipole.MultipoleTensorTest.permanentEnergy;
58 import static ffx.numerics.multipole.MultipoleTensorTest.polarGradICoulomb;
59 import static ffx.numerics.multipole.MultipoleTensorTest.polarTorqueICoulomb;
60 import static ffx.numerics.multipole.MultipoleTensorTest.polarTorqueKCoulomb;
61 import static ffx.numerics.multipole.MultipoleTensorTest.polarizationEnergyCoulomb;
62 import static ffx.numerics.multipole.MultipoleTensorTest.scaleMutual;
63 import static java.lang.String.format;
64 import static org.junit.Assert.assertEquals;
65
66
67
68
69
70
71
72 @RunWith(Parameterized.class)
73 public class MultipoleTensorQISIMDTest {
74 private final double tolerance = 1.0e-13;
75 private final double fdTolerance = 1.0e-6;
76 private final double[] r = new double[3];
77 private DoubleVector xVector;
78 private DoubleVector yVector;
79 private DoubleVector zVector;
80 private DoubleVector[] rVectors;
81 private DoubleVector[] Fi;
82 private DoubleVector[] Fk;
83 private DoubleVector[] Ti;
84 private DoubleVector[] Tk;
85 private double[][] qI;
86 private double[][] uI;
87 private double[][] qK;
88 private double[][] uK;
89 private final int order;
90 private final String info;
91 private final Operator operator;
92
93 public MultipoleTensorQISIMDTest(String info, int order, Operator operator) {
94 this.info = info;
95 this.order = order;
96 r[0] = MultipoleTensorTest.xyz[0];
97 r[1] = MultipoleTensorTest.xyz[1];
98 r[2] = MultipoleTensorTest.xyz[2];
99 xVector = DoubleVector.broadcast(DoubleVector.SPECIES_PREFERRED, r[0]);
100 yVector = DoubleVector.broadcast(DoubleVector.SPECIES_PREFERRED, r[1]);
101 zVector = DoubleVector.broadcast(DoubleVector.SPECIES_PREFERRED, r[2]);
102 rVectors = new DoubleVector[]{xVector, yVector, zVector};
103 int tensorCount = MultipoleUtilities.tensorCount(order);
104 this.operator = operator;
105
106
107 qI = new double[10][DoubleVector.SPECIES_PREFERRED.length()];
108 uI = new double[3][DoubleVector.SPECIES_PREFERRED.length()];
109 qK = new double[10][DoubleVector.SPECIES_PREFERRED.length()];
110 uK = new double[3][DoubleVector.SPECIES_PREFERRED.length()];
111 for (int i = 0; i < DoubleVector.SPECIES_PREFERRED.length(); i++) {
112 for (int j = 0; j < 10; j++) {
113 qI[j][i] = Qi[j];
114 qK[j][i] = Qk[j];
115 }
116 for (int j = 0; j < 3; j++) {
117 uI[j][i] = Ui[j];
118 uK[j][i] = Uk[j];
119 }
120 }
121 Fi = new DoubleVector[3];
122 Fk = new DoubleVector[3];
123 Ti = new DoubleVector[3];
124 Tk = new DoubleVector[3];
125 for (int i = 0; i < 3; i++) {
126 Fi[i] = DoubleVector.broadcast(DoubleVector.SPECIES_PREFERRED, 0.0);
127 Fk[i] = DoubleVector.broadcast(DoubleVector.SPECIES_PREFERRED, 0.0);
128 Ti[i] = DoubleVector.broadcast(DoubleVector.SPECIES_PREFERRED, 0.0);
129 Tk[i] = DoubleVector.broadcast(DoubleVector.SPECIES_PREFERRED, 0.0);
130 }
131 }
132
133 @Parameters
134 public static Collection<Object[]> data() {
135 return Arrays.asList(new Object[][]{{"Order 5 Coulomb", (Object) 5, Operator.COULOMB}});
136 }
137
138
139 @Test
140 public void permanentMultipoleEnergyAndGradTest() {
141
142 if (operator == Operator.THOLE_FIELD) {
143 return;
144 }
145
146 MultipoleTensorSIMD multipoleTensor = new CoulombTensorQISIMD(order);
147 PolarizableMultipoleSIMD mI = new PolarizableMultipoleSIMD(qI, uI, uI);
148 PolarizableMultipoleSIMD mK = new PolarizableMultipoleSIMD(qK, uK, uK);
149 QIFrameSIMD qiFrame = new QIFrameSIMD();
150 qiFrame.setAndRotate(rVectors, mI, mK);
151 multipoleTensor.setR(xVector, yVector, zVector);
152 multipoleTensor.generateTensor();
153 DoubleVector e = multipoleTensor.multipoleEnergyAndGradient(mI, mK, Fi, Fk, Ti, Tk);
154 qiFrame.toGlobal(Fi);
155 qiFrame.toGlobal(Fk);
156 qiFrame.toGlobal(Ti);
157 qiFrame.toGlobal(Tk);
158 double[] eArray = new double[DoubleVector.SPECIES_PREFERRED.length()];
159 double[][] torI = new double[3][DoubleVector.SPECIES_PREFERRED.length()];
160 double[][] torK = new double[3][DoubleVector.SPECIES_PREFERRED.length()];
161 for (int i = 0; i < 3; i++) {
162 torI[i] = Ti[i].toArray();
163 torK[i] = Tk[i].toArray();
164 }
165 e.intoArray(eArray, 0);
166 for (int i = 0; i < DoubleVector.SPECIES_PREFERRED.length(); i++) {
167 assertEquals(format("Permanent energy %s %d", info, i), permanentEnergy, eArray[i], tolerance);
168 assertEquals(info + " Cartesian Multipole x-Torque I - Lane " + i, permTorqueI[0], torI[0][i], tolerance);
169 assertEquals(info + " Cartesian Multipole x-Torque K - Lane " + i, permTorqueK[0], torK[0][i], tolerance);
170 assertEquals(info + " Cartesian Multipole y-Torque I - Lane " + i, permTorqueI[1], torI[1][i], tolerance);
171 assertEquals(info + " Cartesian Multipole y-Torque K - Lane " + i, permTorqueK[1], torK[1][i], tolerance);
172 assertEquals(info + " Cartesian Multipole z-Torque I - Lane " + i, permTorqueI[2], torI[2][i], tolerance);
173 assertEquals(info + " Cartesian Multipole z-Torque K - Lane " + i, permTorqueK[2], torK[2][i], tolerance);
174 }
175
176
177 DoubleVector aX = Fk[0];
178 DoubleVector aY = Fk[1];
179 DoubleVector aZ = Fk[2];
180
181
182
183 double delta = 1.0e-5;
184 double delta2 = 2.0 * 1.0e-5;
185 xVector = xVector.add(delta);
186 multipoleTensor.setR(xVector, yVector, zVector);
187 multipoleTensor.generateTensor();
188 mI = new PolarizableMultipoleSIMD(qI, uI, uI);
189 mK = new PolarizableMultipoleSIMD(qK, uK, uK);
190 rVectors = new DoubleVector[]{xVector, yVector, zVector};
191 qiFrame.setAndRotate(rVectors, mI, mK);
192 DoubleVector posX = multipoleTensor.multipoleEnergy(mI, mK);
193 xVector = xVector.sub(delta2);
194 multipoleTensor.setR(xVector, yVector, zVector);
195 multipoleTensor.generateTensor();
196 mI = new PolarizableMultipoleSIMD(qI, uI, uI);
197 mK = new PolarizableMultipoleSIMD(qK, uK, uK);
198 rVectors = new DoubleVector[]{xVector, yVector, zVector};
199 qiFrame.setAndRotate(rVectors, mI, mK);
200 DoubleVector negX = multipoleTensor.multipoleEnergy(mI, mK);
201 xVector = xVector.add(delta);
202
203
204 yVector = yVector.add(delta);
205 multipoleTensor.setR(xVector, yVector, zVector);
206 multipoleTensor.generateTensor();
207 mI = new PolarizableMultipoleSIMD(qI, uI, uI);
208 mK = new PolarizableMultipoleSIMD(qK, uK, uK);
209 rVectors = new DoubleVector[]{xVector, yVector, zVector};
210 qiFrame.setAndRotate(rVectors, mI, mK);
211 DoubleVector posY = multipoleTensor.multipoleEnergy(mI, mK);
212 yVector = yVector.sub(delta2);
213 multipoleTensor.setR(xVector, yVector, zVector);
214 multipoleTensor.generateTensor();
215 mI = new PolarizableMultipoleSIMD(qI, uI, uI);
216 mK = new PolarizableMultipoleSIMD(qK, uK, uK);
217 rVectors = new DoubleVector[]{xVector, yVector, zVector};
218 qiFrame.setAndRotate(rVectors, mI, mK);
219 DoubleVector negY = multipoleTensor.multipoleEnergy(mI, mK);
220 yVector = yVector.add(delta);
221
222
223 zVector = zVector.add(delta);
224 multipoleTensor.setR(xVector, yVector, zVector);
225 multipoleTensor.generateTensor();
226 mI = new PolarizableMultipoleSIMD(qI, uI, uI);
227 mK = new PolarizableMultipoleSIMD(qK, uK, uK);
228 rVectors = new DoubleVector[]{xVector, yVector, zVector};
229 qiFrame.setAndRotate(rVectors, mI, mK);
230 DoubleVector posZ = multipoleTensor.multipoleEnergy(mI, mK);
231 zVector = zVector.sub(delta2);
232 multipoleTensor.setR(xVector, yVector, zVector);
233 multipoleTensor.generateTensor();
234 mI = new PolarizableMultipoleSIMD(qI, uI, uI);
235 mK = new PolarizableMultipoleSIMD(qK, uK, uK);
236 rVectors = new DoubleVector[]{xVector, yVector, zVector};
237 qiFrame.setAndRotate(rVectors, mI, mK);
238 DoubleVector negZ = multipoleTensor.multipoleEnergy(mI, mK);
239 zVector = zVector.add(delta);
240
241 double[] expectArray = aX.toArray();
242 DoubleVector actual = posX.sub(negX).div(delta2);
243 double[] actualArray = actual.toArray();
244 double[] expectYArray = aY.toArray();
245 DoubleVector actualY = posY.sub(negY).div(delta2);
246 double[] actualYArray = actualY.toArray();
247 double[] expectZArray = aZ.toArray();
248 DoubleVector actualZ = posZ.sub(negZ).div(delta2);
249 double[] actualZArray = actualZ.toArray();
250
251 for (int i = 0; i < DoubleVector.SPECIES_PREFERRED.length(); i++) {
252 assertEquals(format("Permanent Grad %s - Lane %d", info, i), expectArray[i], actualArray[i], fdTolerance);
253 assertEquals(format("Permanent Grad %s - Lane %d", info, i), expectYArray[i], actualYArray[i], fdTolerance);
254 assertEquals(format("Permanent Grad %s - Lane %d", info, i), expectZArray[i], actualZArray[i], fdTolerance);
255 }
256 }
257
258 @Test
259 public void polarizationEnergyAndGradientTest() {
260 MultipoleTensorSIMD multipoleTensor = new CoulombTensorQISIMD(order);
261
262 PolarizableMultipoleSIMD mI = new PolarizableMultipoleSIMD(qI, uI, uI);
263 PolarizableMultipoleSIMD mK = new PolarizableMultipoleSIMD(qK, uK, uK);
264 double[][] uIEwald = new double[3][DoubleVector.SPECIES_PREFERRED.length()];
265 double[][] uKEwald = new double[3][DoubleVector.SPECIES_PREFERRED.length()];
266 for (int i = 0; i < 3; i++) {
267 for (int j = 0; j < DoubleVector.SPECIES_PREFERRED.length(); j++) {
268 uIEwald[i][j] = UiEwald[i];
269 uKEwald[i][j] = UkEwald[i];
270 }
271 }
272 if (operator == Operator.SCREENED_COULOMB) {
273 mI.setInducedDipole(uIEwald, uIEwald);
274 mK.setInducedDipole(uKEwald, uKEwald);
275 }
276
277 QIFrameSIMD qiFrame = new QIFrameSIMD();
278 qiFrame.setAndRotate(rVectors, mI, mK);
279 multipoleTensor.setR(xVector, yVector, zVector);
280 multipoleTensor.generateTensor();
281 DoubleVector e = multipoleTensor.polarizationEnergyAndGradient(mI, mK,
282 DoubleVector.broadcast(DoubleVector.SPECIES_PREFERRED, 1.0),
283 DoubleVector.broadcast(DoubleVector.SPECIES_PREFERRED, 1.0),
284 DoubleVector.broadcast(DoubleVector.SPECIES_PREFERRED, scaleMutual),
285 Fi, Ti, Tk);
286 qiFrame.toGlobal(Fi);
287 qiFrame.toGlobal(Ti);
288 qiFrame.toGlobal(Tk);
289
290 double[] eArray = e.toArray();
291
292 double[] aX = Fi[0].toArray();
293 double[] aY = Fi[1].toArray();
294 double[] aZ = Fi[2].toArray();
295 double[] axTorqueI = Ti[0].toArray();
296 double[] axTorqueK = Tk[0].toArray();
297 double[] ayTorqueI = Ti[1].toArray();
298 double[] ayTorqueK = Tk[1].toArray();
299 double[] azTorqueI = Ti[2].toArray();
300 double[] azTorqueK = Tk[2].toArray();
301
302 double[] gradI = polarGradICoulomb;
303 double[] torqueI = polarTorqueICoulomb;
304 double[] torqueK = polarTorqueKCoulomb;
305 for (int i = 0; i < DoubleVector.SPECIES_PREFERRED.length(); i++) {
306 assertEquals(info + " Polarization Energy", polarizationEnergyCoulomb, eArray[i], tolerance);
307 assertEquals(info + " Polarization GradX", gradI[0], aX[i], tolerance);
308 assertEquals(info + " Polarization GradY", gradI[1], aY[i], tolerance);
309 assertEquals(info + " Polarization GradZ", gradI[2], aZ[i], tolerance);
310 assertEquals(info + " Polarization TorqueX I", torqueI[0], axTorqueI[i], tolerance);
311 assertEquals(info + " Polarization TorqueY I", torqueI[1], ayTorqueI[i], tolerance);
312 assertEquals(info + " Polarization TorqueZ I", torqueI[2], azTorqueI[i], tolerance);
313 assertEquals(info + " Polarization TorqueX K", torqueK[0], axTorqueK[i], tolerance);
314 assertEquals(info + " Polarization TorqueY K", torqueK[1], ayTorqueK[i], tolerance);
315 assertEquals(info + " Polarization TorqueZ K", torqueK[2], azTorqueK[i], tolerance);
316 }
317
318
319 }
320
321 }