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-2021.
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 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   * The Torsion class represents a torsional angle formed between four bonded atoms.
58   *
59   * @author Michael J. Schnieders
60   * @since 1.0
61   */
62  public class Torsion extends BondedTerm implements LambdaInterface {
63  
64    private static final Logger logger = Logger.getLogger(Torsion.class.getName());
65    /** The force field Torsion type in use. */
66    public TorsionType torsionType = null;
67    /** Unit conversion. */
68    public double units = 0.5;
69    /** Value of lambda. */
70    private double lambda = 1.0;
71    /** Value of dE/dL. */
72    private double dEdL = 0.0;
73    /** Flag to indicate lambda dependence. */
74    private boolean lambdaTerm = false;
75    /** Maps global lambda to either itself or 1 - global lambda. */
76    private final DoubleUnaryOperator lambdaMapper;
77  
78    /**
79     * Torsion constructor.
80     *
81     * @param an1 Angle that combines to form the Torsional Angle
82     * @param an2 Angle that combines to form the Torsional Angle
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     * Torsion constructor.
96     *
97     * @param a Angle that has one Atom in common with Bond b
98     * @param b Bond that has one Atom in common with Angle A
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     // See if bond 2 or bond 3 is the middle bond
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    * Create a Torsion from 3 connected bonds (no error checking)
119    *
120    * @param b1 Bond
121    * @param b2 Bond
122    * @param b3 Bond
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    * Torsion Constructor.
136    *
137    * @param n Torsion id
138    */
139   public Torsion(String n) {
140     super(n);
141     lambdaMapper = (double d) -> d;
142   }
143 
144   /**
145    * Log that no TorsionType exists.
146    *
147    * @param a0 Atom 0.
148    * @param a1 Atom 1.
149    * @param a2 Atom 2.
150    * @param a3 Atom 3.
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         // Match both inner atom classes.
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               // Match one distal atom class.
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    * Attempt to create a new Torsion based on the supplied bonds. There is no error checking to
218    * enforce that the bonds make up a linear series of 4 bonded atoms.
219    *
220    * @param bond1 the first Bond.
221    * @param middleBond the middle Bond.
222    * @param bond3 the last Bond.
223    * @param forceField the ForceField parameters to apply.
224    * @return a new Torsion, or null.
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     // No torsion type found.
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    * Set the torsion type.
250    *
251    * @param torsionType The TorsionType.
252    */
253   public void setTorsionType(TorsionType torsionType) {
254     this.torsionType = torsionType;
255   }
256 
257   /**
258    * compare
259    *
260    * @param a0 a {@link ffx.potential.bonded.Atom} object.
261    * @param a1 a {@link ffx.potential.bonded.Atom} object.
262    * @param a2 a {@link ffx.potential.bonded.Atom} object.
263    * @param a3 a {@link ffx.potential.bonded.Atom} object.
264    * @return a boolean.
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    * Compute the torsional angle in degrees.
275    *
276    * @return The torsion in degrees.
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    * {@inheritDoc}
289    *
290    * <p>Evaluate the Torsional Angle energy.
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    * If the specified atom is not a central atom of <b>this</b> torsion, the atom at the opposite end
377    * is returned. These atoms are said to be 1-4 to each other.
378    *
379    * @param a Atom
380    * @return Atom
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   /** {@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       lambdaTerm = true;
403     }
404     this.lambda = lambdaTerm ? lambdaMapper.applyAsDouble(lambda) : 1.0;
405   }
406 
407   /** {@inheritDoc} */
408   @Override
409   public double getd2EdL2() {
410     return 0.0;
411   }
412 
413   /** {@inheritDoc} */
414   @Override
415   public double getdEdL() {
416     if (lambdaTerm) {
417       return dEdL;
418     } else {
419       return 0.0;
420     }
421   }
422 
423   /** {@inheritDoc} */
424   @Override
425   public void getdEdXdL(double[] gradient) {
426     // The chain rule term is zero.
427   }
428 
429   /** Log details for this Torsional Angle energy term. */
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    * {@inheritDoc}
449    *
450    * <p>Overidden toString Method returns the Term's id.
451    */
452   @Override
453   public String toString() {
454     return format("%s  (%7.1f,%7.2f)", id, value, energy);
455   }
456 
457   /** Initialization */
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 }