View Javadoc
1   // ******************************************************************************
2   //
3   // Title:       Force Field X.
4   // Description: Force Field X - Software for Molecular Biophysics.
5   // Copyright:   Copyright (c) Michael J. Schnieders 2001-2025.
6   //
7   // This file is part of Force Field X.
8   //
9   // Force Field X is free software; you can redistribute it and/or modify it
10  // under the terms of the GNU General Public License version 3 as published by
11  // the Free Software Foundation.
12  //
13  // Force Field X is distributed in the hope that it will be useful, but WITHOUT
14  // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
15  // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
16  // details.
17  //
18  // You should have received a copy of the GNU General Public License along with
19  // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
20  // Place, Suite 330, Boston, MA 02111-1307 USA
21  //
22  // Linking this library statically or dynamically with other modules is making a
23  // combined work based on this library. Thus, the terms and conditions of the
24  // GNU General Public License cover the whole combination.
25  //
26  // As a special exception, the copyright holders of this library give you
27  // permission to link this library with independent modules to produce an
28  // executable, regardless of the license terms of these independent modules, and
29  // to copy and distribute the resulting executable under terms of your choice,
30  // provided that you also meet, for each linked independent module, the terms
31  // and conditions of the license of that module. An independent module is a
32  // module which is not derived from or based on this library. If you modify this
33  // library, you may extend this exception to your version of the library, but
34  // you are not obligated to do so. If you do not wish to do so, delete this
35  // exception statement from your version.
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   * The PolarizeType class defines an isotropic atomic polarizability.
68   * <p>
69   * If this modifier is not present, then charge penetration values will be used for polarization damping, as in the HIPPO polarization model.
70   *
71   * @author Michael J. Schnieders
72   * @since 1.0
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    * Thole damping factor.
105    */
106   public final double thole;
107   /**
108    * Value of polarizability scale factor.
109    */
110   public double pdamp;
111   /**
112    * Direct polarization damping.
113    */
114   public final double ddp;
115   /**
116    * Isotropic polarizability in units of Angstroms^3.
117    */
118   public final double polarizability;
119   /**
120    * Atom type number.
121    */
122   public int type;
123   /**
124    * Connected types in the polarization group of each atom (can be null).
125    */
126   public int[] polarizationGroup;
127 
128   /**
129    * PolarizeType Constructor.
130    *
131    * @param atomType          The atom type.
132    * @param polarizability    The polarizability.
133    * @param thole             The Thole damping constant.
134    * @param polarizationGroup The atom types in the polarization group.
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    * Construct a PolarizeType from a reference type and updated polarizability.
153    *
154    * @param polarizeType   The reference PolarizeType.
155    * @param polarizability The updated polarizability.
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     // Keep pdamp fixed during titration.
161     pdamp = polarizeType.pdamp;
162   }
163 
164   /**
165    * Assign polarization groups to atoms based on their connectivity. This method will
166    * assign 1-1, 1-2, and 1-3 polarization groups to the provided atoms.
167    * The first dimension of the 2D arrays should already be allocated to the number of atoms.
168    *
169    * @param atoms an array of {@link ffx.potential.bonded.Atom} objects.
170    * @param ip11  the 1-1 polarization groups as 2D int array.
171    * @param ip12  the 1-2 polarization groups as 2D int array.
172    * @param ip13  the 1-3 polarization groups as 2D int array.
173    */
174   public static void assignPolarizationGroups(Atom[] atoms, int[][] ip11, int[][] ip12,
175                                               int[][] ip13) {
176 
177     // Find directly connected group members for each atom.
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     // Find 1-2 group relationships.
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     // Find 1-3 group relationships.
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    * Average two PolarizeType instances. The atom types to include in the new polarizationGroup must
287    * be supplied.
288    *
289    * @param polarizeType1     The first PolarizeType.
290    * @param polarizeType2     The second PolarizeType.
291    * @param atomType          The atom type to use for the new PolarizeType.
292    * @param polarizationGroup The atom types to include in the new polarizationGroup.
293    * @return The averaged PolarizeType.
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    * Construct a PolarizeType from an input string.
308    *
309    * @param input  The overall input String.
310    * @param tokens The input String tokenized.
311    * @return a PolarizeType instance.
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         // There might be a direct polarization damping value.
323         double ddp = 0.0;
324         int nextToken = 4;
325         try {
326           // Try to parse the next token as an Int;
327           // If there is an exception then it's a direct damping value.
328           parseInt(tokens[nextToken]);
329         } catch (NumberFormatException e) {
330           ddp = parseDouble(tokens[nextToken]);
331           nextToken = 5;
332         } catch (ArrayIndexOutOfBoundsException e) {
333           // No direct damping value.
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    * A recursive method that checks all atoms bonded to the seed atom for inclusion in the
354    * polarization group. The method is called on each newly found group member.
355    *
356    * @param polarizationGroup Atom types that should be included in the group.
357    * @param group             XYZ indices of current group members.
358    * @param seed              The bonds of the seed atom are queried for inclusion in the group.
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    * A recursive method that checks all atoms bonded to the seed atom for inclusion in the
390    * polarization group. The method is called on each newly found group member.
391    *
392    * @param group XYZ indices of current group members.
393    * @param seed  The bonds of the seed atom are queried for inclusion in the group.
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    * Add an atom type to the polarization group.
420    *
421    * @param atomType The atom type to add.
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    * {@inheritDoc}
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    * {@inheritDoc}
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    * {@inheritDoc}
463    */
464   @Override
465   public int hashCode() {
466     return Objects.hash(type);
467   }
468 
469   /**
470    * {@inheritDoc}
471    *
472    * <p>Nicely formatted polarization type.
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    * Add constant attributes to the AmoebaMultipoleForce
491    *
492    * @param node       the AmoebaMultipoleForce element.
493    * @param forceField the ForceField for collecting constants.
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     // polar-12-intra and polar-13-intra are zero by default. Other values are not supported by FFX (or OpenMM).
509     // polar-15-intra is 1.0 by default. Other values are not supported by FFX or OpenMM.
510     node.setAttribute("polar14Intra", valueOf(forceField.getDouble("polar-14-intra", DEFAULT_POLAR_14_INTRA)));
511   }
512 
513   /**
514    * Write PolarizeType to OpenMM XML format.
515    */
516   public Element toXML(Document doc) {
517     Element node = doc.createElement("Polarize");
518     node.setAttribute("type", format("%d", type));
519     // Convert Ang^3 to nm^3
520     node.setAttribute("polarizability", format("%f", polarizability * ANG_TO_NM * ANG_TO_NM * ANG_TO_NM));
521     node.setAttribute("thole", format("%f", thole));
522     // TODO: OMM does not have support for AMOEBA+ and its use of direct polarization damping (ddp)
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    * Both the type and polarization group are incremented by the same amount.
535    *
536    * @param increment The increment to add to the type.
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    * Add mapped known types to the polarization group of a new patch.
550    *
551    * @param typeMap a lookup between new atom types and known atom types.
552    * @return a boolean.
553    */
554   boolean patchTypes(HashMap<AtomType, AtomType> typeMap) {
555     if (polarizationGroup == null) {
556       return false;
557     }
558 
559     // Append known mapped types.
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 }