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