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-2025.
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.openmm;
39  
40  import ffx.openmm.Force;
41  import ffx.openmm.PeriodicTorsionForce;
42  import ffx.potential.ForceFieldEnergy;
43  import ffx.potential.bonded.Torsion;
44  import ffx.potential.parameters.ForceField;
45  import ffx.potential.parameters.TorsionType;
46  import ffx.potential.terms.TorsionPotentialEnergy;
47  
48  import java.util.logging.Logger;
49  
50  import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_KJPerKcal;
51  import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_RadiansPerDegree;
52  import static java.lang.String.format;
53  
54  /**
55   * Torsion Force.
56   *
57   * @author Michael J. Schnieders
58   * @since 1.0
59   */
60  public class TorsionForce extends PeriodicTorsionForce {
61  
62    private static final Logger logger = Logger.getLogger(TorsionForce.class.getName());
63  
64    private final boolean manyBodyTitration;
65  
66    /**
67     * Torsion Force constructor.
68     *
69     * @param torsionPotentialEnergy The TorsionPotentialEnergy instance that contains the torsions.
70     * @param openMMEnergy           The OpenMM Energy instance that contains the torsions.
71     */
72    public TorsionForce(TorsionPotentialEnergy torsionPotentialEnergy, OpenMMEnergy openMMEnergy) {
73      ForceField forceField = openMMEnergy.getMolecularAssembly().getForceField();
74      manyBodyTitration = forceField.getBoolean("MANYBODY_TITRATION", false);
75      Torsion[] torsions = torsionPotentialEnergy.getTorsionArray();
76      for (Torsion torsion : torsions) {
77        int a1 = torsion.getAtom(0).getArrayIndex();
78        int a2 = torsion.getAtom(1).getArrayIndex();
79        int a3 = torsion.getAtom(2).getArrayIndex();
80        int a4 = torsion.getAtom(3).getArrayIndex();
81        TorsionType torsionType = torsion.torsionType;
82        int nTerms = torsionType.phase.length;
83        for (int j = 0; j < nTerms; j++) {
84          double k = torsion.getTorsionScale() * torsionType.torsionUnit * torsionType.amplitude[j];
85          addTorsion(a1, a2, a3, a4, j + 1,
86              torsionType.phase[j] * OpenMM_RadiansPerDegree, OpenMM_KJPerKcal * k);
87        }
88        // Enforce 6-fold torsions since TorsionType instances can have different lengths
89        // when side-chain protonation changes.
90        if (manyBodyTitration) {
91          for (int j = nTerms; j < 6; j++) {
92            addTorsion(a1, a2, a3, a4, j + 1, 0.0, 0.0);
93          }
94        }
95      }
96  
97      int forceGroup = torsionPotentialEnergy.getForceGroup();
98      setForceGroup(forceGroup);
99      logger.info(format("  Torsions:                          %10d", torsions.length));
100     logger.fine(format("   Force Group:                      %10d", forceGroup));
101   }
102 
103   /**
104    * Torsion Force constructor for Dual Topology.
105    *
106    * @param torsionPotentialEnergy   The TorsionPotentialEnergy instance that contains the torsions.
107    * @param topology                 The topology index for the OpenMM System.
108    * @param openMMDualTopologyEnergy The OpenMMDualTopologyEnergy instance.
109    */
110   public TorsionForce(TorsionPotentialEnergy torsionPotentialEnergy,
111                       int topology, OpenMMDualTopologyEnergy openMMDualTopologyEnergy) {
112     ForceFieldEnergy forceFieldEnergy = openMMDualTopologyEnergy.getForceFieldEnergy(topology);
113     ForceField forceField = forceFieldEnergy.getMolecularAssembly().getForceField();
114     manyBodyTitration = forceField.getBoolean("MANYBODY_TITRATION", false);
115     if (manyBodyTitration) {
116       logger.severe("OpenMM Dual Topology does not suppport many body titration.");
117     }
118     Torsion[] torsions = torsionPotentialEnergy.getTorsionArray();
119     double scale = openMMDualTopologyEnergy.getTopologyScale(topology);
120 
121     for (Torsion torsion : torsions) {
122       int a1 = torsion.getAtom(0).getArrayIndex();
123       int a2 = torsion.getAtom(1).getArrayIndex();
124       int a3 = torsion.getAtom(2).getArrayIndex();
125       int a4 = torsion.getAtom(3).getArrayIndex();
126       a1 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, a1);
127       a2 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, a2);
128       a3 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, a3);
129       a4 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, a4);
130       TorsionType torsionType = torsion.torsionType;
131       int nTerms = torsionType.phase.length;
132       for (int j = 0; j < nTerms; j++) {
133         double k = torsion.getTorsionScale() * torsionType.torsionUnit * torsionType.amplitude[j];
134         // Don't apply lambda scale to alchemical torsion
135         if (!torsion.applyLambda()) {
136           k = k * scale;
137         }
138         addTorsion(a1, a2, a3, a4, j + 1,
139             torsionType.phase[j] * OpenMM_RadiansPerDegree, OpenMM_KJPerKcal * k);
140       }
141     }
142 
143     int forceGroup = torsionPotentialEnergy.getForceGroup();
144     setForceGroup(forceGroup);
145     logger.info(format("  Torsions:                          %10d", torsions.length));
146     logger.fine(format("   Force Group:                      %10d", forceGroup));
147   }
148 
149   /**
150    * Convenience method to construct an OpenMM Torsion Force.
151    *
152    * @param openMMEnergy The OpenMM Energy instance that contains the torsions.
153    * @return A Torsion Force, or null if there are no torsions.
154    */
155   public static Force constructForce(OpenMMEnergy openMMEnergy) {
156     TorsionPotentialEnergy torsionPotentialEnergy = openMMEnergy.getTorsionPotentialEnergy();
157     if (torsionPotentialEnergy == null) {
158       return null;
159     }
160     return new TorsionForce(torsionPotentialEnergy, openMMEnergy);
161   }
162 
163   /**
164    * Convenience method to construct a Dual-Topology OpenMM Torsion Force.
165    *
166    * @param topology                 The topology index for the OpenMM System.
167    * @param openMMDualTopologyEnergy The OpenMMDualTopologyEnergy instance.
168    * @return A Torsion Force, or null if there are no torsions.
169    */
170   public static Force constructForce(int topology, OpenMMDualTopologyEnergy openMMDualTopologyEnergy) {
171     ForceFieldEnergy forceFieldEnergy = openMMDualTopologyEnergy.getForceFieldEnergy(topology);
172     TorsionPotentialEnergy torsionPotentialEnergy = forceFieldEnergy.getTorsionPotentialEnergy();
173     if (torsionPotentialEnergy == null) {
174       return null;
175     }
176     return new TorsionForce(torsionPotentialEnergy, topology, openMMDualTopologyEnergy);
177   }
178 
179   /**
180    * Update the Torsion force.
181    *
182    * @param openMMEnergy The OpenMM Energy that contains the torsions.
183    */
184   public void updateForce(OpenMMEnergy openMMEnergy) {
185     TorsionPotentialEnergy torsionPotentialEnergy = openMMEnergy.getTorsionPotentialEnergy();
186     if (torsionPotentialEnergy == null) {
187       return;
188     }
189     Torsion[] torsions = torsionPotentialEnergy.getTorsionArray();
190 
191     int index = 0;
192     for (Torsion torsion : torsions) {
193       TorsionType torsionType = torsion.torsionType;
194       int nTerms = torsionType.phase.length;
195       int a1 = torsion.getAtom(0).getArrayIndex();
196       int a2 = torsion.getAtom(1).getArrayIndex();
197       int a3 = torsion.getAtom(2).getArrayIndex();
198       int a4 = torsion.getAtom(3).getArrayIndex();
199       for (int j = 0; j < nTerms; j++) {
200         double k = torsion.getTorsionScale() * torsionType.torsionUnit * torsionType.amplitude[j];
201         setTorsionParameters(index++, a1, a2, a3, a4, j + 1,
202             torsionType.phase[j] * OpenMM_RadiansPerDegree, OpenMM_KJPerKcal * k);
203       }
204       // Enforce 6-fold torsions since TorsionType instances can have different lengths
205       // when side-chain protonation changes.
206       if (manyBodyTitration) {
207         for (int j = nTerms; j < 6; j++) {
208           setTorsionParameters(index++, a1, a2, a3, a4, j + 1, 0.0, 0.0);
209         }
210       }
211     }
212     updateParametersInContext(openMMEnergy.getContext());
213   }
214 
215   /**
216    * Update the Torsion force.
217    *
218    * @param topology                 The topology index for the OpenMM System.
219    * @param openMMDualTopologyEnergy The OpenMMDualTopologyEnergy instance.
220    */
221   public void updateForce(int topology, OpenMMDualTopologyEnergy openMMDualTopologyEnergy) {
222     ForceFieldEnergy forceFieldEnergy = openMMDualTopologyEnergy.getForceFieldEnergy(topology);
223     TorsionPotentialEnergy torsionPotentialEnergy = forceFieldEnergy.getTorsionPotentialEnergy();
224     if (torsionPotentialEnergy == null) {
225       return;
226     }
227     // Check if this system has torsions.
228     Torsion[] torsions = torsionPotentialEnergy.getTorsionArray();
229 
230     double scale = openMMDualTopologyEnergy.getTopologyScale(topology);
231 
232     int index = 0;
233     for (Torsion torsion : torsions) {
234       TorsionType torsionType = torsion.torsionType;
235       int nTerms = torsionType.phase.length;
236       int a1 = torsion.getAtom(0).getArrayIndex();
237       int a2 = torsion.getAtom(1).getArrayIndex();
238       int a3 = torsion.getAtom(2).getArrayIndex();
239       int a4 = torsion.getAtom(3).getArrayIndex();
240       a1 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, a1);
241       a2 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, a2);
242       a3 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, a3);
243       a4 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, a4);
244       for (int j = 0; j < nTerms; j++) {
245         double k = torsion.getTorsionScale() * torsionType.torsionUnit * torsionType.amplitude[j];
246         // Don't apply lambda scale to alchemical torsion
247         if (!torsion.applyLambda()) {
248           k = k * scale;
249         }
250         setTorsionParameters(index++, a1, a2, a3, a4, j + 1,
251             torsionType.phase[j] * OpenMM_RadiansPerDegree, OpenMM_KJPerKcal * k);
252       }
253     }
254     updateParametersInContext(openMMDualTopologyEnergy.getContext());
255   }
256 
257 }