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.CustomCompoundBondForce;
41 import ffx.openmm.DoubleArray;
42 import ffx.openmm.Force;
43 import ffx.openmm.IntArray;
44 import ffx.potential.ForceFieldEnergy;
45 import ffx.potential.bonded.AngleTorsion;
46 import ffx.potential.bonded.Atom;
47 import ffx.potential.terms.AnglePotentialEnergy;
48 import ffx.potential.terms.AngleTorsionPotentialEnergy;
49
50 import java.util.logging.Logger;
51
52 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_KJPerKcal;
53 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_RadiansPerDegree;
54 import static java.lang.String.format;
55
56
57
58
59 public class AngleTorsionForce extends CustomCompoundBondForce {
60
61 private static final Logger logger = Logger.getLogger(AngleTorsionForce.class.getName());
62
63
64
65
66
67
68 public AngleTorsionForce(AngleTorsionPotentialEnergy angleTorsionPotentialEnergy) {
69 super(4, AngleTorsion.angleTorsionForm());
70 AngleTorsion[] angleTorsions = angleTorsionPotentialEnergy.getAngleTorsionArray();
71 addGlobalParameter("phi1", 0);
72 addGlobalParameter("phi2", Math.PI);
73 addGlobalParameter("phi3", 0);
74 for (int m = 1; m < 3; m++) {
75 for (int n = 1; n < 4; n++) {
76 addPerBondParameter(format("k%d%d", m, n));
77 }
78 }
79 for (int m = 1; m < 3; m++) {
80 addPerBondParameter(format("a%d", m));
81 }
82 for (AngleTorsion angleTorsion : angleTorsions) {
83 double[] constants = angleTorsion.getConstants();
84 DoubleArray parameters = new DoubleArray(0);
85 for (int m = 0; m < 2; m++) {
86 for (int n = 0; n < 3; n++) {
87 int index = (3 * m) + n;
88 parameters.append(constants[index] * OpenMM_KJPerKcal);
89 }
90 }
91 Atom[] atoms = angleTorsion.getAtomArray(true);
92 parameters.append(angleTorsion.angleType1.angle[0] * OpenMM_RadiansPerDegree);
93 parameters.append(angleTorsion.angleType2.angle[0] * OpenMM_RadiansPerDegree);
94
95 IntArray particles = new IntArray(0);
96 for (int i = 0; i < 4; i++) {
97 particles.append(atoms[i].getArrayIndex());
98 }
99
100 addBond(particles, parameters);
101 parameters.destroy();
102 particles.destroy();
103 }
104
105 int forceGroup = angleTorsionPotentialEnergy.getForceGroup();
106 setForceGroup(forceGroup);
107 logger.info(format(" Angle-Torsions: %10d", angleTorsions.length));
108 logger.fine(format(" Force Group: %10d", forceGroup));
109 }
110
111
112
113
114
115
116
117
118 public AngleTorsionForce(AngleTorsionPotentialEnergy angleTorsionPotentialEnergy,
119 int topology, OpenMMDualTopologyEnergy openMMDualTopologyEnergy) {
120 super(4, AngleTorsion.angleTorsionForm());
121 AngleTorsion[] angleTorsions = angleTorsionPotentialEnergy.getAngleTorsionArray();
122 addGlobalParameter("phi1", 0);
123 addGlobalParameter("phi2", Math.PI);
124 addGlobalParameter("phi3", 0);
125 for (int m = 1; m < 3; m++) {
126 for (int n = 1; n < 4; n++) {
127 addPerBondParameter(format("k%d%d", m, n));
128 }
129 }
130 for (int m = 1; m < 3; m++) {
131 addPerBondParameter(format("a%d", m));
132 }
133
134 double scaleDT = openMMDualTopologyEnergy.getTopologyScale(topology);
135
136 for (AngleTorsion angleTorsion : angleTorsions) {
137 double scale = 1.0;
138
139 if (!angleTorsion.applyLambda()) {
140 scale = scaleDT;
141 }
142 double[] constants = angleTorsion.getConstants();
143 DoubleArray parameters = new DoubleArray(0);
144 for (int m = 0; m < 2; m++) {
145 for (int n = 0; n < 3; n++) {
146 int index = (3 * m) + n;
147 parameters.append(constants[index] * OpenMM_KJPerKcal * scale);
148 }
149 }
150 Atom[] atoms = angleTorsion.getAtomArray(true);
151 parameters.append(angleTorsion.angleType1.angle[0] * OpenMM_RadiansPerDegree);
152 parameters.append(angleTorsion.angleType2.angle[0] * OpenMM_RadiansPerDegree);
153
154 IntArray particles = new IntArray(0);
155 for (int i = 0; i < 4; i++) {
156 int atomIndex = atoms[i].getArrayIndex();
157 atomIndex = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, atomIndex);
158 particles.append(atomIndex);
159 }
160
161 addBond(particles, parameters);
162 parameters.destroy();
163 particles.destroy();
164 }
165
166 int forceGroup = angleTorsionPotentialEnergy.getForceGroup();
167 setForceGroup(forceGroup);
168 logger.info(format(" Angle-Torsions: %10d", angleTorsions.length));
169 logger.fine(format(" Force Group: %10d", forceGroup));
170 }
171
172
173
174
175
176
177
178 public static Force constructForce(OpenMMEnergy openMMEnergy) {
179 AngleTorsionPotentialEnergy angleTorsionPotentialEnergy = openMMEnergy.getAngleTorsionPotentialEnergy();
180 if (angleTorsionPotentialEnergy == null) {
181 return null;
182 }
183 return new AngleTorsionForce(angleTorsionPotentialEnergy);
184 }
185
186
187
188
189
190
191
192
193 public static Force constructForce(int topology, OpenMMDualTopologyEnergy openMMDualTopologyEnergy) {
194 ForceFieldEnergy forceFieldEnergy = openMMDualTopologyEnergy.getForceFieldEnergy(topology);
195 AngleTorsionPotentialEnergy angleTorsionPotentialEnergy = forceFieldEnergy.getAngleTorsionPotentialEnergy();
196 if (angleTorsionPotentialEnergy == null) {
197 return null;
198 }
199 return new AngleTorsionForce(angleTorsionPotentialEnergy, topology, openMMDualTopologyEnergy);
200 }
201
202
203
204
205
206
207
208 public void updateForce(int topology, OpenMMDualTopologyEnergy openMMDualTopologyEnergy) {
209 ForceFieldEnergy forceFieldEnergy = openMMDualTopologyEnergy.getForceFieldEnergy(topology);
210 AngleTorsionPotentialEnergy angleTorsionPotentialEnergy = forceFieldEnergy.getAngleTorsionPotentialEnergy();
211
212 if (angleTorsionPotentialEnergy == null) {
213 return;
214 }
215 AngleTorsion[] angleTorsions = angleTorsionPotentialEnergy.getAngleTorsionArray();
216
217 double scaleDT = openMMDualTopologyEnergy.getTopologyScale(topology);
218
219 int atIndex = 0;
220 for (AngleTorsion angleTorsion : angleTorsions) {
221 double scale = 1.0;
222
223 if (!angleTorsion.applyLambda()) {
224 scale = scaleDT;
225 }
226 double[] constants = angleTorsion.getConstants();
227 DoubleArray parameters = new DoubleArray(0);
228 for (int m = 0; m < 2; m++) {
229 for (int n = 0; n < 3; n++) {
230 int index = (3 * m) + n;
231 parameters.append(constants[index] * OpenMM_KJPerKcal * scale);
232 }
233 }
234 Atom[] atoms = angleTorsion.getAtomArray(true);
235 parameters.append(angleTorsion.angleType1.angle[0] * OpenMM_RadiansPerDegree);
236 parameters.append(angleTorsion.angleType2.angle[0] * OpenMM_RadiansPerDegree);
237
238 IntArray particles = new IntArray(0);
239 for (int i = 0; i < 4; i++) {
240 int atomIndex = atoms[i].getArrayIndex();
241 atomIndex = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, atomIndex);
242 particles.append(atomIndex);
243 }
244
245 setBondParameters(atIndex++, particles, parameters);
246 parameters.destroy();
247 particles.destroy();
248 }
249
250 updateParametersInContext(openMMDualTopologyEnergy.getContext());
251 }
252 }