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