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.BOND;
41 import static ffx.utilities.PropertyGroup.EnergyUnitConversion;
42 import static ffx.utilities.PropertyGroup.LocalGeometryFunctionalForm;
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
48 import ffx.utilities.FFXProperty;
49
50 import java.util.Arrays;
51 import java.util.Comparator;
52 import java.util.HashMap;
53 import java.util.logging.Level;
54 import java.util.logging.Logger;
55
56
57
58
59
60
61
62 @FFXProperty(name = "bond", clazz = String.class, propertyGroup = PotentialFunctionParameter, description = """
63 [2 integers and 2 reals]
64 Provides the values for a single bond stretching parameter.
65 The integer modifiers give the atom class numbers for the two kinds of atoms involved in the bond which is to be defined.
66 The real number modifiers give the force constant value for the bond and the ideal bond length in Angstroms.
67 The default value of 1.0 is used, if the bondunit keyword is not given in the force field parameter file or the keyfile.
68 """)
69 public final class BondType extends BaseType implements Comparator<String> {
70
71 public static final double DEFAULT_BOND_UNIT = 1.0;
72
73
74
75
76 public static final double DEFAULT_BOND_CUBIC = 0.0;
77
78
79
80 public static final double DEFAULT_BOND_QUARTIC = 0.0;
81
82
83
84
85 @FFXProperty(name = "bondunit", propertyGroup = EnergyUnitConversion, defaultValue = "1.0", description = """
86 Sets the scale factor needed to convert the energy value computed by the bond stretching potential into units of kcal/mole.
87 The correct value is force field dependent and typically provided in the header of the master force field parameter file.
88 """)
89 public double bondUnit = DEFAULT_BOND_UNIT;
90
91
92
93
94 @FFXProperty(name = "bond-cubic", propertyGroup = LocalGeometryFunctionalForm, defaultValue = "0.0", description = """
95 Sets the value of the cubic term in the Taylor series expansion form of the bond stretching potential energy.
96 The real number modifier gives the value of the coefficient as a multiple of the quadratic coefficient.
97 This term multiplied by the bond stretching energy unit conversion factor, the force constant,
98 and the cube of the deviation of the bond length from its ideal value gives the cubic contribution to the bond stretching energy.
99 The default value in the absence of the bond-cubic keyword is zero; i.e., the cubic bond stretching term is omitted.
100 """)
101 public double cubic = DEFAULT_BOND_CUBIC;
102
103
104
105
106 @FFXProperty(name = "bond-quartic", propertyGroup = LocalGeometryFunctionalForm, defaultValue = "0.0", description = """
107 Sets the value of the quartic term in the Taylor series expansion form of the bond stretching potential energy.
108 The real number modifier gives the value of the coefficient as a multiple of the quadratic coefficient.
109 This term multiplied by the bond stretching energy unit conversion factor, the force constant,
110 and the forth power of the deviation of the bond length from its ideal value gives the quartic contribution to the bond stretching energy.
111 The default value in the absence of the bond-quartic keyword is zero; i.e., the quartic bond stretching term is omitted.
112 """)
113 public double quartic = DEFAULT_BOND_QUARTIC;
114
115
116
117
118 private static final Logger logger = Logger.getLogger(BondType.class.getName());
119
120
121
122 public final int[] atomClasses;
123
124
125
126 public final double forceConstant;
127
128
129
130 public final double distance;
131
132
133
134
135 public final double flatBottomRadius;
136
137
138
139 public BondFunction bondFunction;
140
141
142
143
144
145
146
147
148 public BondType(int[] atomClasses, double forceConstant, double distance) {
149 this(atomClasses, forceConstant, distance, BondFunction.QUARTIC, 0.0);
150 }
151
152
153
154
155
156
157
158
159
160 public BondType(int[] atomClasses, double forceConstant, double distance,
161 BondFunction bondFunction) {
162 this(atomClasses, forceConstant, distance, bondFunction, 0.0);
163 }
164
165
166
167
168
169
170
171
172
173
174 public BondType(int[] atomClasses, double forceConstant, double distance, BondFunction bondFunction,
175 double flatBottomRadius) {
176 super(BOND, sortKey(atomClasses));
177 this.atomClasses = atomClasses;
178 this.forceConstant = forceConstant;
179 this.distance = distance;
180 this.bondFunction = bondFunction;
181 this.flatBottomRadius = flatBottomRadius;
182 assert (flatBottomRadius == 0 || bondFunction == BondFunction.FLAT_BOTTOM_HARMONIC);
183 }
184
185
186
187
188
189
190
191
192
193 public static BondType average(BondType bondType1, BondType bondType2, int[] atomClasses) {
194 if (bondType1 == null || bondType2 == null || atomClasses == null) {
195 return null;
196 }
197 BondType.BondFunction bondFunction = bondType1.bondFunction;
198 if (bondFunction != bondType2.bondFunction) {
199 return null;
200 }
201
202 double forceConstant = (bondType1.forceConstant + bondType2.forceConstant) / 2.0;
203 double distance = (bondType1.distance + bondType2.distance) / 2.0;
204
205 return new BondType(atomClasses, forceConstant, distance, bondFunction);
206 }
207
208
209
210
211
212
213
214
215 public static BondType parse(String input, String[] tokens) {
216 if (tokens.length < 5) {
217 logger.log(Level.WARNING, "Invalid BOND type:\n{0}", input);
218 } else {
219 try {
220 int[] atomClasses = new int[2];
221 atomClasses[0] = parseInt(tokens[1]);
222 atomClasses[1] = parseInt(tokens[2]);
223 double forceConstant = parseDouble(tokens[3]);
224 double distance = parseDouble(tokens[4]);
225 return new BondType(atomClasses, forceConstant, distance);
226 } catch (NumberFormatException e) {
227 String message = "Exception parsing BOND type:\n" + input + "\n";
228 logger.log(Level.SEVERE, message, e);
229 }
230 }
231 return null;
232 }
233
234
235
236
237
238
239
240 public static String sortKey(int[] c) {
241 if (c == null || c.length != 2) {
242 return null;
243 }
244
245 int temp;
246 if (c[1] <= c[0]) {
247 temp = c[1];
248 c[1] = c[0];
249 c[0] = temp;
250 }
251
252 return c[0] + " " + c[1];
253 }
254
255
256
257
258 @Override
259 public int compare(String key1, String key2) {
260 String[] keys1 = key1.split(" ");
261 String[] keys2 = key2.split(" ");
262 int[] c1 = new int[2];
263 int[] c2 = new int[2];
264 for (int i = 0; i < 2; i++) {
265 c1[i] = parseInt(keys1[i]);
266 c2[i] = parseInt(keys2[i]);
267 }
268
269 if (c1[0] < c2[0]) {
270 return -1;
271 } else if (c1[0] > c2[0]) {
272 return 1;
273 } else if (c1[1] < c2[1]) {
274 return -1;
275 } else if (c1[1] > c2[1]) {
276 return 1;
277 }
278
279 return 0;
280 }
281
282
283
284
285 @Override
286 public boolean equals(Object o) {
287 if (this == o) {
288 return true;
289 }
290 if (o == null || getClass() != o.getClass()) {
291 return false;
292 }
293 BondType bondType = (BondType) o;
294 return Arrays.equals(atomClasses, bondType.atomClasses);
295 }
296
297
298
299
300 @Override
301 public int hashCode() {
302 return Arrays.hashCode(atomClasses);
303 }
304
305
306
307
308
309
310 public void incrementClasses(int increment) {
311 for (int i = 0; i < atomClasses.length; i++) {
312 atomClasses[i] += increment;
313 }
314 setKey(sortKey(atomClasses));
315 }
316
317
318
319
320
321
322
323 public BondType patchClasses(HashMap<AtomType, AtomType> typeMap) {
324
325 int count = 0;
326 int len = atomClasses.length;
327
328
329 for (AtomType newType : typeMap.keySet()) {
330
331 for (int atomClass : atomClasses) {
332 if (atomClass == newType.atomClass) {
333 count++;
334 }
335 }
336 }
337
338
339 if (count == 1) {
340 int[] newClasses = Arrays.copyOf(atomClasses, len);
341 for (AtomType newType : typeMap.keySet()) {
342 for (int i = 0; i < len; i++) {
343 if (atomClasses[i] == newType.atomClass) {
344 AtomType knownType = typeMap.get(newType);
345 newClasses[i] = knownType.atomClass;
346 }
347 }
348 }
349
350 return new BondType(newClasses, forceConstant, distance, bondFunction);
351 }
352 return null;
353 }
354
355 public void setBondFunction(BondFunction bondFunction) {
356 this.bondFunction = bondFunction;
357 }
358
359
360
361
362
363
364 @Override
365 public String toString() {
366 return format("bond %5d %5d %7.2f %7.4f", atomClasses[0], atomClasses[1], forceConstant,
367 distance);
368 }
369
370
371
372
373
374 public enum BondFunction {
375
376
377 HARMONIC("0.5*k*(r-r0)^2", false),
378
379
380 QUARTIC("0.5*k*dv^2*((1+cubic)*dv+(1+quartic)*dv^2);dv=r-r0", false),
381
382
383 FLAT_BOTTOM_HARMONIC(
384 "0.5*k*dv^2;dv=step(dv)*step(dv-fb)*(dv-fb)" + "+step(-dv)*step(-dv-fb)*(-dv-fb);dv=r-r0",
385 true),
386
387
388 FLAT_BOTTOM_QUARTIC("0.5*k*dv^2*((1+cubic)*dv+(1+quartic)*dv^2);"
389 + "dv=step(dv)*step(dv-fb)*(dv-fb)+step(-dv)*step(-dv-fb)*(-dv-fb);dv=r-r0", true);
390
391
392
393
394 private final String mathematicalForm;
395
396
397
398
399 private final boolean hasFlatBottom;
400
401
402
403
404
405
406
407 BondFunction(String mathForm, boolean flatBottom) {
408 this.mathematicalForm = mathForm;
409 this.hasFlatBottom = flatBottom;
410 }
411
412
413
414
415
416
417 public boolean hasFlatBottom() {
418 return hasFlatBottom;
419 }
420
421
422
423
424
425
426 public String toMathematicalForm() {
427 return mathematicalForm;
428 }
429 }
430 }