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