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.ANGLE;
53 import static ffx.potential.parameters.ForceField.ForceFieldType.ANGLEP;
54 import static ffx.utilities.Constants.DEGREES_PER_RADIAN;
55 import static ffx.utilities.Constants.KCAL_TO_KJ;
56 import static ffx.utilities.PropertyGroup.EnergyUnitConversion;
57 import static ffx.utilities.PropertyGroup.LocalGeometryFunctionalForm;
58 import static ffx.utilities.PropertyGroup.PotentialFunctionParameter;
59 import static java.lang.Double.parseDouble;
60 import static java.lang.Integer.parseInt;
61 import static java.lang.String.format;
62 import static java.lang.String.valueOf;
63 import static java.lang.System.arraycopy;
64 import static org.apache.commons.math3.util.FastMath.PI;
65 import static org.apache.commons.math3.util.FastMath.pow;
66
67
68
69
70
71
72
73 @FFXProperty(name = "angle", clazz = String.class, propertyGroup = PotentialFunctionParameter, description = """
74 [3 integers and 4 reals]
75 Provides the values for a single bond angle bending parameter.
76 The integer modifiers give the atom class numbers for the three kinds of atoms involved in the angle which is to be defined.
77 The real number modifiers give the force constant value for the angle and up to three ideal bond angles in degrees.
78 In most cases only one ideal bond angle is given, and that value is used for all occurrences of the specified bond angle.
79 If all three ideal angles are given, the values apply when the central atom of the angle is attached to 0, 1 or 2 additional hydrogen atoms, respectively.
80 This "hydrogen environment" option is provided to implement the corresponding feature of the AMOEBA force field.
81 The default units for the force constant are kcal/mole/radian^2, but this can be controlled via the angleunit keyword.
82 """)
83 @FFXProperty(name = "anglep", clazz = String.class, propertyGroup = PotentialFunctionParameter, description = """
84 [3 integers and 3 reals]
85 Provides the values for a single projected in-plane bond angle bending parameter.
86 The integer modifiers give the atom class numbers for the three kinds of atoms involved in the angle which is to be defined.
87 The real number modifiers give the force constant value for the angle and up to two ideal bond angles in degrees.
88 In most cases only one ideal bond angle is given, and that value is used for all occurrences of the specified bond angle.
89 If all two ideal angles are given, the values apply when the central atom of the angle is attached to 0 or 1 additional hydrogen atoms, respectively.
90 This "hydrogen environment" option is provided to implement the corresponding feature of the AMOEBA force field.
91 The default units for the force constant are kcal/mole/radian^2, but this can be controlled via the angleunit keyword.
92 """)
93 public final class AngleType extends BaseType implements Comparator<String> {
94
95
96
97
98 public static final double DEFAULT_ANGLE_UNIT = pow(PI / 180.0, 2.0);
99
100
101
102 public static final double DEFAULT_ANGLE_CUBIC = 0.0;
103
104
105
106 public static final double DEFAULT_ANGLE_QUARTIC = 0.0;
107
108
109
110 public static final double DEFAULT_ANGLE_PENTIC = 0.0;
111
112
113
114 public static final double DEFAULT_ANGLE_SEXTIC = 0.0;
115
116
117
118
119 @FFXProperty(name = "angleunit", propertyGroup = EnergyUnitConversion, defaultValue = "(Pi/180)^2", description = """
120 Sets the scale factor needed to convert the energy value computed by the bond angle bending potential into units of kcal/mole.
121 The correct value is force field dependent and typically provided in the header of the master force field parameter file.
122 """)
123 public double angleUnit = DEFAULT_ANGLE_UNIT;
124
125
126
127
128 @FFXProperty(name = "angle-cubic", propertyGroup = LocalGeometryFunctionalForm, defaultValue = "0.0", description = """
129 Sets the value of the cubic term in the Taylor series expansion form of the bond angle bending potential energy.
130 The real number modifier gives the value of the coefficient as a multiple of the quadratic coefficient.
131 This term multiplied by the angle bending energy unit conversion factor, the force constant,
132 and the cube of the deviation of the bond angle from its ideal value gives the cubic contribution to the angle bending energy.
133 The default value in the absence of the angle-cubic keyword is zero; i.e., the cubic angle bending term is omitted.
134 """)
135 public double cubic = DEFAULT_ANGLE_CUBIC;
136
137
138
139
140 @FFXProperty(name = "angle-quartic", propertyGroup = LocalGeometryFunctionalForm, defaultValue = "0.0", description = """
141 Sets the value of the quartic term in the Taylor series expansion form of the bond angle bending potential energy.
142 The real number modifier gives the value of the coefficient as a multiple of the quadratic coefficient.
143 This term multiplied by the angle bending energy unit conversion factor, the force constant,
144 and the forth power of the deviation of the bond angle from its ideal value gives the quartic contribution to the angle bending energy.
145 The default value in the absence of the angle-quartic keyword is zero; i.e., the quartic angle bending term is omitted.
146 """)
147 public double quartic = DEFAULT_ANGLE_QUARTIC;
148
149
150
151
152 @FFXProperty(name = "angle-pentic", propertyGroup = LocalGeometryFunctionalForm, defaultValue = "0.0", description = """
153 Sets the value of the fifth power term in the Taylor series expansion form of the bond angle bending potential energy.
154 The real number modifier gives the value of the coefficient as a multiple of the quadratic coefficient.
155 This term multiplied by the angle bending energy unit conversion factor, the force constant,
156 and the fifth power of the deviation of the bond angle from its ideal value gives the pentic contribution to the angle bending energy.
157 The default value in the absence of the angle-pentic keyword is zero; i.e., the pentic angle bending term is omitted.
158 """)
159 public double pentic = DEFAULT_ANGLE_PENTIC;
160
161
162
163
164 @FFXProperty(name = "angle-sextic", propertyGroup = LocalGeometryFunctionalForm, defaultValue = "0.0", description = """
165 Sets the value of the sixth power term in the Taylor series expansion form of the bond angle bending potential energy.
166 The real number modifier gives the value of the coefficient as a multiple of the quadratic coefficient.
167 This term multiplied by the angle bending energy unit conversion factor, the force constant,
168 and the sixth power of the deviation of the bond angle from its ideal value gives the sextic contribution to the angle bending energy.
169 The default value in the absence of the angle-sextic keyword is zero; i.e., the sextic angle bending term is omitted.
170 """)
171 public double sextic = DEFAULT_ANGLE_SEXTIC;
172
173
174
175
176 private static final Logger logger = Logger.getLogger(AngleType.class.getName());
177
178
179
180 public final int[] atomClasses;
181
182
183
184 public final double forceConstant;
185
186
187
188
189 public final double[] angle;
190
191
192
193 public final AngleMode angleMode;
194
195
196
197 public AngleFunction angleFunction;
198
199
200
201
202
203
204
205
206 public AngleType(int[] atomClasses, double forceConstant, double[] angle) {
207 this(atomClasses, forceConstant, angle, AngleMode.NORMAL);
208 }
209
210
211
212
213
214
215
216
217
218 public AngleType(int[] atomClasses, double forceConstant, double[] angle, AngleMode angleMode) {
219 this(atomClasses, forceConstant, angle, angleMode, AngleFunction.SEXTIC);
220 }
221
222
223
224
225
226
227
228
229
230
231 public AngleType(int[] atomClasses, double forceConstant, double[] angle, AngleMode angleMode,
232 AngleFunction angleFunction) {
233 super(ANGLE, sortKey(atomClasses));
234 this.atomClasses = atomClasses;
235 this.forceConstant = forceConstant;
236 this.angle = angle;
237 this.angleMode = angleMode;
238 this.angleFunction = angleFunction;
239 if (angleMode == AngleMode.IN_PLANE) {
240 forceFieldType = ANGLEP;
241 }
242 }
243
244
245
246
247
248
249
250
251
252
253 public static AngleType average(@Nullable AngleType angleType1,
254 @Nullable AngleType angleType2,
255 @Nullable int[] atomClasses) {
256 if (angleType1 == null || angleType2 == null || atomClasses == null) {
257 return null;
258 }
259 AngleMode angleMode = angleType1.angleMode;
260 if (angleMode != angleType2.angleMode) {
261 return null;
262 }
263 AngleFunction angleFunction = angleType1.angleFunction;
264 if (angleFunction != angleType2.angleFunction) {
265 return null;
266 }
267 int len = angleType1.angle.length;
268 if (len != angleType2.angle.length) {
269 return null;
270 }
271 double forceConstant = (angleType1.forceConstant + angleType2.forceConstant) / 2.0;
272 double[] angle = new double[len];
273 for (int i = 0; i < len; i++) {
274 angle[i] = (angleType1.angle[i] + angleType2.angle[i]) / 2.0;
275 }
276
277 return new AngleType(atomClasses, forceConstant, angle, angleMode, angleFunction);
278 }
279
280
281
282
283
284
285
286
287 public static AngleType parse(String input, String[] tokens) {
288 if (tokens.length < 6) {
289 logger.log(Level.WARNING, "Invalid ANGLE type:\n{0}", input);
290 } else {
291 int[] atomClasses = new int[3];
292 double forceConstant;
293 int angles;
294 double[] bondAngle;
295 try {
296 atomClasses[0] = parseInt(tokens[1]);
297 atomClasses[1] = parseInt(tokens[2]);
298 atomClasses[2] = parseInt(tokens[3]);
299 forceConstant = parseDouble(tokens[4]);
300 angles = tokens.length - 5;
301 bondAngle = new double[angles];
302 for (int i = 0; i < angles; i++) {
303 bondAngle[i] = parseDouble(tokens[5 + i]);
304 }
305 } catch (NumberFormatException e) {
306 String message = "Exception parsing ANGLE type:\n" + input + "\n";
307 logger.log(Level.SEVERE, message, e);
308 return null;
309 }
310 double[] newBondAngle = new double[angles];
311 arraycopy(bondAngle, 0, newBondAngle, 0, angles);
312 return new AngleType(atomClasses, forceConstant, newBondAngle);
313 }
314 return null;
315 }
316
317
318
319
320
321
322
323
324 public static AngleType parseInPlane(String input, String[] tokens) {
325 if (tokens.length < 6) {
326 logger.log(Level.WARNING, "Invalid ANGLEP type:\n{0}", input);
327 } else {
328 int[] atomClasses = new int[3];
329 double forceConstant;
330 int angles;
331 double[] bondAngle;
332 try {
333 atomClasses[0] = parseInt(tokens[1]);
334 atomClasses[1] = parseInt(tokens[2]);
335 atomClasses[2] = parseInt(tokens[3]);
336 forceConstant = parseDouble(tokens[4]);
337 angles = tokens.length - 5;
338 bondAngle = new double[angles];
339 for (int i = 0; i < angles; i++) {
340 bondAngle[i] = parseDouble(tokens[5 + i]);
341 }
342 } catch (NumberFormatException e) {
343 String message = "Exception parsing ANGLEP type:\n" + input + "\n";
344 logger.log(Level.SEVERE, message, e);
345 return null;
346 }
347 double[] newBondAngle = new double[angles];
348 arraycopy(bondAngle, 0, newBondAngle, 0, angles);
349 return new AngleType(atomClasses, forceConstant, newBondAngle, AngleMode.IN_PLANE);
350 }
351 return null;
352 }
353
354
355
356
357
358
359
360 public static String sortKey(int[] c) {
361 if (c == null || c.length != 3) {
362 return null;
363 }
364 if (c[0] > c[2]) {
365 int temp = c[0];
366 c[0] = c[2];
367 c[2] = temp;
368 }
369 return c[0] + " " + c[1] + " " + c[2];
370 }
371
372
373
374
375 @Override
376 public int compare(String key1, String key2) {
377 String[] keys1 = key1.split(" ");
378 String[] keys2 = key2.split(" ");
379 int[] c1 = new int[3];
380 int[] c2 = new int[3];
381 for (int i = 0; i < 3; i++) {
382 c1[i] = Integer.parseInt(keys1[i]);
383 c2[i] = Integer.parseInt(keys2[i]);
384 }
385
386 if (c1[1] < c2[1]) {
387 return -1;
388 } else if (c1[1] > c2[1]) {
389 return 1;
390 } else if (c1[0] < c2[0]) {
391 return -1;
392 } else if (c1[0] > c2[0]) {
393 return 1;
394 } else if (c1[2] < c2[2]) {
395 return -1;
396 } else if (c1[2] > c2[2]) {
397 return 1;
398 }
399
400 return 0;
401 }
402
403
404
405
406 @Override
407 public boolean equals(Object o) {
408 if (this == o) {
409 return true;
410 }
411 if (o == null || getClass() != o.getClass()) {
412 return false;
413 }
414 AngleType angleType = (AngleType) o;
415 return Arrays.equals(atomClasses, angleType.atomClasses);
416 }
417
418
419
420
421 @Override
422 public int hashCode() {
423 return Arrays.hashCode(atomClasses);
424 }
425
426
427
428
429
430
431 public void incrementClasses(int increment) {
432 for (int i = 0; i < atomClasses.length; i++) {
433 atomClasses[i] += increment;
434 }
435 setKey(sortKey(atomClasses));
436 }
437
438
439
440
441
442
443
444 public AngleType patchClasses(HashMap<AtomType, AtomType> typeMap) {
445 int count = 0;
446 int len = atomClasses.length;
447
448
449 for (AtomType newType : typeMap.keySet()) {
450 for (int atomClass : atomClasses) {
451 if (atomClass == newType.atomClass) {
452 count++;
453 }
454 }
455 }
456
457
458 if (count == 1 || count == 2) {
459 int[] newClasses = Arrays.copyOf(atomClasses, len);
460 for (AtomType newType : typeMap.keySet()) {
461 for (int i = 0; i < len; i++) {
462 if (atomClasses[i] == newType.atomClass) {
463 AtomType knownType = typeMap.get(newType);
464 newClasses[i] = knownType.atomClass;
465 }
466 }
467 }
468 return new AngleType(newClasses, forceConstant, angle, angleMode, angleFunction);
469 }
470
471 return null;
472 }
473
474
475
476
477
478
479 public void setAngleFunction(AngleFunction angleFunction) {
480 this.angleFunction = angleFunction;
481 }
482
483
484
485
486
487
488 @Override
489 public String toString() {
490 StringBuilder angleString;
491 if (angleMode == AngleMode.NORMAL) {
492 angleString = new StringBuilder(
493 format("angle %5d %5d %5d %6.2f", atomClasses[0], atomClasses[1], atomClasses[2],
494 forceConstant));
495 } else {
496 angleString = new StringBuilder(
497 format("anglep %5d %5d %5d %6.2f", atomClasses[0], atomClasses[1], atomClasses[2],
498 forceConstant));
499 }
500 for (double eq : angle) {
501 angleString.append(format(" %6.2f", eq));
502 }
503 return angleString.toString();
504 }
505
506
507
508
509
510
511
512
513 public static Element getXMLForce(Document doc, ForceField forceField) {
514 Map<String, AngleType> angMap = forceField.getAngleTypes();
515 angMap.putAll(forceField.getAnglepTypes());
516 if (!angMap.values().isEmpty()) {
517 Element node = doc.createElement("AmoebaAngleForce");
518 node.setAttribute("angle-cubic", valueOf(forceField.getDouble("angle-cubic", DEFAULT_ANGLE_CUBIC)));
519 node.setAttribute("angle-quartic", valueOf(forceField.getDouble("angle-quartic", DEFAULT_ANGLE_QUARTIC)));
520 node.setAttribute("angle-pentic", valueOf(forceField.getDouble("angle-pentic", DEFAULT_ANGLE_PENTIC)));
521 node.setAttribute("angle-sextic", valueOf(forceField.getDouble("angle-sextic", DEFAULT_ANGLE_SEXTIC)));
522 for (AngleType angleType : angMap.values()) {
523 node.appendChild(angleType.toXML(doc));
524 }
525 return node;
526 }
527 return null;
528 }
529
530
531
532
533
534
535
536 public Element toXML(Document doc) {
537 Element node = doc.createElement("Angle");
538 node.setAttribute("class1", format("%d", atomClasses[0]));
539 node.setAttribute("class2", format("%d", atomClasses[1]));
540 node.setAttribute("class3", format("%d", atomClasses[2]));
541
542 node.setAttribute("k", format("%.17f", forceConstant * KCAL_TO_KJ / (DEGREES_PER_RADIAN * DEGREES_PER_RADIAN)));
543 int i = 1;
544 for (double eq : angle) {
545 node.setAttribute(format("angle%d", i), format("%f", eq));
546 i++;
547 }
548 if (angleMode == AngleMode.NORMAL) {
549 node.setAttribute("inPlane", "False");
550 } else {
551 node.setAttribute("inPlane", "True");
552 }
553 return node;
554 }
555
556
557
558
559 public enum AngleFunction {
560 HARMONIC, SEXTIC
561 }
562
563
564
565
566 public enum AngleMode {
567 NORMAL, IN_PLANE
568 }
569 }