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