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.bonded;
39
40 import static ffx.numerics.math.DoubleMath.dihedralAngle;
41 import static java.lang.String.format;
42 import static java.lang.System.arraycopy;
43 import static java.util.Arrays.copyOf;
44 import static org.apache.commons.math3.util.FastMath.acos;
45 import static org.apache.commons.math3.util.FastMath.max;
46 import static org.apache.commons.math3.util.FastMath.min;
47 import static org.apache.commons.math3.util.FastMath.sqrt;
48 import static org.apache.commons.math3.util.FastMath.toDegrees;
49
50 import ffx.numerics.atomic.AtomicDoubleArray3D;
51 import ffx.potential.parameters.AngleTorsionType;
52 import ffx.potential.parameters.AngleType;
53 import ffx.potential.parameters.ForceField;
54 import ffx.potential.parameters.TorsionType;
55
56 import java.io.Serial;
57 import java.util.logging.Logger;
58
59
60
61
62
63
64
65 public class AngleTorsion extends BondedTerm implements LambdaInterface {
66
67 @Serial
68 private static final long serialVersionUID = 1L;
69
70 private static final Logger logger = Logger.getLogger(AngleTorsion.class.getName());
71
72 private static final String mathForm;
73
74 static {
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91 StringBuilder mathFormBuilder = new StringBuilder();
92 for (int m = 1; m < 3; m++) {
93 for (int n = 1; n < 4; n++) {
94
95 mathFormBuilder.append(
96 format("k%d%d*(aVal%d-a%d)*(1+cos(%d*tVal+phi%d))+", m, n, m, m, n, n));
97 }
98 }
99 int lenStr = mathFormBuilder.length();
100 mathFormBuilder.replace(lenStr - 1, lenStr, ";");
101 for (int m = 1; m < 3; m++) {
102 mathFormBuilder.append(format("aVal%d=angle(p%d,p%d,p%d);", m, m, (m + 1), (m + 2)));
103 }
104 mathFormBuilder.append("tVal=dihedral(p1,p2,p3,p4)");
105 mathForm = mathFormBuilder.toString();
106 }
107
108
109
110
111
112 private final double[] constants = new double[6];
113
114 private final double[] tsin = new double[] {0.0, 0.0, 0.0};
115 private final double[] tcos = new double[] {1.0, 1.0, 1.0};
116
117 public AngleType angleType1 = null;
118
119 public AngleType angleType2 = null;
120
121 private AngleTorsionType angleTorsionType = null;
122
123 private TorsionType torsionType = null;
124
125 private double lambda = 1.0;
126
127 private double dEdL = 0.0;
128
129 private boolean lambdaTerm = false;
130
131
132
133
134
135
136
137 public AngleTorsion(Angle an1, Angle an2) {
138 super();
139 bonds = new Bond[3];
140 bonds[1] = an1.getCommonBond(an2);
141 bonds[0] = an1.getOtherBond(bonds[1]);
142 bonds[2] = an2.getOtherBond(bonds[1]);
143 initialize();
144 }
145
146
147
148
149
150
151
152
153 public AngleTorsion(Bond b1, Bond b2, Bond b3) {
154 super();
155 bonds = new Bond[3];
156 bonds[0] = b1;
157 bonds[1] = b2;
158 bonds[2] = b3;
159 initialize();
160 }
161
162
163
164
165
166
167 public AngleTorsion(String n) {
168 super(n);
169 }
170
171
172
173
174
175
176 public static String angleTorsionForm() {
177 return mathForm;
178 }
179
180
181
182
183
184
185
186
187 static AngleTorsion angleTorsionFactory(Torsion torsion, ForceField forceField) {
188 TorsionType torsionType = torsion.torsionType;
189 String key = torsionType.getKey();
190 AngleTorsionType angleTorsionType = forceField.getAngleTorsionType(key);
191
192 if (angleTorsionType != null) {
193 Bond bond1 = torsion.bonds[0];
194 Bond middleBond = torsion.bonds[1];
195 Bond bond3 = torsion.bonds[2];
196
197 AngleTorsion angleTorsion = new AngleTorsion(bond1, middleBond, bond3);
198 angleTorsion.angleTorsionType = angleTorsionType;
199 angleTorsion.torsionType = torsion.torsionType;
200 Atom atom1 = torsion.atoms[0];
201 Atom atom2 = torsion.atoms[1];
202 Atom atom3 = torsion.atoms[2];
203 Atom atom4 = torsion.atoms[3];
204
205 Angle angle1 = atom1.getAngle(atom2, atom3);
206 Angle angle2 = atom2.getAngle(atom3, atom4);
207 angleTorsion.angleType1 = angle1.angleType;
208 angleTorsion.angleType2 = angle2.angleType;
209
210 angleTorsion.setFlipped(atom1.getAtomType().atomClass != angleTorsionType.atomClasses[0]
211 || atom2.getAtomType().atomClass != angleTorsionType.atomClasses[1]
212 || atom3.getAtomType().atomClass != angleTorsionType.atomClasses[2]
213 || atom4.getAtomType().atomClass != angleTorsionType.atomClasses[3]);
214
215 return angleTorsion;
216 }
217 return null;
218 }
219
220
221
222
223
224
225
226
227
228
229 public boolean compare(Atom a0, Atom a1, Atom a2, Atom a3) {
230 if (a0 == atoms[0] && a1 == atoms[1] && a2 == atoms[2] && a3 == atoms[3]) {
231 return true;
232 }
233
234 return (a0 == atoms[3] && a1 == atoms[2] && a2 == atoms[1] && a3 == atoms[0]);
235 }
236
237
238
239
240
241
242 @Override
243 public double energy(boolean gradient, int threadID, AtomicDoubleArray3D grad, AtomicDoubleArray3D lambdaGrad) {
244 energy = 0.0;
245 value = 0.0;
246 dEdL = 0.0;
247
248 if (!getUse()) {
249 return energy;
250 }
251 var atomA = atoms[0];
252 var atomB = atoms[1];
253 var atomC = atoms[2];
254 var atomD = atoms[3];
255 var va = atomA.getXYZ();
256 var vb = atomB.getXYZ();
257 var vc = atomC.getXYZ();
258 var vd = atomD.getXYZ();
259 var vba = vb.sub(va);
260 var vcb = vc.sub(vb);
261 var vdc = vd.sub(vc);
262 var rba2 = vba.length2();
263 var rcb2 = vcb.length2();
264 var rdc2 = vdc.length2();
265 if (min(min(rba2, rcb2), rdc2) == 0.0) {
266 return 0.0;
267 }
268 var rcb = sqrt(rcb2);
269 var vt = vba.X(vcb);
270 var vu = vcb.X(vdc);
271 var rt2 = max(vt.length2(), 0.000001);
272 var ru2 = max(vu.length2(), 0.000001);
273 var rtru = sqrt(rt2 * ru2);
274 var vca = vc.sub(va);
275 var vdb = vd.sub(vb);
276 var cosine = vt.dot(vu) / rtru;
277 var sine = vcb.dot(vt.X(vu)) / (rcb * rtru);
278 value = toDegrees(acos(cosine));
279 if (sine < 0.0) {
280 value = -value;
281 }
282
283
284 var cosine2 = cosine * cosine - sine * sine;
285 var sine2 = 2.0 * cosine * sine;
286 var cosine3 = cosine * cosine2 - sine * sine2;
287 var sine3 = cosine * sine2 + sine * cosine2;
288 var phi1 = 1.0 + (cosine * tcos[0] + sine * tsin[0]);
289 var dphi1 = (cosine * tsin[0] - sine * tcos[0]);
290 var phi2 = 1.0 + (cosine2 * tcos[1] + sine2 * tsin[1]);
291 var dphi2 = 2.0 * (cosine2 * tsin[1] - sine2 * tcos[1]);
292 var phi3 = 1.0 + (cosine3 * tcos[2] + sine3 * tsin[2]);
293 var dphi3 = 3.0 * (cosine3 * tsin[2] - sine3 * tcos[2]);
294
295
296 var c1 = constants[0];
297 var c2 = constants[1];
298 var c3 = constants[2];
299 var angle1 = toDegrees(acos(-vba.dot(vcb) / sqrt(rba2 * rcb2)));
300 var dt1 = angle1 - angleType1.angle[0];
301 var s1 = c1 * phi1 + c2 * phi2 + c3 * phi3;
302 var e1 = angleTorsionType.angtorunit * dt1 * s1;
303
304
305 var c4 = constants[3];
306 var c5 = constants[4];
307 var c6 = constants[5];
308 var angle2 = toDegrees(acos(-vcb.dot(vdc) / sqrt(rcb2 * rdc2)));
309 var dt2 = angle2 - angleType2.angle[0];
310 var s2 = c4 * phi1 + c5 * phi2 + c6 * phi3;
311 var e2 = angleTorsionType.angtorunit * dt2 * s2;
312 energy = e1 + e2;
313 energy = energy * lambda;
314 dEdL = energy;
315 if (gradient || lambdaTerm) {
316
317 var dedphi = angleTorsionType.angtorunit * dt1 * (c1 * dphi1 + c2 * dphi2 + c3 * dphi3);
318 var ddt = angleTorsionType.angtorunit * toDegrees(s1);
319 var vdt = vt.X(vcb).scaleI(dedphi / (rt2 * rcb));
320 var vdu = vcb.X(vu).scaleI(dedphi / (ru2 * rcb));
321
322
323 var rt = sqrt(rt2);
324 var sa = -ddt / (rba2 * rt);
325 var sc = ddt / (rcb2 * rt);
326 var ga = vt.X(vba).scaleI(sa).addI(vdt.X(vcb));
327 var gb = vba.X(vt).scaleI(sa).addI(vt.X(vcb).scaleI(sc)).addI(vca.X(vdt)).addI(vdu.X(vdc));
328 var gc = vcb.X(vt).scaleI(sc).addI(vdt.X(vba)).addI(vdb.X(vdu));
329 var gd = vdu.X(vcb);
330
331
332 dedphi = angleTorsionType.angtorunit * dt2 * (c4 * dphi1 + c5 * dphi2 + c6 * dphi3);
333 ddt = angleTorsionType.angtorunit * toDegrees(s2);
334 vdt = vt.X(vcb).scaleI(dedphi / (rt2 * rcb));
335 vdu = vcb.X(vu).scaleI(dedphi / (ru2 * rcb));
336
337
338 var ur = sqrt(ru2);
339 var sb = -ddt / (rcb2 * ur);
340 var sd = ddt / (rdc2 * ur);
341 ga.addI(vdt.X(vcb));
342 gb.addI(vu.X(vcb).scaleI(sb)).addI(vca.X(vdt)).addI(vdu.X(vdc));
343 gc.addI(vcb.X(vu).scaleI(sb)).addI(vu.X(vdc).scaleI(sd)).addI(vdt.X(vba)).addI(vdb.X(vdu));
344 gd.addI(vdc.X(vu).scaleI(sd)).addI(vdu.X(vcb));
345 var ia = atomA.getIndex() - 1;
346 var ib = atomB.getIndex() - 1;
347 var ic = atomC.getIndex() - 1;
348 var id = atomD.getIndex() - 1;
349 if (lambdaTerm) {
350 lambdaGrad.add(threadID, ia, ga);
351 lambdaGrad.add(threadID, ib, gb);
352 lambdaGrad.add(threadID, ic, gc);
353 lambdaGrad.add(threadID, id, gd);
354 }
355 if (gradient) {
356 grad.add(threadID, ia, ga.scaleI(lambda));
357 grad.add(threadID, ib, gb.scaleI(lambda));
358 grad.add(threadID, ic, gc.scaleI(lambda));
359 grad.add(threadID, id, gd.scaleI(lambda));
360 }
361 }
362
363 return energy;
364 }
365
366
367
368
369
370
371
372
373 public Atom get1_4(ffx.potential.bonded.Atom a) {
374 if (a == atoms[0]) {
375 return atoms[3];
376 }
377 if (a == atoms[3]) {
378 return atoms[0];
379 }
380 return null;
381 }
382
383
384
385
386
387
388 public double[] getConstants() {
389 return copyOf(constants, constants.length);
390 }
391
392
393 @Override
394 public double getLambda() {
395 return lambda;
396 }
397
398
399 @Override
400 public void setLambda(double lambda) {
401 if (applyAllLambda()) {
402 this.lambda = lambda;
403 lambdaTerm = true;
404 } else {
405 this.lambda = 1.0;
406 }
407 }
408
409
410 @Override
411 public double getd2EdL2() {
412 return 0.0;
413 }
414
415
416 @Override
417 public double getdEdL() {
418 if (lambdaTerm) {
419 return dEdL;
420 } else {
421 return 0.0;
422 }
423 }
424
425
426 @Override
427 public void getdEdXdL(double[] gradient) {
428
429 }
430
431
432 public void log() {
433 logger.info(format(" %-8s %6d-%s %6d-%s %6d-%s %6d-%s %10.4f %10.4f", "Angle-Torsion",
434 atoms[0].getIndex(), atoms[0].getAtomType().name,
435 atoms[1].getIndex(), atoms[1].getAtomType().name,
436 atoms[2].getIndex(), atoms[2].getAtomType().name,
437 atoms[3].getIndex(), atoms[3].getAtomType().name,
438 value, energy));
439 }
440
441
442
443
444
445
446 @Override
447 public String toString() {
448 return format("%s (%7.1f,%7.2f)", id, value, energy);
449 }
450
451
452 private void initialize() {
453 atoms = new Atom[4];
454 atoms[1] = bonds[0].getCommonAtom(bonds[1]);
455 atoms[0] = bonds[0].get1_2(atoms[1]);
456 atoms[2] = bonds[1].get1_2(atoms[1]);
457 atoms[3] = bonds[2].get1_2(atoms[2]);
458 setID_Key(false);
459 value = dihedralAngle(atoms[0].getXYZ(null), atoms[1].getXYZ(null),
460 atoms[2].getXYZ(null), atoms[3].getXYZ(null));
461 }
462
463
464
465
466
467
468 private void setFlipped(boolean flipped) {
469 if (flipped) {
470 constants[0] = angleTorsionType.forceConstants[3];
471 constants[1] = angleTorsionType.forceConstants[4];
472 constants[2] = angleTorsionType.forceConstants[5];
473 constants[3] = angleTorsionType.forceConstants[0];
474 constants[4] = angleTorsionType.forceConstants[1];
475 constants[5] = angleTorsionType.forceConstants[2];
476 } else {
477 arraycopy(angleTorsionType.forceConstants, 0, constants, 0, 6);
478 }
479
480 arraycopy(torsionType.sine, 0, tsin, 0, min(torsionType.sine.length, 3));
481 arraycopy(torsionType.cosine, 0, tcos, 0, min(torsionType.cosine.length, 3));
482 }
483 }