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.numerics.math.DoubleMath.add;
41  import static ffx.numerics.math.DoubleMath.dot;
42  import static ffx.numerics.math.DoubleMath.normalize;
43  import static ffx.numerics.math.DoubleMath.sub;
44  import static ffx.potential.parameters.ForceField.ELEC_FORM.FIXED_CHARGE;
45  import static ffx.potential.parameters.ForceField.ForceFieldType.MULTIPOLE;
46  import static ffx.utilities.Constants.BOHR;
47  import static ffx.utilities.Constants.BOHR2;
48  import static ffx.utilities.KeywordGroup.PotentialFunctionParameter;
49  import static java.lang.Double.parseDouble;
50  import static java.lang.Integer.parseInt;
51  import static java.lang.String.format;
52  import static java.lang.System.arraycopy;
53  import static java.util.Arrays.fill;
54  import static org.apache.commons.math3.util.FastMath.abs;
55  
56  import ffx.numerics.math.DoubleMath;
57  import ffx.potential.bonded.Atom;
58  import ffx.potential.parameters.ForceField.ELEC_FORM;
59  import ffx.utilities.FFXKeyword;
60  import java.io.BufferedReader;
61  import java.util.ArrayList;
62  import java.util.Arrays;
63  import java.util.Comparator;
64  import java.util.HashMap;
65  import java.util.List;
66  import java.util.logging.Level;
67  import java.util.logging.Logger;
68  
69  /**
70   * The MultipoleType class defines a multipole in its local frame.
71   *
72   * @author Michael J. Schnieders
73   * @since 1.0
74   */
75  @FFXKeyword(name = "multipole", clazz = String[].class, keywordGroup = PotentialFunctionParameter,
76      description = "[5 lines with: 3 or 4 integers and 1 real; 3 reals; 1 real; 2 reals; 3 reals] "
77          + "Provides the values for a set of atomic multipole parameters at a single site. "
78          + "A complete keyword entry consists of five consecutive lines, the first line containing the multipole keyword and the four following lines. "
79          + "The first line contains three integers which define the atom type on which the multipoles are centered, and the Z-axis and X-axis defining atom types for this center. "
80          + "The optional fourth integer contains the Y-axis defining atom type, and is only required for locally chiral multipole sites. "
81          + "The real number on the first line gives the monopole (atomic charge) in electrons."
82          + "The second line contains three real numbers which give the X-, Y- and Z-components of the atomic dipole in electron-Ang. "
83          + "The final three lines, consisting of one, two and three real numbers give the upper triangle of the traceless atomic quadrupole tensor in electron-Ang^2.")
84  public final class MultipoleType extends BaseType implements Comparator<String> {
85  
86    /** Constant <code>zeroM</code> */
87    public static final double[] zeroM =
88        new double[] {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
89    /** Constant <code>zeroD</code> */
90    public static final double[] zeroD = new double[] {0.0, 0.0, 0.0};
91    /** Constant <code>zeroD</code> */
92    public static final double[][] zeroQ
93        = new double[][] {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
94  
95    /** Constant <code>chrg=t000</code> */
96    public static final int t000 = 0;
97    /** Constant <code>t100=1</code> */
98    public static final int t100 = 1;
99    /** Constant <code>t010=2</code> */
100   public static final int t010 = 2;
101   /** Constant <code>t001=3</code> */
102   public static final int t001 = 3;
103   /** Constant <code>t200=4</code> */
104   public static final int t200 = 4;
105   /** Constant <code>t020=5</code> */
106   public static final int t020 = 5;
107   /** Constant <code>t002=6</code> */
108   public static final int t002 = 6;
109   /** Constant <code>t110=7</code> */
110   public static final int t110 = 7;
111   /** Constant <code>t101=8</code> */
112   public static final int t101 = 8;
113   /** Constant <code>t011=9</code> */
114   public static final int t011 = 9;
115   /** Constant <code>t300=10</code> */
116   public static final int t300 = 10;
117   /** Constant <code>t030=11</code> */
118   public static final int t030 = 11;
119   /** Constant <code>t003=12</code> */
120   public static final int t003 = 12;
121   /** Constant <code>t210=13</code> */
122   public static final int t210 = 13;
123   /** Constant <code>t201=14</code> */
124   public static final int t201 = 14;
125   /** Constant <code>t120=15</code> */
126   public static final int t120 = 15;
127   /** Constant <code>t021=16</code> */
128   public static final int t021 = 16;
129   /** Constant <code>t102=17</code> */
130   public static final int t102 = 17;
131   /** Constant <code>t012=18</code> */
132   public static final int t012 = 18;
133   /** Constant <code>t111=19</code> */
134   public static final int t111 = 19;
135 
136   private static final Logger logger = Logger.getLogger(MultipoleType.class.getName());
137   /** Partial atomic charge (e). */
138   public final double charge;
139   /** Atomic dipole. 1 x 3 (e Angstroms). */
140   public final double[] dipole;
141   /** Atomic quadrupole. 3 x 3 (e Angstroms^2). */
142   public final double[][] quadrupole;
143   /** Local frame definition method. */
144   public final MultipoleFrameDefinition frameDefinition;
145   /** Atom types that define the local frame of this multipole. */
146   public final int[] frameAtomTypes;
147   /**
148    * Charge, dipole, and quadrupole packed into tensor notation: c, dx, dy, dz, qxx, qyy, qzz, qxy,
149    * qxz, qyz
150    */
151   private final double[] multipole;
152 
153   /**
154    * Multipole Constructor. Conversion to electron Angstroms should be requested only when reading
155    * multipole values from the force field file.
156    *
157    * @param multipole an array of {@link double} objects.
158    * @param frameAtomTypes an array of {@link int} objects.
159    * @param frameDefinition a
160    *     {@link ffx.potential.parameters.MultipoleType.MultipoleFrameDefinition} object.
161    * @param convertFromBohr a boolean.
162    */
163   public MultipoleType(
164       double[] multipole,
165       int[] frameAtomTypes,
166       MultipoleFrameDefinition frameDefinition,
167       boolean convertFromBohr) {
168     super(MULTIPOLE, frameAtomTypes);
169     this.multipole = (convertFromBohr) ? bohrToElectronAngstroms(multipole) : multipole;
170     this.frameAtomTypes = frameAtomTypes;
171     this.frameDefinition = frameDefinition;
172     charge = multipole[t000];
173     dipole = unpackDipole(multipole);
174     quadrupole = unpackQuad(multipole);
175     checkMultipole();
176   }
177 
178   /**
179    * Construct a MultipoleType from a reference type and updated coefficients.
180    *
181    * @param multipoleType Reference MultipoleType.
182    * @param multipole The updated multipole parameters.
183    */
184   public MultipoleType(MultipoleType multipoleType, double[] multipole) {
185     this(multipole, Arrays.copyOf(multipoleType.frameAtomTypes, multipoleType.frameAtomTypes.length),
186         multipoleType.frameDefinition, false);
187   }
188 
189   /**
190    * Constructor for MultipoleType.
191    *
192    * @param charge a double.
193    * @param dipole an array of {@link double} objects.
194    * @param quadrupole an array of {@link double} objects.
195    * @param frameAtomTypes an array of {@link int} objects.
196    * @param frameDefinition a
197    *     {@link ffx.potential.parameters.MultipoleType.MultipoleFrameDefinition} object.
198    * @param convertFromBohr a boolean.
199    */
200   public MultipoleType(
201       double charge,
202       double[] dipole,
203       double[][] quadrupole,
204       int[] frameAtomTypes,
205       MultipoleFrameDefinition frameDefinition,
206       boolean convertFromBohr) {
207     super(MULTIPOLE, frameAtomTypes);
208     this.charge = charge;
209     if (convertFromBohr) {
210       this.multipole = bohrToElectronAngstroms(pack(charge, dipole, quadrupole));
211       this.dipole = unpackDipole(multipole);
212       this.quadrupole = unpackQuad(multipole);
213     } else {
214       this.multipole = pack(charge, dipole, quadrupole);
215       this.dipole = dipole;
216       this.quadrupole = quadrupole;
217     }
218     this.frameAtomTypes = frameAtomTypes;
219     this.frameDefinition = frameDefinition;
220     checkMultipole();
221   }
222 
223   /**
224    * assignMultipole.
225    *
226    * @param elecForm a {@link ffx.potential.parameters.ForceField.ELEC_FORM} object.
227    * @param atom a {@link ffx.potential.bonded.Atom} object.
228    * @param forceField a {@link ffx.potential.parameters.ForceField} object.
229    * @param multipole an array of {@link double} objects.
230    * @param i a int.
231    * @param axisAtom an array of {@link int} objects.
232    * @param frame an array of
233    *     {@link ffx.potential.parameters.MultipoleType.MultipoleFrameDefinition} objects.
234    * @return a boolean.
235    */
236   public static boolean assignMultipole(
237       ELEC_FORM elecForm,
238       Atom atom,
239       ForceField forceField,
240       double[] multipole,
241       int i,
242       int[][] axisAtom,
243       MultipoleFrameDefinition[] frame) {
244     MultipoleType type = multipoleTypeFactory(elecForm, atom, forceField);
245     if (type == null) {
246       return false;
247     }
248     arraycopy(type.getMultipole(), 0, multipole, 0, 10);
249     axisAtom[i] = atom.getAxisAtomIndices();
250     frame[i] = atom.getMultipoleType().frameDefinition;
251     return true;
252   }
253 
254   /**
255    * Average two MultipoleType instances. The atom types that define the frame of the new type must
256    * be supplied.
257    *
258    * @param multipoleType1 a {@link ffx.potential.parameters.MultipoleType} object.
259    * @param multipoleType2 a {@link ffx.potential.parameters.MultipoleType} object.
260    * @param multipoleFrameTypes an array of {@link int} objects.
261    * @return a {@link ffx.potential.parameters.MultipoleType} object.
262    */
263   public static MultipoleType averageTypes(
264       MultipoleType multipoleType1, MultipoleType multipoleType2, int[] multipoleFrameTypes) {
265     if (multipoleType1.frameDefinition != multipoleType2.frameDefinition) {
266       return null;
267     }
268     MultipoleType[] types = {multipoleType1, multipoleType2};
269     double[] weights = {0.5, 0.5};
270     double[] averagedMultipole = weightMultipole(types, weights);
271     if (averagedMultipole == null) {
272       return null;
273     }
274     return new MultipoleType(
275         averagedMultipole, multipoleFrameTypes, multipoleType1.frameDefinition, false);
276   }
277 
278   /**
279    * checkMultipoleChirality.
280    *
281    * @param frame a {@link ffx.potential.parameters.MultipoleType.MultipoleFrameDefinition}
282    *     object.
283    * @param frameCoords an array of {@link double} objects.
284    * @param localOrigin an array of {@link double} objects.
285    * @return Whether this multipole underwent chiral inversion.
286    */
287   public static boolean checkMultipoleChirality(
288       MultipoleFrameDefinition frame,
289       double[] localOrigin,
290       double[][] frameCoords) {
291     if (frame != MultipoleFrameDefinition.ZTHENX) {
292       return false;
293     }
294     double[] zAxis = new double[3];
295     double[] xAxis = new double[3];
296     double[] yAxis = new double[3];
297     double[] yMinOrigin = new double[3];
298     zAxis[0] = frameCoords[0][0];
299     zAxis[1] = frameCoords[0][1];
300     zAxis[2] = frameCoords[0][2];
301     xAxis[0] = frameCoords[1][0];
302     xAxis[1] = frameCoords[1][1];
303     xAxis[2] = frameCoords[1][2];
304     yAxis[0] = frameCoords[2][0];
305     yAxis[1] = frameCoords[2][1];
306     yAxis[2] = frameCoords[2][2];
307     sub(localOrigin, yAxis, yMinOrigin);
308     sub(zAxis, yAxis, zAxis);
309     sub(xAxis, yAxis, xAxis);
310     double c1 = zAxis[1] * xAxis[2] - zAxis[2] * xAxis[1];
311     double c2 = xAxis[1] * yMinOrigin[2] - xAxis[2] * yMinOrigin[1];
312     double c3 = yMinOrigin[1] * zAxis[2] - yMinOrigin[2] * zAxis[1];
313     double vol = yMinOrigin[0] * c1 + zAxis[0] * c2 + xAxis[0] * c3;
314     return (vol < 0.0);
315   }
316 
317   /**
318    * Return the rotation matrix for the local to lab frame.
319    *
320    * @param frame the multipole frame definition
321    * @param frameCoords the coordinates of the frame atoms
322    * @param localOrigin the local origin of the frame
323    * @return the rotation matrix
324    */
325   public static double[][] getRotationMatrix(
326       MultipoleFrameDefinition frame, double[] localOrigin, double[][] frameCoords) {
327     double[][] rotmat = new double[3][3];
328     getRotationMatrix(frame, localOrigin, frameCoords, rotmat);
329     return rotmat;
330   }
331 
332   /**
333    * Return the rotation matrix for the local to lab frame.
334    *
335    * @param frame the multipole frame definition
336    * @param frameCoords the coordinates of the frame atoms
337    * @param localOrigin the local origin of the frame
338    * @param rotmat the rotation matrix.
339    */
340   public static void getRotationMatrix(
341       MultipoleFrameDefinition frame, double[] localOrigin, double[][] frameCoords,
342       double[][] rotmat) {
343     double[] zAxis = new double[3];
344     double[] xAxis = new double[3];
345     double[] yAxis = new double[3];
346 
347     // Use the identity matrix as the default rotation matrix
348     rotmat[0][0] = 1.0;
349     rotmat[1][0] = 0.0;
350     rotmat[2][0] = 0.0;
351     rotmat[0][1] = 0.0;
352     rotmat[1][1] = 1.0;
353     rotmat[2][1] = 0.0;
354     rotmat[0][2] = 0.0;
355     rotmat[1][2] = 0.0;
356     rotmat[2][2] = 1.0;
357 
358     switch (frame) {
359       case NONE:
360         return;
361       case BISECTOR:
362         zAxis[0] = frameCoords[0][0];
363         zAxis[1] = frameCoords[0][1];
364         zAxis[2] = frameCoords[0][2];
365         xAxis[0] = frameCoords[1][0];
366         xAxis[1] = frameCoords[1][1];
367         xAxis[2] = frameCoords[1][2];
368         sub(zAxis, localOrigin, zAxis);
369         normalize(zAxis, zAxis);
370         sub(xAxis, localOrigin, xAxis);
371         normalize(xAxis, xAxis);
372         add(xAxis, zAxis, zAxis);
373         normalize(zAxis, zAxis);
374         rotmat[0][2] = zAxis[0];
375         rotmat[1][2] = zAxis[1];
376         rotmat[2][2] = zAxis[2];
377         double dot = dot(xAxis, zAxis);
378         DoubleMath.scale(zAxis, dot, zAxis);
379         sub(xAxis, zAxis, xAxis);
380         normalize(xAxis, xAxis);
381         break;
382       case ZTHENBISECTOR:
383         zAxis[0] = frameCoords[0][0];
384         zAxis[1] = frameCoords[0][1];
385         zAxis[2] = frameCoords[0][2];
386         xAxis[0] = frameCoords[1][0];
387         xAxis[1] = frameCoords[1][1];
388         xAxis[2] = frameCoords[1][2];
389         yAxis[0] = frameCoords[2][0];
390         yAxis[1] = frameCoords[2][1];
391         yAxis[2] = frameCoords[2][2];
392         // Z-Axis
393         sub(zAxis, localOrigin, zAxis);
394         normalize(zAxis, zAxis);
395         rotmat[0][2] = zAxis[0];
396         rotmat[1][2] = zAxis[1];
397         rotmat[2][2] = zAxis[2];
398         sub(xAxis, localOrigin, xAxis);
399         normalize(xAxis, xAxis);
400         sub(yAxis, localOrigin, yAxis);
401         normalize(yAxis, yAxis);
402         // Sum the normalized vectors to the bisector atoms.
403         add(xAxis, yAxis, xAxis);
404         normalize(xAxis, xAxis);
405         dot = dot(xAxis, zAxis);
406         DoubleMath.scale(zAxis, dot, zAxis);
407         sub(xAxis, zAxis, xAxis);
408         normalize(xAxis, xAxis);
409         break;
410       case ZONLY:
411         zAxis[0] = frameCoords[0][0];
412         zAxis[1] = frameCoords[0][1];
413         zAxis[2] = frameCoords[0][2];
414         sub(zAxis, localOrigin, zAxis);
415         normalize(zAxis, zAxis);
416         rotmat[0][2] = zAxis[0];
417         rotmat[1][2] = zAxis[1];
418         rotmat[2][2] = zAxis[2];
419         // X-Axis: initially assume its along the global X-axis.
420         xAxis[0] = 1.0;
421         xAxis[1] = 0.0;
422         xAxis[2] = 0.0;
423         // If the Z-axis is close to the global X-axis,
424         // assume a X-axis along the global Y-axis.
425         dot = rotmat[0][2];
426         if (abs(dot) > 0.866) {
427           xAxis[0] = 0.0;
428           xAxis[1] = 1.0;
429           dot = rotmat[1][2];
430         }
431         DoubleMath.scale(zAxis, dot, zAxis);
432         sub(xAxis, zAxis, xAxis);
433         normalize(xAxis, xAxis);
434         break;
435       // 3-Fold frame rotation matrix elements for Z- and X-axes.
436       case THREEFOLD:
437         zAxis[0] = frameCoords[0][0];
438         zAxis[1] = frameCoords[0][1];
439         zAxis[2] = frameCoords[0][2];
440         sub(zAxis, localOrigin, zAxis);
441         normalize(zAxis, zAxis);
442         xAxis[0] = frameCoords[1][0];
443         xAxis[1] = frameCoords[1][1];
444         xAxis[2] = frameCoords[1][2];
445         sub(xAxis, localOrigin, xAxis);
446         normalize(xAxis, xAxis);
447         yAxis[0] = frameCoords[2][0];
448         yAxis[1] = frameCoords[2][1];
449         yAxis[2] = frameCoords[2][2];
450         sub(yAxis, localOrigin, yAxis);
451         normalize(yAxis, yAxis);
452         add(zAxis, xAxis, zAxis);
453         add(zAxis, yAxis, zAxis);
454         normalize(zAxis, zAxis);
455         rotmat[0][2] = zAxis[0];
456         rotmat[1][2] = zAxis[1];
457         rotmat[2][2] = zAxis[2];
458         dot = dot(xAxis, zAxis);
459         DoubleMath.scale(zAxis, dot, zAxis);
460         sub(xAxis, zAxis, xAxis);
461         normalize(xAxis, xAxis);
462         break;
463       case ZTHENX:
464         zAxis[0] = frameCoords[0][0];
465         zAxis[1] = frameCoords[0][1];
466         zAxis[2] = frameCoords[0][2];
467         xAxis[0] = frameCoords[1][0];
468         xAxis[1] = frameCoords[1][1];
469         xAxis[2] = frameCoords[1][2];
470         sub(zAxis, localOrigin, zAxis);
471         normalize(zAxis, zAxis);
472         rotmat[0][2] = zAxis[0];
473         rotmat[1][2] = zAxis[1];
474         rotmat[2][2] = zAxis[2];
475         sub(xAxis, localOrigin, xAxis);
476         dot = dot(xAxis, zAxis);
477         DoubleMath.scale(zAxis, dot, zAxis);
478         sub(xAxis, zAxis, xAxis);
479         normalize(xAxis, xAxis);
480         break;
481       default:
482         throw new IllegalStateException("Unexpected value: " + frame);
483     }
484     // Set the X elements.
485     rotmat[0][0] = xAxis[0];
486     rotmat[1][0] = xAxis[1];
487     rotmat[2][0] = xAxis[2];
488     // Set the Y elements.
489     rotmat[0][1] = rotmat[2][0] * rotmat[1][2] - rotmat[1][0] * rotmat[2][2];
490     rotmat[1][1] = rotmat[0][0] * rotmat[2][2] - rotmat[2][0] * rotmat[0][2];
491     rotmat[2][1] = rotmat[1][0] * rotmat[0][2] - rotmat[0][0] * rotmat[1][2];
492   }
493 
494   /**
495    * multipoleTypeFactory.
496    *
497    * @param elecForm The electrostatics form being used.
498    * @param atom a {@link ffx.potential.bonded.Atom} object.
499    * @param forceField a {@link ffx.potential.parameters.ForceField} object.
500    * @return a {@link ffx.potential.parameters.MultipoleType} object.
501    */
502   public static MultipoleType multipoleTypeFactory(
503       ELEC_FORM elecForm, Atom atom, ForceField forceField) {
504     AtomType atomType = atom.getAtomType();
505     if (atomType == null) {
506       String message = " Multipoles can only be assigned to atoms that have been typed.";
507       logger.severe(message);
508       return null;
509     }
510 
511     if (elecForm != FIXED_CHARGE) {
512       PolarizeType polarizeType = forceField.getPolarizeType(atomType.getKey());
513       if (polarizeType != null) {
514         atom.setPolarizeType(polarizeType);
515       } else {
516         String message = " No polarization type was found for " + atom;
517         logger.info(message);
518         double polarizability = 0.0;
519         double thole = 0.0;
520         double ddp = 0.0;
521         int[] polarizationGroup = null;
522         polarizeType = new PolarizeType(atomType.type, polarizability, thole, ddp,
523             polarizationGroup);
524         forceField.addForceFieldType(polarizeType);
525         atom.setPolarizeType(polarizeType);
526       }
527     }
528 
529     // No reference atoms.
530     String key1 = atomType.getKey();
531     MultipoleType multipoleType = forceField.getMultipoleType(key1);
532     if (multipoleType != null) {
533       atom.setMultipoleType(multipoleType);
534       atom.setAxisAtoms((Atom[]) null);
535       return multipoleType;
536     }
537 
538     // List of sorted of 1-2 iteractions.
539     List<Atom> n12 = atom.get12List();
540 
541     // No bonds.
542     if (n12 == null || n12.size() < 1) {
543       String message = "Multipoles can only be assigned after bonded relationships are defined.\n";
544       logger.severe(message);
545       return null;
546     }
547 
548     // Only 1-2 connected atoms: 1 reference atom.
549     for (Atom atom2 : n12) {
550       String key = key1 + " " + atom2.getAtomType().getKey();
551       multipoleType = forceField.getMultipoleType(key);
552       if (multipoleType != null) {
553         atom.setMultipoleType(multipoleType);
554         atom.setAxisAtoms(atom2);
555         return multipoleType;
556       }
557     }
558 
559     // Only 1-2 connected atoms: 2 reference atoms.
560     for (Atom atom2 : n12) {
561       String key2 = atom2.getAtomType().getKey();
562       for (Atom atom3 : n12) {
563         if (atom2 == atom3) {
564           continue;
565         }
566         String key3 = atom3.getAtomType().getKey();
567         String key = key1 + " " + key2 + " " + key3;
568         multipoleType = forceField.getMultipoleType(key);
569         if (multipoleType != null) {
570           atom.setMultipoleType(multipoleType);
571           atom.setAxisAtoms(atom2, atom3);
572           return multipoleType;
573         }
574       }
575     }
576 
577     // Only 1-2 connected atoms: 3 reference atoms.
578     for (Atom atom2 : n12) {
579       String key2 = atom2.getAtomType().getKey();
580       for (Atom atom3 : n12) {
581         if (atom2 == atom3) {
582           continue;
583         }
584         String key3 = atom3.getAtomType().getKey();
585         for (Atom atom4 : n12) {
586           if (atom2 == atom4 || atom3 == atom4) {
587             continue;
588           }
589           String key4 = atom4.getAtomType().getKey();
590           String key = key1 + " " + key2 + " " + key3 + " " + key4;
591           multipoleType = forceField.getMultipoleType(key);
592           if (multipoleType != null) {
593             atom.setMultipoleType(multipoleType);
594             atom.setAxisAtoms(atom2, atom3, atom4);
595             return multipoleType;
596           }
597         }
598       }
599     }
600 
601     List<Atom> n13 = atom.get13List();
602 
603     // Revert to a reference atom definition that may include a 1-3 site. For example a hydrogen on
604     // water.
605     for (Atom atom2 : n12) {
606       String key2 = atom2.getAtomType().getKey();
607       for (Atom atom3 : n13) {
608         String key3 = atom3.getAtomType().getKey();
609         String key = key1 + " " + key2 + " " + key3;
610         multipoleType = forceField.getMultipoleType(key);
611         if (multipoleType != null) {
612           atom.setMultipoleType(multipoleType);
613           atom.setAxisAtoms(atom2, atom3);
614           return multipoleType;
615         }
616         for (Atom atom4 : n13) {
617           if (atom4 != null && atom4 != atom3) {
618             String key4 = atom4.getAtomType().getKey();
619             key = key1 + " " + key2 + " " + key3 + " " + key4;
620             multipoleType = forceField.getMultipoleType(key);
621             if (multipoleType != null) {
622               atom.setMultipoleType(multipoleType);
623               atom.setAxisAtoms(atom2, atom3, atom4);
624               return multipoleType;
625             }
626           }
627         }
628       }
629     }
630 
631     return null;
632   }
633 
634   /**
635    * Assign local multipole frame defining atoms.
636    *
637    * @param atom Assign local multipole frame defining atoms for this atom.
638    */
639   public static void assignAxisAtoms(Atom atom) {
640     MultipoleType multipoleType = atom.getMultipoleType();
641     String mutipoleKey = multipoleType.getKey();
642     String atomKey = atom.getAtomType().getKey();
643     String[] types = mutipoleKey.split(" +");
644     int numAxisAtoms = types.length - 1;
645 
646     if (numAxisAtoms == 0) {
647       // No reference atoms.
648       atom.setAxisAtoms((Atom[]) null);
649       return;
650     }
651 
652     // List sorted of 1-2 iteractions.
653     List<Atom> n12 = atom.get12List();
654     if (n12 == null || n12.size() < 1) {
655       String message = "Multipoles can only be assigned after bonded relationships are defined.\n";
656       logger.severe(message);
657       return;
658     }
659 
660     // Only 1-2 connected atoms: 1 reference atom.
661     for (Atom atom2 : n12) {
662       String key = atomKey + " " + atom2.getAtomType().getKey();
663       if (key.equalsIgnoreCase(mutipoleKey)) {
664         atom.setAxisAtoms(atom2);
665         return;
666       }
667     }
668 
669     // Only 1-2 connected atoms: 2 reference atoms.
670     for (Atom atom2 : n12) {
671       String key2 = atom2.getAtomType().getKey();
672       for (Atom atom3 : n12) {
673         if (atom2 == atom3) {
674           continue;
675         }
676         String key3 = atom3.getAtomType().getKey();
677         String key = atomKey + " " + key2 + " " + key3;
678         if (key.equalsIgnoreCase(mutipoleKey)) {
679           atom.setAxisAtoms(atom2, atom3);
680           return;
681         }
682       }
683     }
684 
685     // Only 1-2 connected atoms: 3 reference atoms.
686     for (Atom atom2 : n12) {
687       String key2 = atom2.getAtomType().getKey();
688       for (Atom atom3 : n12) {
689         if (atom2 == atom3) {
690           continue;
691         }
692         String key3 = atom3.getAtomType().getKey();
693         for (Atom atom4 : n12) {
694           if (atom2 == atom4 || atom3 == atom4) {
695             continue;
696           }
697           String key4 = atom4.getAtomType().getKey();
698           String key = atomKey + " " + key2 + " " + key3 + " " + key4;
699           if (key.equalsIgnoreCase(mutipoleKey)) {
700             atom.setAxisAtoms(atom2, atom3);
701             return;
702           }
703         }
704       }
705     }
706 
707     List<Atom> n13 = atom.get13List();
708     // Revert to a reference atom definition that may include a 1-3 site.
709     // For example a hydrogen on water.
710     for (Atom atom2 : n12) {
711       String key2 = atom2.getAtomType().getKey();
712       for (Atom atom3 : n13) {
713         String key3 = atom3.getAtomType().getKey();
714         String key = atomKey + " " + key2 + " " + key3;
715         if (key.equalsIgnoreCase(mutipoleKey)) {
716           atom.setAxisAtoms(atom2, atom3);
717           return;
718         }
719         for (Atom atom4 : n13) {
720           if (atom4 != null && atom4 != atom3) {
721             String key4 = atom4.getAtomType().getKey();
722             key = atomKey + " " + key2 + " " + key3 + " " + key4;
723             if (key.equalsIgnoreCase(mutipoleKey)) {
724               atom.setAxisAtoms(atom2, atom3, atom4);
725               return;
726             }
727           }
728         }
729       }
730     }
731 
732     String message = format(" Assignment of axis atoms failed for %s  %s.", atom, multipoleType);
733     logger.severe(message);
734   }
735 
736 
737   /**
738    * Parse a single line multipole.
739    *
740    * @param input Input String.
741    * @param tokens Input tokens.
742    * @param br A BufferedReader instance.
743    * @return a MultipoleType instance.
744    * @since 1.0
745    */
746   public static MultipoleType parse(String input, String[] tokens, BufferedReader br) {
747     if (tokens == null || tokens.length < 3 || tokens.length > 6) {
748       logger.log(Level.WARNING, "Invalid MULTIPOLE type:\n{0}", input);
749       return null;
750     }
751 
752     try {
753       int nTokens = tokens.length;
754       ArrayList<Integer> frameAtoms = new ArrayList<>();
755       // Loop over integer tokens (i.e. except for the initial 'multipole' keyword and final charge
756       // value).
757       for (int i = 1; i < nTokens - 1; i++) {
758         int frameType = parseInt(tokens[i]);
759         // Ignore atom types of '0'.
760         if (frameType != 0) {
761           frameAtoms.add(frameType);
762         }
763       }
764       int nAtomTypes = frameAtoms.size();
765       int[] atomTypes = new int[nAtomTypes];
766       for (int i = 0; i < nAtomTypes; i++) {
767         atomTypes[i] = frameAtoms.get(i);
768       }
769       // Last token is the monopole.
770       double charge = parseDouble(tokens[nTokens - 1]);
771 
772       MultipoleFrameDefinition frameDefinition = null;
773       if (nAtomTypes == 1) {
774         frameDefinition = MultipoleFrameDefinition.NONE;
775       } else if (nAtomTypes == 2) {
776         frameDefinition = MultipoleFrameDefinition.ZONLY;
777       } else if (nAtomTypes == 3) {
778         // ZTHENX or BISECTOR
779         if (atomTypes[1] < 0 || atomTypes[2] < 0) {
780           frameDefinition = MultipoleFrameDefinition.BISECTOR;
781         } else {
782           frameDefinition = MultipoleFrameDefinition.ZTHENX;
783         }
784       } else if (nAtomTypes == 4) {
785         // ZTHENBISECTOR or THREEFOLD
786         if (atomTypes[2] < 0 && atomTypes[3] < 0) {
787           frameDefinition = MultipoleFrameDefinition.ZTHENBISECTOR;
788           if (atomTypes[1] < 0) {
789             frameDefinition = MultipoleFrameDefinition.THREEFOLD;
790           }
791         }
792       }
793 
794       // Warn if the frame definition is ambiguous.
795       if (frameDefinition == null) {
796         logger.log(Level.FINE, "Ambiguous MULTIPOLE type:\n{0}", input);
797         frameDefinition = MultipoleFrameDefinition.ZTHENX;
798       }
799 
800       for (int i = 0; i < nAtomTypes; i++) {
801         atomTypes[i] = abs(atomTypes[i]);
802       }
803 
804       input = br.readLine().split("#")[0];
805       tokens = input.trim().split(" +");
806       if (tokens.length != 3) {
807         logger.log(Level.WARNING, "Invalid MULTIPOLE type:\n{0}", input);
808         return null;
809       }
810       double[] dipole = new double[3];
811       dipole[0] = parseDouble(tokens[0]);
812       dipole[1] = parseDouble(tokens[1]);
813       dipole[2] = parseDouble(tokens[2]);
814       input = br.readLine().split("#")[0];
815       tokens = input.trim().split(" +");
816       if (tokens.length != 1) {
817         logger.log(Level.WARNING, "Invalid MULTIPOLE type:\n{0}", input);
818         return null;
819       }
820       double[][] quadrupole = new double[3][3];
821       quadrupole[0][0] = parseDouble(tokens[0]);
822       input = br.readLine().split("#")[0];
823       tokens = input.trim().split(" +");
824       if (tokens.length != 2) {
825         logger.log(Level.WARNING, "Invalid MULTIPOLE type:\n{0}", input);
826         return null;
827       }
828       quadrupole[1][0] = parseDouble(tokens[0]);
829       quadrupole[1][1] = parseDouble(tokens[1]);
830       input = br.readLine().split("#")[0];
831       tokens = input.trim().split(" +");
832       if (tokens.length != 3) {
833         logger.log(Level.WARNING, "Invalid MULTIPOLE type:\n{0}", input);
834         return null;
835       }
836       quadrupole[2][0] = parseDouble(tokens[0]);
837       quadrupole[2][1] = parseDouble(tokens[1]);
838       quadrupole[2][2] = parseDouble(tokens[2]);
839       // Fill in symmetric components.
840       quadrupole[0][1] = quadrupole[1][0];
841       quadrupole[0][2] = quadrupole[2][0];
842       quadrupole[1][2] = quadrupole[2][1];
843       return new MultipoleType(charge, dipole, quadrupole, atomTypes, frameDefinition, true);
844     } catch (Exception e) {
845       String message = "Exception parsing MULTIPOLE type:\n" + input + "\n";
846       logger.log(Level.SEVERE, message, e);
847     }
848     return null;
849   }
850 
851   /**
852    * Parse a single line multipole.
853    *
854    * @param input Input String.
855    * @param tokens Input tokens.
856    * @return a MultipoleType instance.
857    * @since 1.0
858    */
859   public static MultipoleType parse(String input, String[] tokens) {
860     if (tokens == null || tokens.length < 12 || tokens.length > 15) {
861       logger.log(Level.WARNING, "Invalid MULTIPOLE type:\n{0}", input);
862       return null;
863     }
864 
865     try {
866       int nTokens = tokens.length;
867       ArrayList<Integer> frameAtoms = new ArrayList<>();
868       // Loop over integer tokens (i.e. except for the initial 'multipole' keyword and final 10
869       // multipole values).
870       for (int i = 1; i < nTokens - 10; i++) {
871         int frameType = parseInt(tokens[i]);
872         // Ignore atom types of '0'.
873         if (frameType != 0) {
874           frameAtoms.add(frameType);
875         }
876       }
877       int nAtomTypes = frameAtoms.size();
878       int[] atomTypes = new int[nAtomTypes];
879       for (int i = 0; i < nAtomTypes; i++) {
880         atomTypes[i] = frameAtoms.get(i);
881       }
882 
883       MultipoleFrameDefinition frameDefinition = null;
884       if (nAtomTypes == 1) {
885         frameDefinition = MultipoleFrameDefinition.NONE;
886       } else if (nAtomTypes == 2) {
887         frameDefinition = MultipoleFrameDefinition.ZONLY;
888       } else if (nAtomTypes == 3) {
889         // ZTHENX or BISECTOR
890         if (atomTypes[1] < 0 || atomTypes[2] < 0) {
891           frameDefinition = MultipoleFrameDefinition.BISECTOR;
892         } else {
893           frameDefinition = MultipoleFrameDefinition.ZTHENX;
894         }
895       } else if (nAtomTypes == 4) {
896         // ZTHENBISECTOR or THREEFOLD
897         if (atomTypes[2] < 0 && atomTypes[3] < 0) {
898           frameDefinition = MultipoleFrameDefinition.ZTHENBISECTOR;
899           if (atomTypes[1] < 0) {
900             frameDefinition = MultipoleFrameDefinition.THREEFOLD;
901           }
902         }
903       }
904 
905       // Warn if the frame definition is ambiguous.
906       if (frameDefinition == null) {
907         logger.log(Level.FINE, "Ambiguous MULTIPOLE type:\n{0}", input);
908         frameDefinition = MultipoleFrameDefinition.ZTHENX;
909       }
910 
911       for (int i = 0; i < nAtomTypes; i++) {
912         atomTypes[i] = abs(atomTypes[i]);
913       }
914 
915       double[] dipole = new double[3];
916       double[][] quadrupole = new double[3][3];
917       double charge = parseDouble(tokens[nTokens - 10]);
918       dipole[0] = parseDouble(tokens[nTokens - 9]);
919       dipole[1] = parseDouble(tokens[nTokens - 8]);
920       dipole[2] = parseDouble(tokens[nTokens - 7]);
921       quadrupole[0][0] = parseDouble(tokens[nTokens - 6]);
922       quadrupole[1][0] = parseDouble(tokens[nTokens - 5]);
923       quadrupole[1][1] = parseDouble(tokens[nTokens - 4]);
924       quadrupole[2][0] = parseDouble(tokens[nTokens - 3]);
925       quadrupole[2][1] = parseDouble(tokens[nTokens - 2]);
926       quadrupole[2][2] = parseDouble(tokens[nTokens - 1]);
927       // Fill in symmetric components.
928       quadrupole[0][1] = quadrupole[1][0];
929       quadrupole[0][2] = quadrupole[2][0];
930       quadrupole[1][2] = quadrupole[2][1];
931       return new MultipoleType(charge, dipole, quadrupole, atomTypes, frameDefinition, true);
932     } catch (Exception e) {
933       String message = "Exception parsing MULTIPOLE type:\n" + input + "\n";
934       logger.log(Level.SEVERE, message, e);
935     }
936     return null;
937   }
938 
939   /**
940    * Map charge parameters to a Multipole instance.
941    *
942    * @param input Input string.
943    * @param tokens Input string tokens.
944    * @return a MultipoleType instance
945    */
946   public static MultipoleType parseChargeType(String input, String[] tokens) {
947     if (tokens.length < 3) {
948       logger.log(Level.WARNING, "Invalid CHARGE type:\n{0}", input);
949     } else {
950       try {
951         int[] atomTypes = new int[] {parseInt(tokens[1]), 0, 0};
952         double charge = parseDouble(tokens[2]);
953         double[] dipole = new double[3];
954         double[][] quadrupole = new double[3][3];
955         MultipoleFrameDefinition frameDefinition = MultipoleFrameDefinition.NONE;
956         return new MultipoleType(charge, dipole, quadrupole, atomTypes, frameDefinition, true);
957       } catch (NumberFormatException e) {
958         String message = "Exception parsing CHARGE type:\n" + input + "\n";
959         logger.log(Level.SEVERE, message, e);
960       }
961     }
962     return null;
963   }
964 
965   /**
966    * rotateDipole.
967    *
968    * @param rotmat an array of {@link double} objects.
969    * @param dipole an array of {@link double} objects.
970    * @param rotatedDipole an array of {@link double} objects.
971    */
972   public static void rotateDipole(double[][] rotmat, double[] dipole, double[] rotatedDipole) {
973     for (int i = 0; i < 3; i++) {
974       double[] rotmati = rotmat[i];
975       for (int j = 0; j < 3; j++) {
976         rotatedDipole[i] += rotmati[j] * dipole[j];
977       }
978     }
979   }
980 
981   /**
982    * rotateMultipole.
983    *
984    * @param rotmat an array of {@link double} objects.
985    * @param dipole an array of {@link double} objects.
986    * @param quadrupole an array of {@link double} objects.
987    * @param rotatedDipole an array of {@link double} objects.
988    * @param rotatedQuadrupole an array of {@link double} objects.
989    */
990   public static void rotateMultipole(
991       double[][] rotmat,
992       double[] dipole,
993       double[][] quadrupole,
994       double[] rotatedDipole,
995       double[][] rotatedQuadrupole) {
996 
997     // Initialize the rotated multipole to zero.
998     fill(rotatedDipole, 0.0);
999     for (int i = 0; i < 3; i++) {
1000       fill(rotatedQuadrupole[i], 0.0);
1001     }
1002 
1003     // Compute the rotated multipole.
1004     for (int i = 0; i < 3; i++) {
1005       double[] rotmati = rotmat[i];
1006       double[] quadrupolei = rotatedQuadrupole[i];
1007       rotatedDipole[i] = 0.0;
1008       for (int j = 0; j < 3; j++) {
1009         double[] rotmatj = rotmat[j];
1010         rotatedDipole[i] += rotmati[j] * dipole[j];
1011         if (j < i) {
1012           quadrupolei[j] = rotatedQuadrupole[j][i];
1013         } else {
1014           for (int k = 0; k < 3; k++) {
1015             double[] localQuadrupolek = quadrupole[k];
1016             quadrupolei[j] +=
1017                 rotmati[k]
1018                     * (rotmatj[0] * localQuadrupolek[0]
1019                     + rotmatj[1] * localQuadrupolek[1]
1020                     + rotmatj[2] * localQuadrupolek[2]);
1021           }
1022         }
1023       }
1024     }
1025   }
1026 
1027   /**
1028    * scale.
1029    *
1030    * @param type a {@link ffx.potential.parameters.MultipoleType} object.
1031    * @param cdtScales an array of {@link double} objects.
1032    * @return an array of {@link double} objects.
1033    */
1034   public static double[] scale(MultipoleType type, double[] cdtScales) {
1035     return scale(type.getMultipole(), cdtScales);
1036   }
1037 
1038   /**
1039    * scale.
1040    *
1041    * @param multipole an array of {@link double} objects.
1042    * @param cdtScales an array of {@link double} objects.
1043    * @return an array of {@link double} objects.
1044    */
1045   public static double[] scale(double[] multipole, double[] cdtScales) {
1046     double chargeScale = cdtScales[0];
1047     double dipoleScale = cdtScales[1];
1048     double quadScale = cdtScales[2];
1049     return new double[] {
1050         multipole[t000] * chargeScale,
1051         multipole[t100] * dipoleScale,
1052         multipole[t010] * dipoleScale,
1053         multipole[t001] * dipoleScale,
1054         multipole[t200] * quadScale,
1055         multipole[t020] * quadScale,
1056         multipole[t002] * quadScale,
1057         multipole[t110] * quadScale,
1058         multipole[t101] * quadScale,
1059         multipole[t011] * quadScale
1060     };
1061   }
1062 
1063   /* Indices into a 1D tensor array based on compressed tensor notation. This makes multipole code easier to read. */
1064 
1065   /**
1066    * weightMultipoleTypes.
1067    *
1068    * @param types an array of {@link ffx.potential.parameters.MultipoleType} objects.
1069    * @param weights an array of {@link double} objects.
1070    * @param frameAtomTypes an array of {@link int} objects.
1071    * @return a {@link ffx.potential.parameters.MultipoleType} object.
1072    */
1073   public static MultipoleType weightMultipoleTypes(
1074       MultipoleType[] types, double[] weights, int[] frameAtomTypes) {
1075     double[] weightedMultipole = weightMultipole(types, weights);
1076     if (weightedMultipole == null) {
1077       return null;
1078     }
1079     return new MultipoleType(weightedMultipole, frameAtomTypes, types[0].frameDefinition, false);
1080   }
1081 
1082   private static double[] bohrToElectronAngstroms(double[] multipole) {
1083     return new double[] {
1084         multipole[t000],
1085         multipole[t100] *= BOHR,
1086         multipole[t010] *= BOHR,
1087         multipole[t001] *= BOHR,
1088         multipole[t200] *= BOHR2,
1089         multipole[t020] *= BOHR2,
1090         multipole[t002] *= BOHR2,
1091         multipole[t110] *= BOHR2,
1092         multipole[t101] *= BOHR2,
1093         multipole[t011] *= BOHR2
1094     };
1095   }
1096 
1097   /**
1098    * Pack charge, dipole, quad into 1d tensor array form.
1099    *
1100    * @param charge a double.
1101    * @param dipole an array of {@link double} objects.
1102    * @param quad an array of {@link double} objects.
1103    * @return an array of {@link double} objects.
1104    */
1105   private static double[] pack(double charge, double[] dipole, double[][] quad) {
1106     return new double[] {
1107         charge,
1108         dipole[0],
1109         dipole[1],
1110         dipole[2],
1111         quad[0][0],
1112         quad[1][1],
1113         quad[2][2],
1114         quad[0][1],
1115         quad[0][2],
1116         quad[1][2]
1117     };
1118   }
1119 
1120   /** Unpack dipole from 1d tensor-form multipole. */
1121   private static double[] unpackDipole(double[] mpole) {
1122     return new double[] {mpole[t100], mpole[t010], mpole[t001]};
1123   }
1124 
1125   /** Unpack quadrupole from 1d tensor-form multipole. */
1126   private static double[][] unpackQuad(double[] mpole) {
1127     return new double[][] {
1128         {mpole[t200], mpole[t110], mpole[t101]},
1129         {mpole[t110], mpole[t020], mpole[t011]},
1130         {mpole[t101], mpole[t011], mpole[t002]}
1131     };
1132   }
1133 
1134   /**
1135    * Create a new multipole representing a weighted average.
1136    *
1137    * @param types an array of {@link ffx.potential.parameters.MultipoleType} objects.
1138    * @param weights an array of {@link double} objects.
1139    * @return an array of {@link double} objects.
1140    */
1141   private static double[] weightMultipole(MultipoleType[] types, double[] weights) {
1142     if (types == null || weights == null || types.length != weights.length) {
1143       throw new IllegalArgumentException();
1144     }
1145     if (Arrays.asList(types).contains(null)) {
1146       // Multipoles have not yet been assigned.
1147       return null;
1148     }
1149     for (MultipoleType type : types) {
1150       if (type.frameDefinition != types[0].frameDefinition) {
1151         logger.warning(
1152             format(
1153                 "Multipole frame definition mismatch during weighting:\n\t%s->%s,\n\t%s->%s",
1154                 types[0],
1155                 types[0].frameDefinition.toString(),
1156                 type,
1157                 type.frameDefinition.toString()));
1158         throw new IllegalArgumentException();
1159       }
1160     }
1161     double[] weightedMultipole = new double[10];
1162     fill(weightedMultipole, 0.0);
1163     for (int idx = 0; idx < types.length; idx++) {
1164       double[] multipole = types[idx].getMultipole();
1165       for (int comp = 0; comp < 10; comp++) {
1166         weightedMultipole[comp] += weights[idx] * multipole[comp];
1167       }
1168     }
1169     return weightedMultipole;
1170   }
1171 
1172   /** {@inheritDoc} */
1173   @Override
1174   public int compare(String s1, String s2) {
1175     String[] keys1 = s1.split(" ");
1176     String[] keys2 = s2.split(" ");
1177 
1178     int len = keys1.length;
1179     if (keys1.length > keys2.length) {
1180       len = keys2.length;
1181     }
1182     int[] c1 = new int[len];
1183     int[] c2 = new int[len];
1184     for (int i = 0; i < len; i++) {
1185       c1[i] = abs(parseInt(keys1[i]));
1186       c2[i] = abs(parseInt(keys2[i]));
1187       if (c1[i] < c2[i]) {
1188         return -1;
1189       } else if (c1[i] > c2[i]) {
1190         return 1;
1191       }
1192     }
1193 
1194     if (keys1.length < keys2.length) {
1195       return -1;
1196     } else if (keys1.length > keys2.length) {
1197       return 1;
1198     }
1199 
1200     return 0;
1201   }
1202 
1203   /** {@inheritDoc} */
1204   @Override
1205   public boolean equals(Object o) {
1206     if (this == o) {
1207       return true;
1208     }
1209     if (o == null || getClass() != o.getClass()) {
1210       return false;
1211     }
1212     MultipoleType multipoleType = (MultipoleType) o;
1213     return Arrays.equals(frameAtomTypes, multipoleType.frameAtomTypes);
1214   }
1215 
1216   /**
1217    * Getter for the field <code>charge</code>.
1218    *
1219    * @return An uneditable copy of this type's charge. To make changes, use getMultipoleReference().
1220    */
1221   public double getCharge() {
1222     return multipole[t000];
1223   }
1224 
1225   /**
1226    * Getter for the field <code>dipole</code>.
1227    *
1228    * @return An uneditable copy of this type's dipole. To make changes, use getMultipoleReference().
1229    */
1230   public double[] getDipole() {
1231     return new double[] {multipole[t100], multipole[t010], multipole[t001]};
1232   }
1233 
1234   /**
1235    * Getter for the field <code>multipole</code>.
1236    *
1237    * @return An uneditable copy of this type's multipole. To make changes, use
1238    *     getMultipoleReference().
1239    */
1240   public double[] getMultipole() {
1241     return new double[] {
1242         multipole[t000],
1243         multipole[t100],
1244         multipole[t010],
1245         multipole[t001],
1246         multipole[t200],
1247         multipole[t020],
1248         multipole[t002],
1249         multipole[t110],
1250         multipole[t101],
1251         multipole[t011]
1252     };
1253   }
1254 
1255   /**
1256    * Getter for the field <code>quadrupole</code>.
1257    *
1258    * @return An uneditable copy of this type's quadrupole. To make changes, use
1259    *     getMultipoleReference().
1260    */
1261   public double[][] getQuadrupole() {
1262     return new double[][] {
1263         {multipole[t200], multipole[t110], multipole[t101]},
1264         {multipole[t110], multipole[t020], multipole[t011]},
1265         {multipole[t101], multipole[t011], multipole[t002]}
1266     };
1267   }
1268 
1269   /** {@inheritDoc} */
1270   @Override
1271   public int hashCode() {
1272     return Arrays.hashCode(frameAtomTypes);
1273   }
1274 
1275   /**
1276    * Nicely formatted multipole string. Dipole and qaudrupole are in electron-Bohrs and
1277    * electron-Bohrs^2, respectively.
1278    *
1279    * @return String
1280    */
1281   public String toCompactBohrString() {
1282     StringBuilder multipoleBuffer = new StringBuilder("mpol ");
1283     switch (frameDefinition) {
1284       case NONE:
1285         multipoleBuffer.append(format("(%3d):", frameAtomTypes[0]));
1286         break;
1287       case ZONLY:
1288         multipoleBuffer.append(format("(%3d,%3d):", frameAtomTypes[0], frameAtomTypes[1]));
1289         break;
1290       case ZTHENX:
1291         multipoleBuffer.append(
1292             format("(%3d,%3d,%3d):", frameAtomTypes[0], frameAtomTypes[1], frameAtomTypes[2]));
1293         break;
1294       case BISECTOR:
1295         multipoleBuffer.append(
1296             format("(%3d,%3d,%3d):", frameAtomTypes[0], -frameAtomTypes[1], -frameAtomTypes[2]));
1297         break;
1298       case ZTHENBISECTOR:
1299         multipoleBuffer.append(
1300             format(
1301                 "(%3d,%3d,%3d,%3d):",
1302                 frameAtomTypes[0], frameAtomTypes[1], -frameAtomTypes[2], -frameAtomTypes[3]));
1303         break;
1304       case THREEFOLD:
1305         multipoleBuffer.append(
1306             format(
1307                 "(%3d,%3d,%3d,%3d):",
1308                 frameAtomTypes[0], -frameAtomTypes[1], -frameAtomTypes[2], -frameAtomTypes[3]));
1309     }
1310     multipoleBuffer.append(
1311         format(
1312             "[%6.3f / %6.3f %6.3f %6.3f / %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f]",
1313             charge,
1314             dipole[0] / BOHR,
1315             dipole[1] / BOHR,
1316             dipole[2] / BOHR,
1317             quadrupole[0][0] / BOHR2,
1318             quadrupole[1][0] / BOHR2,
1319             quadrupole[1][1] / BOHR2,
1320             quadrupole[2][0] / BOHR2,
1321             quadrupole[2][1] / BOHR2,
1322             quadrupole[2][2] / BOHR2));
1323     return multipoleBuffer.toString();
1324   }
1325 
1326   /** {@inheritDoc} */
1327   @Override
1328   public String toString() {
1329     return toBohrString();
1330   }
1331 
1332   /**
1333    * incrementType
1334    *
1335    * @param increment a int.
1336    */
1337   void incrementType(int increment) {
1338     for (int i = 0; i < frameAtomTypes.length; i++) {
1339       // Frame atom types of 0 are unchanged.
1340       if (frameAtomTypes[i] > 0) {
1341         frameAtomTypes[i] += increment;
1342       } else if (frameAtomTypes[i] < 0) {
1343         frameAtomTypes[i] -= increment;
1344       }
1345     }
1346     setKey(frameAtomTypes);
1347   }
1348 
1349   /**
1350    * Remap new atom types to known internal ones.
1351    *
1352    * @param typeMap a lookup between new atom types and known atom types.
1353    * @return a {@link ffx.potential.parameters.MultipoleType} object.
1354    */
1355   MultipoleType patchTypes(HashMap<AtomType, AtomType> typeMap) {
1356     int count = 0;
1357     int len = frameAtomTypes.length;
1358 
1359     // Look for a MultipoleType that contain a mapped atom class.
1360     for (AtomType newType : typeMap.keySet()) {
1361       for (int frameAtomType : frameAtomTypes) {
1362         if (frameAtomType == newType.type || frameAtomType == 0) {
1363           count++;
1364         }
1365       }
1366     }
1367 
1368     // If found, create a new MultipoleType that bridges to known classes.
1369     if (count > 0 && count < len) {
1370       int[] newFrame = Arrays.copyOf(frameAtomTypes, len);
1371       for (AtomType newType : typeMap.keySet()) {
1372         for (int i = 0; i < len; i++) {
1373           if (frameAtomTypes[i] == newType.type) {
1374             AtomType knownType = typeMap.get(newType);
1375             newFrame[i] = knownType.type;
1376           }
1377         }
1378       }
1379       return new MultipoleType(multipole, newFrame, frameDefinition, false);
1380     }
1381     return null;
1382   }
1383 
1384   private void checkMultipole() {
1385     double[][] quadrupole = unpackQuad(multipole);
1386     // Check symmetry.
1387     if (abs(quadrupole[0][1] - quadrupole[1][0]) > 1.0e-6) {
1388       logger.warning("Multipole component Qxy != Qyx");
1389       logger.info(this.toString());
1390     }
1391     if (abs(quadrupole[0][2] - quadrupole[2][0]) > 1.0e-6) {
1392       logger.warning("Multipole component Qxz != Qzx");
1393       logger.info(this.toString());
1394     }
1395     if (abs(quadrupole[1][2] - quadrupole[2][1]) > 1.0e-6) {
1396       logger.warning("Multipole component Qyz != Qzy");
1397       logger.info(this.toString());
1398     }
1399     // Warn if the multipole is not traceless.
1400     if (abs(quadrupole[0][0] + quadrupole[1][1] + quadrupole[2][2]) > 1.0e-5) {
1401       logger.log(
1402           Level.WARNING,
1403           format(
1404               "Multipole is not traceless: %12.8f",
1405               abs(quadrupole[0][0] + quadrupole[1][1] + quadrupole[2][2])));
1406       logger.info(this.toString());
1407     }
1408   }
1409 
1410   /**
1411    * Nicely formatted multipole string. Dipole and qaudrupole are in electron-Bohrs and
1412    * electron-Bohrs^2, respectively.
1413    *
1414    * @return String
1415    */
1416   private String toBohrString() {
1417     StringBuilder multipoleBuffer = new StringBuilder("multipole");
1418     switch (frameDefinition) {
1419       case NONE:
1420         multipoleBuffer.append(format("  %5d  %5s  %5s  %5s", frameAtomTypes[0], "", "", ""));
1421         break;
1422       case ZONLY:
1423         multipoleBuffer.append(
1424             format("  %5d  %5d  %5s  %5s", frameAtomTypes[0], frameAtomTypes[1], "", ""));
1425         break;
1426       case ZTHENX:
1427         multipoleBuffer.append(
1428             format(
1429                 "  %5d  %5d  %5d  %5s",
1430                 frameAtomTypes[0], frameAtomTypes[1], frameAtomTypes[2], ""));
1431         break;
1432       case BISECTOR:
1433         multipoleBuffer.append(
1434             format(
1435                 "  %5d  %5d  %5d  %5s",
1436                 frameAtomTypes[0], -frameAtomTypes[1], -frameAtomTypes[2], ""));
1437         break;
1438       case ZTHENBISECTOR:
1439         multipoleBuffer.append(
1440             format(
1441                 "  %5d  %5d  %5d  %5d",
1442                 frameAtomTypes[0], frameAtomTypes[1], -frameAtomTypes[2], -frameAtomTypes[3]));
1443         break;
1444       case THREEFOLD:
1445         multipoleBuffer.append(
1446             format(
1447                 "  %5d  %5d  %5d  %5d",
1448                 frameAtomTypes[0], -frameAtomTypes[1], -frameAtomTypes[2], -frameAtomTypes[3]));
1449     }
1450     if (frameAtomTypes.length == 3) {
1451       multipoleBuffer.append("       ");
1452     }
1453     multipoleBuffer.append(
1454         format(
1455             "  % 7.5f \\\n"
1456                 + "%11$s % 7.5f % 7.5f % 7.5f \\\n"
1457                 + "%11$s % 7.5f \\\n"
1458                 + "%11$s % 7.5f % 7.5f \\\n"
1459                 + "%11$s % 7.5f % 7.5f % 7.5f",
1460             multipole[t000],
1461             multipole[t100] / BOHR,
1462             multipole[t010] / BOHR,
1463             multipole[t001] / BOHR,
1464             multipole[t200] / BOHR2,
1465             multipole[t110] / BOHR2,
1466             multipole[t020] / BOHR2,
1467             multipole[t101] / BOHR2,
1468             multipole[t011] / BOHR2,
1469             multipole[t002] / BOHR2,
1470             "                                      "));
1471     return multipoleBuffer.toString();
1472   }
1473 
1474   /** The local multipole frame is defined by the Z-then-X or Bisector convention. */
1475   public enum MultipoleFrameDefinition {
1476     NONE,
1477     ZONLY,
1478     ZTHENX,
1479     BISECTOR,
1480     ZTHENBISECTOR,
1481     THREEFOLD,
1482   }
1483 }