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.potential.parameters.AngleType.AngleMode.IN_PLANE;
41  import static org.apache.commons.math3.util.FastMath.acos;
42  import static org.apache.commons.math3.util.FastMath.max;
43  import static org.apache.commons.math3.util.FastMath.min;
44  import static org.apache.commons.math3.util.FastMath.signum;
45  import static org.apache.commons.math3.util.FastMath.sqrt;
46  import static org.apache.commons.math3.util.FastMath.toDegrees;
47  
48  import ffx.numerics.atomic.AtomicDoubleArray3D;
49  import ffx.potential.parameters.ForceField;
50  import ffx.potential.parameters.OutOfPlaneBendType;
51  
52  import java.io.Serial;
53  import java.util.logging.Logger;
54  
55  /**
56   * The OutOfPlaneBend class represents an Out-Of-Plane Bend.
57   *
58   * @author Michael J. Schnieders
59   * @since 1.0
60   */
61  public class OutOfPlaneBend extends BondedTerm {
62  
63    @Serial
64    private static final long serialVersionUID = 1L;
65  
66    private static final Logger logger = Logger.getLogger(OutOfPlaneBend.class.getName());
67    /** Force field parameters to compute the Out-of-Plane Bend energy. */
68    public OutOfPlaneBendType outOfPlaneBendType = null;
69  
70    /**
71     * OutOfPlaneBend constructor.
72     *
73     * @param angle Angle that contains 3 of 4 OutOfPlaneBend atoms.
74     * @param atom The 4th atom of the trigonal center.
75     */
76    public OutOfPlaneBend(Angle angle, Atom atom) {
77      super();
78      atoms = new Atom[4];
79      atoms[0] = angle.atoms[0];
80      atoms[1] = angle.atoms[1];
81      atoms[2] = angle.atoms[2];
82      atoms[3] = atom;
83      bonds = new Bond[3];
84      bonds[0] = angle.bonds[0];
85      bonds[1] = angle.bonds[1];
86      bonds[2] = atoms[1].getBond(atom);
87      setID_Key(false);
88    }
89  
90    /**
91     * The atom of this out-of-plane bend that was not part of the Angle.
92     *
93     * @return Fourth atom.
94     */
95    public Atom getFourthAtom() {
96      return atoms[3];
97    }
98  
99    /**
100    * Get the trigonal atom of this out-of-plane bend (central atom of the Angle).
101    *
102    * @return Fourth atom.
103    */
104   public Atom getTrigonalAtom() {
105     return atoms[1];
106   }
107 
108   /**
109    * Get the first atom of the Angle.
110    *
111    * @return Fourth atom.
112    */
113   public Atom getFirstAngleAtom() {
114     return atoms[0];
115   }
116 
117   /**
118    * Get the first atom of the Angle.
119    *
120    * @return Fourth atom.
121    */
122   public Atom getLastAngleAtom() {
123     return atoms[2];
124   }
125 
126   /**
127    * Attempt to create a new OutOfPlaneBend instance for a given Angle and Force Field.
128    *
129    * @param angle the Angle to create an OutOfPlaneBend around.
130    * @param forceField the ForceField parameters to use.
131    * @return a new OutOfPlaneBend if the central atom of the angle is trigonal and a force field type
132    *     exists.
133    */
134   public static OutOfPlaneBend outOfPlaneBendFactory(Angle angle, ForceField forceField) {
135     Atom centralAtom = angle.atoms[1];
136     if (centralAtom.isTrigonal()) {
137       Atom fourthAtom = angle.getFourthAtomOfTrigonalCenter();
138       Atom[] atoms = angle.atoms;
139 
140       OutOfPlaneBendType outOfPlaneBendType = forceField.getOutOfPlaneBendType(
141           fourthAtom.getAtomType(), atoms[0].getAtomType(), atoms[1].getAtomType(),
142           atoms[2].getAtomType());
143 
144       if (outOfPlaneBendType != null) {
145         if (angle.getAngleMode() == IN_PLANE) {
146           angle.setInPlaneAtom(fourthAtom);
147         }
148         OutOfPlaneBend outOfPlaneBend = new OutOfPlaneBend(angle, fourthAtom);
149         outOfPlaneBend.setOutOfPlaneBendType(outOfPlaneBendType);
150         return outOfPlaneBend;
151       }
152     }
153     return null;
154   }
155 
156   /** {@inheritDoc} */
157   @Override
158   public int compareTo(BondedTerm o) {
159     if (o == null) {
160       throw new NullPointerException();
161     }
162     if (o == this) {
163       return 0;
164     }
165     if (!o.getClass().isInstance(this)) {
166       return super.compareTo(o);
167     }
168     int this1 = atoms[1].getIndex();
169     int a1 = o.atoms[1].getIndex();
170     if (this1 < a1) {
171       return -1;
172     }
173     if (this1 > a1) {
174       return 1;
175     }
176     int this3 = atoms[3].getIndex();
177     int a3 = o.atoms[3].getIndex();
178     return Integer.compare(this3, a3);
179   }
180 
181   /**
182    * {@inheritDoc}
183    *
184    * <p>Evaluate this Out-of-Plane Bend energy.
185    */
186   @Override
187   public double energy(boolean gradient, int threadID, AtomicDoubleArray3D grad, AtomicDoubleArray3D lambdaGrad) {
188     energy = 0.0;
189     value = 0.0;
190     // Only compute this term if at least one atom is being used.
191     if (!getUse()) {
192       return energy;
193     }
194     var atomA = atoms[0];
195     var atomB = atoms[1];
196     var atomC = atoms[2];
197     var atomD = atoms[3];
198     var va = atomA.getXYZ();
199     var vb = atomB.getXYZ();
200     var vc = atomC.getXYZ();
201     var vd = atomD.getXYZ();
202     var vab = va.sub(vb);
203     var vcb = vc.sub(vb);
204     var vdb = vd.sub(vb);
205     var vad = va.sub(vd);
206     var vcd = vc.sub(vd);
207     var rdb2 = vdb.length2();
208     var rad2 = vad.length2();
209     var rcd2 = vcd.length2();
210     var vp = vcb.X(vdb);
211     var ee = vab.dot(vp);
212     var rac2 = vad.dot(vcd);
213     var cc = rad2 * rcd2 - rac2 * rac2;
214     if (rdb2 != 0.0 && cc != 0.0) {
215       var bkk2 = rdb2 - ee * ee / cc;
216       var cosine = min(1.0, max(-1.0, sqrt(bkk2 / rdb2)));
217       value = toDegrees(acos(cosine));
218       var dv = value;
219       var dv2 = dv * dv;
220       var dv3 = dv2 * dv;
221       var dv4 = dv2 * dv2;
222       energy = outOfPlaneBendType.opBendUnit * outOfPlaneBendType.forceConstant * dv2 * (1.0
223           + outOfPlaneBendType.cubic * dv + outOfPlaneBendType.quartic * dv2
224           + outOfPlaneBendType.pentic * dv3 + outOfPlaneBendType.sextic * dv4);
225       if (gradient) {
226         var deddt = outOfPlaneBendType.opBendUnit * outOfPlaneBendType.forceConstant * dv * toDegrees(
227                 2.0 + 3.0 * outOfPlaneBendType.cubic * dv + 4.0 * outOfPlaneBendType.quartic * dv2
228                     + 5.0 * outOfPlaneBendType.pentic * dv3 + 6.0 * outOfPlaneBendType.sextic * dv4);
229         var dedcos = 0.0;
230         if (ee != 0.0) {
231           dedcos = -deddt * signum(ee) / sqrt(cc * bkk2);
232         } else {
233           dedcos = -deddt / sqrt(cc * bkk2);
234         }
235         var term = ee / cc;
236 
237         // Chain rule terms for first derivative components.
238         var svad = vad.scale(rcd2);
239         var svcd = vcd.scale(rac2);
240         var dcda = svad.sub(svcd).scaleI(term);
241         svad = vad.scale(rac2);
242         svcd = vcd.scale(rad2);
243         var dadc = svcd.sub(svad).scaleI(term);
244         var dcdd = dcda.add(dadc).scaleI(-1.0);
245         var deda = vdb.X(vcb);
246         var dedc = vab.X(vdb);
247         var dedd = vcb.X(vab);
248         dedd.addI(vdb.scaleI(ee / rdb2));
249 
250         // Atomic gradient.
251         var ga = dcda.add(deda).scaleI(dedcos);
252         var gc = dadc.add(dedc).scaleI(dedcos);
253         var gd = dcdd.add(dedd).scaleI(dedcos);
254         var gb = ga.add(gc).addI(gd).scaleI(-1.0);
255         var ia = atomA.getIndex() - 1;
256         var ib = atomB.getIndex() - 1;
257         var ic = atomC.getIndex() - 1;
258         var id = atomD.getIndex() - 1;
259         grad.add(threadID, ia, ga);
260         grad.add(threadID, ib, gb);
261         grad.add(threadID, ic, gc);
262         grad.add(threadID, id, gd);
263       }
264     }
265     return energy;
266   }
267 
268   /** Log details for this Out-of-Plane Bend energy term. */
269   public void log() {
270     logger.info(
271         String.format(" %s %6d-%s %6d-%s %6.4f %10.4f", "Out-of-Plane Bend", atoms[1].getIndex(),
272             atoms[1].getAtomType().name, atoms[3].getIndex(), atoms[3].getAtomType().name, value,
273             energy));
274   }
275 
276   /**
277    * Set a reference to the force field parameters for <b>this</b> Angle.
278    *
279    * @param a a {@link ffx.potential.parameters.OutOfPlaneBendType} object.
280    */
281   public void setOutOfPlaneBendType(OutOfPlaneBendType a) {
282     outOfPlaneBendType = a;
283   }
284 
285   /**
286    * {@inheritDoc}
287    *
288    * <p>Overridden toString Method returns the Term's id.
289    */
290   @Override
291   public String toString() {
292     return String.format("%s  (%7.1f,%7.2f)", id, value, energy);
293   }
294 }