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