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