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.parameters;
39
40 import ffx.utilities.FFXProperty;
41 import org.w3c.dom.Document;
42 import org.w3c.dom.Element;
43
44 import javax.annotation.Nullable;
45 import java.util.Arrays;
46 import java.util.Comparator;
47 import java.util.HashMap;
48 import java.util.Map;
49 import java.util.logging.Level;
50 import java.util.logging.Logger;
51
52 import static ffx.potential.parameters.ForceField.ForceFieldType.OPBEND;
53 import static ffx.utilities.Constants.DEGREES_PER_RADIAN;
54 import static ffx.utilities.Constants.KCAL_TO_KJ;
55 import static ffx.utilities.PropertyGroup.EnergyUnitConversion;
56 import static ffx.utilities.PropertyGroup.LocalGeometryFunctionalForm;
57 import static ffx.utilities.PropertyGroup.PotentialFunctionParameter;
58 import static java.lang.Double.parseDouble;
59 import static java.lang.Integer.parseInt;
60 import static java.lang.String.format;
61 import static java.lang.String.valueOf;
62 import static java.util.Arrays.copyOf;
63 import static org.apache.commons.math3.util.FastMath.PI;
64 import static org.apache.commons.math3.util.FastMath.pow;
65
66
67
68
69
70
71
72 @FFXProperty(name = "opbend", clazz = String.class, propertyGroup = PotentialFunctionParameter, description = """
73 [4 integers and 1 real]
74 Provides the values for a single out-of-plane bending potential parameter.
75 The first integer modifier is the atom class of the out-of-plane atom and the second integer is the atom class of the central trigonal atom.
76 The third and fourth integers give the atom classes of the two remaining atoms attached to the trigonal atom.
77 Values of zero for the third and fourth integers are treated as wildcards, and can represent any atom type.
78 The real number modifier gives the force constant value for the out-of-plane angle.
79 The default units for the force constant are kcal/mole/radian^2, but this can be controlled via the opbendunit keyword.
80 """)
81 public final class OutOfPlaneBendType extends BaseType implements Comparator<String> {
82
83
84
85
86 public static final double DEFAULT_OPBEND_CUBIC = 0.0;
87
88
89
90 public static final double DEFAULT_OPBEND_QUARTIC = 0.0;
91
92
93
94 public static final double DEFAULT_OPBEND_PENTIC = 0.0;
95
96
97
98 public static final double DEFAULT_OPBEND_SEXTIC = 0.0;
99
100
101
102
103 @FFXProperty(name = "opbend-cubic", propertyGroup = LocalGeometryFunctionalForm, defaultValue = "0.0", description = """
104 Sets the value of the cubic term in the Taylor series expansion form of the out-of-plane bending potential energy.
105 The real number modifier gives the value of the coefficient as a multiple of the quadratic coefficient.
106 This term multiplied by the out-of-plane bending energy unit conversion factor, the force constant,
107 and the cube of the deviation of the out-of-plane angle from zero gives the cubic contribution to the out-of-plane bending energy.
108 The default value in the absence of the opbend-cubic keyword is zero; i.e., the cubic out-of-plane bending term is omitted.
109 """)
110 public double cubic = DEFAULT_OPBEND_CUBIC;
111
112
113
114
115 @FFXProperty(name = "opbend-quartic", propertyGroup = LocalGeometryFunctionalForm, defaultValue = "0.0", description = """
116 Sets the value of the quartic term in the Taylor series expansion form of the out-of-plane bending potential energy.
117 The real number modifier gives the value of the coefficient as a multiple of the quadratic coefficient.
118 This term multiplied by the out-of-plane bending energy unit conversion factor, the force constant,
119 and the forth power of the deviation of the out-of-plane angle from zero gives the quartic contribution to the out-of-plane bending energy.
120 The default value in the absence of the opbend-quartic keyword is zero; i.e., the quartic out-of-plane bending term is omitted.
121 """)
122 public double quartic = DEFAULT_OPBEND_QUARTIC;
123
124
125
126
127 @FFXProperty(name = "opbend-pentic", propertyGroup = LocalGeometryFunctionalForm, defaultValue = "0.0", description = """
128 Sets the value of the fifth power term in the Taylor series expansion form of the out-of-plane bending potential energy.
129 The real number modifier gives the value of the coefficient as a multiple of the quadratic coefficient.
130 This term multiplied by the out-of-plane bending energy unit conversion factor, the force constant,
131 and the fifth power of the deviation of the out-of-plane angle from zero gives the pentic contribution to the out-of-plane bending energy.
132 The default value in the absence of the opbend-pentic keyword is zero; i.e., the pentic out-of-plane bending term is omitted.
133 """)
134 public double pentic = DEFAULT_OPBEND_PENTIC;
135
136
137
138
139 @FFXProperty(name = "opbend-sextic", propertyGroup = LocalGeometryFunctionalForm, defaultValue = "0.0", description = """
140 Sets the value of the sixth power term in the Taylor series expansion form of the out-of-plane bending potential energy.
141 The real number modifier gives the value of the coefficient as a multiple of the quadratic coefficient.
142 This term multiplied by the out-of-plane bending energy unit conversion factor, the force constant,
143 and the sixth power of the deviation of the out-of-plane angle from zero gives the sextic contribution to the out-of-plane bending energy.
144 The default value in the absence of the opbend-sextic keyword is zero; i.e., the sextic out-of-plane bending term is omitted.
145 """)
146 public double sextic = DEFAULT_OPBEND_SEXTIC;
147
148
149
150
151
152
153
154
155
156
157 @FFXProperty(name = "opbendunit", propertyGroup = EnergyUnitConversion, defaultValue = "(Pi/180)^2", description = """
158 Sets the scale factor needed to convert the energy value computed by the out-of-plane bending potential into units of kcal/mole. "
159 The correct value is force field dependent and typically provided in the header of the master force field parameter file.
160 """)
161 public double opBendUnit = DEFAULT_OPBEND_UNIT;
162
163 public static final double DEFAULT_OPBEND_UNIT = pow(PI / 180.0, 2);
164
165
166
167 private static final Logger logger = Logger.getLogger(OutOfPlaneBendType.class.getName());
168
169
170
171 public final int[] atomClasses;
172
173
174
175 public final double forceConstant;
176
177
178
179
180
181
182
183 public OutOfPlaneBendType(int[] atomClasses, double forceConstant) {
184 super(OPBEND, sortKey(atomClasses));
185 this.atomClasses = atomClasses;
186 this.forceConstant = forceConstant;
187 }
188
189
190
191
192
193
194
195
196
197
198 public static OutOfPlaneBendType average(@Nullable OutOfPlaneBendType outOfPlaneBendType1,
199 @Nullable OutOfPlaneBendType outOfPlaneBendType2,
200 @Nullable int[] atomClasses) {
201 if (outOfPlaneBendType1 == null || outOfPlaneBendType2 == null || atomClasses == null) {
202 return null;
203 }
204
205 double forceConstant =
206 (outOfPlaneBendType1.forceConstant + outOfPlaneBendType2.forceConstant) / 2.0;
207
208 return new OutOfPlaneBendType(atomClasses, forceConstant);
209 }
210
211
212
213
214
215
216
217
218 public static OutOfPlaneBendType parse(String input, String[] tokens) {
219 if (tokens.length < 6) {
220 logger.log(Level.WARNING, "Invalid OPBEND type:\n{0}", input);
221 } else {
222 try {
223 int[] atomClasses = new int[4];
224 atomClasses[0] = parseInt(tokens[1]);
225 atomClasses[1] = parseInt(tokens[2]);
226 atomClasses[2] = parseInt(tokens[3]);
227 atomClasses[3] = parseInt(tokens[4]);
228 double forceConstant = parseDouble(tokens[5]);
229 return new OutOfPlaneBendType(atomClasses, forceConstant);
230 } catch (NumberFormatException e) {
231 String message = "Exception parsing OPBEND type:\n" + input + "\n";
232 logger.log(Level.SEVERE, message, e);
233 }
234 }
235 return null;
236 }
237
238
239
240
241
242
243
244 public static String sortKey(int[] c) {
245 if (c == null || c.length != 4) {
246 return null;
247 }
248 return c[0] + " " + c[1] + " " + c[2] + " " + c[3];
249 }
250
251
252
253
254 @Override
255 public int compare(String s1, String s2) {
256 String[] keys1 = s1.split(" ");
257 String[] keys2 = s2.split(" ");
258
259 for (int i = 0; i < 4; i++) {
260 int c1 = parseInt(keys1[i]);
261 int c2 = parseInt(keys2[i]);
262 if (c1 < c2) {
263 return -1;
264 } else if (c1 > c2) {
265 return 1;
266 }
267 }
268 return 0;
269 }
270
271
272
273
274 @Override
275 public boolean equals(Object o) {
276 if (this == o) {
277 return true;
278 }
279 if (o == null || getClass() != o.getClass()) {
280 return false;
281 }
282 OutOfPlaneBendType outOfPlaneBendType = (OutOfPlaneBendType) o;
283 return Arrays.equals(atomClasses, outOfPlaneBendType.atomClasses);
284 }
285
286
287
288
289 @Override
290 public int hashCode() {
291 return Arrays.hashCode(atomClasses);
292 }
293
294
295
296
297
298
299 public void incrementClasses(int increment) {
300 for (int i = 0; i < atomClasses.length; i++) {
301 if (atomClasses[i] > 0) {
302 atomClasses[i] += increment;
303 }
304 }
305 setKey(sortKey(atomClasses));
306 }
307
308
309
310
311
312
313
314 public OutOfPlaneBendType patchClasses(HashMap<AtomType, AtomType> typeMap) {
315 int count = 0;
316 int len = atomClasses.length;
317
318
319 for (AtomType newType : typeMap.keySet()) {
320
321 for (int atomClass : atomClasses) {
322 if (atomClass == newType.atomClass) {
323 count++;
324 }
325 }
326 }
327
328
329 if (count == 1) {
330 int[] newClasses = copyOf(atomClasses, len);
331 for (AtomType newType : typeMap.keySet()) {
332 for (int i = 0; i < len; i++) {
333 if (atomClasses[i] == newType.atomClass) {
334 AtomType knownType = typeMap.get(newType);
335 newClasses[i] = knownType.atomClass;
336 }
337 }
338 }
339 return new OutOfPlaneBendType(newClasses, forceConstant);
340 }
341 return null;
342 }
343
344
345
346
347
348
349 @Override
350 public String toString() {
351 return String.format("opbend %5d %5d %5d %5d %6.2f", atomClasses[0], atomClasses[1],
352 atomClasses[2], atomClasses[3], forceConstant);
353 }
354
355
356
357
358
359
360
361 public static Element getXMLForce(Document doc, ForceField forceField) {
362 Map<String, OutOfPlaneBendType> types = forceField.getOutOfPlaneBendTypes();
363 if (!types.values().isEmpty()) {
364 Element node = doc.createElement("AmoebaOutOfPlaneBendForce");
365 node.setAttribute("type", forceField.getString("opbendtype", "ALLINGER"));
366 node.setAttribute("opbend-cubic", valueOf(forceField.getDouble("opbend-cubic", DEFAULT_OPBEND_CUBIC)));
367 node.setAttribute("opbend-quartic", valueOf(forceField.getDouble("opbend-quartic", DEFAULT_OPBEND_QUARTIC)));
368 node.setAttribute("opbend-pentic", valueOf(forceField.getDouble("opbend-pentic", DEFAULT_OPBEND_PENTIC)));
369 node.setAttribute("opbend-sextic", valueOf(forceField.getDouble("opbend-sextic", DEFAULT_OPBEND_SEXTIC)));
370 for (OutOfPlaneBendType outOfPlaneBendType : types.values()) {
371 node.appendChild(outOfPlaneBendType.toXML(doc));
372 }
373 return node;
374 }
375 return null;
376 }
377
378
379
380
381 public Element toXML(Document doc) {
382 Element node = doc.createElement("Angle");
383 int i = 1;
384 for (int ac : atomClasses) {
385 if (ac == 0) {
386 node.setAttribute(format("class%d", i), "");
387 } else {
388 node.setAttribute(format("class%d", i), format("%d", ac));
389 }
390 i++;
391 }
392
393 node.setAttribute("k", format("%f", forceConstant * KCAL_TO_KJ / (DEGREES_PER_RADIAN * DEGREES_PER_RADIAN)));
394 return node;
395 }
396 }