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