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  
39  package ffx.potential.extended;
40  
41  import edu.rit.pj.reduction.SharedDouble;
42  import ffx.numerics.Constraint;
43  import ffx.numerics.Potential;
44  import ffx.potential.ForceFieldEnergy;
45  import ffx.potential.MolecularAssembly;
46  import ffx.potential.SystemState;
47  import ffx.potential.bonded.AminoAcidUtils;
48  import ffx.potential.bonded.AminoAcidUtils.AminoAcid3;
49  import ffx.potential.bonded.Atom;
50  import ffx.potential.bonded.Residue;
51  import ffx.potential.constraint.ShakeChargeConstraint;
52  import ffx.potential.parameters.ForceField;
53  import ffx.potential.parameters.PolarizeType;
54  import ffx.potential.parameters.TitrationUtils;
55  import ffx.potential.parsers.ESVFilter;
56  import ffx.utilities.Constants;
57  import ffx.utilities.FileUtils;
58  import org.apache.commons.configuration2.CompositeConfiguration;
59  import org.apache.commons.io.FilenameUtils;
60  
61  import java.io.File;
62  import java.util.*;
63  import java.util.logging.Level;
64  import java.util.logging.Logger;
65  
66  import static ffx.potential.bonded.BondedUtils.hasAttachedAtom;
67  import static ffx.utilities.Constants.kB;
68  import static java.lang.String.format;
69  import static org.apache.commons.math3.util.FastMath.*;
70  
71  /**
72   * ExtendedSystem class.
73   *
74   * @author Andrew Thiel
75   * @since 1.0
76   */
77  public class ExtendedSystem implements Potential {
78  
79      private static final double DISCR_BIAS = 1.0; // kcal/mol
80      private static final double LOG10 = log(10.0);
81      private static final double THETA_FRICTION = 0.5; // psec^-1
82      private static final double THETA_MASS = 5.0; //Atomic Mass Units
83      private static final int dDiscr_dTautIndex = 6;
84      private static final int dDiscr_dTitrIndex = 3;
85      private static final int dModel_dTautIndex = 8;
86      private static final int dModel_dTitrIndex = 5;
87      private static final int dPh_dTautIndex = 7;
88      private static final int dPh_dTitrIndex = 4;
89      private static final int discrBiasIndex = 0;
90      private static final Logger logger = Logger.getLogger(ExtendedSystem.class.getName());
91      private static final int modelBiasIndex = 2;
92      private static final int pHBiasIndex = 1;
93      public final boolean guessTitrState;
94      /**
95       * Array of AminoAcid3 initialized  to match the number of atoms in the system.
96       * Used to know how to apply vdW or electrostatic ESV terms for the atom.
97       */
98      public final AminoAcid3[] residueNames;
99      /**
100      * Array of ints that is initialized to match the number of atoms in the molecular assembly.
101      * 1 indicates that the tautomer lambda direction is normal.
102      * -1 indicates that the tautomer lambda direction is reversed (1-x).
103      * 0 indicates that the atom is not a tautomerizing atom.
104      */
105     public final int[] tautomerDirections;
106     private final double[] ASH1toASH2 = new double[4];
107     /**
108      * Coefficients that define the per residue type Fmod polynomial
109      * quadratic * titrationLambda^2 + linear * titrationLambda
110      */
111     private final double[] ASHFmod = new double[4];
112     private final double ASHrestraintConstant;
113     /**
114      * Descritizer Bias Magnitude. Default is 1 kcal/mol.
115      */
116     private final double ASHtautBiasMag;
117     private final double ASHtitrBiasMag;
118     private final double[] CYSFmod = new double[4];
119     private final double CYSrestraintConstant;
120     /**
121      * Descritizer Bias Magnitude. Default is 1 kcal/mol.
122      */
123     private final double CYStitrBiasMag;
124     private final double[] GLH1toGLH2 = new double[4];
125     private final double[] GLHFmod = new double[4];
126     private final double GLHrestraintConstant;
127     /**
128      * Descritizer Bias Magnitude. Default is 1 kcal/mol.
129      */
130     private final double GLHtautBiasMag;
131     private final double GLHtitrBiasMag;
132     /**
133      * Coefficients that define the tautomer component of the bivariate Histidine Fmod
134      * quadratic * tautomerLambda^2 + linear * tautomerLambda
135      */
136     private final double[] HIDFmod = new double[4];
137     private final double[] HIDtoHIEFmod = new double[4];
138     private final double[] HIEFmod = new double[4];
139     private final double HISrestraintConstant;
140     /**
141      * Descritizer Bias Magnitude. Default is 1 kcal/mol.
142      */
143     private final double HIStautBiasMag;
144     private final double HIStitrBiasMag;
145     private final double[] LYSFmod = new double[4];
146     private final double LYSrestraintConstant;
147     /**
148      * Descritizer Bias Magnitude. Default is 1 kcal/mol.
149      */
150     private final double LYStitrBiasMag;
151     private final List<Constraint> constraints;
152     private final boolean doBias;
153     private final boolean doElectrostatics;
154     private final boolean doPolarization;
155     // Controls for turning of certain terms for testing.
156     private final boolean doVDW;
157     /**
158      * 3D array to store the titration and tautomer population states for each ESV
159      */
160     final private int[][][] esvHistogram;
161     private final SharedDouble[] esvIndElecDerivs;
162     private final SharedDouble[] esvPermElecDerivs;
163     /**
164      * Array of ints that is initialized to match the number of atoms in the molecular assembly.
165      * Elements correspond to residue index in the tautomerizingResidueList. Only set for tautomerizing residues, -1 otherwise.
166      */
167     private final SystemState esvState;
168     /**
169      * Shared double that is initialized to match the number of ESVs in the system.
170      * Once reduced, will equal either dU_Titr/dLambda or dU_Taut/dLambda for specific ESV
171      */
172     private final SharedDouble[] esvVdwDerivs;
173     /**
174      * Extended Atoms
175      */
176     private final Atom[] extendedAtoms;
177     /**
178      * Array of lambda values that matches residues in extendedResidueList
179      */
180     private final double[] extendedLambdas;
181     /**
182      * Extended Molecules
183      */
184     private final int[] extendedMolecules;
185     /**
186      * Concatenated list of titrating residues + tautomerizing residues.
187      */
188     private final List<Residue> extendedResidueList;
189     /**
190      * ForceField Energy Instance. This instance is only used for Potential implementations for grabbing the energy components.
191      */
192     private final ForceFieldEnergy forceFieldEnergy;
193     /**
194      * Array of booleans that is initialized to match the number of atoms in the molecular assembly
195      * noting whether the atom is tautomerizing. Note that any side chain atom that belongs to a tautomerizing residue
196      * will be flagged as tautomerizing for purposes of scaling electrostatic parameters.
197      */
198     private final boolean[] isTautomerizing;
199     /**
200      * Array of booleans that is initialized to match the number of atoms in the molecular assembly
201      * noting whether the atom is titrating. Note that any side chain atom that belongs to a titrating residue
202      * will be flagged as titrating for purposes of scaling electrostatic parameters.
203      */
204     private final boolean[] isTitrating;
205     /**
206      * Array of booleans that is initialized to match the number of atoms in the assembly to note whether an atom is
207      * specifically a heavy atom with changing polarizability (CYS SG, ASP OD1/OD2, GLU OE1/OE2).
208      */
209     private final boolean[] isTitratingHeavy;
210     /**
211      * Array of booleans that is initialized to match the number of atoms in the assembly to note whether an atom is
212      * specifically a titrating hydrogen.
213      */
214     private final boolean[] isTitratingHydrogen;
215     /**
216      * Boolean similar to fixTitrationState/fixTautomerState but is even more restrictive in that set methods
217      * are not allowed to change lambda values from their initialized values.
218      * Mainly used when evaluating archive snapshots at different initialized lambda values.
219      * If not set to true the archive and esv files set the lambdas automatically.
220      */
221     private final boolean lockStates;
222     /**
223      * Number of atoms in the molecular assembly. Since all protons are instantiated at start, this int will not change.
224      */
225     private final int nAtoms;
226     /**
227      * Number of ESVs attached to the molecular assembly. This number is the sum of tautomerizing + titrating ESVs.
228      */
229     private final int nESVs;
230     /**
231      * Number of titrating ESVs attached to the molecular assembly.
232      */
233     private final int nTitr;
234     private final ArrayList<Double> specialInitTautomer;
235     private final ArrayList<Double> specialInitTitration;
236     private final ArrayList<Double> specialResiduePKAs;
237     private final ArrayList<Double> specialResidues;
238     private final int[] tautomerIndexMap;
239     /**
240      * Array of doubles that is initialized to match the number of atoms in the molecular assembly.
241      * Only elements that match a tautomerizing atom will have their lambda updated.
242      */
243     private final double[] tautomerLambdas;
244     /**
245      * List of tautomerizing residues.
246      */
247     private final List<Residue> tautomerizingResidueList;
248     /**
249      * Friction for the ESV system
250      */
251     private final double thetaFriction;
252     /**
253      * The system defined theta mass of the fictional particle. Used to fill theta mass array.
254      */
255     private final double thetaMass;
256     /**
257      * List of titrating residues.
258      */
259     private final List<Residue> titratingResidueList;
260     /**
261      * Array of ints that is initialized to match the number of atoms in the molecular assembly.
262      * Elements correspond to residue index in the titratingResidueList. Only set for titrating residues, -1 otherwise.
263      */
264     private final int[] titrationIndexMap;
265     /**
266      * Array of doubles that is initialized to match the number of atoms in the molecular assembly.
267      * Only elements that match a titrating atom will have their lambda updated.
268      */
269     private final double[] titrationLambdas;
270     /**
271      * Titration Utils instance. This instance is the master copy that will be distributed to electrostatics classes
272      * when an Extended System is attached.
273      */
274     private final TitrationUtils titrationUtils;
275     private final boolean useChargeConstraint;
276     /**
277      * Dynamics restart file.
278      */
279     File restartFile;
280     /**
281      * Boolean to keep the lambdas from updating over the course of dynamics. Useful for running dynamics
282      * with extended system variables at fixed windows (i.e. BAR)
283      */
284     private boolean fixTautomerState;
285     private boolean fixTitrationState;
286     /**
287      * Filter to parse the dynamics restart file.
288      */
289     private ESVFilter esvFilter = null;
290     /**
291      * System PH.
292      */
293     private double constantSystemPh = 7.4;
294     /**
295      * Target system temperature.
296      */
297     private double currentTemperature = Constants.ROOM_TEMPERATURE;
298 
299     /**
300      * Construct extended system with the provided configuration.
301      *
302      * @param mola a {@link MolecularAssembly} object.
303      */
304     public ExtendedSystem(MolecularAssembly mola, double pH, final File esvFile) {
305         extendedAtoms = mola.getAtomArray();
306         extendedMolecules = mola.getMoleculeNumbers();
307         setConstantPh(pH);
308         ForceField forceField = mola.getForceField();
309         forceFieldEnergy = mola.getPotentialEnergy();
310         if (forceFieldEnergy == null) {
311             logger.severe("No potential energy found?");
312         }
313         CompositeConfiguration properties = mola.getProperties();
314         titrationUtils = new TitrationUtils(forceField);
315 
316         doVDW = properties.getBoolean("esv.vdW", true);
317         doElectrostatics = properties.getBoolean("esv.elec", true);
318         doBias = properties.getBoolean("esv.bias", true);
319         doPolarization = properties.getBoolean("esv.polarization", true);
320         thetaFriction = properties.getDouble("esv.friction", ExtendedSystem.THETA_FRICTION);
321         thetaMass = properties.getDouble("esv.mass", ExtendedSystem.THETA_MASS);
322 
323         lockStates = properties.getBoolean("lock.esv.states", false); // Prevents setTitrationLambda/setTautomerLambda
324         double initialTitrationLambda = properties.getDouble("lambda.titration.initial", 0.5);
325         double initialTautomerLambda = properties.getDouble("lambda.tautomer.initial", 0.5);
326         guessTitrState = properties.getBoolean("guess.titration.state", false);
327         specialResidues = getPropertyList(properties, "esv.special.residues");
328         specialResiduePKAs = getPropertyList(properties, "esv.special.residues.pka");
329         specialInitTautomer = getPropertyList(properties, "esv.special.residues.tautomer");
330         specialInitTitration = getPropertyList(properties, "esv.special.residues.titration");
331         initSpecialResidues(specialResidues, specialResiduePKAs, specialInitTautomer, specialInitTitration, mola);
332 
333         fixTitrationState = properties.getBoolean("fix.titration.lambda", false);
334         fixTautomerState = properties.getBoolean("fix.tautomer.lambda", false);
335         useChargeConstraint = properties.getBoolean("esv.charge.constraint", false);
336         int totalCharge = properties.getInt("esv.charge.constraint.value", 0);
337 
338 
339         ASHFmod[3] = properties.getDouble("ASH.cubic", TitrationUtils.Titration.ASHtoASP.cubic);
340         ASHFmod[2] = properties.getDouble("ASH.quadratic", TitrationUtils.Titration.ASHtoASP.quadratic);
341         ASHFmod[1] = properties.getDouble("ASH.linear", TitrationUtils.Titration.ASHtoASP.linear);
342         ASHFmod[0] = TitrationUtils.Titration.ASHtoASP.offset;
343         ASH1toASH2[3] = properties.getDouble("ASH1toASH2.cubic", TitrationUtils.Titration.ASH1toASH2.cubic);
344         ASH1toASH2[2] = properties.getDouble("ASH1toASH2.quadratic", TitrationUtils.Titration.ASH1toASH2.quadratic);
345         ASH1toASH2[1] = properties.getDouble("ASH1toASH2.linear", TitrationUtils.Titration.ASH1toASH2.linear);
346         ASH1toASH2[0] = TitrationUtils.Titration.ASH1toASH2.offset;
347         ASHtitrBiasMag = properties.getDouble("ASH.titration.bias.magnitude", DISCR_BIAS);
348         ASHtautBiasMag = properties.getDouble("ASH.tautomer.bias.magnitude", DISCR_BIAS);
349         ASHrestraintConstant = properties.getDouble("ASH.restraint.constant", 0.0);
350 
351         GLHFmod[3] = properties.getDouble("GLH.cubic", TitrationUtils.Titration.GLHtoGLU.cubic);
352         GLHFmod[2] = properties.getDouble("GLH.quadratic", TitrationUtils.Titration.GLHtoGLU.quadratic);
353         GLHFmod[1] = properties.getDouble("GLH.linear", TitrationUtils.Titration.GLHtoGLU.linear);
354         GLHFmod[0] = TitrationUtils.Titration.GLHtoGLU.offset;
355         GLH1toGLH2[3] = properties.getDouble("GLH1toGLH2.cubic", TitrationUtils.Titration.GLH1toGLH2.cubic);
356         GLH1toGLH2[2] = properties.getDouble("GLH1toGLH2.quadratic", TitrationUtils.Titration.GLH1toGLH2.quadratic);
357         GLH1toGLH2[1] = properties.getDouble("GLH1toGLH2.linear", TitrationUtils.Titration.GLH1toGLH2.linear);
358         GLH1toGLH2[0] = TitrationUtils.Titration.GLH1toGLH2.offset;
359         GLHtitrBiasMag = properties.getDouble("GLH.titration.bias.magnitude", DISCR_BIAS);
360         GLHtautBiasMag = properties.getDouble("GLH.tautomer.bias.magnitude", DISCR_BIAS);
361         GLHrestraintConstant = properties.getDouble("GLH.restraint.constant", 0.0);
362 
363         LYSFmod[3] = properties.getDouble("LYS.cubic", TitrationUtils.Titration.LYStoLYD.cubic);
364         LYSFmod[2] = properties.getDouble("LYS.quadratic", TitrationUtils.Titration.LYStoLYD.quadratic);
365         LYSFmod[1] = properties.getDouble("LYS.linear", TitrationUtils.Titration.LYStoLYD.linear);
366         LYSFmod[0] = TitrationUtils.Titration.LYStoLYD.offset;
367         LYStitrBiasMag = properties.getDouble("LYS.titration.bias.magnitude", DISCR_BIAS);
368         LYSrestraintConstant = properties.getDouble("LYS.restraint.constant", 0.0);
369 
370         CYSFmod[3] = properties.getDouble("CYS.cubic", TitrationUtils.Titration.CYStoCYD.cubic);
371         CYSFmod[2] = properties.getDouble("CYS.quadratic", TitrationUtils.Titration.CYStoCYD.quadratic);
372         CYSFmod[1] = properties.getDouble("CYS.linear", TitrationUtils.Titration.CYStoCYD.linear);
373         CYSFmod[0] = TitrationUtils.Titration.CYStoCYD.offset;
374         CYStitrBiasMag = properties.getDouble("CYS.titration.bias.magnitude", DISCR_BIAS);
375         CYSrestraintConstant = properties.getDouble("CYS.restraint.constant", 0.0);
376 
377         HIDFmod[3] = properties.getDouble("HID.cubic", TitrationUtils.Titration.HIStoHID.cubic);
378         HIDFmod[2] = properties.getDouble("HID.quadratic", TitrationUtils.Titration.HIStoHID.quadratic);
379         HIDFmod[1] = properties.getDouble("HID.linear", TitrationUtils.Titration.HIStoHID.linear);
380         HIDFmod[0] = TitrationUtils.Titration.HIStoHID.offset;
381         HIEFmod[3] = properties.getDouble("HIE.cubic", TitrationUtils.Titration.HIStoHIE.cubic);
382         HIEFmod[2] = properties.getDouble("HIE.quadratic", TitrationUtils.Titration.HIStoHIE.quadratic);
383         HIEFmod[1] = properties.getDouble("HIE.linear", TitrationUtils.Titration.HIStoHIE.linear);
384         HIEFmod[0] = TitrationUtils.Titration.HIStoHIE.offset;
385         HIDtoHIEFmod[3] = properties.getDouble("HIDtoHIE.cubic", TitrationUtils.Titration.HIDtoHIE.cubic);
386         HIDtoHIEFmod[2] = properties.getDouble("HIDtoHIE.quadratic", TitrationUtils.Titration.HIDtoHIE.quadratic);
387         HIDtoHIEFmod[1] = properties.getDouble("HIDtoHIE.linear", TitrationUtils.Titration.HIDtoHIE.linear);
388         HIDtoHIEFmod[0] = TitrationUtils.Titration.HIDtoHIE.offset;
389         HIStitrBiasMag = properties.getDouble("HIS.titration.bias.magnitude", DISCR_BIAS);
390         HIStautBiasMag = properties.getDouble("HIS.tautomer.bias.magnitude", DISCR_BIAS);
391         HISrestraintConstant = properties.getDouble("HIS.restraint.constant", 0.0);
392 
393         // Log all of the titration bias magnitudes for each titratable residue.
394         logger.info("\n Titration bias magnitudes:");
395         logger.info(" Lysine titration bias magnitude: " + LYStitrBiasMag);
396         logger.info(" Cysteine titration bias magnitude: " + CYStitrBiasMag);
397         logger.info(" Histidine titration bias magnitude: " + HIStitrBiasMag);
398         logger.info(" Glutamic acid titration bias magnitude: " + GLHtitrBiasMag);
399         logger.info(" Aspartic acid titration bias magnitude: " + ASHtitrBiasMag);
400 
401         // Log all of the tautomer bias magnitudes for each tautomerizable residue.
402         logger.info("\n Tautomer bias magnitudes:");
403         logger.info(" Aspartic acid tautomer bias magnitude: " + ASHtautBiasMag);
404         logger.info(" Glutamic acid tautomer bias magnitude: " + GLHtautBiasMag);
405         logger.info(" Histidine tautomer bias magnitude: " + HIStautBiasMag);
406 
407         // Log all of  th cubic, quadratic, and linear terms for all of the amino acid terms
408         logger.info("\n Titration bias terms:");
409         logger.info(" Lysine cubic term: " + LYSFmod[3]);
410         logger.info(" Lysine quadratic term: " + LYSFmod[2]);
411         logger.info(" Lysine linear term: " + LYSFmod[1]);
412         logger.info(" Cysteine cubic term: " + CYSFmod[3]);
413         logger.info(" Cysteine quadratic term: " + CYSFmod[2]);
414         logger.info(" Cysteine linear term: " + CYSFmod[1]);
415         logger.info(" Histidine cubic term: " + HIDFmod[3]);
416         logger.info(" HID quadratic term: " + HIDFmod[2]);
417         logger.info(" HID linear term: " + HIDFmod[1]);
418         logger.info(" HIE cubic term: " + HIEFmod[3]);
419         logger.info(" HIE quadratic term: " + HIEFmod[2]);
420         logger.info(" HIE linear term: " + HIEFmod[1]);
421         logger.info(" HID-HIE cubic term: " + HIDtoHIEFmod[3]);
422         logger.info(" HID-HIE quadratic term: " + HIDtoHIEFmod[2]);
423         logger.info(" HID-HIE linear term: " + HIDtoHIEFmod[1]);
424         logger.info(" Glutamic acid cubic term: " + GLHFmod[3]);
425         logger.info(" Glutamic acid quadratic term: " + GLHFmod[2]);
426         logger.info(" Glutamic acid linear term: " + GLHFmod[1]);
427         logger.info(" Aspartic acid cubic term: " + ASHFmod[3]);
428         logger.info(" Aspartic acid quadratic term: " + ASHFmod[2]);
429         logger.info(" Aspartic acid linear term: " + ASHFmod[1]);
430 
431         titratingResidueList = new ArrayList<>();
432         tautomerizingResidueList = new ArrayList<>();
433         extendedResidueList = new ArrayList<>();
434         // Initialize atom arrays with the existing assembly.
435         Atom[] atoms = mola.getAtomArray();
436         nAtoms = atoms.length;
437         isTitrating = new boolean[nAtoms];
438         isTitratingHydrogen = new boolean[nAtoms];
439         isTitratingHeavy = new boolean[nAtoms];
440         isTautomerizing = new boolean[nAtoms];
441         titrationLambdas = new double[nAtoms];
442         tautomerLambdas = new double[nAtoms];
443         titrationIndexMap = new int[nAtoms];
444         tautomerIndexMap = new int[nAtoms];
445         tautomerDirections = new int[nAtoms];
446         residueNames = new AminoAcid3[nAtoms];
447 
448         Arrays.fill(isTitrating, false);
449         Arrays.fill(isTitratingHydrogen, false);
450         Arrays.fill(isTitratingHeavy, false);
451         Arrays.fill(isTautomerizing, false);
452         Arrays.fill(titrationLambdas, 1.0);
453         Arrays.fill(tautomerLambdas, 0.0);
454         Arrays.fill(titrationIndexMap, -1);
455         Arrays.fill(tautomerIndexMap, -1);
456         Arrays.fill(tautomerDirections, 0);
457         Arrays.fill(residueNames, AminoAcid3.UNK);
458 
459         // Cycle through each residue to determine if it is titratable or tautomerizing.
460         // If a residue is one of these, add to titration or tautomer lists.
461         // Next, loop through all atoms and check to see if the atom belongs to this residue.
462         // If the atom does belong to this residue, set all corresponding variables in the respective titration or tautomer array (size = numAtoms).
463         // Store the index of the residue in the respective list into a map array (size = numAtoms).
464         List<Residue> residueList = mola.getResidueList();
465         //logger.info(residueList.toString());
466         List<Residue> preprocessList = new ArrayList<>(residueList);
467         for (Residue residue : preprocessList) {
468             List<Atom> atomList = residue.getSideChainAtoms();
469             for (Atom atom : atomList) {
470                 //Detect disulfide sulfurs, so we can exclude these when setting up titrating residues.
471                 if (atom.getAtomicNumber() == 16) {
472                     if (hasAttachedAtom(atom, 16)) {
473                         residueList.remove(residue);
474                     }
475                 }
476             }
477         }
478         //logger.info(residueList.toString());
479         // Use only a list that contains AminoAcid residues so remove Nucleic Acid residues
480         residueList.removeIf(residue -> (residue.getResidueType() == Residue.ResidueType.NA));
481         for (Residue residue : residueList) {
482             if (isTitratable(residue)) {
483                 titratingResidueList.add(residue);
484                 List<Atom> atomList = residue.getSideChainAtoms();
485                 for (Atom atom : atomList) {
486                     int atomIndex = atom.getArrayIndex();
487                     residueNames[atomIndex] = residue.getAminoAcid3();
488                     isTitrating[atomIndex] = true;
489                     titrationLambdas[atomIndex] = initialTitrationState(residue, initialTitrationLambda, guessTitrState);
490                     int titrationIndex = titratingResidueList.indexOf(residue);
491                     titrationIndexMap[atomIndex] = titrationIndex;
492                     isTitratingHydrogen[atomIndex] = TitrationUtils.isTitratingHydrogen(residue.getAminoAcid3(), atom);
493                     isTitratingHeavy[atomIndex] = TitrationUtils.isTitratingHeavy(residue.getAminoAcid3(), atom);
494                     // Average out pdamp values of the atoms with changing polarizability which will then be used as fixed values throughout simulation.
495                     // When testing end state energies don't average pdamp.
496                     // Default pdamp is set from protonated polarizability, so it must be changed when testing deprotonated end state (Titration lambda = 0.0.)
497                     if (isTitratingHeavy(atomIndex)) {
498                         //If polarization is turned off atom.getPolarizeType() will return null
499                         if (atom.getPolarizeType() != null) {
500                             double deprotPolar = titrationUtils.getPolarizability(atom, 0.0, 0.0, atom.getPolarizeType().polarizability);
501                             double protPolar = titrationUtils.getPolarizability(atom, 1.0, 1.0, atom.getPolarizeType().polarizability);
502                             double avgPolar = 0.5 * deprotPolar + 0.5 * protPolar;
503                             PolarizeType esvPolarizingHeavyAtom = new PolarizeType(atom.getType(), avgPolar, atom.getPolarizeType().thole,
504                                     atom.getPolarizeType().ddp, atom.getPolarizeType().polarizationGroup);
505                             atom.setPolarizeType(esvPolarizingHeavyAtom);
506                         }
507                     }
508                 }
509                 // If is a tautomer, it must also be titrating.
510                 if (isTautomer(residue)) {
511                     tautomerizingResidueList.add(residue);
512                     for (Atom atom : atomList) {
513                         int atomIndex = atom.getArrayIndex();
514                         if (isTitratingHydrogen[atomIndex]) {
515                             logger.info("Residue: " + residue + " Atom: " + atom + " atomType: " +
516                                     atom.getAtomType().type + " " + atom.getAtomType().atomClass);
517                         }
518                         isTautomerizing[atomIndex] = true;
519                         double resNum = residue.getResidueNumber();
520                         tautomerLambdas[atomIndex] = specialResidues.contains(resNum) && !specialInitTautomer.isEmpty() &&
521                                 Math.abs(specialInitTitration.get(specialResidues.indexOf(resNum)) + 1) > 1e-4
522                                 ? specialInitTitration.get(specialResidues.indexOf(resNum)) : initialTautomerLambda;
523                         int tautomerIndex = tautomerizingResidueList.indexOf(residue);
524                         tautomerIndexMap[atomIndex] = tautomerIndex;
525                         tautomerDirections[atomIndex] = TitrationUtils.getTitratingHydrogenDirection(residue.getAminoAcid3(), atom);
526                     }
527                 }
528                 // Test the multipole frames during unit testing.
529                 assert (titrationUtils.testResidueTypes(residue));
530             }
531         }
532 
533         //Concatenate titratingResidueList and tautomerizingResidueList
534         extendedResidueList.addAll(titratingResidueList);
535         extendedResidueList.addAll(tautomerizingResidueList);
536         //Arrays that are sent to integrator are based on extendedResidueList size
537         nESVs = extendedResidueList.size();
538         nTitr = titratingResidueList.size();
539         extendedLambdas = new double[nESVs];
540         esvState = new SystemState(nESVs);
541         //thetaPosition = new double[nESVs];
542         //thetaVelocity = new double[nESVs];
543         //thetaAccel = new double[nESVs];
544         //thetaMassArray = new double[nESVs];
545         esvHistogram = new int[nTitr][10][10];
546         esvVdwDerivs = new SharedDouble[nESVs];
547         esvPermElecDerivs = new SharedDouble[nESVs];
548         esvIndElecDerivs = new SharedDouble[nESVs];
549         for (int i = 0; i < nESVs; i++) {
550             esvVdwDerivs[i] = new SharedDouble(0.0);
551             esvPermElecDerivs[i] = new SharedDouble(0.0);
552             esvIndElecDerivs[i] = new SharedDouble(0.0);
553         }
554         if (useChargeConstraint) {
555             constraints = new ArrayList<>();
556             ShakeChargeConstraint chargeConstraint = new ShakeChargeConstraint(nTitr, totalCharge, 0.001);
557             constraints.add(chargeConstraint);
558         } else {
559             constraints = Collections.emptyList();
560         }
561 
562         //Theta masses should always be the same for each ESV
563         Arrays.fill(esvState.getMass(), thetaMass);
564 
565         for (int i = 0; i < nESVs; i++) {
566             if (i < nTitr) {
567                 Residue residue = extendedResidueList.get(i);
568                 double initialTitrLambda = initialTitrationState(residue, initialTitrationLambda, guessTitrState);
569                 initializeThetaArrays(i, initialTitrLambda);
570             } else {
571                 double resNum = extendedResidueList.get(i).getResidueNumber();
572                 initialTautomerLambda = specialResidues.contains(resNum) && !specialInitTautomer.isEmpty() &&
573                         Math.abs(specialInitTitration.get(specialResidues.indexOf(resNum)) + 1) > 1e-4
574                         ? specialInitTitration.get(specialResidues.indexOf(resNum)) : initialTautomerLambda;
575                 initializeThetaArrays(i, initialTautomerLambda);
576             }
577         }
578         if (esvFilter == null) {
579             esvFilter = new ESVFilter(mola.getName());
580         }
581         if (esvFile == null) {
582             String firstFileName = FilenameUtils.removeExtension(mola.getFile().getAbsolutePath());
583             restartFile = new File(firstFileName + ".esv");
584         } else {
585             double[] thetaPosition = esvState.x();
586             double[] thetaVelocity = esvState.v();
587             double[] thetaAccel = esvState.a();
588             if (!esvFilter.readESV(esvFile, thetaPosition, thetaVelocity, thetaAccel, esvHistogram)) {
589                 String message = " Could not load the restart file - dynamics terminated.";
590                 logger.log(Level.WARNING, message);
591                 throw new IllegalStateException(message);
592             } else {
593                 restartFile = esvFile;
594                 updateLambdas();
595             }
596         }
597     }
598 
599     /**
600      * Initialize special residues specified in the key file
601      */
602     private void initSpecialResidues(ArrayList<Double> specialResidues, ArrayList<Double> specialResiduePKAs, ArrayList<Double> specialInitTautomer, ArrayList<Double> specialInitTitration, MolecularAssembly mola) {
603         if ((specialResidues.size() != specialResiduePKAs.size() && !specialResiduePKAs.isEmpty())
604                 || (specialResidues.size() != specialInitTautomer.size() && !specialInitTautomer.isEmpty())
605                 || (specialResidues.size() != specialInitTitration.size() && !specialInitTitration.isEmpty())
606         ) {
607             logger.severe("The number of special residues and their associated values do not match.");
608         } else if (!specialResidues.isEmpty()) {
609             logger.info("\nSpecial residues and their associated values:");
610             for (int i = 0; i < specialResidues.size(); i++) {
611                 int resNum = (int) (double) specialResidues.get(i) - mola.getResidueList().get(0).getResidueNumber(); // Shift pdb index by first residue number
612                 if (!specialResiduePKAs.isEmpty()) {
613                     logger.info("Residue: " + specialResidues.get(i) + "-" +
614                             mola.getResidueList().get(resNum).getName()
615                             + " Pka: " + specialResiduePKAs.get(i));
616                 }
617                 if (!specialInitTautomer.isEmpty()) {
618                     logger.info("Residue: " + specialResidues.get(i) + "-" +
619                             mola.getResidueList().get(resNum).getName()
620                             + " Tautomer: " + specialInitTautomer.get(i));
621                 }
622                 if (!specialInitTitration.isEmpty()) {
623                     logger.info("Residue: " + specialResidues.get(i) + "-" +
624                             mola.getResidueList().get(resNum).getName()
625                             + " Titration: " + specialInitTitration.get(i));
626                 }
627             }
628             logger.info(" ");
629         }
630         logger.info("Special residues: " + specialResidues);
631         logger.info("Special residues pKa: " + specialResiduePKAs);
632         logger.info("Special residues titration: " + specialInitTitration);
633         logger.info("Special residues tautomer: " + specialInitTautomer);
634         for (Residue res : mola.getResidueList()) {
635             if (!isTitratable(res) && specialResidues.contains((double) res.getResidueNumber())) {
636                 logger.severe("Given special residue: " + res + " is not titratable.");
637             }
638         }
639     }
640 
641     /**
642      * Returns the titratibility of the passed residue
643      */
644     public boolean isTitratable(Residue residue) {
645         if (residue.getResidueType() == Residue.ResidueType.NA) {
646             return false;
647         }
648         AminoAcidUtils.AminoAcid3 AA3 = residue.getAminoAcid3();
649         return AA3.isConstantPhTitratable;
650     }
651 
652     // Method that takes in properties, a string.
653     private ArrayList<Double> getPropertyList(CompositeConfiguration properties, String s) {
654         ArrayList<Double> list = new ArrayList<>();
655         String[] split = properties.getString(s, "").trim()
656                 .replace("[", "")
657                 .replace("]", "")
658                 .replace(",", " ")
659                 .split(" ");
660         for (String s1 : split) {
661             if (s1.isEmpty()) {
662                 continue;
663             }
664             list.add(Double.parseDouble(s1));
665         }
666         return list;
667     }
668 
669     /**
670      * During constructor, initialize arrays that will hold theta positions, velocities, and accelerations.
671      * Positions determined from starting lambda.
672      * Velocities randomly set according to Maxwell Boltzmann Distribution based on temperature.
673      * Accelerations determined from initial forces.
674      *
675      * @param index  index of ExtendedResidueList for which to set these values.
676      * @param lambda starting lambda value for each ESV.
677      */
678     private void initializeThetaArrays(int index, double lambda) {
679         extendedLambdas[index] = lambda;
680         double[] thetaPosition = esvState.x();
681         double[] thetaVelocity = esvState.v();
682         double[] thetaAccel = esvState.a();
683         thetaPosition[index] = Math.asin(Math.sqrt(lambda));
684         //Perform unit analysis carefully
685         Random random = new Random();
686         thetaVelocity[index] = random.nextGaussian() * sqrt(kB * 298.15 / thetaMass);
687         double dUdL = getDerivatives()[index];
688         double dUdTheta = dUdL * sin(2 * thetaPosition[index]);
689         thetaAccel[index] = -Constants.KCAL_TO_GRAM_ANG2_PER_PS2 * dUdTheta / thetaMass;
690         //logger.info(format("Index: %d, dU/dL: %6.8f, dU/dTheta: %6.8f Theta Accel: %6.8f, Theta Velocity: %6.8f", index, -Constants.KCAL_TO_GRAM_ANG2_PER_PS2*dUdL, -Constants.KCAL_TO_GRAM_ANG2_PER_PS2*dUdTheta, thetaAccel[index], thetaVelocity[index]));
691     }
692 
693     /**
694      * get array of dU/dL for each titrating residue
695      *
696      * @return esvDeriv a double[]
697      */
698     public double[] getDerivatives() {
699         double[] esvDeriv = new double[nESVs];
700         double[] biasDerivComponents = new double[9];
701 
702         for (Residue residue : titratingResidueList) {
703             int resTitrIndex = titratingResidueList.indexOf(residue);
704             //Bias Terms
705             if (doBias) {
706                 getBiasTerms(residue, biasDerivComponents);
707                 //Sum up titration bias derivs
708                 esvDeriv[resTitrIndex] += biasDerivComponents[dDiscr_dTitrIndex] + biasDerivComponents[dPh_dTitrIndex] + biasDerivComponents[dModel_dTitrIndex];
709                 //Sum up tautomer bias derivs
710                 if (isTautomer(residue)) {
711                     int resTautIndex = tautomerizingResidueList.indexOf(residue) + nTitr;
712                     esvDeriv[resTautIndex] += biasDerivComponents[dDiscr_dTautIndex] + biasDerivComponents[dPh_dTautIndex] + biasDerivComponents[dModel_dTautIndex];
713                 }
714             }
715 
716             if (doVDW) {
717                 esvDeriv[resTitrIndex] += getVdwDeriv(resTitrIndex);
718                 //Sum up tautomer bias derivs
719                 if (isTautomer(residue)) {
720                     int resTautIndex = tautomerizingResidueList.indexOf(residue) + nTitr;
721                     esvDeriv[resTautIndex] += getVdwDeriv(resTautIndex);
722                 }
723             }
724 
725             if (doElectrostatics) {
726                 esvDeriv[resTitrIndex] += getPermElecDeriv(resTitrIndex);
727                 //Sum up tautomer bias derivs
728                 if (isTautomer(residue)) {
729                     int resTautIndex = tautomerizingResidueList.indexOf(residue) + nTitr;
730                     esvDeriv[resTautIndex] += getPermElecDeriv(resTautIndex);
731                 }
732                 if (doPolarization) {
733                     esvDeriv[resTitrIndex] += getIndElecDeriv(resTitrIndex);
734                     if (isTautomer(residue)) {
735                         int resTautIndex = tautomerizingResidueList.indexOf(residue) + nTitr;
736                         esvDeriv[resTautIndex] += getIndElecDeriv(resTautIndex);
737                     }
738                 }
739             }
740         }
741         return esvDeriv;
742     }
743 
744     /**
745      * Collect respective pH, model, and discr bias terms and their derivatives for each titrating residue.
746      */
747     private void getBiasTerms(Residue residue, double[] biasEnergyAndDerivs) {
748         AminoAcidUtils.AminoAcid3 AA3 = residue.getAminoAcid3();
749         double titrationLambda = getTitrationLambda(residue);
750         double titrationLambdaSquared = titrationLambda * titrationLambda;
751         double titrationLambdaCubed = titrationLambdaSquared * titrationLambda;
752         double discrBias;
753         double pHBias;
754         double modelBias;
755         double dDiscr_dTitr;
756         double dDiscr_dTaut;
757         double dPh_dTitr;
758         double dPh_dTaut;
759         double dMod_dTitr;
760         double dMod_dTaut;
761         //If bias terms shouldn't be computed, set AA3 to UNK so that default case executes and all terms are set to zero.
762         if (!doBias) {
763             AA3 = AminoAcid3.UNK;
764         }
765         switch (AA3) {
766             case ASD:
767             case ASH:
768             case ASP:
769                 // Discr Bias & Derivs
770                 double tautomerLambda = getTautomerLambda(residue);
771                 discrBias = -4.0 * ASHtitrBiasMag * (titrationLambda - 0.5) * (titrationLambda - 0.5);
772                 discrBias += -4.0 * ASHtautBiasMag * (tautomerLambda - 0.5) * (tautomerLambda - 0.5);
773                 dDiscr_dTitr = -8.0 * ASHtitrBiasMag * (titrationLambda - 0.5);
774                 dDiscr_dTaut = -8.0 * ASHtautBiasMag * (tautomerLambda - 0.5);
775 
776                 // pH Bias & Derivs
777                 double pKa1 = TitrationUtils.Titration.ASHtoASP.pKa;
778                 double pKa2 = pKa1;
779                 pHBias = LOG10 * Constants.R * currentTemperature * (1.0 - titrationLambda)
780                         * (tautomerLambda * (pKa1 - constantSystemPh) + (1.0 - tautomerLambda) * (pKa2 - constantSystemPh));
781                 dPh_dTitr = LOG10 * Constants.R * currentTemperature * -1.0
782                         * (tautomerLambda * (pKa1 - constantSystemPh) + (1.0 - tautomerLambda) * (pKa2 - constantSystemPh));
783                 dPh_dTaut = LOG10 * Constants.R * currentTemperature * (1.0 - titrationLambda)
784                         * ((pKa1 - constantSystemPh) - (pKa2 - constantSystemPh));
785 
786                 // pH restraint to add to pH Bias
787                 double restraint = ASHrestraintConstant;
788                 double pHSignedSquared = (pKa1 - constantSystemPh) * Math.abs(pKa1 - constantSystemPh);
789                 // x*abs(x) is quadratic-like but sign is preserved
790                 pHBias += restraint * (1.0 - titrationLambda) * pHSignedSquared;
791                 dPh_dTitr += restraint * -1.0 * pHSignedSquared;
792 
793                 // Model Bias & Derivs
794                 double coeffA0 = ASH1toASH2[3];
795                 double coeffA1 = ASH1toASH2[2];
796                 double coeffA2 = ASH1toASH2[1];
797                 double coeffB0 = ASHFmod[3];
798                 double coeffB1 = ASHFmod[2];
799                 double coeffB2 = ASHFmod[1];
800                 double coeffC0 = ASHFmod[3];
801                 double coeffC1 = ASHFmod[2];
802                 double coeffC2 = ASHFmod[1];
803                 double tautomerLambdaSquared = tautomerLambda * tautomerLambda;
804                 double tautomerLambdaCubed = tautomerLambdaSquared * tautomerLambda;
805                 double oneMinusTitrationLambda = (1.0 - titrationLambda);
806                 double oneMinusTautomerLambda = (1.0 - tautomerLambda);
807                 double coeffBSum = coeffB0 + coeffB1 + coeffB2;
808                 double coeffCSum = coeffC0 + coeffC1 + coeffC2;
809                 modelBias = titrationLambda * (coeffA0 * tautomerLambdaCubed + coeffA1 * tautomerLambdaSquared + coeffA2 * tautomerLambda) +
810                         tautomerLambda * (coeffB0 * titrationLambdaCubed + coeffB1 * titrationLambdaSquared + coeffB2 * titrationLambda) +
811                         oneMinusTautomerLambda * (coeffC0 * titrationLambdaCubed + coeffC1 * titrationLambdaSquared + coeffC2 * titrationLambda) +
812                         //Enforce that HIS(titration==1) state is equal energy no matter tautomer value
813                         oneMinusTitrationLambda * (coeffCSum - coeffBSum) * tautomerLambda;
814                 dMod_dTitr = (coeffA0 * tautomerLambdaCubed + coeffA1 * tautomerLambdaSquared + coeffA2 * tautomerLambda) +
815                         tautomerLambda * (3.0 * coeffB0 * titrationLambdaSquared + 2.0 * coeffB1 * titrationLambda + coeffB2) +
816                         oneMinusTautomerLambda * (3.0 * coeffC0 * titrationLambdaSquared + 2.0 * coeffC1 * titrationLambda + coeffC2) +
817                         -tautomerLambda * (coeffCSum - coeffBSum);
818                 dMod_dTaut = titrationLambda * (3.0 * coeffA0 * tautomerLambdaSquared + 2.0 * coeffA1 * tautomerLambda + coeffA2) +
819                         (coeffB0 * titrationLambdaCubed + coeffB1 * titrationLambdaSquared + coeffB2 * titrationLambda) +
820                         -1.0 * (coeffC0 * titrationLambdaCubed + coeffC1 * titrationLambdaSquared + coeffC2 * titrationLambda) +
821                         oneMinusTautomerLambda * (coeffCSum - coeffBSum);
822                 break;
823             case GLD:
824             case GLH:
825             case GLU:
826                 // Discr Bias & Derivs
827                 tautomerLambda = getTautomerLambda(residue);
828                 discrBias = -4.0 * GLHtitrBiasMag * (titrationLambda - 0.5) * (titrationLambda - 0.5);
829                 discrBias += -4.0 * GLHtautBiasMag * (tautomerLambda - 0.5) * (tautomerLambda - 0.5);
830                 dDiscr_dTitr = -8.0 * GLHtitrBiasMag * (titrationLambda - 0.5);
831                 dDiscr_dTaut = -8.0 * GLHtautBiasMag * (tautomerLambda - 0.5);
832 
833                 // pH Bias & Derivs
834                 pKa1 = TitrationUtils.Titration.GLHtoGLU.pKa;
835                 pKa2 = pKa1;
836                 pHBias = LOG10 * Constants.R * currentTemperature * (1.0 - titrationLambda)
837                         * (tautomerLambda * (pKa1 - constantSystemPh) + (1.0 - tautomerLambda) * (pKa2 - constantSystemPh));
838                 dPh_dTitr = LOG10 * Constants.R * currentTemperature * -1.0
839                         * (tautomerLambda * (pKa1 - constantSystemPh) + (1.0 - tautomerLambda) * (pKa2 - constantSystemPh));
840                 dPh_dTaut = LOG10 * Constants.R * currentTemperature * (1.0 - titrationLambda)
841                         * ((pKa1 - constantSystemPh) - (pKa2 - constantSystemPh));
842 
843                 // pH restraint to add to pH Bias
844                 restraint = GLHrestraintConstant;
845                 pHSignedSquared = (pKa1 - constantSystemPh) * Math.abs(pKa1 - constantSystemPh);
846                 // x*abs(x) is quadratic-like but sign is preserved
847                 pHBias += restraint * (1.0 - titrationLambda) * pHSignedSquared;
848                 dPh_dTitr += restraint * -1.0 * pHSignedSquared;
849 
850                 // Model Bias & Derivs
851                 coeffA0 = GLH1toGLH2[3];
852                 coeffA1 = GLH1toGLH2[2];
853                 coeffA2 = GLH1toGLH2[1];
854                 coeffB0 = GLHFmod[3];
855                 coeffB1 = GLHFmod[2];
856                 coeffB2 = GLHFmod[1];
857                 coeffC0 = GLHFmod[3];
858                 coeffC1 = GLHFmod[2];
859                 coeffC2 = GLHFmod[1];
860                 tautomerLambdaSquared = tautomerLambda * tautomerLambda;
861                 tautomerLambdaCubed = tautomerLambdaSquared * tautomerLambda;
862                 oneMinusTitrationLambda = (1.0 - titrationLambda);
863                 oneMinusTautomerLambda = (1.0 - tautomerLambda);
864                 coeffBSum = coeffB0 + coeffB1 + coeffB2;
865                 coeffCSum = coeffC0 + coeffC1 + coeffC2;
866                 modelBias = titrationLambda * (coeffA0 * tautomerLambdaCubed + coeffA1 * tautomerLambdaSquared + coeffA2 * tautomerLambda) +
867                         tautomerLambda * (coeffB0 * titrationLambdaCubed + coeffB1 * titrationLambdaSquared + coeffB2 * titrationLambda) +
868                         oneMinusTautomerLambda * (coeffC0 * titrationLambdaCubed + coeffC1 * titrationLambdaSquared + coeffC2 * titrationLambda) +
869                         //Enforce that HIS(titration==1) state is equal energy no matter tautomer value
870                         oneMinusTitrationLambda * (coeffCSum - coeffBSum) * tautomerLambda;
871                 dMod_dTitr = (coeffA0 * tautomerLambdaCubed + coeffA1 * tautomerLambdaSquared + coeffA2 * tautomerLambda) +
872                         tautomerLambda * (3.0 * coeffB0 * titrationLambdaSquared + 2.0 * coeffB1 * titrationLambda + coeffB2) +
873                         oneMinusTautomerLambda * (3.0 * coeffC0 * titrationLambdaSquared + 2.0 * coeffC1 * titrationLambda + coeffC2) +
874                         -tautomerLambda * (coeffCSum - coeffBSum);
875                 dMod_dTaut = titrationLambda * (3.0 * coeffA0 * tautomerLambdaSquared + 2.0 * coeffA1 * tautomerLambda + coeffA2) +
876                         (coeffB0 * titrationLambdaCubed + coeffB1 * titrationLambdaSquared + coeffB2 * titrationLambda) +
877                         -1.0 * (coeffC0 * titrationLambdaCubed + coeffC1 * titrationLambdaSquared + coeffC2 * titrationLambda) +
878                         oneMinusTautomerLambda * (coeffCSum - coeffBSum);
879                 break;
880             case HIS:
881             case HID:
882             case HIE:
883                 // Discr Bias & Derivs
884                 tautomerLambda = getTautomerLambda(residue);
885 
886                 discrBias = -4.0 * HIStitrBiasMag * (titrationLambda - 0.5) * (titrationLambda - 0.5);
887                 discrBias += -4.0 * HIStautBiasMag * (tautomerLambda - 0.5) * (tautomerLambda - 0.5);
888                 dDiscr_dTitr = -8.0 * HIStitrBiasMag * (titrationLambda - 0.5);
889                 dDiscr_dTaut = -8.0 * HIStautBiasMag * (tautomerLambda - 0.5);
890 
891                 // pH Bias & Derivs
892                 // At tautomerLambda=1 HIE is fully on.
893                 pKa1 = TitrationUtils.Titration.HIStoHIE.pKa;
894                 // At tautomerLambda=0 HID is fully on.
895                 pKa2 = TitrationUtils.Titration.HIStoHID.pKa;
896                 pHBias = LOG10 * Constants.R * currentTemperature * (1.0 - titrationLambda)
897                         * (tautomerLambda * (pKa1 - constantSystemPh) + (1.0 - tautomerLambda) * (pKa2 - constantSystemPh));
898                 dPh_dTitr = LOG10 * Constants.R * currentTemperature * -1.0
899                         * (tautomerLambda * (pKa1 - constantSystemPh) + (1.0 - tautomerLambda) * (pKa2 - constantSystemPh));
900                 dPh_dTaut = LOG10 * Constants.R * currentTemperature * (1.0 - titrationLambda)
901                         * ((pKa1 - constantSystemPh) - (pKa2 - constantSystemPh));
902 
903                 // pH restraint to add to pH Bias
904                 restraint = HISrestraintConstant;
905                 double pH1SignedSquared = (pKa1 - constantSystemPh) * Math.abs(pKa1 - constantSystemPh);
906                 double pH2SignedSquared = (pKa2 - constantSystemPh) * Math.abs(pKa2 - constantSystemPh);
907                 // x*abs(x) is quadratic-like but sign is preserved
908                 pHBias += restraint * (1.0 - titrationLambda)
909                         * (tautomerLambda * pH1SignedSquared + (1.0 - tautomerLambda) * pH2SignedSquared);
910                 dPh_dTitr += restraint * -1.0
911                         * (tautomerLambda * pH1SignedSquared + (1.0 - tautomerLambda) * pH2SignedSquared);
912                 dPh_dTaut += restraint * (1.0 - titrationLambda)
913                         * (pH1SignedSquared - pH2SignedSquared);
914                 // Model Bias & Derivs
915 
916                 coeffA0 = HIDtoHIEFmod[3];
917                 coeffA1 = HIDtoHIEFmod[2];
918                 coeffA2 = HIDtoHIEFmod[1];
919                 coeffB0 = HIEFmod[3];
920                 coeffB1 = HIEFmod[2];
921                 coeffB2 = HIEFmod[1];
922                 coeffC0 = HIDFmod[3];
923                 coeffC1 = HIDFmod[2];
924                 coeffC2 = HIDFmod[1];
925                 tautomerLambdaSquared = tautomerLambda * tautomerLambda;
926                 tautomerLambdaCubed = tautomerLambdaSquared * tautomerLambda;
927                 oneMinusTitrationLambda = (1.0 - titrationLambda);
928                 oneMinusTautomerLambda = (1.0 - tautomerLambda);
929                 coeffBSum = coeffB0 + coeffB1 + coeffB2;
930                 coeffCSum = coeffC0 + coeffC1 + coeffC2;
931                 modelBias = oneMinusTitrationLambda * (coeffA0 * tautomerLambdaCubed + coeffA1 * tautomerLambdaSquared + coeffA2 * tautomerLambda) +
932                         tautomerLambda * (coeffB0 * titrationLambdaCubed + coeffB1 * titrationLambdaSquared + coeffB2 * titrationLambda) +
933                         oneMinusTautomerLambda * (coeffC0 * titrationLambdaCubed + coeffC1 * titrationLambdaSquared + coeffC2 * titrationLambda) +
934                         //Enforce that HIS(titration==1) state is equal energy no matter tautomer value
935                         titrationLambda * (coeffCSum - coeffBSum) * tautomerLambda;
936                 dMod_dTitr = -(coeffA0 * tautomerLambdaCubed + coeffA1 * tautomerLambdaSquared + coeffA2 * tautomerLambda) +
937                         tautomerLambda * (3.0 * coeffB0 * titrationLambdaSquared + 2.0 * coeffB1 * titrationLambda + coeffB2) +
938                         oneMinusTautomerLambda * (3.0 * coeffC0 * titrationLambdaSquared + 2.0 * coeffC1 * titrationLambda + coeffC2) +
939                         tautomerLambda * (coeffCSum - coeffBSum);
940                 dMod_dTaut = oneMinusTitrationLambda * (3.0 * coeffA0 * tautomerLambdaSquared + 2.0 * coeffA1 * tautomerLambda + coeffA2) +
941                         (coeffB0 * titrationLambdaCubed + coeffB1 * titrationLambdaSquared + coeffB2 * titrationLambda) +
942                         -1.0 * (coeffC0 * titrationLambdaCubed + coeffC1 * titrationLambdaSquared + coeffC2 * titrationLambda) +
943                         titrationLambda * (coeffCSum - coeffBSum);
944 
945                 break;
946             case LYS:
947             case LYD:
948                 // Discr Bias & Derivs
949                 discrBias = -4.0 * LYStitrBiasMag * (titrationLambda - 0.5) * (titrationLambda - 0.5);
950                 dDiscr_dTitr = -8.0 * LYStitrBiasMag * (titrationLambda - 0.5);
951                 dDiscr_dTaut = 0.0;
952 
953                 // pH Bias & Derivs
954                 pKa1 = TitrationUtils.Titration.LYStoLYD.pKa;
955                 pHBias = LOG10 * Constants.R * currentTemperature * (1.0 - titrationLambda) * (pKa1 - constantSystemPh);
956                 dPh_dTitr = LOG10 * Constants.R * currentTemperature * -1.0 * (pKa1 - constantSystemPh);
957                 dPh_dTaut = 0.0;
958 
959                 // pH restraint to add to pH Bias
960                 restraint = LYSrestraintConstant;
961                 pHSignedSquared = (pKa1 - constantSystemPh) * Math.abs(pKa1 - constantSystemPh);
962                 // x*abs(x) is quadratic-like but sign is preserved
963                 pHBias += restraint * (1.0 - titrationLambda) * pHSignedSquared;
964                 dPh_dTitr += restraint * -1.0 * pHSignedSquared;
965 
966                 // Model Bias & Derivs
967                 double cubic = LYSFmod[3];
968                 double quadratic = LYSFmod[2];
969                 double linear = LYSFmod[1];
970                 modelBias = cubic * titrationLambdaCubed + quadratic * titrationLambdaSquared + linear * titrationLambda;
971                 dMod_dTitr = 3 * cubic * titrationLambdaSquared + 2 * quadratic * titrationLambda + linear;
972                 dMod_dTaut = 0.0;
973                 break;
974             case CYS:
975             case CYD:
976                 // Discr Bias & Derivs
977                 discrBias = -4.0 * CYStitrBiasMag * (titrationLambda - 0.5) * (titrationLambda - 0.5);
978                 dDiscr_dTitr = -8.0 * CYStitrBiasMag * (titrationLambda - 0.5);
979                 dDiscr_dTaut = 0.0;
980 
981                 // pH Bias & Derivs
982                 pKa1 = TitrationUtils.Titration.CYStoCYD.pKa;
983                 pHBias = LOG10 * Constants.R * currentTemperature * (1.0 - titrationLambda) * (pKa1 - constantSystemPh);
984                 dPh_dTitr = LOG10 * Constants.R * currentTemperature * -1.0 * (pKa1 - constantSystemPh);
985                 dPh_dTaut = 0.0;
986 
987                 // pH restraint to add to pH Bias
988                 restraint = CYSrestraintConstant;
989                 pHSignedSquared = (pKa1 - constantSystemPh) * Math.abs(pKa1 - constantSystemPh);
990                 // x*abs(x) is quadratic-like but sign is preserved
991                 pHBias += restraint * (1.0 - titrationLambda) * pHSignedSquared;
992                 dPh_dTitr += restraint * -1.0 * pHSignedSquared;
993 
994 
995                 // Model Bias & Derivs
996                 cubic = CYSFmod[3];
997                 quadratic = CYSFmod[2];
998                 linear = CYSFmod[1];
999                 modelBias = cubic * titrationLambdaCubed + quadratic * titrationLambdaSquared + linear * titrationLambda;
1000                 dMod_dTitr = 3 * cubic * titrationLambdaSquared + 2 * quadratic * titrationLambda + linear;
1001                 dMod_dTaut = 0.0;
1002                 break;
1003             default:
1004                 discrBias = 0.0;
1005                 pHBias = 0.0;
1006                 modelBias = 0.0;
1007                 dDiscr_dTitr = 0.0;
1008                 dDiscr_dTaut = 0.0;
1009                 dPh_dTitr = 0.0;
1010                 dPh_dTaut = 0.0;
1011                 dMod_dTitr = 0.0;
1012                 dMod_dTaut = 0.0;
1013                 break;
1014         }
1015         biasEnergyAndDerivs[0] = discrBias;
1016         biasEnergyAndDerivs[1] = pHBias;
1017         biasEnergyAndDerivs[2] = -modelBias;
1018         biasEnergyAndDerivs[3] = dDiscr_dTitr;
1019         biasEnergyAndDerivs[4] = dPh_dTitr;
1020         biasEnergyAndDerivs[5] = -dMod_dTitr;
1021         biasEnergyAndDerivs[6] = dDiscr_dTaut;
1022         biasEnergyAndDerivs[7] = dPh_dTaut;
1023         biasEnergyAndDerivs[8] = -dMod_dTaut;
1024     }
1025 
1026     /**
1027      * Gets the titration lambda for the input residue if the residue is titrating
1028      *
1029      * @param residue a titrating residue
1030      * @return the titration lambda for the residue
1031      */
1032     public double getTitrationLambda(Residue residue) {
1033         if (titratingResidueList.contains(residue)) {
1034             int resIndex = titratingResidueList.indexOf(residue);
1035             return extendedLambdas[resIndex];
1036         } else {
1037             return 1.0;
1038         }
1039     }
1040 
1041     /**
1042      * Gets the tautomer lambda for the input residue if the residue is tautomerizing
1043      *
1044      * @param residue a tautomer residue
1045      * @return the tautomer lambda for the residue
1046      */
1047     public double getTautomerLambda(Residue residue) {
1048         if (tautomerizingResidueList.contains(residue)) {
1049             int resIndex = tautomerizingResidueList.indexOf(residue);
1050             return extendedLambdas[nTitr + resIndex];
1051         } else {
1052             return 1.0;
1053         }
1054     }
1055 
1056     /**
1057      * get total dUvdw/dL for the selected extended system variable
1058      */
1059     private double getVdwDeriv(int esvID) {
1060         return esvVdwDerivs[esvID].get();
1061     }
1062 
1063     /**
1064      * get total dUpermElec/dL for the selected extended system variable.
1065      */
1066     private double getPermElecDeriv(int esvID) {
1067         return esvPermElecDerivs[esvID].get();
1068     }
1069 
1070     /**
1071      * get total dUindElec/dL for the selected extended system variable
1072      */
1073     private double getIndElecDeriv(int esvID) {
1074         return esvIndElecDerivs[esvID].get();
1075     }
1076 
1077     /**
1078      * Returns the tautomerizibility of a residue
1079      */
1080     public boolean isTautomer(Residue residue) {
1081         if (residue.getResidueType() == Residue.ResidueType.NA) {
1082             return false;
1083         }
1084         AminoAcidUtils.AminoAcid3 AA3 = residue.getAminoAcid3();
1085         return AA3.isConstantPhTautomer;
1086     }
1087 
1088     /**
1089      * Update all theta (lambda) positions after each move from the Stochastic integrator
1090      */
1091     private void updateLambdas() {
1092         //If lockStates is true, then the titration and tautomer states are permanently locked.
1093         if (lockStates) {
1094             return;
1095         }
1096         double[] thetaPosition = esvState.x();
1097         //This will prevent recalculating multiple sinTheta*sinTheta that are the same number.
1098         for (int i = 0; i < nESVs; i++) {
1099             //Check to see if titration/tautomer lambdas are to be fixed
1100             if ((!fixTitrationState && i < nTitr) || (!fixTautomerState && i >= nTitr)) {
1101                 double sinTheta = Math.sin(thetaPosition[i]);
1102                 extendedLambdas[i] = sinTheta * sinTheta;
1103             }
1104         }
1105         for (int i = 0; i < nAtoms; i++) {
1106             int mappedTitrationIndex = titrationIndexMap[i];
1107             int mappedTautomerIndex = tautomerIndexMap[i] + nTitr;
1108             if (isTitrating(i) && mappedTitrationIndex != -1) {
1109                 titrationLambdas[i] = extendedLambdas[mappedTitrationIndex];
1110             }
1111             if (isTautomerizing(i) && mappedTautomerIndex >= nTitr) {
1112                 tautomerLambdas[i] = extendedLambdas[mappedTautomerIndex];
1113             }
1114         }
1115         setESVHistogram();
1116     }
1117 
1118     /**
1119      * Returns whether an atom is titrating
1120      */
1121     public boolean isTitrating(int atomIndex) {
1122         return isTitrating[atomIndex];
1123     }
1124 
1125     /**
1126      * Returns whether an atom is tautomerizing
1127      */
1128     public boolean isTautomerizing(int atomIndex) {
1129         return isTautomerizing[atomIndex];
1130     }
1131 
1132     /**
1133      * Goes through residues updates the ESV histogram based on the residues current state
1134      */
1135     private void setESVHistogram() {
1136         for (Residue residue : titratingResidueList) {
1137             int index = titratingResidueList.indexOf(residue);
1138             if (residue.getAminoAcid3().equals(AminoAcid3.LYS)) { // TODO: Add support for CYS? by adding "|| residue.getAminoAcid3().equals(AminoAcid3.CYS))"
1139                 double titrLambda = getTitrationLambda(residue);
1140                 esvHistogram(index, titrLambda);
1141             } else {
1142                 double titrLambda = getTitrationLambda(residue);
1143                 double tautLambda = getTautomerLambda(residue);
1144                 esvHistogram(index, titrLambda, tautLambda);
1145             }
1146         }
1147     }
1148 
1149     /**
1150      * Updates the ESV histogram of the passed residue at the given lambda
1151      *
1152      * @param esv    the index of the esv to be updated
1153      * @param lambda the lambda value to be updated
1154      */
1155     private void esvHistogram(int esv, double lambda) {
1156         int value = (int) (lambda * 10.0);
1157         //Cover the case where lambda could be exactly 1.0
1158         if (value == 10) {
1159             value = 9;
1160         }
1161         esvHistogram[esv][value][0]++;
1162     }
1163 
1164     /**
1165      * Updates the ESV histogram of the passed residue at the given titr and taut state
1166      *
1167      * @param esv        the index of the esv to be updated
1168      * @param titrLambda the titration lambda coordinate to be updated
1169      * @param tautLambda the tautomer lambda coordinate to be updated
1170      */
1171     private void esvHistogram(int esv, double titrLambda, double tautLambda) {
1172         int titrValue = (int) (titrLambda * 10.0);
1173         //Cover the case where lambda could be exactly 1.0
1174         if (titrValue == 10) {
1175             titrValue = 9;
1176         }
1177         int tautValue = (int) (tautLambda * 10.0);
1178         if (tautValue == 10) {
1179             tautValue = 9;
1180         }
1181         esvHistogram[esv][titrValue][tautValue]++;
1182     }
1183 
1184     /**
1185      * Naive guess as to what the best starting state should be based purely on the acidostat term.
1186      */
1187     private double initialTitrationState(Residue residue, double initialLambda, boolean guessTitrState) {
1188         AminoAcid3 AA3 = residue.getAminoAcid3();
1189         double residueNumber = residue.getResidueNumber();
1190         double initialTitrationLambda = 0.0;
1191         if (specialResidues.contains(residueNumber)) { // If we set a special residue for this state
1192             if (!specialResiduePKAs.isEmpty()) {
1193                 initialTitrationLambda =
1194                         (constantSystemPh < specialResiduePKAs.get(specialResidues.indexOf(residueNumber))) ? 1.0 : 0.0;
1195             }
1196             if (!specialInitTitration.isEmpty()) { // Override the pKa value if we have a special initial titration state
1197                 double value = specialInitTitration.get(specialResidues.indexOf(residueNumber));
1198                 if (Math.abs(value + 1) > 1e-4) { // If the value is not -1
1199                     initialTitrationLambda = value;
1200                 }
1201             }
1202         } else if (!guessTitrState) { // If we do not want to guess the titration state but instead use the value
1203             // passed in via the initialLambda parameter.
1204             initialTitrationLambda = initialLambda;
1205         } else { // Guess titration state
1206             initialTitrationLambda = switch (AA3) {
1207                 case ASD -> (constantSystemPh < TitrationUtils.Titration.ASHtoASP.pKa) ? 1.0 : 0.0;
1208                 case GLD -> (constantSystemPh < TitrationUtils.Titration.GLHtoGLU.pKa) ? 1.0 : 0.0;
1209                 case HIS -> (constantSystemPh < TitrationUtils.Titration.HIStoHID.pKa) ? 1.0 : 0.0;
1210                 case LYS -> (constantSystemPh < TitrationUtils.Titration.LYStoLYD.pKa) ? 1.0 : 0.0;
1211                 case CYS -> (constantSystemPh < TitrationUtils.Titration.CYStoCYD.pKa) ? 1.0 : 0.0;
1212                 default -> initialLambda;
1213             };
1214         }
1215         return initialTitrationLambda;
1216     }
1217 
1218     /**
1219      * Questions whether the current non-hydrogen atom's polarizability is changing in response to lambda being updated.
1220      * Only affects carboxylic oxygen and sulfur.
1221      *
1222      * @param atomIndex
1223      * @return
1224      */
1225     public boolean isTitratingHeavy(int atomIndex) {
1226         return isTitratingHeavy[atomIndex];
1227     }
1228 
1229     // Getter for specialResidues
1230     public ArrayList<Double> getSpecialResidueList() {
1231         return specialResidues;
1232     }
1233 
1234     /**
1235      * Does not allow for changes to the tautomer states of tautomerizing residues
1236      */
1237     public void setFixedTautomerState(boolean fixTautomerState) {
1238         this.fixTautomerState = fixTautomerState;
1239     }
1240 
1241     /**
1242      * Does not allow for changes to the tautomer states of titrating residues
1243      */
1244     public void setFixedTitrationState(boolean fixTitrationState) {
1245         this.fixTitrationState = fixTitrationState;
1246     }
1247 
1248     /**
1249      * Reset initialized lambdas to a naive guess based on the model pKa for each extended residue
1250      */
1251     public void reGuessLambdas() {
1252         logger.info(" Reinitializing lambdas to match RepEx window pH");
1253         for (Residue residue : titratingResidueList) {
1254             double lambda = initialTitrationState(residue, 1.0, guessTitrState);
1255             setTitrationLambda(residue, lambda);
1256             double resNum = residue.getResidueNumber();
1257             double tautomerLambda = specialResidues.contains(resNum) && !specialInitTautomer.isEmpty() &&
1258                     Math.abs(specialInitTitration.get(specialResidues.indexOf(resNum)) + 1) > 1e-4
1259                     ? specialInitTitration.get(specialResidues.indexOf(resNum)) : (int) Math.round(random());
1260             setTautomerLambda(residue, tautomerLambda);
1261         }
1262     }
1263 
1264     /**
1265      * Set the tautomer lambda of a residue and update corresponding theta
1266      *
1267      * @param residue
1268      * @param lambda
1269      */
1270     public void setTautomerLambda(Residue residue, double lambda) {
1271         setTautomerLambda(residue, lambda, true);
1272     }
1273 
1274     /**
1275      * Set the tautomer lambda of a residue and update corresponding theta if desired
1276      *
1277      * @param residue      residue to set the lambda of
1278      * @param lambda       value to set the residue to
1279      * @param changeThetas whether to change the theta positions ~comes with information loss~
1280      */
1281     public void setTautomerLambda(Residue residue, double lambda, boolean changeThetas) {
1282         double[] thetaPosition = esvState.x();
1283         if (tautomerizingResidueList.contains(residue) && !lockStates) {
1284             // The correct index in the theta arrays for tautomer coordinates is after the titration list.
1285             // So titrationList.size() + tautomerIndex should match with appropriate spot in thetaPosition, etc.
1286             int index = tautomerizingResidueList.indexOf(residue) + nTitr;
1287             extendedLambdas[index] = lambda;
1288             if (changeThetas) {
1289                 thetaPosition[index] = Math.asin(Math.sqrt(lambda));
1290             }
1291             List<Atom> currentAtomList = residue.getSideChainAtoms();
1292             for (Atom atom : currentAtomList) {
1293                 int atomIndex = atom.getArrayIndex();
1294                 tautomerLambdas[atomIndex] = lambda;
1295             }
1296         } /*else {
1297             logger.warning(format("This residue %s does not have any titrating tautomers.", residue.getName()));
1298         }*/
1299     }
1300 
1301     /**
1302      * Set the titration lambda of a residue and update corresponding theta
1303      *
1304      * @param residue residue to set the lambda of
1305      * @param lambda  value to set the residue to
1306      */
1307     public void setTitrationLambda(Residue residue, double lambda) {
1308         setTitrationLambda(residue, lambda, true);
1309     }
1310 
1311     /**
1312      * Set the titration lambda of a residue and update corresponding theta if desired
1313      *
1314      * @param residue      residue to set the lambda of
1315      * @param lambda       value to set the residue to
1316      * @param changeThetas whether to change the theta positions ~comes with information loss~
1317      */
1318     public void setTitrationLambda(Residue residue, double lambda, boolean changeThetas) {
1319         double[] thetaPosition = esvState.x();
1320         if (titratingResidueList.contains(residue) && !lockStates) {
1321             int index = titratingResidueList.indexOf(residue);
1322             extendedLambdas[index] = lambda;
1323             if (changeThetas) {
1324                 thetaPosition[index] = Math.asin(Math.sqrt(lambda));
1325             }
1326             List<Atom> currentAtomList = residue.getSideChainAtoms();
1327             for (Atom atom : currentAtomList) {
1328                 int atomIndex = atom.getArrayIndex();
1329                 titrationLambdas[atomIndex] = lambda;
1330             }
1331         }
1332     }
1333 
1334     /**
1335      * Overwrites the histogram passed into it and returns the new one out ~output never used?~
1336      *
1337      * @param histogram 2D histogram list with the tautomer and titration states compressed to a 1D array
1338      * @return another compressed histogram
1339      */
1340     public int[][] getESVHistogram(int[][] histogram) {
1341         for (int i = 0; i < titratingResidueList.size(); i++) {
1342             int h = 0;
1343             for (int j = 0; j < 10; j++) {
1344                 for (int k = 0; k < 10; k++) {
1345                     histogram[i][h++] = esvHistogram[i][j][k];
1346                 }
1347             }
1348         }
1349         return histogram;
1350     }
1351 
1352     /**
1353      * Changes this ESV's histogram to equal the one passed
1354      *
1355      * @param histogram histogram to set this ESV histogram to
1356      */
1357     public void copyESVHistogramTo(int[][] histogram) {
1358         for (int i = 0; i < titratingResidueList.size(); i++) {
1359             int h = 0;
1360             for (int j = 0; j < 10; j++) {
1361                 for (int k = 0; k < 10; k++) {
1362                     esvHistogram[i][j][k] = histogram[i][h++];
1363                 }
1364             }
1365         }
1366     }
1367 
1368     /**
1369      * Zero out each array element of the vdW ESV array
1370      */
1371     public void initEsvVdw() {
1372         for (int i = 0; i < nESVs; i++) {
1373             esvVdwDerivs[i].set(0.0);
1374         }
1375     }
1376 
1377     /**
1378      * Zero out each array element of the permElec ESV array
1379      */
1380     public void initEsvPermElec() {
1381         for (int i = 0; i < nESVs; i++) {
1382             esvPermElecDerivs[i].set(0.0);
1383         }
1384     }
1385 
1386     /**
1387      * Zero out each array element of the indElec ESV array
1388      */
1389     public void initEsvIndElec() {
1390         for (int i = 0; i < nESVs; i++) {
1391             esvIndElecDerivs[i].set(0.0);
1392         }
1393     }
1394 
1395     /**
1396      * get Titration Lambda for an extended atom
1397      *
1398      * @param atomIndex
1399      * @return titrationLambdas[atomIndex]
1400      */
1401     public double getTitrationLambda(int atomIndex) {
1402         return titrationLambdas[atomIndex];
1403     }
1404 
1405     /**
1406      * get the index of the extended residue list that corresponds to this atom
1407      *
1408      * @param i
1409      * @return titrationIndexMap[i]
1410      */
1411     public int getTitrationESVIndex(int i) {
1412         return titrationIndexMap[i];
1413     }
1414 
1415     /**
1416      * get Tautomer Lambda for an extended atom
1417      *
1418      * @param atomIndex
1419      * @return tautomerLambdas[atomIndex]
1420      */
1421     public double getTautomerLambda(int atomIndex) {
1422         return tautomerLambdas[atomIndex];
1423     }
1424 
1425     public int getTautomerESVIndex(int i) {
1426         return tautomerIndexMap[i];
1427     }
1428 
1429     /**
1430      * Return the List of Titrating Residues
1431      *
1432      * @return titratingResidueList
1433      */
1434     public List<Residue> getTitratingResidueList() {
1435         return titratingResidueList;
1436     }
1437 
1438     /**
1439      * Return the List of Tautomerizing Residues
1440      *
1441      * @return tautomerizingResidueList
1442      */
1443     public List<Residue> getTautomerizingResidueList() {
1444         return tautomerizingResidueList;
1445     }
1446 
1447     /**
1448      * Return the List of Extended Residues which = TitratingResidueList + TautomerizingResidueList
1449      *
1450      * @return extendedResidueList
1451      */
1452     public List<Residue> getExtendedResidueList() {
1453         return extendedResidueList;
1454     }
1455 
1456     public double getThetaMass() {
1457         return thetaMass;
1458     }
1459 
1460     public double getThetaFriction() {
1461         return thetaFriction;
1462     }
1463 
1464     /**
1465      * Gets a copy of the array of doubles that matches the nESVs correspoding to each titration and tautomer lambda
1466      *
1467      * @return double array of length nESVs
1468      */
1469     public double[] getExtendedLambdas() {
1470         double[] lambdas = new double[nESVs];
1471         System.arraycopy(extendedLambdas, 0, lambdas, 0, lambdas.length);
1472         return lambdas;
1473     }
1474 
1475     /**
1476      * getLambdaList.
1477      *
1478      * @return a {@link String} object.
1479      */
1480     public String getLambdaList() {
1481         if (nESVs < 1) {
1482             return "";
1483         }
1484         StringBuilder sb = new StringBuilder();
1485         for (int i = 0; i < nESVs; i++) {
1486             if (i == 0) {
1487                 sb.append("\n  Titration Lambdas: ");
1488             }
1489             if (i > 0) {
1490                 sb.append(", ");
1491             }
1492             if (i == nTitr) {
1493                 sb.append("\n  Tautomer Lambdas: ");
1494             }
1495             sb.append(format("%6.4f", extendedLambdas[i]));
1496         }
1497         return sb.toString();
1498     }
1499 
1500     /**
1501      * All atoms of the fully-protonated system (not just those affected by this system).
1502      *
1503      * @return an array of {@link Atom} objects.
1504      */
1505     public Atom[] getExtendedAtoms() {
1506         return extendedAtoms;
1507     }
1508 
1509     /**
1510      * Companion to getExtendedAtoms() for vdw::setAtoms and pme::setAtoms.
1511      *
1512      * @return an array of {@link int} objects.
1513      */
1514     public int[] getExtendedMolecule() {
1515         return extendedMolecules;
1516     }
1517 
1518     /**
1519      * Sum up total bias (Ubias = UpH + Udiscr - Umod)
1520      *
1521      * @return totalBiasEnergy
1522      */
1523     public double getBiasEnergy() {
1524         double totalBiasEnergy = 0.0;
1525         double[] biasEnergyComponents = new double[9];
1526         //Do not double count residues in tautomer list.
1527         for (Residue residue : titratingResidueList) {
1528             getBiasTerms(residue, biasEnergyComponents);
1529             double biasEnergy = biasEnergyComponents[discrBiasIndex] + biasEnergyComponents[pHBiasIndex] + biasEnergyComponents[modelBiasIndex];
1530             totalBiasEnergy += biasEnergy;
1531         }
1532         return totalBiasEnergy;
1533     }
1534 
1535     /**
1536      * getBiasDecomposition.
1537      *
1538      * @return a {@link String} object.
1539      */
1540     public String getBiasDecomposition() {
1541         double discrBias = 0.0;
1542         double phBias = 0.0;
1543         double modelBias = 0.0;
1544         double[] biasEnergyAndDerivs = new double[9];
1545         for (Residue residue : titratingResidueList) {
1546             getBiasTerms(residue, biasEnergyAndDerivs);
1547             discrBias += biasEnergyAndDerivs[discrBiasIndex];
1548             phBias += biasEnergyAndDerivs[pHBiasIndex];
1549             //Reminder: Ubias = UpH + Udiscr - Umod but the Umod terms already have their sign reversed
1550             modelBias += biasEnergyAndDerivs[modelBiasIndex];
1551         }
1552         return format("    %-16s %16.8f\n", "Discretizer", discrBias)
1553                 + format("    %-16s %16.8f\n", "Acidostat", phBias)
1554                 + format("    %-16s %16.8f\n", "Fmod", modelBias);
1555     }
1556 
1557     /**
1558      * Calculate prefactor for scaling the van der Waals based on titration/tautomer state if titrating proton
1559      */
1560     public void getVdwPrefactor(int atomIndex, double[] vdwPrefactorAndDerivs) {
1561         double prefactor = 1.0;
1562         double titrationDeriv = 0.0;
1563         double tautomerDeriv = 0.0;
1564         if (!isTitratingHydrogen(atomIndex) || !doVDW) {
1565             vdwPrefactorAndDerivs[0] = prefactor;
1566             vdwPrefactorAndDerivs[1] = titrationDeriv;
1567             vdwPrefactorAndDerivs[2] = tautomerDeriv;
1568             return;
1569         }
1570         AminoAcid3 AA3 = residueNames[atomIndex];
1571         switch (AA3) {
1572             case ASD:
1573             case GLD:
1574                 if (tautomerDirections[atomIndex] == 1) {
1575                     prefactor = titrationLambdas[atomIndex] * tautomerLambdas[atomIndex];
1576                     titrationDeriv = tautomerLambdas[atomIndex];
1577                     tautomerDeriv = titrationLambdas[atomIndex];
1578                 } else if (tautomerDirections[atomIndex] == -1) {
1579                     prefactor = titrationLambdas[atomIndex] * (1.0 - tautomerLambdas[atomIndex]);
1580                     titrationDeriv = (1.0 - tautomerLambdas[atomIndex]);
1581                     tautomerDeriv = -titrationLambdas[atomIndex];
1582                 }
1583                 break;
1584             case HIS:
1585                 if (tautomerDirections[atomIndex] == 1) {
1586                     prefactor = (1.0 - titrationLambdas[atomIndex]) * tautomerLambdas[atomIndex] + titrationLambdas[atomIndex];
1587                     titrationDeriv = -tautomerLambdas[atomIndex] + 1.0;
1588                     tautomerDeriv = (1 - titrationLambdas[atomIndex]);
1589                 } else if (tautomerDirections[atomIndex] == -1) {
1590                     prefactor = (1.0 - titrationLambdas[atomIndex]) * (1.0 - tautomerLambdas[atomIndex]) + titrationLambdas[atomIndex];
1591                     titrationDeriv = tautomerLambdas[atomIndex];
1592                     tautomerDeriv = -(1.0 - titrationLambdas[atomIndex]);
1593                 }
1594                 break;
1595             case LYS:
1596             case CYS:
1597                 prefactor = titrationLambdas[atomIndex];
1598                 titrationDeriv = 1.0;
1599                 tautomerDeriv = 0.0;
1600                 break;
1601         }
1602         vdwPrefactorAndDerivs[0] = prefactor;
1603         vdwPrefactorAndDerivs[1] = titrationDeriv;
1604         vdwPrefactorAndDerivs[2] = tautomerDeriv;
1605     }
1606 
1607     /**
1608      * Questions whether the current hydrogen's polarizability is changing in response to lambda being updated.
1609      *
1610      * @param atomIndex
1611      * @return
1612      */
1613     public boolean isTitratingHydrogen(int atomIndex) {
1614         return isTitratingHydrogen[atomIndex];
1615     }
1616 
1617     /**
1618      * Add van der Waals deriv to appropriate dU/dL term given the atom index and its contributions.
1619      *
1620      * @param atomI
1621      * @param vdwEnergy
1622      * @param vdwPrefactorAndDerivI
1623      * @param vdwPrefactorJ
1624      */
1625     public void addVdwDeriv(int atomI, double vdwEnergy, double[] vdwPrefactorAndDerivI, double vdwPrefactorJ) {
1626         if (!isTitratingHydrogen(atomI)) {
1627             return;
1628         }
1629 
1630         //Sum up dU/dL for titration ESV if atom i is titrating hydrogen
1631         //Sum up dU/dL for tautomer ESV if atom i is titrating hydrogen
1632         int titrationEsvIndex = titrationIndexMap[atomI];
1633         int tautomerEsvIndex = tautomerIndexMap[atomI] + nTitr;// tautomerIndexMap[atomI] = -1 for non tautomerizing residues
1634         double dTitr_dLambda;
1635         double dTaut_dLambda;
1636 
1637         dTitr_dLambda = vdwPrefactorAndDerivI[1] * vdwPrefactorJ * vdwEnergy;
1638         dTaut_dLambda = vdwPrefactorAndDerivI[2] * vdwPrefactorJ * vdwEnergy;
1639 
1640         esvVdwDerivs[titrationEsvIndex].addAndGet(dTitr_dLambda);
1641         if (tautomerEsvIndex >= nTitr) {
1642             esvVdwDerivs[tautomerEsvIndex].addAndGet(dTaut_dLambda);
1643         }
1644     }
1645 
1646     /**
1647      * Add Perm Elec deriv to appropriate dU/dL term given the atom index and its contributions.
1648      *
1649      * @param atomI
1650      * @param titrationEnergy
1651      * @param tautomerEnergy
1652      */
1653     public void addPermElecDeriv(int atomI, double titrationEnergy, double tautomerEnergy) {
1654         int titrationEsvIndex = titrationIndexMap[atomI];
1655         int tautomerEsvIndex = tautomerIndexMap[atomI] + nTitr;// tautomerIndexMap[atomI] = -1 for non tautomerizing residues
1656         esvPermElecDerivs[titrationEsvIndex].addAndGet(titrationEnergy);
1657         if (tautomerEsvIndex >= nTitr) {
1658             esvPermElecDerivs[tautomerEsvIndex].addAndGet(tautomerEnergy);
1659         }
1660     }
1661 
1662     /**
1663      * Add Induced Elec deriv to appropriate dU/dL term given the atom index and its contributions.
1664      *
1665      * @param atomI
1666      * @param titrationEnergy
1667      * @param tautomerEnergy
1668      */
1669     public void addIndElecDeriv(int atomI, double titrationEnergy, double tautomerEnergy) {
1670         int titrationEsvIndex = titrationIndexMap[atomI];
1671         int tautomerEsvIndex = tautomerIndexMap[atomI] + nTitr; // tautomerIndexMap[atomI] = -1 for non tautomerizing residues
1672         esvIndElecDerivs[titrationEsvIndex].addAndGet(titrationEnergy);
1673         if (tautomerEsvIndex >= nTitr) {
1674             esvIndElecDerivs[tautomerEsvIndex].addAndGet(tautomerEnergy);
1675         }
1676     }
1677 
1678     /**
1679      * getConstantPh.
1680      *
1681      * @return a double.
1682      */
1683     public double getConstantPh() {
1684         return constantSystemPh;
1685     }
1686 
1687     /**
1688      * setConstantPh.
1689      *
1690      * @param pH a double.
1691      */
1692     public void setConstantPh(double pH) {
1693         constantSystemPh = pH;
1694     }
1695 
1696     /**
1697      * @return titrationUtils master copy from ExtendedSystem.
1698      */
1699     public TitrationUtils getTitrationUtils() {
1700         return titrationUtils;
1701     }
1702 
1703     /**
1704      * Sets the restartFile field of this extended system to the passed file. This does not read the file, it only
1705      * determines where the writeRestartFile() will write to.
1706      *
1707      * @param esvFile
1708      */
1709     public void setRestartFile(File esvFile) {
1710         restartFile = esvFile;
1711     }
1712 
1713     /**
1714      * Writes out the current state of the extended system to the specified file without setting the file to that location.
1715      *
1716      * @param esvFile file to be written to
1717      * @return whether the read was successful or not
1718      */
1719     public boolean writeESVInfoTo(File esvFile) {
1720         logger.info("Writing pH Dynamics out to: " + esvFile.getParentFile().getName() + File.separator + esvFile.getName());
1721         double[] thetaPosition = esvState.x();
1722         double[] thetaVelocity = esvState.v();
1723         double[] thetaAccel = esvState.a();
1724         return esvFilter.writeESV(esvFile, thetaPosition, thetaVelocity, thetaAccel, titratingResidueList, esvHistogram, constantSystemPh);
1725     }
1726 
1727     /**
1728      * Method overwrites whatever is in the extended system at the time with the read data.
1729      * <p>
1730      * CAUTION: If the old data is not written out to file before this is called, the data will be lost.
1731      *
1732      * @param esvFile esvFile to read
1733      * @return whether the read was successful or not
1734      */
1735     public boolean readESVInfoFrom(File esvFile) {
1736         double[] thetaPosition = esvState.x();
1737         double[] thetaVelocity = esvState.v();
1738         double[] thetaAccel = esvState.a();
1739         return esvFilter.readESV(esvFile, thetaPosition, thetaVelocity, thetaAccel, esvHistogram);
1740     }
1741 
1742     /**
1743      * setTemperature.
1744      *
1745      * @param set a double.
1746      */
1747     public void setTemperature(double set) {
1748         currentTemperature = set;
1749     }
1750 
1751     /**
1752      * Processes lambda values based on propagation of theta value from Stochastic integrator in Molecular dynamics
1753      */
1754     public void preForce() {
1755         updateLambdas();
1756     }
1757 
1758     /**
1759      * Execute writeESV from esvFilter that is contained within ExtendedSystem
1760      */
1761     public void writeRestart() {
1762         String esvName = FileUtils.relativePathTo(restartFile).toString();
1763         double[] thetaPosition = esvState.x();
1764         double[] thetaVelocity = esvState.v();
1765         double[] thetaAccel = esvState.a();
1766         if (esvFilter.writeESV(restartFile, thetaPosition, thetaVelocity, thetaAccel, titratingResidueList, esvHistogram, constantSystemPh)) {
1767             logger.info(" Wrote PhDynamics restart file to " + esvName);
1768         } else {
1769             logger.info(" Writing PhDynamics restart file to " + esvName + " failed");
1770         }
1771     }
1772 
1773     /**
1774      * Execute writeLambdaHistogrm from esvFilter that is contained within ExtendedSystem
1775      * Prints the ESV histogram for each titrating residue
1776      */
1777     public void writeLambdaHistogram(boolean printHistograms) {
1778         printProtonationRatios();
1779         if (printHistograms) {
1780             logger.info(esvFilter.getLambdaHistogram(titratingResidueList, esvHistogram, constantSystemPh));
1781         }
1782     }
1783 
1784     /**
1785      * Calculate the Deprotonation Fraction from the ESV histogram
1786      */
1787     public void printProtonationRatios() {
1788         for (int i = 0; i < esvHistogram.length; i++) {
1789             int[] rowSums = new int[esvHistogram[i].length];
1790             for (int j = 0; j < esvHistogram[i].length; j++) {
1791                 for (int k = 0; k < esvHistogram[i][j].length; k++) {
1792                     rowSums[j] += esvHistogram[i][j][k];
1793                 }
1794             }
1795             int i1 = rowSums[0] + rowSums[rowSums.length - 1];
1796             double buf = i1 == 0 ? 0.0 : .001;
1797             logger.info(" " + extendedResidueList.get(i).toString() + " Deprotonation Fraction at pH " + constantSystemPh + ": " + (rowSums[0] / (i1 + buf)));
1798             if (buf == 0.0) {
1799                 logger.info(" Buffer required to avoid division by 0");
1800             }
1801         }
1802     }
1803 
1804     @Override
1805     public double[] getAcceleration(double[] acceleration) {
1806         return getThetaAccel();
1807     }
1808 
1809     public double[] getThetaAccel() {
1810         return esvState.a();
1811     }
1812 
1813     @Override
1814     public double energyAndGradient(double[] x, double[] g) {
1815         double[] thetaPosition = esvState.x();
1816         System.arraycopy(x, 0, thetaPosition, 0, thetaPosition.length);
1817         updateLambdas();
1818         fillESVgradient(g);
1819         return energy(x);
1820     }
1821 
1822     @Override
1823     public List<Constraint> getConstraints() {
1824         return constraints.isEmpty() ? Collections.emptyList() : new ArrayList<>(constraints);
1825     }
1826 
1827     @Override
1828     public STATE getEnergyTermState() {
1829         return null;
1830     }
1831 
1832     private double[] fillESVgradient(double[] g) {
1833         double[] gradESV = postForce();
1834         System.arraycopy(gradESV, 0, g, 0, g.length);
1835         return g;
1836     }
1837 
1838     @Override
1839     public void setEnergyTermState(STATE state) {
1840 
1841     }
1842 
1843     @Override
1844     public double[] getMass() {
1845         return getThetaMassArray();
1846     }
1847 
1848     /**
1849      * Applies a chain rule term to the derivative to account for taking a derivative of lambda = sin(theta)^2
1850      *
1851      * @return dE/dL a double[]
1852      */
1853     public double[] postForce() {
1854         double[] thetaPosition = esvState.x();
1855         double[] dEdL = ExtendedSystem.this.getDerivatives();
1856         double[] dEdTheta = new double[dEdL.length];
1857         for (int i = 0; i < nESVs; i++) {
1858             dEdTheta[i] = dEdL[i] * sin(2 * thetaPosition[i]);
1859             //logger.info("dEdL["+i+"]: "+dEdL[i]);
1860         }
1861         return dEdTheta;
1862     }
1863 
1864     public double[] getThetaMassArray() {
1865         return esvState.getMass();
1866     }
1867 
1868     @Override
1869     public double[] getPreviousAcceleration(double[] previousAcceleration) {
1870         return new double[0];
1871     }
1872     @Override
1873     public double energy(double[] x) {
1874         double[] coords = new double[forceFieldEnergy.getNumberOfVariables()];
1875         forceFieldEnergy.getCoordinates(coords);
1876         return forceFieldEnergy.energy(coords);
1877     }
1878 
1879     @Override
1880     public VARIABLE_TYPE[] getVariableTypes() {
1881         return new VARIABLE_TYPE[0];
1882     }
1883 
1884     @Override
1885     public double[] getVelocity(double[] velocity) {
1886         return getThetaVelocity();
1887     }
1888 
1889     public double[] getThetaVelocity() {
1890         return esvState.v();
1891     }
1892 
1893     @Override
1894     public void setAcceleration(double[] acceleration) {
1895 
1896     }
1897     @Override
1898     public double[] getCoordinates(double[] parameters) {
1899         return getThetaPosition();
1900     }
1901 
1902     @Override
1903     public void setPreviousAcceleration(double[] previousAcceleration) {
1904 
1905     }
1906 
1907     @Override
1908     public void setVelocity(double[] velocity) {
1909 
1910     }
1911 
1912     public SystemState getState() {
1913         return esvState;
1914     }
1915 
1916     public double[] getThetaPosition() {
1917         return esvState.x();
1918     }
1919 
1920     @Override
1921     public int getNumberOfVariables() {
1922         return nESVs;
1923     }
1924 
1925     @Override
1926     public double[] getScaling() {
1927         double[] scaling = new double[nESVs];
1928         Arrays.fill(scaling, 1.0);
1929         return scaling;
1930     }
1931 
1932     @Override
1933     public void setScaling(double[] scaling) {
1934 
1935     }
1936 
1937     @Override
1938     public double getTotalEnergy() {
1939         return 0;
1940     }
1941 
1942 }