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.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
56
57
58
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
68
69
70
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
89
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
105
106
107
108
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
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
151
152
153
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
165
166
167
168
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
181
182
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
205
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
217
218
219
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
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
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 }