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.potential.nonbonded.pme;
39
40 import ffx.potential.parameters.MultipoleType.MultipoleFrameDefinition;
41
42 import java.util.logging.Level;
43 import java.util.logging.Logger;
44
45 import static ffx.numerics.math.DoubleMath.X;
46 import static ffx.numerics.math.DoubleMath.add;
47 import static ffx.numerics.math.DoubleMath.dot;
48 import static ffx.numerics.math.DoubleMath.length;
49 import static ffx.numerics.math.DoubleMath.normalize;
50 import static ffx.numerics.math.DoubleMath.scale;
51 import static ffx.numerics.math.DoubleMath.sub;
52 import static ffx.potential.parameters.MultipoleType.MultipoleFrameDefinition.NONE;
53 import static ffx.potential.parameters.MultipoleType.MultipoleFrameDefinition.THREEFOLD;
54 import static ffx.potential.parameters.MultipoleType.MultipoleFrameDefinition.ZONLY;
55 import static ffx.potential.parameters.MultipoleType.MultipoleFrameDefinition.ZTHENBISECTOR;
56 import static org.apache.commons.math3.util.FastMath.abs;
57 import static org.apache.commons.math3.util.FastMath.sqrt;
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77 public class Torque {
78
79 private static final Logger logger = Logger.getLogger(Torque.class.getName());
80
81 private final double[] vecZ = new double[3];
82 private final double[] vecX = new double[3];
83 private final double[] vecY = new double[3];
84 private final double[] sumXY = new double[3];
85 private final double[] sumZXY = new double[3];
86 private final double[] r = new double[3];
87 private final double[] del = new double[3];
88 private final double[] eps = new double[3];
89 private final double[] xZXY = new double[3];
90 private final double[] xXZ = new double[3];
91 private final double[] xYZ = new double[3];
92 private final double[] xYX = new double[3];
93 private final double[] xXYZ = new double[3];
94 private final double[] xxZXYZ = new double[3];
95 private final double[] xxZXYX = new double[3];
96 private final double[] xxZXYY = new double[3];
97 private final double[] t1 = new double[3];
98 private final double[] t2 = new double[3];
99 private final double[] localOrigin = new double[3];
100
101 private int[][] axisAtom;
102 private MultipoleFrameDefinition[] frame;
103 private double[][][] coordinates;
104
105 public void init(int[][] axisAtom, MultipoleFrameDefinition[] frame, double[][][] coordinates) {
106 this.axisAtom = axisAtom;
107 this.frame = frame;
108 this.coordinates = coordinates;
109 }
110
111 public void torque(int i, int iSymm, double[] trq, int[] frameIndex, double[][] g) {
112 final int[] ax = axisAtom[i];
113
114 if (frame[i] == NONE) {
115 return;
116 }
117
118
119 int nAxisAtoms = ax.length;
120 final int iz = ax[0];
121 int ix = -1;
122 int iy = -1;
123 if (nAxisAtoms > 1) {
124 ix = ax[1];
125 if (nAxisAtoms > 2) {
126 iy = ax[2];
127 }
128 }
129 frameIndex[0] = iz;
130 frameIndex[1] = ix;
131 frameIndex[2] = iy;
132 frameIndex[3] = i;
133 double[] x = coordinates[iSymm][0];
134 double[] y = coordinates[iSymm][1];
135 double[] z = coordinates[iSymm][2];
136
137 localOrigin[0] = x[i];
138 localOrigin[1] = y[i];
139 localOrigin[2] = z[i];
140 vecZ[0] = x[iz];
141 vecZ[1] = y[iz];
142 vecZ[2] = z[iz];
143 sub(vecZ, localOrigin, vecZ);
144 double rZ = length(vecZ);
145 if (frame[i] != ZONLY) {
146 vecX[0] = x[ix];
147 vecX[1] = y[ix];
148 vecX[2] = z[ix];
149 sub(vecX, localOrigin, vecX);
150 } else {
151 vecX[0] = 1.0;
152 vecX[1] = 0.0;
153 vecX[2] = 0.0;
154 double dot = vecZ[0] / rZ;
155 if (abs(dot) > 0.866) {
156 vecX[0] = 0.0;
157 vecX[1] = 1.0;
158 }
159 }
160 double rX = length(vecX);
161 if (frame[i] == ZTHENBISECTOR || frame[i] == THREEFOLD) {
162 vecY[0] = x[iy];
163 vecY[1] = y[iy];
164 vecY[2] = z[iy];
165 sub(vecY, localOrigin, vecY);
166 } else {
167 X(vecZ, vecX, vecY);
168 }
169 double rY = length(vecY);
170 scale(vecZ, 1.0 / rZ, vecZ);
171 scale(vecX, 1.0 / rX, vecX);
172 scale(vecY, 1.0 / rY, vecY);
173
174 X(vecX, vecZ, xXZ);
175 X(vecY, vecZ, xYZ);
176 X(vecY, vecX, xYX);
177 normalize(xXZ, xXZ);
178 normalize(xYZ, xYZ);
179 normalize(xYX, xYX);
180
181 double cosZX = dot(vecZ, vecX);
182 double sinZX = sqrt(1.0 - cosZX * cosZX);
183
184
185
186
187 double dPhidZ = -(trq[0] * vecZ[0] + trq[1] * vecZ[1] + trq[2] * vecZ[2]);
188 double dPhidX = -(trq[0] * vecX[0] + trq[1] * vecX[1] + trq[2] * vecX[2]);
189 double dPhidY = -(trq[0] * vecY[0] + trq[1] * vecY[1] + trq[2] * vecY[2]);
190 switch (frame[i]) {
191 case ZONLY:
192 for (int j = 0; j < 3; j++) {
193 double dZ = xXZ[j] * dPhidX / (rZ * sinZX) + xYZ[j] * dPhidY / rZ;
194 g[0][j] = dZ;
195 g[3][j] = -dZ;
196 }
197 break;
198 case ZTHENX:
199 for (int j = 0; j < 3; j++) {
200 double dZ = xXZ[j] * dPhidX / (rZ * sinZX) + xYZ[j] * dPhidY / rZ;
201 double dX = -xXZ[j] * dPhidZ / (rX * sinZX);
202 g[0][j] = dZ;
203 g[1][j] = dX;
204 g[3][j] = -(dZ + dX);
205 }
206 break;
207 case BISECTOR:
208 for (int j = 0; j < 3; j++) {
209 double dZ = xXZ[j] * dPhidX / (rZ * sinZX) + 0.5 * xYZ[j] * dPhidY / rZ;
210 double dX = -xXZ[j] * dPhidZ / (rX * sinZX) + 0.5 * xYX[j] * dPhidY / rX;
211 g[0][j] = dZ;
212 g[1][j] = dX;
213 g[3][j] = -(dZ + dX);
214 }
215 break;
216 case ZTHENBISECTOR:
217
218 add(vecX, vecY, sumXY);
219 X(vecZ, sumXY, xZXY);
220 normalize(sumXY, sumXY);
221 normalize(xZXY, xZXY);
222
223 X(sumXY, vecZ, xXYZ);
224 X(xZXY, vecZ, xxZXYZ);
225 X(xZXY, vecX, xxZXYX);
226 X(xZXY, vecY, xxZXYY);
227 normalize(xXYZ, xXYZ);
228 normalize(xxZXYZ, xxZXYZ);
229 normalize(xxZXYX, xxZXYX);
230 normalize(xxZXYY, xxZXYY);
231
232 double cosZXY = dot(vecZ, sumXY);
233 double sinZXY = sqrt(1.0 - cosZXY * cosZXY);
234 double cosXZXY = dot(vecX, xZXY);
235 double sinXZXY = sqrt(1.0 - cosXZXY * cosXZXY);
236 double cosYZXY = dot(vecY, xZXY);
237 double sinYZXY = sqrt(1.0 - cosYZXY * cosYZXY);
238
239 scale(xZXY, -cosXZXY, t1);
240 scale(xZXY, -cosYZXY, t2);
241 add(vecX, t1, t1);
242 add(vecY, t2, t2);
243 normalize(t1, t1);
244 normalize(t2, t2);
245 double ut1cos = dot(vecZ, t1);
246 double ut1sin = sqrt(1.0 - ut1cos * ut1cos);
247 double ut2cos = dot(vecZ, t2);
248 double ut2sin = sqrt(1.0 - ut2cos * ut2cos);
249 double dPhidR = -(trq[0] * sumXY[0] + trq[1] * sumXY[1] + trq[2] * sumXY[2]);
250 double dPhidZXY = -(trq[0] * xZXY[0] + trq[1] * xZXY[1] + trq[2] * xZXY[2]);
251 for (int j = 0; j < 3; j++) {
252 double dZ = xXYZ[j] * dPhidR / (rZ * sinZXY) + xxZXYZ[j] * dPhidZXY / rZ;
253 double dX = (sinXZXY * xZXY[j] - cosXZXY * t1[j]) * dPhidZ / (rX * (ut1sin + ut2sin));
254 double dY = (sinYZXY * xZXY[j] - cosYZXY * t2[j]) * dPhidZ / (rY * (ut1sin + ut2sin));
255 g[0][j] = dZ;
256 g[1][j] = dX;
257 g[2][j] = dY;
258 g[3][j] = -(dZ + dX + dY);
259 }
260 break;
261 case THREEFOLD:
262 add(vecZ, vecX, sumZXY);
263 add(vecY, sumZXY, sumZXY);
264 double rZXY = length(sumZXY);
265 scale(sumZXY, 1.0 / rZXY, sumZXY);
266 double cosYSum = dot(vecY, sumZXY);
267 double cosZSum = dot(vecZ, sumZXY);
268 double cosXSum = dot(vecX, sumZXY);
269 add(vecX, vecY, r);
270 normalize(r, r);
271 double cosRZ = dot(r, vecZ);
272 double sinRZ = sqrt(1.0 - cosRZ * cosRZ);
273 dPhidR = -trq[0] * r[0] - trq[1] * r[1] - trq[2] * r[2];
274 X(r, vecZ, del);
275 normalize(del, del);
276 double dPhidDel = -trq[0] * del[0] - trq[1] * del[1] - trq[2] * del[2];
277 X(del, vecZ, eps);
278 for (int j = 0; j < 3; j++) {
279 double dZ = del[j] * dPhidR / (rZ * sinRZ) + eps[j] * dPhidDel * cosZSum / (rZ * rZXY);
280 g[0][j] = dZ;
281 g[3][j] = -dZ;
282 }
283 add(vecZ, vecY, r);
284 normalize(r, r);
285 double cosRX = dot(r, vecX);
286 double sinRX = sqrt(1.0 - cosRX * cosRX);
287 dPhidR = -trq[0] * r[0] - trq[1] * r[1] - trq[2] * r[2];
288 X(r, vecX, del);
289 normalize(del, del);
290 dPhidDel = -trq[0] * del[0] - trq[1] * del[1] - trq[2] * del[2];
291 X(del, vecX, eps);
292 for (int j = 0; j < 3; j++) {
293 double dX = del[j] * dPhidR / (rX * sinRX) + eps[j] * dPhidDel * cosXSum / (rX * rZXY);
294 g[1][j] = dX;
295 g[3][j] -= dX;
296 }
297 add(vecZ, vecX, r);
298 normalize(r, r);
299 double cosRY = dot(r, vecY);
300 double sinRY = sqrt(1.0 - cosRY * cosRY);
301 dPhidR = -trq[0] * r[0] - trq[1] * r[1] - trq[2] * r[2];
302 X(r, vecY, del);
303 normalize(del, del);
304 dPhidDel = -trq[0] * del[0] - trq[1] * del[1] - trq[2] * del[2];
305 X(del, vecY, eps);
306 for (int j = 0; j < 3; j++) {
307 double dY = del[j] * dPhidR / (rY * sinRY) + eps[j] * dPhidDel * cosYSum / (rY * rZXY);
308 g[2][j] = dY;
309 g[3][j] -= dY;
310 }
311 break;
312 default:
313 String message = "Fatal exception: Unsupported frame definition: " + frame[i] + "\n";
314 logger.log(Level.SEVERE, message);
315 }
316 }
317 }