View Javadoc
1   // ******************************************************************************
2   //
3   // Title:       Force Field X.
4   // Description: Force Field X - Software for Molecular Biophysics.
5   // Copyright:   Copyright (c) Michael J. Schnieders 2001-2024.
6   //
7   // This file is part of Force Field X.
8   //
9   // Force Field X is free software; you can redistribute it and/or modify it
10  // under the terms of the GNU General Public License version 3 as published by
11  // the Free Software Foundation.
12  //
13  // Force Field X is distributed in the hope that it will be useful, but WITHOUT
14  // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
15  // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
16  // details.
17  //
18  // You should have received a copy of the GNU General Public License along with
19  // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
20  // Place, Suite 330, Boston, MA 02111-1307 USA
21  //
22  // Linking this library statically or dynamically with other modules is making a
23  // combined work based on this library. Thus, the terms and conditions of the
24  // GNU General Public License cover the whole combination.
25  //
26  // As a special exception, the copyright holders of this library give you
27  // permission to link this library with independent modules to produce an
28  // executable, regardless of the license terms of these independent modules, and
29  // to copy and distribute the resulting executable under terms of your choice,
30  // provided that you also meet, for each linked independent module, the terms
31  // and conditions of the license of that module. An independent module is a
32  // module which is not derived from or based on this library. If you modify this
33  // library, you may extend this exception to your version of the library, but
34  // you are not obligated to do so. If you do not wish to do so, delete this
35  // exception statement from your version.
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   * The AngleTorsion class represents an angle torsion coupling between four bonded atoms.
61   *
62   * @author Michael J. Schnieders
63   * @since 1.0
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    /** Functional form for OpenMM. */
72    private static final String mathForm;
73  
74    static {
75      /*
76       Defined constants:
77       p1-p4 are particles 1-4.
78       m is an angle number, from 1-2, representing angles p1-p2-p3, p2-p3-p4.
79       n is a periodicity, from 1-3.
80  
81       k[m][n] is a set of 6 energy constants defined in the parameter file for this angle-torsion.
82  
83       aVal[m] is the current value for angle m.
84       a[m] is the equilibrium value for angle m.
85  
86       tVal is the current value of the 1-2-3-4 dihedral angle.
87  
88       phi[m] is a phase offset constant; phi1 = phi3 = 0, phi2 = pi.
89      */
90  
91      StringBuilder mathFormBuilder = new StringBuilder();
92      for (int m = 1; m < 3; m++) {
93        for (int n = 1; n < 4; n++) {
94          // kmn * (am - am0) * (1 + cos(n*tors + phi(n)))
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    * Angle Torsion force constants (may be reversed compared to storage in the AngleTorsionType
110    * instance).
111    */
112   private final double[] constants = new double[6];
113   /** The AngleTorsion may use more sine and cosine terms than are defined in the TorsionType. */
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   /** First angle force field type. */
117   public AngleType angleType1 = null;
118   /** Second angle force field type. */
119   public AngleType angleType2 = null;
120   /** Angle Torsion force field type. */
121   private AngleTorsionType angleTorsionType = null;
122   /** Torsion force field type. */
123   private TorsionType torsionType = null;
124   /** Value of lambda. */
125   private double lambda = 1.0;
126   /** Value of dE/dL. */
127   private double dEdL = 0.0;
128   /** Flag to indicate lambda dependence. */
129   private boolean lambdaTerm = false;
130 
131   /**
132    * AngleTorsion constructor.
133    *
134    * @param an1 Angle that combines to form the Torsional Angle
135    * @param an2 Angle that combines to form the Torsional Angle
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    * Create a AngleTorsion from 3 connected bonds (no error checking)
148    *
149    * @param b1 Bond
150    * @param b2 Bond
151    * @param b3 Bond
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    * AngleTorsion Constructor.
164    *
165    * @param n Torsion id
166    */
167   public AngleTorsion(String n) {
168     super(n);
169   }
170 
171   /**
172    * Returns the mathematical form of an angle-torsion as an OpenMM-parsable String.
173    *
174    * @return Mathematical form of the angle-torsion coupling.
175    */
176   public static String angleTorsionForm() {
177     return mathForm;
178   }
179 
180   /**
181    * Attempt to create a new AngleTorsion based on the supplied torsion.
182    *
183    * @param torsion the Torsion.
184    * @param forceField the ForceField parameters to apply.
185    * @return a new Torsion, or null.
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    * compare
222    *
223    * @param a0 a {@link ffx.potential.bonded.Atom} object.
224    * @param a1 a {@link ffx.potential.bonded.Atom} object.
225    * @param a2 a {@link ffx.potential.bonded.Atom} object.
226    * @param a3 a {@link ffx.potential.bonded.Atom} object.
227    * @return a boolean.
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    * {@inheritDoc}
239    *
240    * <p>Evaluate the Angle-Torsion energy.
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     // Only compute this term if at least one atom is being used.
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     // Compute multiple angle trigonometry and phase terms
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     // Set the angle-torsion parameters for the first angle
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     // Set the angle-torsion values for the second angle
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       // Compute derivative components for this interaction.
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       // Determine chain rule components for the first angle.
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       // Compute derivative components for the 2nd angle.
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       // Increment chain rule components for the 2nd angle.
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    * If the specified atom is not a central atom of <b>this</b> torsion, the atom at the opposite end
368    * is returned. These atoms are said to be 1-4 to each other.
369    *
370    * @param a Atom
371    * @return Atom
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    * Returns the array of stretch-torsion constants, in units of kcal/mol/degree.
385    *
386    * @return Stretch-torsion constants.
387    */
388   public double[] getConstants() {
389     return copyOf(constants, constants.length);
390   }
391 
392   /** {@inheritDoc} */
393   @Override
394   public double getLambda() {
395     return lambda;
396   }
397 
398   /** {@inheritDoc} */
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   /** {@inheritDoc} */
410   @Override
411   public double getd2EdL2() {
412     return 0.0;
413   }
414 
415   /** {@inheritDoc} */
416   @Override
417   public double getdEdL() {
418     if (lambdaTerm) {
419       return dEdL;
420     } else {
421       return 0.0;
422     }
423   }
424 
425   /** {@inheritDoc} */
426   @Override
427   public void getdEdXdL(double[] gradient) {
428     // The dEdXdL terms are zero.
429   }
430 
431   /** Log details for this Angle-Torsion energy term. */
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    * {@inheritDoc}
443    *
444    * <p>Overridden toString Method returns the Term's id.
445    */
446   @Override
447   public String toString() {
448     return format("%s  (%7.1f,%7.2f)", id, value, energy);
449   }
450 
451   /** Initialization */
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    * setFlipped.
465    *
466    * @param flipped a boolean.
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 }