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.potential.bonded.Atom;
41 import ffx.potential.bonded.Bond;
42 import ffx.utilities.FFXProperty;
43 import org.w3c.dom.Document;
44 import org.w3c.dom.Element;
45
46 import java.util.ArrayList;
47 import java.util.Arrays;
48 import java.util.Collections;
49 import java.util.Comparator;
50 import java.util.HashMap;
51 import java.util.List;
52 import java.util.Objects;
53 import java.util.logging.Level;
54 import java.util.logging.Logger;
55
56 import static ffx.potential.parameters.ForceField.ForceFieldType.POLARIZE;
57 import static ffx.utilities.Constants.ANG_TO_NM;
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.pow;
65
66
67
68
69
70
71
72
73
74 @FFXProperty(name = "polarize", clazz = String.class, propertyGroup = PotentialFunctionParameter, description = """
75 [1 integer, up to 3 reals and up to 8 integers]
76 Provides the values for a single atomic dipole polarizability parameter.
77 The initial integer modifier, if positive, gives the atom type number for which a polarizability parameter is to be defined.
78 If the first integer modifier is negative, then the parameter value to follow applies only to the specific atom whose atom number is the negative of the modifier.
79 The first real number modifier gives the value of the dipole polarizability in Ang^3.
80 The second real number modifier, if present, gives the Thole damping value.
81 A Thole value of zero implies undamped polarization.
82 The third real modifier, if present, gives a direct field damping value only used with the AMOEBA+ polarization model.
83 The remaining integer modifiers list the atom type numbers of atoms directly bonded to the current atom and which will be considered to be part of the current atom’s polarization group.
84 If the parameter is for a specific atom, then the integers defining the polarization group are ignored.
85 """)
86 public final class PolarizeType extends BaseType implements Comparator<String> {
87
88 private static final Logger logger = Logger.getLogger(PolarizeType.class.getName());
89
90 private static final double sixth = 1.0 / 6.0;
91 public static final double DEFAULT_DIRECT_11_SCALE = 0.0;
92 public static final double DEFAULT_DIRECT_12_SCALE = 1.0;
93 public static final double DEFAULT_DIRECT_13_SCALE = 1.0;
94 public static final double DEFAULT_DIRECT_14_SCALE = 1.0;
95 public static final double DEFAULT_POLAR_12_SCALE = 0.0;
96 public static final double DEFAULT_POLAR_13_SCALE = 0.0;
97 public static final double DEFAULT_POLAR_14_SCALE = 1.0;
98 public static final double DEFAULT_POLAR_15_SCALE = 1.0;
99 public static final double DEFAULT_POLAR_12_INTRA = 0.0;
100 public static final double DEFAULT_POLAR_13_INTRA = 0.0;
101 public static final double DEFAULT_POLAR_14_INTRA = 0.5;
102 public static final double DEFAULT_POLAR_15_INTRA = 1.0;
103
104
105
106 public final double thole;
107
108
109
110 public double pdamp;
111
112
113
114 public final double ddp;
115
116
117
118 public final double polarizability;
119
120
121
122 public int type;
123
124
125
126 public int[] polarizationGroup;
127
128
129
130
131
132
133
134
135
136 public PolarizeType(int atomType, double polarizability, double thole, double ddp,
137 int[] polarizationGroup) {
138 super(POLARIZE, Integer.toString(atomType));
139 this.type = atomType;
140 this.thole = thole;
141 this.polarizability = polarizability;
142 this.ddp = 0.0;
143 this.polarizationGroup = polarizationGroup;
144 if (thole == 0.0) {
145 pdamp = 0.0;
146 } else {
147 pdamp = pow(polarizability, sixth);
148 }
149 }
150
151
152
153
154
155
156
157 public PolarizeType(PolarizeType polarizeType, double polarizability) {
158 this(polarizeType.type, polarizability, polarizeType.thole, polarizeType.ddp,
159 Arrays.copyOf(polarizeType.polarizationGroup, polarizeType.polarizationGroup.length));
160
161 pdamp = polarizeType.pdamp;
162 }
163
164
165
166
167
168
169
170
171
172
173
174 public static void assignPolarizationGroups(Atom[] atoms, int[][] ip11, int[][] ip12,
175 int[][] ip13) {
176
177
178 List<Integer> group = new ArrayList<>();
179 List<Integer> polarizationGroup = new ArrayList<>();
180 for (Atom ai : atoms) {
181 group.clear();
182 polarizationGroup.clear();
183 int index = ai.getIndex() - 1;
184 group.add(index);
185 polarizationGroup.add(ai.getType());
186 PolarizeType polarizeType = ai.getPolarizeType();
187 if (polarizeType != null) {
188 if (polarizeType.polarizationGroup != null) {
189 for (int i : polarizeType.polarizationGroup) {
190 if (!polarizationGroup.contains(i)) {
191 polarizationGroup.add(i);
192 }
193 }
194 growGroup(polarizationGroup, group, ai);
195 Collections.sort(group);
196 ip11[index] = new int[group.size()];
197 int j = 0;
198 for (int k : group) {
199 ip11[index][j++] = k;
200 }
201 } else {
202 ip11[index] = new int[group.size()];
203 int j = 0;
204 for (int k : group) {
205 ip11[index][j++] = k;
206 }
207 }
208 } else {
209 String message = "The polarize keyword was not found for atom " + (index + 1) + " with type "
210 + ai.getType();
211 logger.severe(message);
212 }
213 }
214
215
216 int nAtoms = atoms.length;
217 int[] mask = new int[nAtoms];
218 List<Integer> list = new ArrayList<>();
219 List<Integer> keep = new ArrayList<>();
220 for (int i = 0; i < nAtoms; i++) {
221 mask[i] = -1;
222 }
223 for (int i = 0; i < nAtoms; i++) {
224 list.clear();
225 for (int j : ip11[i]) {
226 list.add(j);
227 mask[j] = i;
228 }
229 keep.clear();
230 for (int j : list) {
231 Atom aj = atoms[j];
232 List<Bond> bonds = aj.getBonds();
233 for (Bond b : bonds) {
234 Atom ak = b.get1_2(aj);
235 int k = ak.getIndex() - 1;
236 if (mask[k] != i) {
237 keep.add(k);
238 }
239 }
240 }
241 list.clear();
242 for (int j : keep) {
243 for (int k : ip11[j]) {
244 list.add(k);
245 }
246 }
247 Collections.sort(list);
248 ip12[i] = new int[list.size()];
249 int j = 0;
250 for (int k : list) {
251 ip12[i][j++] = k;
252 }
253 }
254
255
256 for (int i = 0; i < nAtoms; i++) {
257 mask[i] = -1;
258 }
259 for (int i = 0; i < nAtoms; i++) {
260 for (int j : ip11[i]) {
261 mask[j] = i;
262 }
263 for (int j : ip12[i]) {
264 mask[j] = i;
265 }
266 list.clear();
267 for (int j : ip12[i]) {
268 for (int k : ip12[j]) {
269 if (mask[k] != i) {
270 if (!list.contains(k)) {
271 list.add(k);
272 }
273 }
274 }
275 }
276 ip13[i] = new int[list.size()];
277 Collections.sort(list);
278 int j = 0;
279 for (int k : list) {
280 ip13[i][j++] = k;
281 }
282 }
283 }
284
285
286
287
288
289
290
291
292
293
294
295 public static PolarizeType average(PolarizeType polarizeType1, PolarizeType polarizeType2,
296 int atomType, int[] polarizationGroup) {
297 if (polarizeType1 == null || polarizeType2 == null) {
298 return null;
299 }
300 double thole = (polarizeType1.thole + polarizeType2.thole) / 2.0;
301 double polarizability = (polarizeType1.polarizability + polarizeType2.polarizability) / 2.0;
302 double ddp = (polarizeType1.ddp + polarizeType2.ddp) / 2.0;
303 return new PolarizeType(atomType, polarizability, thole, ddp, polarizationGroup);
304 }
305
306
307
308
309
310
311
312
313 public static PolarizeType parse(String input, String[] tokens) {
314 if (tokens.length < 4) {
315 logger.log(Level.WARNING, "Invalid POLARIZE type:\n{0}", input);
316 } else {
317 try {
318 int atomType = parseInt(tokens[1]);
319 double polarizability = parseDouble(tokens[2]);
320 double thole = parseDouble(tokens[3]);
321
322
323 double ddp = 0.0;
324 int nextToken = 4;
325 try {
326
327
328 parseInt(tokens[nextToken]);
329 } catch (NumberFormatException e) {
330 ddp = parseDouble(tokens[nextToken]);
331 nextToken = 5;
332 } catch (ArrayIndexOutOfBoundsException e) {
333
334 }
335 int entries = tokens.length - nextToken;
336 int[] polarizationGroup = null;
337 if (entries > 0) {
338 polarizationGroup = new int[entries];
339 for (int i = nextToken; i < tokens.length; i++) {
340 polarizationGroup[i - nextToken] = parseInt(tokens[i]);
341 }
342 }
343 return new PolarizeType(atomType, polarizability, thole, ddp, polarizationGroup);
344 } catch (NumberFormatException e) {
345 String message = "Exception parsing POLARIZE type:\n" + input + "\n";
346 logger.log(Level.SEVERE, message, e);
347 }
348 }
349 return null;
350 }
351
352
353
354
355
356
357
358
359
360 private static void growGroup(List<Integer> polarizationGroup, List<Integer> group, Atom seed) {
361 List<Bond> bonds = seed.getBonds();
362 for (Bond bi : bonds) {
363 Atom aj = bi.get1_2(seed);
364 int tj = aj.getType();
365 boolean added = false;
366 for (int g : polarizationGroup) {
367 if (g == tj) {
368 Integer index = aj.getIndex() - 1;
369 if (!group.contains(index)) {
370 group.add(index);
371 added = true;
372 break;
373 }
374 }
375 }
376 if (added) {
377 PolarizeType polarizeType = aj.getPolarizeType();
378 for (int i : polarizeType.polarizationGroup) {
379 if (!polarizationGroup.contains(i)) {
380 polarizationGroup.add(i);
381 }
382 }
383 growGroup(polarizationGroup, group, aj);
384 }
385 }
386 }
387
388
389
390
391
392
393
394
395 public static void growGroup(List<Integer> group, Atom seed) {
396 List<Bond> bonds = seed.getBonds();
397 for (Bond bond : bonds) {
398 Atom atom = bond.get1_2(seed);
399 int tj = atom.getType();
400 boolean added = false;
401 PolarizeType polarizeType = seed.getPolarizeType();
402 for (int type : polarizeType.polarizationGroup) {
403 if (type == tj) {
404 Integer index = atom.getIndex() - 1;
405 if (!group.contains(index)) {
406 group.add(index);
407 added = true;
408 break;
409 }
410 }
411 }
412 if (added) {
413 growGroup(group, atom);
414 }
415 }
416 }
417
418
419
420
421
422
423 public void add(int atomType) {
424 for (int i : polarizationGroup) {
425 if (atomType == i) {
426 return;
427 }
428 }
429 int len = polarizationGroup.length;
430 int[] newGroup = new int[len + 1];
431 arraycopy(polarizationGroup, 0, newGroup, 0, len);
432 newGroup[len] = atomType;
433 polarizationGroup = newGroup;
434 }
435
436
437
438
439 @Override
440 public int compare(String s1, String s2) {
441 int t1 = parseInt(s1);
442 int t2 = parseInt(s2);
443 return Integer.compare(t1, t2);
444 }
445
446
447
448
449 @Override
450 public boolean equals(Object o) {
451 if (this == o) {
452 return true;
453 }
454 if (o == null || getClass() != o.getClass()) {
455 return false;
456 }
457 PolarizeType polarizeType = (PolarizeType) o;
458 return polarizeType.type == this.type;
459 }
460
461
462
463
464 @Override
465 public int hashCode() {
466 return Objects.hash(type);
467 }
468
469
470
471
472
473
474 @Override
475 public String toString() {
476 StringBuilder polarizeString = new StringBuilder(
477 format("polarize %5d %8.5f %8.5f", type, polarizability, thole));
478 if (ddp != 0.0) {
479 polarizeString.append(format(" %8.5f", ddp));
480 }
481 if (polarizationGroup != null) {
482 for (int a : polarizationGroup) {
483 polarizeString.append(format(" %5d", a));
484 }
485 }
486 return polarizeString.toString();
487 }
488
489
490
491
492
493
494
495 public static void addXMLAttributes(Element node, ForceField forceField) {
496 node.setAttribute("direct11Scale", valueOf(forceField.getDouble("direct-11-scale", DEFAULT_DIRECT_11_SCALE)));
497 node.setAttribute("direct12Scale", valueOf(forceField.getDouble("direct-12-scale", DEFAULT_DIRECT_12_SCALE)));
498 node.setAttribute("direct13Scale", valueOf(forceField.getDouble("direct-13-scale", DEFAULT_DIRECT_13_SCALE)));
499 node.setAttribute("direct14Scale", valueOf(forceField.getDouble("direct-14-scale", DEFAULT_DIRECT_14_SCALE)));
500 node.setAttribute("mutual11Scale", valueOf(forceField.getDouble("mutual-11-scale", 1.0)));
501 node.setAttribute("mutual12Scale", valueOf(forceField.getDouble("mutual-12-scale", 1.0)));
502 node.setAttribute("mutual13Scale", valueOf(forceField.getDouble("mutual-13-scale", 1.0)));
503 node.setAttribute("mutual14Scale", valueOf(forceField.getDouble("mutual-14-scale", 1.0)));
504 node.setAttribute("polar12Scale", valueOf(forceField.getDouble("polar-12-scale", DEFAULT_POLAR_12_SCALE)));
505 node.setAttribute("polar13Scale", valueOf(forceField.getDouble("polar-13-scale", DEFAULT_POLAR_13_SCALE)));
506 node.setAttribute("polar14Scale", valueOf(forceField.getDouble("polar-14-scale", DEFAULT_POLAR_14_SCALE)));
507 node.setAttribute("polar15Scale", valueOf(forceField.getDouble("polar-15-scale", DEFAULT_POLAR_15_SCALE)));
508
509
510 node.setAttribute("polar14Intra", valueOf(forceField.getDouble("polar-14-intra", DEFAULT_POLAR_14_INTRA)));
511 }
512
513
514
515
516 public Element toXML(Document doc) {
517 Element node = doc.createElement("Polarize");
518 node.setAttribute("type", format("%d", type));
519
520 node.setAttribute("polarizability", format("%f", polarizability * ANG_TO_NM * ANG_TO_NM * ANG_TO_NM));
521 node.setAttribute("thole", format("%f", thole));
522
523 int i = 1;
524 if (polarizationGroup != null) {
525 for (int a : polarizationGroup) {
526 node.setAttribute(format("pgrp%d", i), format("%d", a));
527 i++;
528 }
529 }
530 return node;
531 }
532
533
534
535
536
537
538 void incrementType(int increment) {
539 type += increment;
540 setKey(Integer.toString(type));
541 if (polarizationGroup != null) {
542 for (int i = 0; i < polarizationGroup.length; i++) {
543 polarizationGroup[i] += increment;
544 }
545 }
546 }
547
548
549
550
551
552
553
554 boolean patchTypes(HashMap<AtomType, AtomType> typeMap) {
555 if (polarizationGroup == null) {
556 return false;
557 }
558
559
560 int len = polarizationGroup.length;
561 int added = 0;
562 for (AtomType newType : typeMap.keySet()) {
563 for (int i = 1; i < len; i++) {
564 if (polarizationGroup[i] == newType.type) {
565 AtomType knownType = typeMap.get(newType);
566 added++;
567 polarizationGroup = Arrays.copyOf(polarizationGroup, len + added);
568 polarizationGroup[len + added - 1] = knownType.type;
569 }
570 }
571 }
572
573 return added > 0;
574 }
575 }