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.CustomBondForce;
41 import ffx.openmm.DoubleArray;
42 import ffx.openmm.Force;
43 import ffx.potential.ForceFieldEnergy;
44 import ffx.potential.bonded.Bond;
45 import ffx.potential.parameters.BondType;
46 import ffx.potential.terms.BondPotentialEnergy;
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_NmPerAngstrom;
52 import static java.lang.String.format;
53
54
55
56
57 public class BondForce extends CustomBondForce {
58
59 private static final Logger logger = Logger.getLogger(BondForce.class.getName());
60
61
62
63
64
65
66 public BondForce(BondPotentialEnergy bondPotentialEnergy) {
67 super(bondPotentialEnergy.getBondEnergyString());
68 Bond[] bonds = bondPotentialEnergy.getBondArray();
69 addPerBondParameter("r0");
70 addPerBondParameter("k");
71 setName("AmoebaBond");
72
73 double kParameterConversion = OpenMM_KJPerKcal / (OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom);
74 DoubleArray parameters = new DoubleArray(0);
75 for (Bond bond : bonds) {
76 int i1 = bond.getAtom(0).getArrayIndex();
77 int i2 = bond.getAtom(1).getArrayIndex();
78 BondType bondType = bond.bondType;
79 double r0 = bondType.distance * OpenMM_NmPerAngstrom;
80 double k = kParameterConversion * bondType.forceConstant * bond.bondType.bondUnit;
81 parameters.append(r0);
82 parameters.append(k);
83 addBond(i1, i2, parameters);
84 parameters.resize(0);
85 }
86 parameters.destroy();
87
88 int forceGroup = bondPotentialEnergy.getForceGroup();
89 setForceGroup(forceGroup);
90 logger.info(format(" Bonds: %10d", bonds.length));
91 logger.fine(format(" Force Group: %10d", forceGroup));
92 }
93
94
95
96
97
98
99
100
101 public BondForce(BondPotentialEnergy bondPotentialEnergy,
102 int topology, OpenMMDualTopologyEnergy openMMDualTopologyEnergy) {
103 super(bondPotentialEnergy.getBondEnergyString());
104 Bond[] bonds = bondPotentialEnergy.getBondArray();
105 addPerBondParameter("r0");
106 addPerBondParameter("k");
107 setName("AmoebaBond");
108
109 double scale = openMMDualTopologyEnergy.getTopologyScale(topology);
110
111 double kParameterConversion = OpenMM_KJPerKcal / (OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom);
112 DoubleArray parameters = new DoubleArray(0);
113 for (Bond bond : bonds) {
114 int i1 = bond.getAtom(0).getArrayIndex();
115 int i2 = bond.getAtom(1).getArrayIndex();
116 BondType bondType = bond.bondType;
117 double r0 = bondType.distance * OpenMM_NmPerAngstrom;
118 double k = kParameterConversion * bondType.forceConstant * bond.bondType.bondUnit;
119
120 if (!bond.applyLambda()) {
121 k = k * scale;
122 }
123 parameters.append(r0);
124 parameters.append(k);
125 i1 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, i1);
126 i2 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, i2);
127 addBond(i1, i2, parameters);
128 parameters.resize(0);
129 }
130 parameters.destroy();
131
132 int forceGroup = bondPotentialEnergy.getForceGroup();
133 setForceGroup(forceGroup);
134 logger.info(format(" Bonds: %10d", bonds.length));
135 logger.fine(format(" Force Group: %10d", forceGroup));
136 }
137
138
139
140
141
142
143 public static Force constructForce(OpenMMEnergy openMMEnergy) {
144 BondPotentialEnergy bondPotentialEnergy = openMMEnergy.getBondPotentialEnergy();
145 if (bondPotentialEnergy == null) {
146 return null;
147 }
148 return new BondForce(bondPotentialEnergy);
149 }
150
151
152
153
154
155
156
157 public static Force constructForce(int topology, OpenMMDualTopologyEnergy openMMDualTopologyEnergy) {
158 ForceFieldEnergy forceFieldEnergy = openMMDualTopologyEnergy.getForceFieldEnergy(topology);
159 BondPotentialEnergy bondPotentialEnergy = forceFieldEnergy.getBondPotentialEnergy();
160 if (bondPotentialEnergy == null) {
161 return null;
162 }
163 return new BondForce(bondPotentialEnergy, topology, openMMDualTopologyEnergy);
164 }
165
166
167
168
169 public void updateForce(OpenMMEnergy openMMEnergy) {
170 BondPotentialEnergy bondPotentialEnergy = openMMEnergy.getBondPotentialEnergy();
171 if (bondPotentialEnergy == null) {
172 return;
173 }
174 Bond[] bonds = bondPotentialEnergy.getBondArray();
175
176 double kParameterConversion = OpenMM_KJPerKcal / (OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom);
177 DoubleArray parameters = new DoubleArray(0);
178 int index = 0;
179 for (Bond bond : bonds) {
180 int i1 = bond.getAtom(0).getArrayIndex();
181 int i2 = bond.getAtom(1).getArrayIndex();
182 BondType bondType = bond.bondType;
183 double r0 = bondType.distance * OpenMM_NmPerAngstrom;
184 double k = kParameterConversion * bondType.forceConstant * bondType.bondUnit;
185 parameters.append(r0);
186 parameters.append(k);
187 setBondParameters(index++, i1, i2, parameters);
188 parameters.resize(0);
189 }
190 parameters.destroy();
191 updateParametersInContext(openMMEnergy.getContext());
192 }
193
194
195
196
197
198
199
200 public void updateForce(int topology, OpenMMDualTopologyEnergy openMMDualTopologyEnergy) {
201 ForceFieldEnergy forceFieldEnergy = openMMDualTopologyEnergy.getForceFieldEnergy(topology);
202 BondPotentialEnergy bondPotentialEnergy = forceFieldEnergy.getBondPotentialEnergy();
203 if (bondPotentialEnergy == null) {
204 return;
205 }
206 Bond[] bonds = bondPotentialEnergy.getBondArray();
207
208 double scale = openMMDualTopologyEnergy.getTopologyScale(topology);
209
210 double kParameterConversion = OpenMM_KJPerKcal / (OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom);
211 DoubleArray parameters = new DoubleArray(0);
212 int index = 0;
213 for (Bond bond : bonds) {
214 int i1 = bond.getAtom(0).getArrayIndex();
215 int i2 = bond.getAtom(1).getArrayIndex();
216 BondType bondType = bond.bondType;
217 double r0 = bondType.distance * OpenMM_NmPerAngstrom;
218 double k = kParameterConversion * bondType.forceConstant * bondType.bondUnit;
219
220 if (!bond.applyLambda()) {
221 k = k * scale;
222 }
223 parameters.append(r0);
224 parameters.append(k);
225 i1 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, i1);
226 i2 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, i2);
227 setBondParameters(index++, i1, i2, parameters);
228 parameters.resize(0);
229 }
230 parameters.destroy();
231 updateParametersInContext(openMMDualTopologyEnergy.getContext());
232 }
233
234 }