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