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.PiOrbitalTorsion;
46 import ffx.potential.parameters.PiOrbitalTorsionType;
47 import ffx.potential.terms.PiOrbitalTorsionPotentialEnergy;
48
49 import java.util.logging.Logger;
50
51 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_KJPerKcal;
52 import static java.lang.String.format;
53
54
55
56
57 public class PiOrbitalTorsionForce extends CustomCompoundBondForce {
58
59 private static final Logger logger = Logger.getLogger(PiOrbitalTorsionForce.class.getName());
60
61
62
63
64
65
66 public PiOrbitalTorsionForce(PiOrbitalTorsionPotentialEnergy piOrbitalTorsionPotentialEnergy) {
67 super(6, PiOrbitalTorsionPotentialEnergy.getPiOrbitalTorsionEnergyString());
68 PiOrbitalTorsion[] piOrbitalTorsions = piOrbitalTorsionPotentialEnergy.getPiOrbitalTorsionArray();
69 addPerBondParameter("k");
70 setName("PiOrbitalTorsion");
71
72 IntArray particles = new IntArray(0);
73 DoubleArray parameters = new DoubleArray(0);
74 for (PiOrbitalTorsion piOrbitalTorsion : piOrbitalTorsions) {
75 int a1 = piOrbitalTorsion.getAtom(0).getArrayIndex();
76 int a2 = piOrbitalTorsion.getAtom(1).getArrayIndex();
77 int a3 = piOrbitalTorsion.getAtom(2).getArrayIndex();
78 int a4 = piOrbitalTorsion.getAtom(3).getArrayIndex();
79 int a5 = piOrbitalTorsion.getAtom(4).getArrayIndex();
80 int a6 = piOrbitalTorsion.getAtom(5).getArrayIndex();
81 PiOrbitalTorsionType type = piOrbitalTorsion.piOrbitalTorsionType;
82 double k = OpenMM_KJPerKcal * type.forceConstant * piOrbitalTorsion.piOrbitalTorsionType.piTorsUnit;
83 particles.append(a1);
84 particles.append(a2);
85 particles.append(a3);
86 particles.append(a4);
87 particles.append(a5);
88 particles.append(a6);
89 parameters.append(k);
90 addBond(particles, parameters);
91 particles.resize(0);
92 parameters.resize(0);
93 }
94 particles.destroy();
95 parameters.destroy();
96
97 int forceGroup = piOrbitalTorsionPotentialEnergy.getForceGroup();
98 setForceGroup(forceGroup);
99 logger.info(format(" Pi-Orbital Torsions: %10d", piOrbitalTorsions.length));
100 logger.fine(format(" Force Group: %10d", forceGroup));
101 }
102
103
104
105
106
107
108
109
110
111 public PiOrbitalTorsionForce(PiOrbitalTorsionPotentialEnergy piOrbitalTorsionPotentialEnergy,
112 int topology, OpenMMDualTopologyEnergy openMMDualTopologyEnergy) {
113 super(6, PiOrbitalTorsionPotentialEnergy.getPiOrbitalTorsionEnergyString());
114 PiOrbitalTorsion[] piOrbitalTorsions = piOrbitalTorsionPotentialEnergy.getPiOrbitalTorsionArray();
115 addPerBondParameter("k");
116 setName("PiOrbitalTorsion");
117
118 double scale = openMMDualTopologyEnergy.getTopologyScale(topology);
119
120 IntArray particles = new IntArray(0);
121 DoubleArray parameters = new DoubleArray(0);
122 for (PiOrbitalTorsion piOrbitalTorsion : piOrbitalTorsions) {
123 int a1 = piOrbitalTorsion.getAtom(0).getArrayIndex();
124 int a2 = piOrbitalTorsion.getAtom(1).getArrayIndex();
125 int a3 = piOrbitalTorsion.getAtom(2).getArrayIndex();
126 int a4 = piOrbitalTorsion.getAtom(3).getArrayIndex();
127 int a5 = piOrbitalTorsion.getAtom(4).getArrayIndex();
128 int a6 = piOrbitalTorsion.getAtom(5).getArrayIndex();
129 PiOrbitalTorsionType type = piOrbitalTorsion.piOrbitalTorsionType;
130 double k = OpenMM_KJPerKcal * type.forceConstant * piOrbitalTorsion.piOrbitalTorsionType.piTorsUnit;
131
132 if (!piOrbitalTorsion.applyLambda()) {
133 k = k * scale;
134 }
135 a1 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, a1);
136 a2 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, a2);
137 a3 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, a3);
138 a4 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, a4);
139 a5 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, a5);
140 a6 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, a6);
141 particles.append(a1);
142 particles.append(a2);
143 particles.append(a3);
144 particles.append(a4);
145 particles.append(a5);
146 particles.append(a6);
147 parameters.append(k);
148 addBond(particles, parameters);
149 particles.resize(0);
150 parameters.resize(0);
151 }
152 particles.destroy();
153 parameters.destroy();
154
155 int forceGroup = piOrbitalTorsionPotentialEnergy.getForceGroup();
156 setForceGroup(forceGroup);
157 logger.info(format(" Pi-Orbital Torsions: %10d", piOrbitalTorsions.length));
158 logger.fine(format(" Force Group: %10d", forceGroup));
159 }
160
161
162
163
164
165
166
167 public static Force constructForce(OpenMMEnergy openMMEnergy) {
168 PiOrbitalTorsionPotentialEnergy piOrbitalTorsionPotentialEnergy = openMMEnergy.getPiOrbitalTorsionPotentialEnergy();
169 if (piOrbitalTorsionPotentialEnergy == null) {
170 return null;
171 }
172 return new PiOrbitalTorsionForce(piOrbitalTorsionPotentialEnergy);
173 }
174
175
176
177
178
179
180
181
182 public static Force constructForce(int topology, OpenMMDualTopologyEnergy openMMDualTopologyEnergy) {
183 ForceFieldEnergy forceFieldEnergy = openMMDualTopologyEnergy.getForceFieldEnergy(topology);
184 PiOrbitalTorsionPotentialEnergy piOrbitalTorsionPotentialEnergy = forceFieldEnergy.getPiOrbitalTorsionPotentialEnergy();
185 if (piOrbitalTorsionPotentialEnergy == null) {
186 return null;
187 }
188 return new PiOrbitalTorsionForce(piOrbitalTorsionPotentialEnergy, topology, openMMDualTopologyEnergy);
189 }
190
191
192
193
194
195
196 public void updateForce(OpenMMEnergy openMMEnergy) {
197 PiOrbitalTorsionPotentialEnergy piOrbitalTorsionPotentialEnergy = openMMEnergy.getPiOrbitalTorsionPotentialEnergy();
198 if (piOrbitalTorsionPotentialEnergy == null) {
199 return;
200 }
201 PiOrbitalTorsion[] piOrbitalTorsions = piOrbitalTorsionPotentialEnergy.getPiOrbitalTorsionArray();
202
203 IntArray particles = new IntArray(0);
204 DoubleArray parameters = new DoubleArray(0);
205 int index = 0;
206 for (PiOrbitalTorsion piOrbitalTorsion : piOrbitalTorsions) {
207 int a1 = piOrbitalTorsion.getAtom(0).getArrayIndex();
208 int a2 = piOrbitalTorsion.getAtom(1).getArrayIndex();
209 int a3 = piOrbitalTorsion.getAtom(2).getArrayIndex();
210 int a4 = piOrbitalTorsion.getAtom(3).getArrayIndex();
211 int a5 = piOrbitalTorsion.getAtom(4).getArrayIndex();
212 int a6 = piOrbitalTorsion.getAtom(5).getArrayIndex();
213 PiOrbitalTorsionType type = piOrbitalTorsion.piOrbitalTorsionType;
214 double k = OpenMM_KJPerKcal * type.forceConstant * piOrbitalTorsion.piOrbitalTorsionType.piTorsUnit;
215 particles.append(a1);
216 particles.append(a2);
217 particles.append(a3);
218 particles.append(a4);
219 particles.append(a5);
220 particles.append(a6);
221 parameters.append(k);
222 setBondParameters(index++, particles, parameters);
223 particles.resize(0);
224 parameters.resize(0);
225 }
226 particles.destroy();
227 parameters.destroy();
228 updateParametersInContext(openMMEnergy.getContext());
229 }
230
231
232
233
234
235
236
237 public void updateForce(int topology, OpenMMDualTopologyEnergy openMMDualTopologyEnergy) {
238 ForceFieldEnergy forceFieldEnergy = openMMDualTopologyEnergy.getForceFieldEnergy(topology);
239 PiOrbitalTorsionPotentialEnergy piOrbitalTorsionPotentialEnergy = forceFieldEnergy.getPiOrbitalTorsionPotentialEnergy();
240 if (piOrbitalTorsionPotentialEnergy == null) {
241 return;
242 }
243 PiOrbitalTorsion[] piOrbitalTorsions = piOrbitalTorsionPotentialEnergy.getPiOrbitalTorsionArray();
244
245 double scale = openMMDualTopologyEnergy.getTopologyScale(topology);
246
247 IntArray particles = new IntArray(0);
248 DoubleArray parameters = new DoubleArray(0);
249 int index = 0;
250 for (PiOrbitalTorsion piOrbitalTorsion : piOrbitalTorsions) {
251 int a1 = piOrbitalTorsion.getAtom(0).getArrayIndex();
252 int a2 = piOrbitalTorsion.getAtom(1).getArrayIndex();
253 int a3 = piOrbitalTorsion.getAtom(2).getArrayIndex();
254 int a4 = piOrbitalTorsion.getAtom(3).getArrayIndex();
255 int a5 = piOrbitalTorsion.getAtom(4).getArrayIndex();
256 int a6 = piOrbitalTorsion.getAtom(5).getArrayIndex();
257 PiOrbitalTorsionType type = piOrbitalTorsion.piOrbitalTorsionType;
258 double k = OpenMM_KJPerKcal * type.forceConstant * piOrbitalTorsion.piOrbitalTorsionType.piTorsUnit;
259
260 if (!piOrbitalTorsion.applyLambda()) {
261 k = k * scale;
262 }
263 a1 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, a1);
264 a2 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, a2);
265 a3 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, a3);
266 a4 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, a4);
267 a5 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, a5);
268 a6 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, a6);
269 particles.append(a1);
270 particles.append(a2);
271 particles.append(a3);
272 particles.append(a4);
273 particles.append(a5);
274 particles.append(a6);
275 parameters.append(k);
276 setBondParameters(index++, particles, parameters);
277 particles.resize(0);
278 parameters.resize(0);
279 }
280 particles.destroy();
281 parameters.destroy();
282 updateParametersInContext(openMMDualTopologyEnergy.getContext());
283 }
284 }