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