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 org.apache.commons.math3.util.FastMath.acos;
43 import static org.apache.commons.math3.util.FastMath.sqrt;
44 import static org.apache.commons.math3.util.FastMath.toDegrees;
45
46 import ffx.numerics.atomic.AtomicDoubleArray3D;
47 import ffx.numerics.math.DoubleMath;
48 import ffx.potential.parameters.AtomType;
49 import ffx.potential.parameters.ForceField;
50 import ffx.potential.parameters.TorsionType;
51 import java.util.ArrayList;
52 import java.util.List;
53 import java.util.function.DoubleUnaryOperator;
54 import java.util.logging.Logger;
55
56
57
58
59
60
61
62 public class Torsion extends BondedTerm implements LambdaInterface {
63
64 private static final Logger logger = Logger.getLogger(Torsion.class.getName());
65
66 public TorsionType torsionType = null;
67
68 public double units = 0.5;
69
70 private double lambda = 1.0;
71
72 private double dEdL = 0.0;
73
74 private boolean lambdaTerm = false;
75
76 private final DoubleUnaryOperator lambdaMapper;
77
78
79
80
81
82
83
84 public Torsion(Angle an1, Angle an2) {
85 super();
86 bonds = new Bond[3];
87 bonds[1] = an1.getCommonBond(an2);
88 bonds[0] = an1.getOtherBond(bonds[1]);
89 bonds[2] = an2.getOtherBond(bonds[1]);
90 lambdaMapper = (double d) -> d;
91 initialize();
92 }
93
94
95
96
97
98
99
100 public Torsion(Angle a, Bond b) {
101 super();
102 bonds = new Bond[3];
103 bonds[0] = b;
104 bonds[1] = a.getBond(0);
105 bonds[2] = a.getBond(1);
106
107 Atom atom = bonds[1].getCommonAtom(b);
108 if (atom == null) {
109 Bond temp = bonds[1];
110 bonds[1] = bonds[2];
111 bonds[2] = temp;
112 }
113 lambdaMapper = (double d) -> d;
114 initialize();
115 }
116
117
118
119
120
121
122
123
124 public Torsion(Bond b1, Bond b2, Bond b3) {
125 super();
126 bonds = new Bond[3];
127 bonds[0] = b1;
128 bonds[1] = b2;
129 bonds[2] = b3;
130 lambdaMapper = (double d) -> d;
131 initialize();
132 }
133
134
135
136
137
138
139 public Torsion(String n) {
140 super(n);
141 lambdaMapper = (double d) -> d;
142 }
143
144
145
146
147
148
149
150
151
152 public static void logNoTorsionType(Atom a0, Atom a1, Atom a2, Atom a3, ForceField forceField) {
153 AtomType atomType0 = a0.getAtomType();
154 AtomType atomType1 = a1.getAtomType();
155 AtomType atomType2 = a2.getAtomType();
156 AtomType atomType3 = a3.getAtomType();
157 int[] c = {atomType0.atomClass, atomType1.atomClass, atomType2.atomClass, atomType3.atomClass};
158 String key = TorsionType.sortKey(c);
159 StringBuilder sb = new StringBuilder(
160 format(
161 "No TorsionType for key: %s\n %s -> %s\n %s -> %s\n %s -> %s\n %s -> %s",
162 key, a0, atomType0, a1, atomType1, a2, atomType2, a3, atomType3));
163 if (matchTorsions(a0, a1, a2, a3, forceField, sb, true) <= 0) {
164 matchTorsions(a0, a1, a2, a3, forceField, sb, false);
165 }
166 logger.severe(sb.toString());
167 }
168
169 private static int matchTorsions(Atom a0, Atom a1, Atom a2, Atom a3,
170 ForceField forceField, StringBuilder sb, boolean strict) {
171 AtomType atomType0 = a0.getAtomType();
172 AtomType atomType1 = a1.getAtomType();
173 AtomType atomType2 = a2.getAtomType();
174 AtomType atomType3 = a3.getAtomType();
175 int c0 = atomType0.atomClass;
176 int c1 = atomType1.atomClass;
177 int c2 = atomType2.atomClass;
178 int c3 = atomType3.atomClass;
179 List<AtomType> types0 = forceField.getSimilarAtomTypes(atomType0);
180 List<AtomType> types1 = forceField.getSimilarAtomTypes(atomType1);
181 List<AtomType> types2 = forceField.getSimilarAtomTypes(atomType2);
182 List<AtomType> types3 = forceField.getSimilarAtomTypes(atomType3);
183 List<TorsionType> torsionTypes = new ArrayList<>();
184 boolean match = false;
185 for (AtomType type1 : types1) {
186 for (AtomType type2 : types2) {
187
188 if ((type1.atomClass == c1 && type2.atomClass == c2) ||
189 (type1.atomClass == c2 && type2.atomClass == c1)) {
190 for (AtomType type0 : types0) {
191 for (AtomType type3 : types3) {
192
193 if (strict) {
194 if ((type0.atomClass != c0) && (type0.atomClass != c3) &&
195 (type3.atomClass != c0) && (type3.atomClass != c3)) {
196 continue;
197 }
198 }
199 TorsionType torsionType = forceField.getTorsionType(type0, type1, type2, type3);
200 if (torsionType != null && !torsionTypes.contains(torsionType)) {
201 if (!match) {
202 match = true;
203 sb.append("\n Similar Angle Types:");
204 }
205 torsionTypes.add(torsionType);
206 sb.append(format("\n %s", torsionType));
207 }
208 }
209 }
210 }
211 }
212 }
213 return torsionTypes.size();
214 }
215
216
217
218
219
220
221
222
223
224
225
226 static Torsion torsionFactory(Bond bond1, Bond middleBond, Bond bond3, ForceField forceField) {
227 Atom a0 = bond1.getOtherAtom(middleBond);
228 Atom a1 = middleBond.getAtom(0);
229 Atom a2 = middleBond.getAtom(1);
230 Atom a3 = bond3.getOtherAtom(middleBond);
231
232 TorsionType torsionType = forceField.getTorsionType(a0.getAtomType(), a1.getAtomType(),
233 a2.getAtomType(), a3.getAtomType());
234
235
236 if (torsionType == null) {
237 logNoTorsionType(a0, a1, a2, a3, forceField);
238 return null;
239 }
240
241 Torsion torsion = new Torsion(bond1, middleBond, bond3);
242 torsion.setTorsionType(torsionType);
243 torsion.units = forceField.getDouble("TORSIONUNIT", 1.0);
244
245 return torsion;
246 }
247
248
249
250
251
252
253 public void setTorsionType(TorsionType torsionType) {
254 this.torsionType = torsionType;
255 }
256
257
258
259
260
261
262
263
264
265
266 public boolean compare(Atom a0, Atom a1, Atom a2, Atom a3) {
267 if (a0 == atoms[0] && a1 == atoms[1] && a2 == atoms[2] && a3 == atoms[3]) {
268 return true;
269 }
270 return (a0 == atoms[3] && a1 == atoms[2] && a2 == atoms[1] && a3 == atoms[0]);
271 }
272
273
274
275
276
277
278 public double measure() {
279 double angle = DoubleMath.dihedralAngle(
280 atoms[0].getXYZ(null),
281 atoms[1].getXYZ(null),
282 atoms[2].getXYZ(null),
283 atoms[3].getXYZ(null));
284 return toDegrees(angle);
285 }
286
287
288
289
290
291
292 @Override
293 public double energy(
294 boolean gradient, int threadID, AtomicDoubleArray3D grad, AtomicDoubleArray3D lambdaGrad) {
295 energy = 0.0;
296 value = 0.0;
297 dEdL = 0.0;
298 var atomA = atoms[0];
299 var atomB = atoms[1];
300 var atomC = atoms[2];
301 var atomD = atoms[3];
302 var va = atomA.getXYZ();
303 var vb = atomB.getXYZ();
304 var vc = atomC.getXYZ();
305 var vd = atomD.getXYZ();
306 var vba = vb.sub(va);
307 var vcb = vc.sub(vb);
308 var vdc = vd.sub(vc);
309 var vt = vba.X(vcb);
310 var vu = vcb.X(vdc);
311 var rt2 = vt.length2();
312 var ru2 = vu.length2();
313 var rtru2 = rt2 * ru2;
314 if (rtru2 != 0.0) {
315 var rr = sqrt(rtru2);
316 var rcb = vcb.length();
317 var cosine = vt.dot(vu) / rr;
318 var sine = vcb.dot(vt.X(vu)) / (rcb * rr);
319 value = toDegrees(acos(cosine));
320 if (sine < 0.0) {
321 value = -value;
322 }
323 var amp = torsionType.amplitude;
324 var tsin = torsionType.sine;
325 var tcos = torsionType.cosine;
326 energy = amp[0] * (1.0 + cosine * tcos[0] + sine * tsin[0]);
327 var dedphi = amp[0] * (cosine * tsin[0] - sine * tcos[0]);
328 var cosprev = cosine;
329 var sinprev = sine;
330 var n = torsionType.terms;
331 for (int i = 1; i < n; i++) {
332 var cosn = cosine * cosprev - sine * sinprev;
333 var sinn = sine * cosprev + cosine * sinprev;
334 var phi = 1.0 + cosn * tcos[i] + sinn * tsin[i];
335 var dphi = (1.0 + i) * (cosn * tsin[i] - sinn * tcos[i]);
336 energy = energy + amp[i] * phi;
337 dedphi = dedphi + amp[i] * dphi;
338 cosprev = cosn;
339 sinprev = sinn;
340 }
341 energy = units * energy * lambda;
342 dEdL = units * energy;
343 if (gradient || lambdaTerm) {
344 dedphi = units * dedphi;
345 var vca = vc.sub(va);
346 var vdb = vd.sub(vb);
347 var dedt = vt.X(vcb).scaleI(dedphi / (rt2 * rcb));
348 var dedu = vu.X(vcb).scaleI(-dedphi / (ru2 * rcb));
349 var ga = dedt.X(vcb);
350 var gb = vca.X(dedt).addI(dedu.X(vdc));
351 var gc = dedt.X(vba).addI(vdb.X(dedu));
352 var gd = dedu.X(vcb);
353 int ia = atomA.getIndex() - 1;
354 int ib = atomB.getIndex() - 1;
355 int ic = atomC.getIndex() - 1;
356 int id = atomD.getIndex() - 1;
357 if (lambdaTerm) {
358 lambdaGrad.add(threadID, ia, ga);
359 lambdaGrad.add(threadID, ib, gb);
360 lambdaGrad.add(threadID, ic, gc);
361 lambdaGrad.add(threadID, id, gd);
362 }
363 if (gradient) {
364 grad.add(threadID, ia, ga.scaleI(lambda));
365 grad.add(threadID, ib, gb.scaleI(lambda));
366 grad.add(threadID, ic, gc.scaleI(lambda));
367 grad.add(threadID, id, gd.scaleI(lambda));
368 }
369 }
370 }
371
372 return energy;
373 }
374
375
376
377
378
379
380
381
382 public Atom get1_4(Atom a) {
383 if (a == atoms[0]) {
384 return atoms[3];
385 }
386 if (a == atoms[3]) {
387 return atoms[0];
388 }
389 return null;
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 lambdaTerm = true;
403 }
404 this.lambda = lambdaTerm ? lambdaMapper.applyAsDouble(lambda) : 1.0;
405 }
406
407
408 @Override
409 public double getd2EdL2() {
410 return 0.0;
411 }
412
413
414 @Override
415 public double getdEdL() {
416 if (lambdaTerm) {
417 return dEdL;
418 } else {
419 return 0.0;
420 }
421 }
422
423
424 @Override
425 public void getdEdXdL(double[] gradient) {
426
427 }
428
429
430 public void log() {
431 logger.info(
432 format(
433 " %-8s %6d-%s %6d-%s %6d-%s %6d-%s %10.4f %10.4f",
434 "Torsional-Angle",
435 atoms[0].getIndex(),
436 atoms[0].getAtomType().name,
437 atoms[1].getIndex(),
438 atoms[1].getAtomType().name,
439 atoms[2].getIndex(),
440 atoms[2].getAtomType().name,
441 atoms[3].getIndex(),
442 atoms[3].getAtomType().name,
443 value,
444 energy));
445 }
446
447
448
449
450
451
452 @Override
453 public String toString() {
454 return format("%s (%7.1f,%7.2f)", id, value, energy);
455 }
456
457
458 private void initialize() {
459 atoms = new Atom[4];
460 atoms[1] = bonds[0].getCommonAtom(bonds[1]);
461 atoms[0] = bonds[0].get1_2(atoms[1]);
462 atoms[2] = bonds[1].get1_2(atoms[1]);
463 atoms[3] = bonds[2].get1_2(atoms[2]);
464 atoms[0].setTorsion(this);
465 atoms[1].setTorsion(this);
466 atoms[2].setTorsion(this);
467 atoms[3].setTorsion(this);
468 setID_Key(false);
469 value =
470 dihedralAngle(
471 atoms[0].getXYZ(null),
472 atoms[1].getXYZ(null),
473 atoms[2].getXYZ(null),
474 atoms[3].getXYZ(null));
475 }
476 }