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.IMPROPER;
41 import static ffx.potential.parameters.ForceField.ForceFieldType.TORSION;
42 import static ffx.utilities.PropertyGroup.EnergyUnitConversion;
43 import static ffx.utilities.PropertyGroup.PotentialFunctionParameter;
44 import static java.lang.Double.parseDouble;
45 import static java.lang.Integer.parseInt;
46 import static java.lang.String.format;
47 import static java.util.Arrays.copyOf;
48 import static org.apache.commons.math3.util.FastMath.cos;
49 import static org.apache.commons.math3.util.FastMath.sin;
50 import static org.apache.commons.math3.util.FastMath.toRadians;
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 = "improper", clazz = String.class, propertyGroup = PotentialFunctionParameter, description = """
67 [4 integers and 2 reals]"
68 Provides the values for a single CHARMM-style improper dihedral angle parameter.
69 The integer modifiers give the atom class numbers for the four kinds of atoms involved in the torsion which is to be defined.
70 The real number modifiers give the force constant value for the deviation from the target improper torsional angle, and the target value for the torsional angle, respectively.
71 The default units for the improper force constant are kcal/mole/radian^2, but this can be controlled via the impropunit keyword.
72 """)
73 @FFXProperty(name = "torsion", clazz = String.class, propertyGroup = PotentialFunctionParameter, description = """
74 [4 integers and up to 6 real/real/integer triples]
75 Provides the values for a single torsional angle parameter.
76 The first four integer modifiers give the atom class numbers for the atoms involved in the torsional angle to be defined.
77 Each of the remaining triples of real/real/integer modifiers give the amplitude,
78 phase offset in degrees and periodicity of a particular torsional function term, respectively.
79 Periodicities through 6-fold are allowed for torsional parameters.
80 """)
81 public final class TorsionType extends BaseType implements Comparator<String> {
82
83 private static final Logger logger = Logger.getLogger(TorsionType.class.getName());
84
85
86
87 public final int[] atomClasses;
88
89
90
91 public final int terms;
92
93
94
95 public final double[] amplitude;
96
97
98
99 public final double[] phase;
100
101
102
103 public final double[] cosine;
104
105
106
107 public final double[] sine;
108
109
110
111 public final int[] periodicity;
112
113
114
115 private final TorsionMode torsionMode;
116
117 public static final double DEFAULT_TORSION_UNIT = 1.0;
118
119
120
121 @FFXProperty(name = "torsionunit", propertyGroup = EnergyUnitConversion, defaultValue = "1.0", description = """
122 Sets the scale factor needed to convert the energy value computed by the torsional angle potential into units of kcal/mole.
123 The correct value is force field dependent and typically provided in the header of the master force field parameter file.
124 """)
125 public double torsionUnit = DEFAULT_TORSION_UNIT;
126
127
128
129
130
131
132
133
134
135 public TorsionType(int[] atomClasses, double[] amplitude, double[] phase, int[] periodicity) {
136 this(atomClasses, amplitude, phase, periodicity, TorsionMode.NORMAL);
137 }
138
139
140
141
142
143
144
145
146
147
148 public TorsionType(int[] atomClasses, double[] amplitude, double[] phase, int[] periodicity,
149 TorsionMode torsionMode) {
150 super(TORSION, sortKey(atomClasses));
151 this.atomClasses = atomClasses;
152 int max = 1;
153 for (int i1 : periodicity) {
154 if (i1 > max) {
155 max = i1;
156 }
157 }
158 terms = max;
159 this.amplitude = new double[terms];
160 this.phase = new double[terms];
161 this.periodicity = new int[terms];
162 for (int i = 0; i < terms; i++) {
163 this.periodicity[i] = i + 1;
164 }
165 for (int i = 0; i < amplitude.length; i++) {
166 int j = periodicity[i] - 1;
167 if (j >= 0 && j < terms) {
168 this.amplitude[j] = amplitude[i];
169 this.phase[j] = phase[i];
170 }
171 }
172
173 cosine = new double[terms];
174 sine = new double[terms];
175 for (int i = 0; i < terms; i++) {
176 double angle = toRadians(this.phase[i]);
177 cosine[i] = cos(angle);
178 sine[i] = sin(angle);
179 }
180
181 this.torsionMode = torsionMode;
182 if (torsionMode == TorsionMode.IMPROPER) {
183 forceFieldType = IMPROPER;
184 }
185 }
186
187
188
189
190
191
192
193
194
195 public static TorsionType average(TorsionType torsionType1, TorsionType torsionType2,
196 int[] atomClasses) {
197 if (torsionType1 == null || torsionType2 == null || atomClasses == null) {
198 return null;
199 }
200 int len = torsionType1.amplitude.length;
201 if (len != torsionType2.amplitude.length) {
202 return null;
203 }
204 double[] amplitude = new double[len];
205 double[] phase = new double[len];
206 int[] periodicity = new int[len];
207 for (int i = 0; i < len; i++) {
208 amplitude[i] = (torsionType1.amplitude[i] + torsionType2.amplitude[i]) / 2.0;
209 phase[i] = (torsionType1.phase[i] + torsionType2.phase[i]) / 2.0;
210 periodicity[i] = (torsionType1.periodicity[i] + torsionType2.periodicity[i]) / 2;
211 }
212
213 return new TorsionType(atomClasses, amplitude, phase, periodicity);
214 }
215
216
217
218
219
220
221
222
223 public static TorsionType parse(String input, String[] tokens) {
224 if (tokens.length < 5) {
225 logger.log(Level.WARNING, "Invalid TORSION type:\n{0}", input);
226 } else {
227 try {
228 int[] atomClasses = new int[4];
229 atomClasses[0] = parseInt(tokens[1]);
230 atomClasses[1] = parseInt(tokens[2]);
231 atomClasses[2] = parseInt(tokens[3]);
232 atomClasses[3] = parseInt(tokens[4]);
233 int terms = (tokens.length - 5) / 3;
234 double[] amplitude = new double[terms];
235 double[] phase = new double[terms];
236 int[] periodicity = new int[terms];
237 int index = 5;
238 for (int i = 0; i < terms; i++) {
239 amplitude[i] = parseDouble(tokens[index++]);
240 phase[i] = parseDouble(tokens[index++]);
241 periodicity[i] = parseInt(tokens[index++]);
242 }
243 return new TorsionType(atomClasses, amplitude, phase, periodicity);
244 } catch (NumberFormatException e) {
245 String message = "NumberFormatException parsing TORSION type:\n" + input + "\n";
246 logger.log(Level.SEVERE, message, e);
247 } catch (Exception e) {
248 String message = "Exception parsing TORSION type:\n" + input + "\n";
249 logger.log(Level.SEVERE, message, e);
250 }
251 }
252 return null;
253 }
254
255
256
257
258
259
260
261
262 public static TorsionType parseImproper(String input, String[] tokens) {
263 if (tokens.length < 5) {
264 logger.log(Level.WARNING, "Invalid IMPROPER type:\n{0}", input);
265 } else {
266 try {
267 int[] atomClasses = new int[4];
268 atomClasses[0] = parseInt(tokens[1]);
269 atomClasses[1] = parseInt(tokens[2]);
270 atomClasses[2] = parseInt(tokens[3]);
271 atomClasses[3] = parseInt(tokens[4]);
272 double[] amplitude = new double[1];
273 double[] phase = new double[1];
274 int[] periodicity = new int[1];
275 int index = 5;
276 amplitude[0] = parseDouble(tokens[index++]);
277 phase[0] = parseDouble(tokens[index]);
278 periodicity[0] = 1;
279 return new TorsionType(atomClasses, amplitude, phase, periodicity, TorsionMode.IMPROPER);
280 } catch (NumberFormatException e) {
281 String message = "Exception parsing IMPROPER type:\n" + input + "\n";
282 logger.log(Level.SEVERE, message, e);
283 }
284 }
285 return null;
286 }
287
288
289
290
291
292
293
294
295 public static String sortKey(int[] c) {
296 if (c == null || c.length != 4) {
297 return null;
298 }
299 if (c[1] < c[2]) {
300
301 } else if (c[2] < c[1]) {
302
303 int temp = c[0];
304 c[0] = c[3];
305 c[3] = temp;
306 temp = c[1];
307 c[1] = c[2];
308 c[2] = temp;
309 } else if (c[1] == c[2]) {
310 if (c[0] > c[3]) {
311
312 int temp = c[0];
313 c[0] = c[3];
314 c[3] = temp;
315 }
316 } else if (c[0] <= c[3]) {
317
318 } else {
319
320 int temp = c[0];
321 c[0] = c[3];
322 c[3] = temp;
323 temp = c[1];
324 c[1] = c[2];
325 c[2] = temp;
326 }
327
328 return c[0] + " " + c[1] + " " + c[2] + " " + c[3];
329 }
330
331
332
333
334
335
336 @Override
337 public int compare(String s1, String s2) {
338 String[] keys1 = s1.split(" ");
339 String[] keys2 = s2.split(" ");
340 int[] c1 = new int[4];
341 int[] c2 = new int[4];
342
343 for (int i = 0; i < 4; i++) {
344 c1[i] = parseInt(keys1[i]);
345 c2[i] = parseInt(keys2[i]);
346 }
347
348 if (c1[1] < c2[1]) {
349 return -1;
350 } else if (c1[1] > c2[1]) {
351 return 1;
352 } else if (c1[2] < c2[2]) {
353 return -1;
354 } else if (c1[2] > c2[2]) {
355 return 1;
356 } else if (c1[0] < c2[0]) {
357 return -1;
358 } else if (c1[0] > c2[0]) {
359 return 1;
360 } else if (c1[3] < c2[3]) {
361 return -1;
362 } else if (c1[3] > c2[3]) {
363 return 1;
364 }
365
366 return 0;
367 }
368
369
370
371
372
373
374
375
376 @Override
377 public boolean equals(Object o) {
378 if (this == o) {
379 return true;
380 }
381 if (o == null || getClass() != o.getClass()) {
382 return false;
383 }
384 TorsionType torsionType = (TorsionType) o;
385 return Arrays.equals(atomClasses, torsionType.atomClasses);
386 }
387
388
389
390
391
392
393
394
395 @Override
396 public int hashCode() {
397 return Arrays.hashCode(atomClasses);
398 }
399
400
401
402
403
404
405 public void incrementClasses(int increment) {
406 for (int i = 0; i < atomClasses.length; i++) {
407 if (atomClasses[i] != 0) {
408 atomClasses[i] += increment;
409 }
410 }
411 setKey(sortKey(atomClasses));
412 }
413
414
415
416
417
418
419
420 public TorsionType patchClasses(HashMap<AtomType, AtomType> typeMap) {
421 int count = 0;
422 int len = atomClasses.length;
423
424
425 for (AtomType newType : typeMap.keySet()) {
426 for (int atomClass : atomClasses) {
427 if (atomClass == newType.atomClass) {
428 count++;
429 }
430 }
431 }
432
433
434 if (count == 1 || count == 2 || count == 3) {
435 int[] newClasses = copyOf(atomClasses, len);
436 for (AtomType newType : typeMap.keySet()) {
437 for (int i = 0; i < len; i++) {
438 if (atomClasses[i] == newType.atomClass) {
439 AtomType knownType = typeMap.get(newType);
440 newClasses[i] = knownType.atomClass;
441 }
442 }
443 }
444 return new TorsionType(newClasses, amplitude, phase, periodicity);
445 }
446 return null;
447 }
448
449
450
451
452
453
454
455
456 @Override
457 public String toString() {
458 StringBuilder torsionBuffer;
459 if (torsionMode == TorsionMode.IMPROPER) {
460 torsionBuffer = new StringBuilder("improper");
461 } else {
462 torsionBuffer = new StringBuilder("torsion");
463 }
464 for (int i : atomClasses) {
465 torsionBuffer.append(format(" %5d", i));
466 }
467
468 boolean nonZero = false;
469 for (double v : amplitude) {
470 if (v != 0.0) {
471 nonZero = true;
472 break;
473 }
474 }
475
476 for (int i = 0; i < amplitude.length; i++) {
477 if (amplitude[i] == 0.0 && (i > 0 || nonZero)) {
478 continue;
479 }
480 if (torsionMode == TorsionMode.NORMAL) {
481 torsionBuffer.append(format(" %7.3f %7.3f %1d", amplitude[i], phase[i], periodicity[i]));
482 } else {
483 torsionBuffer.append(format(" %7.3f %7.3f", amplitude[i], phase[i]));
484 }
485 }
486 return torsionBuffer.toString();
487 }
488
489
490
491
492 public enum TorsionMode {
493 NORMAL, IMPROPER
494 }
495 }