// ******************************************************************************
   //
   // Title:       Force Field X.
   // Description: Force Field X - Software for Molecular Biophysics.
   // Copyright:   Copyright (c) Michael J. Schnieders 2001-2025.
   //
   // This file is part of Force Field X.
   //
   // Force Field X is free software; you can redistribute it and/or modify it
  // under the terms of the GNU General Public License version 3 as published by
  // the Free Software Foundation.
  //
  // Force Field X is distributed in the hope that it will be useful, but WITHOUT
  // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  // details.
  //
  // You should have received a copy of the GNU General Public License along with
  // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
  // Place, Suite 330, Boston, MA 02111-1307 USA
  //
  // Linking this library statically or dynamically with other modules is making a
  // combined work based on this library. Thus, the terms and conditions of the
  // GNU General Public License cover the whole combination.
  //
  // As a special exception, the copyright holders of this library give you
  // permission to link this library with independent modules to produce an
  // executable, regardless of the license terms of these independent modules, and
  // to copy and distribute the resulting executable under terms of your choice,
  // provided that you also meet, for each linked independent module, the terms
  // and conditions of the license of that module. An independent module is a
  // module which is not derived from or based on this library. If you modify this
  // library, you may extend this exception to your version of the library, but
  // you are not obligated to do so. If you do not wish to do so, delete this
  // exception statement from your version.
  //
  // ******************************************************************************
  package ffx.algorithms.thermodynamics;
  
  import edu.rit.mp.DoubleBuf;
  import edu.rit.pj.Comm;
  import ffx.algorithms.AlgorithmListener;
  import ffx.algorithms.cli.DynamicsOptions;
  import ffx.algorithms.cli.LambdaParticleOptions;
  import ffx.algorithms.dynamics.Barostat;
  import ffx.algorithms.dynamics.integrators.Stochastic;
  import ffx.algorithms.optimize.Minimize;
  import ffx.crystal.Crystal;
  import ffx.crystal.CrystalPotential;
  import ffx.numerics.Potential;
  import ffx.numerics.integrate.DataSet;
  import ffx.numerics.integrate.DoublesDataSet;
  import ffx.numerics.integrate.Integrate1DNumeric;
  import ffx.numerics.integrate.Integrate1DNumeric.IntegrationType;
  import ffx.potential.MolecularAssembly;
  import ffx.potential.SystemState;
  import ffx.potential.Utilities;
  import ffx.potential.bonded.LambdaInterface;
  import ffx.potential.parsers.PDBFilter;
  import ffx.potential.parsers.SystemFilter;
  import ffx.potential.parsers.XYZFilter;
  import org.apache.commons.configuration2.CompositeConfiguration;
  import org.apache.commons.io.FilenameUtils;
  
  import javax.annotation.Nullable;
  import java.io.File;
  import java.io.PrintWriter;
  import java.util.ArrayList;
  import java.util.Arrays;
  import java.util.List;
  import java.util.Optional;
  import java.util.logging.Level;
  import java.util.logging.Logger;
  
  import static ffx.numerics.integrate.Integrate1DNumeric.IntegrationType.SIMPSONS;
  import static ffx.utilities.Constants.R;
  import static java.lang.String.format;
  import static java.lang.System.arraycopy;
  import static java.util.Arrays.fill;
  import static org.apache.commons.math3.util.FastMath.PI;
  import static org.apache.commons.math3.util.FastMath.abs;
  import static org.apache.commons.math3.util.FastMath.asin;
  import static org.apache.commons.math3.util.FastMath.exp;
  import static org.apache.commons.math3.util.FastMath.floor;
  import static org.apache.commons.math3.util.FastMath.sin;
  import static org.apache.commons.math3.util.FastMath.sqrt;
  
  /**
   * An implementation of the Orthogonal Space Tempering algorithm.
   *
   * <p>This only partially implements the LambdaInterface, since it does not return 2nd lambda
   * derivatives. The 2nd derivatives of the bias require 3rd derivatives of the underlying Hamiltonian
   * (not these derivatives are not needed for OST MD).
   *
   * @author Michael J. Schnieders, James Dama, Wei Yang and Pengyu Ren
   * @since 1.0
   */
  public class OrthogonalSpaceTempering implements CrystalPotential, LambdaInterface {
  
     private static final Logger logger = Logger.getLogger(OrthogonalSpaceTempering.class.getName());
     /**
      * The potential energy of the system.
      */
     protected final CrystalPotential potential;
     /**
      * The AlgorithmListener is called each time a count is added.
      */
     protected final AlgorithmListener algorithmListener;
     /**
      * Print detailed energy information.
      */
     protected final boolean print = false;
     /**
      * Number of variables.
      */
     protected final int nVariables;
     /**
      * A potential energy that implements the LambdaInterface.
      */
     private final LambdaInterface lambdaInterface;
     /**
      * List of additional Histograms this OST can switch to.
      */
     private final List<Histogram> allHistograms = new ArrayList<>();
     /**
      * Parameters to control saving local optimizations.
      */
     private final OptimizationParameters optimizationParameters;
     /**
      * Properties.
      */
     private final CompositeConfiguration properties;
     /**
      * The MolecularAssembly being simulated.
      */
     protected MolecularAssembly molecularAssembly;
     /**
      * Are FAST varying energy terms being computed, SLOW varying energy terms, or BOTH. OST is not
      * active when only FAST varying energy terms are being propagated.
      */
     protected Potential.STATE state = Potential.STATE.BOTH;
     /**
      * Force Field Potential Energy (i.e. with no bias terms added).
      */
     protected double forceFieldEnergy;
     /**
      * Contains counts for the OST bias.
      */
     private Histogram histogram;
     /**
      * Index of the current Histogram.
      */
     private int histogramIndex;
     /**
      * Flag to indicate that the Lambda particle should be propagated.
      */
     private boolean propagateLambda = true;
     /**
      * Mixed second partial derivative with respect to coordinates and lambda.
      */
     private final double[] dUdXdL;
 
     /**
      * Gradient array needed when the OST Energy method is called.
      */
     private final double[] tempGradient;
 
     /**
      * Partial derivative of the force field energy with respect to lambda.
      */
     private double dForceFieldEnergydL;
     /**
      * Magnitude of the 2D orthogonal space bias G(L,dE/dL).
      */
     private double gLdEdL = 0.0;
     /**
      * OST Bias energy.
      */
     private double biasEnergy;
     /**
      * Total system energy.
      */
     private double totalEnergy;
     /**
      * Total partial derivative of the potential (U) being sampled with respect to lambda.
      */
     private double dUdLambda;
     /**
      * Second partial derivative of the potential being sampled with respect to lambda.
      */
     private double d2UdL2;
     /**
      * If true, values of (lambda, dU/dL) that have not been observed are rejected.
      */
     private boolean hardWallConstraint = false;
 
     private final DynamicsOptions dynamicsOptions;
 
     private final LambdaParticleOptions lambdaParticleOptions;
 
     /**
      * OST Constructor.
      *
      * @param lambdaInterface       defines Lambda and dU/dL.
      * @param potential             defines the Potential energy.
      * @param histogramData         contains histogram restart data.
      * @param lambdaData            contains lambda restart data.
      * @param properties            defines System properties.
      * @param dynamicsOptions       defines molecular dynamics parameters.
      * @param lambdaParticleOptions defines lambda particle parameters.
      * @param algorithmListener     the AlgorithmListener to be notified of progress.
      */
     public OrthogonalSpaceTempering(LambdaInterface lambdaInterface, CrystalPotential potential,
                                     HistogramData histogramData, LambdaData lambdaData, CompositeConfiguration properties,
                                     DynamicsOptions dynamicsOptions, LambdaParticleOptions lambdaParticleOptions,
                                     AlgorithmListener algorithmListener) {
         this.lambdaInterface = lambdaInterface;
         this.potential = potential;
         this.properties = properties;
         this.dynamicsOptions = dynamicsOptions;
         this.lambdaParticleOptions = lambdaParticleOptions;
         this.algorithmListener = algorithmListener;
         nVariables = potential.getNumberOfVariables();
 
         if (potential instanceof Barostat) {
             logger.severe(" The Barostat must wrap the OST instance.");
         }
 
         dUdXdL = new double[nVariables];
         tempGradient = new double[nVariables];
 
         // Init the Histogram.
         histogram = new Histogram(properties, histogramData, lambdaData);
         histogramIndex = 0;
         allHistograms.add(histogram);
 
         // Configure optimization parameters.
         optimizationParameters = new OptimizationParameters(properties);
     }
 
     /**
      * Add an alternate Histogram this OST can use.
      *
      * @param histogramData Settings to use for the new Histogram.
      */
     public void addHistogram(HistogramData histogramData, LambdaData lambdaData) {
         Histogram newHisto = new Histogram(properties, histogramData, lambdaData);
         histogramData.asynchronous = this.histogram.hd.asynchronous;
         allHistograms.add(newHisto);
     }
 
     /**
      * {@inheritDoc}
      */
     @Override
     public boolean dEdLZeroAtEnds() {
         return false;
     }
 
     /**
      * {@inheritDoc}
      */
     @Override
     public boolean destroy() {
         // Shut down the CountReceiveThread.
         histogram.destroy();
         return potential.destroy();
     }
 
     /**
      * Compute the force field + bias energy.
      */
     public double energy(double[] x) {
 
         // OST is propagated with the slowly varying terms.
         if (state == Potential.STATE.FAST) {
             forceFieldEnergy = potential.energy(x);
             return forceFieldEnergy;
         }
 
         // Stochastic dynamics pre-force propagation.
         if (propagateLambda) {
             histogram.stochasticPreForce();
         }
 
         // We have to compute the energy and gradient, since we need dU/dL to be computed.
         fill(tempGradient, 0.0);
         forceFieldEnergy = potential.energyAndGradient(x, tempGradient);
 
         dForceFieldEnergydL = lambdaInterface.getdEdL();
         d2UdL2 = lambdaInterface.getd2EdL2();
         int lambdaBin = histogram.indexForLambda();
         dUdLambda = dForceFieldEnergydL;
 
         gLdEdL = 0.0;
         double bias1D;
         // Update the free energy difference.
         histogram.updateFreeEnergyDifference(false, false);
         if (histogram.hd.metaDynamics) {
             bias1D = histogram.energyAndGradientMeta(true);
         } else {
             // Calculate recursion kernel G(L, F_L) and its derivatives with respect to L and F_L.
             if (histogram.hd.biasMag > 0.0) {
                 double[] chainRule = new double[2];
                 gLdEdL = histogram.energyAndGradient2D(dUdLambda, chainRule);
                 double dGdLambda = chainRule[0];
                 double dGdFLambda = chainRule[1];
                 dUdLambda += dGdLambda + dGdFLambda * d2UdL2;
             }
 
             // Compute the energy and gradient for the recursion worker at F(L) using interpolation.
             bias1D = histogram.energyAndGradient1D(true);
         }
 
         // The total bias energy is the sum of the 1D and 2D terms.
         biasEnergy = bias1D + gLdEdL;
 
         if (print) {
             logger.info(format(" Bias Energy        %16.8f", biasEnergy));
             logger.info(format(" %s %16.8f  (Kcal/mole)", "OST Potential    ", forceFieldEnergy + biasEnergy));
         }
 
         if (propagateLambda) {
             histogram.stochasticPostForce();
 
             histogram.ld.stepsTaken++;
             long energyCount = histogram.ld.stepsTaken;
 
             // Log the current Lambda state.
             int printFrequency = dynamicsOptions.getReportFrequency(100);
             if (energyCount % printFrequency == 0) {
                 double dBdL = dUdLambda - dForceFieldEnergydL;
                 double lambda = histogram.ld.lambda;
                 int lambdaBins = histogram.hd.getLambdaBins();
                 if (lambdaBins < 1000) {
                     logger.info(format(" L=%6.4f (%3d) F_LU=%10.4f F_LB=%10.4f F_L=%10.4f V_L=%10.4f", lambda,
                             lambdaBin, dForceFieldEnergydL, dBdL, dUdLambda, histogram.ld.thetaVelocity));
                 } else {
                     logger.info(format(" L=%6.4f (%4d) F_LU=%10.4f F_LB=%10.4f F_L=%10.4f V_L=%10.4f", lambda,
                             lambdaBin, dForceFieldEnergydL, dBdL, dUdLambda, histogram.ld.thetaVelocity));
                 }
             }
 
             // Metadynamics grid counts (every 'countInterval' steps).
             if (energyCount % histogram.hd.countInterval == 0) {
                 histogram.addBias(dForceFieldEnergydL);
 
                 // Locally optimize the current state.
                 if (optimizationParameters.doOptimization) {
                     optimizationParameters.optimize(forceFieldEnergy, x, tempGradient);
                 }
 
                 if (algorithmListener != null) {
                     algorithmListener.algorithmUpdate(molecularAssembly);
                 }
             }
         }
 
         totalEnergy = forceFieldEnergy + biasEnergy;
 
         return totalEnergy;
     }
 
     /**
      * {@inheritDoc}
      */
     @Override
     public double energyAndGradient(double[] x, double[] gradient) {
 
         // Stochastic dynamics pre-force propagation.
         if (state != STATE.FAST && propagateLambda) {
             histogram.stochasticPreForce();
         }
 
         // Compute the force field energy and gradient.
         forceFieldEnergy = potential.energyAndGradient(x, gradient);
 
         // OST is propagated with the slowly varying terms.
         if (state == STATE.FAST) {
             return forceFieldEnergy;
         }
 
         // dU/dL is the partial derivative of the force field energy with respect to lambda.
         dForceFieldEnergydL = lambdaInterface.getdEdL();
         dUdLambda = dForceFieldEnergydL;
         d2UdL2 = lambdaInterface.getd2EdL2();
 
         gLdEdL = 0.0;
         double bias1D;
         // Update the free energy difference.
         histogram.updateFreeEnergyDifference(false, false);
         if (histogram.hd.metaDynamics) {
             bias1D = histogram.energyAndGradientMeta(true);
         } else {
             if (histogram.hd.biasMag > 0) {
                 double[] chainRule = new double[2];
                 gLdEdL = histogram.energyAndGradient2D(dUdLambda, chainRule);
                 double dGdLambda = chainRule[0];
                 double dGdFLambda = chainRule[1];
 
                 // Lambda gradient due to recursion kernel G(L, F_L).
                 dUdLambda += dGdLambda + dGdFLambda * d2UdL2;
 
                 // Cartesian coordinate gradient due to recursion kernel G(L, F_L).
                 fill(dUdXdL, 0.0);
                 lambdaInterface.getdEdXdL(dUdXdL);
                 for (int i = 0; i < nVariables; i++) {
                     gradient[i] += dGdFLambda * dUdXdL[i];
                 }
             }
 
             // Compute the energy and gradient for the recursion worker at F(L) using interpolation.
             bias1D = histogram.energyAndGradient1D(true);
         }
 
         // The total bias is the sum of 1D and 2D terms.
         biasEnergy = bias1D + gLdEdL;
 
         if (print) {
             logger.info(format(" %s %16.8f", "Bias Energy       ", biasEnergy));
             logger.info(format(" %s %16.8f  %s", "OST Potential    ", forceFieldEnergy + biasEnergy, "(Kcal/mole)"));
         }
 
         if (propagateLambda) {
             histogram.stochasticPostForce();
 
             histogram.ld.stepsTaken++;
             long energyCount = histogram.ld.stepsTaken;
 
             // Log the current Lambda state.
             int printFrequency = dynamicsOptions.getReportFrequency(100);
             if (energyCount % printFrequency == 0) {
                 double dBdL = dUdLambda - dForceFieldEnergydL;
                 int lambdaBin = histogram.indexForLambda();
                 double lambda = histogram.ld.lambda;
                 int lambdaBins = histogram.hd.getLambdaBins();
                 if (lambdaBins < 1000) {
                     logger.info(format(" L=%6.4f (%3d) F_LU=%10.4f F_LB=%10.4f F_L=%10.4f V_L=%10.4f", lambda,
                             lambdaBin, dForceFieldEnergydL, dBdL, dUdLambda, histogram.ld.thetaVelocity));
                 } else {
                     logger.info(format(" L=%6.4f (%4d) F_LU=%10.4f F_LB=%10.4f F_L=%10.4f V_L=%10.4f", lambda,
                             lambdaBin, dForceFieldEnergydL, dBdL, dUdLambda, histogram.ld.thetaVelocity));
                 }
             }
 
             // Metadynamics grid counts (every 'countInterval' steps).
             if (energyCount % histogram.hd.countInterval == 0) {
                 histogram.addBias(dForceFieldEnergydL);
 
                 // Locally optimize the current state.
                 if (optimizationParameters.doOptimization) {
                     optimizationParameters.optimize(forceFieldEnergy, x, gradient);
                 }
 
                 if (algorithmListener != null) {
                     algorithmListener.algorithmUpdate(molecularAssembly);
                 }
             }
         }
 
         totalEnergy = forceFieldEnergy + biasEnergy;
 
         return totalEnergy;
     }
 
     /**
      * {@inheritDoc}
      */
     @Override
     public double[] getAcceleration(double[] acceleration) {
         return potential.getAcceleration(acceleration);
     }
 
     public Histogram[] getAllHistograms() {
         int nHisto = allHistograms.size();
         Histogram[] ret = new Histogram[nHisto];
         ret = allHistograms.toArray(ret);
         return ret;
     }
 
     /**
      * {@inheritDoc}
      */
     @Override
     public double[] getCoordinates(double[] doubles) {
         return potential.getCoordinates(doubles);
     }
 
     /**
      * {@inheritDoc}
      */
     @Override
     public void setCoordinates(double[] doubles) {
         potential.setCoordinates(doubles);
     }
 
     /**
      * {@inheritDoc}
      */
     @Override
     public Crystal getCrystal() {
         return potential.getCrystal();
     }
 
     /**
      * {@inheritDoc}
      */
     @Override
     public void setCrystal(Crystal crystal) {
         potential.setCrystal(crystal);
     }
 
     /**
      * Returns the number of energy evaluations performed by this OST, including those picked up in the
      * lambda file.
      *
      * @return Number of energy steps taken by this walker.
      */
     public long getEnergyCount() {
         return histogram.ld.stepsTaken;
     }
 
     /**
      * {@inheritDoc}
      */
     @Override
     public Potential.STATE getEnergyTermState() {
         return state;
     }
 
     /**
      * {@inheritDoc}
      */
     @Override
     public void setEnergyTermState(Potential.STATE state) {
         this.state = state;
         potential.setEnergyTermState(state);
     }
 
     /**
      * Getter for the field <code>forceFieldEnergy</code>.
      *
      * @return a double.
      */
     public double getForceFieldEnergy() {
         return forceFieldEnergy;
     }
 
     /**
      * Return the current 2D Histogram of counts.
      *
      * @return the Histogram.
      */
     public Histogram getHistogram() {
         return histogram;
     }
 
     /**
      * Getter for the field <code>lambda</code>.
      *
      * @return a double.
      */
     public double getLambda() {
         return histogram.ld.lambda;
     }
 
     /**
      * Setter for the field <code>lambda</code>.
      *
      * @param lambda a double.
      */
     public void setLambda(double lambda) {
         if (histogram != null) {
             lambda = histogram.mapLambda(lambda);
         } else {
             logger.warning(" OrthogonalSpaceTempering.setLambda was called before histogram constructed!");
             logger.info(Utilities.stackTraceToString(new RuntimeException()));
         }
         lambdaInterface.setLambda(lambda);
         histogram.ld.lambda = lambda;
         histogram.ld.theta = asin(sqrt(lambda));
     }
 
     /**
      * {@inheritDoc}
      */
     @Override
     public double[] getMass() {
         return potential.getMass();
     }
 
     /**
      * {@inheritDoc}
      */
     @Override
     public int getNumberOfVariables() {
         return potential.getNumberOfVariables();
     }
 
     /**
      * Return the OST optimization information.
      *
      * @return The OST optimization parameters.
      */
     public OptimizationParameters getOptimizationParameters() {
         return optimizationParameters;
     }
 
     /**
      * getPotentialEnergy.
      *
      * @return a {@link ffx.numerics.Potential} object.
      */
     public Potential getPotentialEnergy() {
         return potential;
     }
 
     /**
      * {@inheritDoc}
      */
     @Override
     public double[] getPreviousAcceleration(double[] previousAcceleration) {
         return potential.getPreviousAcceleration(previousAcceleration);
     }
 
     /**
      * {@inheritDoc}
      */
     @Override
     public double[] getScaling() {
         return potential.getScaling();
     }
 
     /**
      * {@inheritDoc}
      */
     @Override
     public void setScaling(double[] scaling) {
         potential.setScaling(scaling);
     }
 
     /**
      * {@inheritDoc}
      */
     @Override
     public double getTotalEnergy() {
         return totalEnergy;
     }
 
     /**
      * getTotaldEdLambda.
      *
      * @return a double.
      */
     public double getTotaldEdLambda() {
         return dUdLambda;
     }
 
     /**
      * {@inheritDoc}
      *
      * <p>Return a reference to each variables type.
      */
     @Override
     public Potential.VARIABLE_TYPE[] getVariableTypes() {
         return potential.getVariableTypes();
     }
 
     /**
      * {@inheritDoc}
      */
     @Override
     public double[] getVelocity(double[] velocity) {
         return potential.getVelocity(velocity);
     }
 
     /**
      * {@inheritDoc}
      */
     @Override
     public double getd2EdL2() {
         throw new UnsupportedOperationException(
                 " Second derivatives of the bias are not implemented, as they require third derivatives of the potential.");
     }
 
     /**
      * {@inheritDoc}
      */
     @Override
     public double getdEdL() {
         return getTotaldEdLambda();
     }
 
     /**
      * {@inheritDoc}
      */
     @Override
     public void getdEdXdL(double[] gradient) {
         throw new UnsupportedOperationException(
                 " Second derivatives of the bias are not implemented, as they require third derivatives of the potential.");
     }
 
     // TODO: Delete method when debugging of RepexOST is done.
     public void logOutputFiles(int index) {
         logger.info(format(" OST: Lambda file %s, histogram %s", histogram.ld.getLambdaFile(),
                 allHistograms.get(index).hd.getHistogramFile()));
     }
 
     /**
      * {@inheritDoc}
      */
     @Override
     public void setAcceleration(double[] acceleration) {
         potential.setAcceleration(acceleration);
     }
 
     /**
      * If this flag is true, (lambda, dU/dL) Monte Carlo samples that have no weight in the Histogram
      * are rejected.
      *
      * @param hardWallConstraint If true, MC samples outside the current range are rejected.
      */
     public void setHardWallConstraint(boolean hardWallConstraint) {
         this.hardWallConstraint = hardWallConstraint;
     }
 
     public void setMolecularAssembly(MolecularAssembly molecularAssembly) {
         this.molecularAssembly = molecularAssembly;
     }
 
     /**
      * {@inheritDoc}
      */
     @Override
     public void setPreviousAcceleration(double[] previousAcceleration) {
         potential.setPreviousAcceleration(previousAcceleration);
     }
 
     /**
      * Indicate if the Lambda extended system particle should be propagated using Langevin dynamics.
      *
      * @param propagateLambda If true, Lambda will be propagated using Langevin dynamics.
      */
     public void setPropagateLambda(boolean propagateLambda) {
         this.propagateLambda = propagateLambda;
     }
 
     /**
      * If true, the Lambda extended system particle is propagated using Langevin dynamics.
      */
     public boolean getPropagateLambda() {
         return propagateLambda;
     }
 
     /**
      * {@inheritDoc}
      */
     @Override
     public void setVelocity(double[] velocity) {
         potential.setVelocity(velocity);
     }
 
     /**
      * Switch to an alternate Histogram.
      *
      * @param index Index of the Histogram to use.
      */
     public void switchHistogram(int index) {
         histogramIndex = index;
         histogram = allHistograms.get(histogramIndex);
         logger.info(" OST switching to histogram " + histogramIndex);
     }
 
     @Override
     public void writeAdditionalRestartInfo(boolean recursive) {
         histogram.writeRestart();
         if (recursive) {
             potential.writeAdditionalRestartInfo(recursive);
         }
     }
 
     double getLambdaWriteOut() {
         return histogram.hd.resetHistogramAtLambda;
     }
 
     /**
      * getForceFielddEdL.
      *
      * @return a double.
      */
     double getForceFielddEdL() {
         return dForceFieldEnergydL;
     }
 
     /**
      * Getter for the field <code>biasEnergy</code>.
      *
      * @return a double.
      */
     double getBiasEnergy() {
         return biasEnergy;
     }
 
     /**
      * If the dUdLHardWall flag is set to true, this method will return false if the (lambda, dU/dL)
      * sample is has not been seen.
      *
      * @param lambda The proposed lambda value.
      * @param dUdL   The proposed dU/dL value.
      * @return Returns false only if the dUdLHardWall flag is true, and the (lambda, dU/dL) sample has
      * not been seen.
      */
     boolean insideHardWallConstraint(double lambda, double dUdL) {
         if (hardWallConstraint) {
             double weight = histogram.getRecursionKernelValue(lambda, dUdL);
             return weight > 0.0;
         }
         return true;
     }
 
     /**
      * Parameters for running local optimizations during OST sampling.
      */
     public class OptimizationParameters {
 
         /**
          * Flag to turn on OST optimization.
          *
          * <p>The default doOptimization = false.
          */
         private boolean doOptimization = false;
         /**
          * Reset unit cell parameters, molecular orientation and translation.
          */
         private final boolean doUnitCellReset;
         /**
          * OST optimization only runs if Lambda is greater than the lambdaCutoff.
          *
          * <p>The default lambdaCutoff = 0.8.
          */
         private final double lambdaCutoff;
         /**
          * The lowest energy found via optimizations.
          *
          * <p>The optimumEnergy is initially set to Double.MAX_VALUE.
          */
         private double optimumEnergy = Double.MAX_VALUE;
         /**
          * The OST optimization frequency
          *
          * <p>The default is once every 10,000 steps.
          */
         private final int frequency;
         /**
          * The OST optimization convergence criteria.
          *
          * <p>The default eps = 0.1.
          */
         private final double eps;
         /**
          * The OST tolerance when checking for equal energy after coordinate reversion.
          *
          * <p>The default is 1.0e-8 kcal/mol.
          */
         private final double tolerance;
         /**
          * The OST optimization energy window.
          *
          * <p>The default is 4.0 kcal/mol, which is convenient for small organic crystals.
          */
         private final double energyWindow;
         /**
          * Holds the lowest potential energy coordinates.
          */
         private double[] optimumCoords;
         /**
          * File instance used for saving optimized structures.
          */
         private File optimizationFile;
         /**
          * SystemFilter used to save optimized structures.
          */
         private SystemFilter optimizationFilter;
 
         /**
          * Empty constructor.
          */
         OptimizationParameters(CompositeConfiguration properties) {
             energyWindow = properties.getDouble("ost-opt-energy-window", 10.0);
             eps = properties.getDouble("ost-opt-eps", 0.1);
             tolerance = properties.getDouble("ost-opt-tolerance", 1.0e-8);
             frequency = properties.getInt("ost-opt-frequency", 10000);
             lambdaCutoff = properties.getDouble("ost-opt-lambda-cutoff", 0.8);
             doUnitCellReset = properties.getBoolean("ost-opt-unitcell-reset", false);
         }
 
         private void log() {
             logger.info("\n Optimization Parameters");
             logger.info(format("  Energy Window:                  %6.3f (kcal/mol)", energyWindow));
             logger.info(format("  EPS:                            %6.4f RMS (kcal/mol/Å)", eps));
             logger.info(format("  Tolerance:                      %6.4f (kcal/mol)", tolerance));
             logger.info(format("  Frequency:                      %6d (steps)", frequency));
             logger.info(format("  Lambda Cutoff:                  %6.4f", lambdaCutoff));
             logger.info(format("  Unit Cell Reset:                %6B", doUnitCellReset));
             logger.info(format("  File:                           %s", optimizationFile.getName()));
         }
 
         /**
          * The coordinates of the lowest energy structure found during optimization.
          *
          * @return the coordinates of the lowest energy structure found during optimization.
          */
         public double[] getOptimumCoordinates() {
             if (optimumEnergy < Double.MAX_VALUE) {
                 return optimumCoords;
             } else {
                 logger.info("Lambda optimization cutoff was not reached. Try increasing the number of timesteps.");
                 return null;
             }
         }
 
         /**
          * getOptimumEnergy.
          *
          * @return a double.
          */
         public double getOptimumEnergy() {
             if (optimumEnergy == Double.MAX_VALUE) {
                 logger.info(
                         "Lambda optimization cutoff was not reached. Try increasing the number of timesteps.");
             }
             return optimumEnergy;
         }
 
         /**
          * Run a local optimization.
          *
          * @param e        Current energy.
          * @param x        Current atomic coordinates.
          * @param gradient Work array for collecting the gradient.
          */
         public void optimize(double e, double[] x, @Nullable double[] gradient) {
 
             // Return if the optimization flag is not set, or if lambda is not beyond the cutoff.
             double lambda = histogram.ld.lambda;
             long stepsTaken = histogram.ld.stepsTaken;
             if (doOptimization && lambda > lambdaCutoff && stepsTaken % frequency == 0) {
                 if (gradient == null) {
                     gradient = new double[x.length];
                 }
             } else {
                 return;
             }
 
             logger.info(format("\n OST Minimization (Step %d)", stepsTaken));
 
             // Set the underlying Potential's Lambda value to 1.0.
             lambdaInterface.setLambda(1.0);
 
             // Use all energy terms.
             potential.setEnergyTermState(Potential.STATE.BOTH);
 
             // Optimize the system.
             try {
                 double startingEnergy = potential.energy(x);
                 Minimize minimize = new Minimize(molecularAssembly, potential, algorithmListener);
                 minimize.minimize(eps);
                 // Collect the minimum energy.
                 double minEnergy = potential.getTotalEnergy();
                 // Check for a new minimum within an energy window of the lowest energy structure found.
                 if (minEnergy < optimumEnergy + energyWindow) {
                     if (minEnergy < optimumEnergy) {
                         optimumEnergy = minEnergy;
                     }
                     int n = potential.getNumberOfVariables();
                     optimumCoords = new double[n];
                     optimumCoords = potential.getCoordinates(optimumCoords);
                     double mass = molecularAssembly.getMass();
                     Crystal crystal = molecularAssembly.getCrystal();
                     double density = crystal.getDensity(mass);
                     optimizationFilter.writeFile(optimizationFile, true);
                     Crystal uc = crystal.getUnitCell();
                     logger.info(format(" Minimum: %12.6f %s (%12.6f g/cc) optimized from %12.6f at step %d.",
                             minEnergy, uc.toShortString(), density, startingEnergy, stepsTaken));
                 }
             } catch (Exception ex) {
                 String message = ex.getMessage();
                 logger.info(
                         format(" Exception minimizing coordinates at lambda=%8.6f\n %s.", lambda, message));
                 logger.info(" Sampling will continue.");
             }
 
             // Set the underlying Potential's Lambda value back to current lambda value.
             lambdaInterface.setLambda(lambda);
 
             // Reset the Potential State
             potential.setEnergyTermState(state);
 
             if (doUnitCellReset) {
                 logger.info("\n Resetting Unit Cell");
                 double mass = molecularAssembly.getMass();
                 double density = molecularAssembly.getCrystal().getDensity(mass);
                 molecularAssembly.applyRandomDensity(density);
                 molecularAssembly.applyRandomSymOp(0.0);
                 lambda = 0.0;
                 lambdaInterface.setLambda(lambda);
             } else {
                 // Revert to the coordinates and gradient prior to optimization.
                 double eCheck = potential.energyAndGradient(x, gradient);
 
                 if (abs(eCheck - e) > tolerance) {
                     logger.warning(format(" Optimization could not revert coordinates %16.8f vs. %16.8f.", e, eCheck));
                 }
             }
         }
 
         /**
          * setOptimization.
          *
          * @param doOptimization    a boolean.
          * @param molecularAssembly a {@link ffx.potential.MolecularAssembly} object.
          */
         public void setOptimization(boolean doOptimization, MolecularAssembly molecularAssembly) {
             this.doOptimization = doOptimization;
             OrthogonalSpaceTempering.this.molecularAssembly = molecularAssembly;
             File file = molecularAssembly.getFile();
             String fileName = FilenameUtils.removeExtension(file.getAbsolutePath());
             String ext = FilenameUtils.getExtension(file.getAbsolutePath());
             if (optimizationFilter == null) {
                 if (ext.toUpperCase().contains("XYZ")) {
                     optimizationFile = new File(fileName + "_opt.arc");
                     optimizationFilter = new XYZFilter(optimizationFile, molecularAssembly,
                             molecularAssembly.getForceField(), molecularAssembly.getProperties());
                 } else {
                     optimizationFile = new File(fileName + "_opt.pdb");
                     PDBFilter pdbFilter = new PDBFilter(optimizationFile, molecularAssembly,
                             molecularAssembly.getForceField(), molecularAssembly.getProperties());
                     int models = pdbFilter.countNumModels();
                     pdbFilter.setModelNumbering(models);
                     optimizationFilter = pdbFilter;
                 }
             }
 
             if (this.doOptimization) {
                 log();
             }
         }
     }
 
     /**
      * Store and operate on a 2D Histogram of (Lambda, dU/dL) observations to produce an OST bias.
      *
      * @author Michael J. Schnieders
      * @since 1.0
      */
     public class Histogram {
 
         /**
          * If a 2D bias is in use (i.e. lambda-bias-cutoff > 0), do not normalize bias height.
          */
         private static final double TWO_D_NORMALIZATION = 1.0;
         /**
          * If a 1D bias is in use (i.e. lambda-bias-cutoff == 0), normalize bias height to deposit the same
          * volume of bias.
          */
         private static final double ONE_D_NORMALIZATION = Math.sqrt(2.0 * Math.PI);
         /**
          * If the real bias magnitude is 0 (for no 2D bias), temporarily set biasMag to this value to calculate the ensemble average dU/dL.
          *
          * <p>Any value that does not overflow/underflow double precision summations of the count matrix
          * will give identical results. For example, values of 1.0e-20, 1.0 and 1.0e20 were tested and
          * found to give identical results.
          *
          * <p>Thus, a value of 1.0 is good choice.
          */
         private static final double PSEUDO_BIAS_MAGNITUDE = 1.0;
         /**
          * Parallel Java world communicator.
          */
         protected final Comm world;
         /**
          * Rank of this process.
          */
         protected final int rank;
         /**
          * 1D PMF with respect to lambda F(L).
          */
         final double[] ensembleAveragedUdL;
         /**
          * This deltaT is used to determine the tempering weight as described below for the
          * temperingWeight variable.
          *
          * <p>deltaT = temperingFactor * kB * T.
          */
         private final double deltaT;
         /**
          * The integration algorithm used for thermodynamic integration.
          */
         private final IntegrationType integrationType;
         /**
          * Send OST counts asynchronously.
          */
         private final SendAsynchronous sendAsynchronous;
         /**
          * Send OST counts synchronously.
          */
         private final SendSynchronous sendSynchronous;
         /**
          * The Dama et al. transition-tempering weight. The initial temperingWeight = 1.0,
          * and more generally temperingWeight = exp(-biasHeight/deltaT)
          */
         private double temperingWeight = 1.0;
         /**
          * If the recursion kernel becomes too large or too small for some combinations of (Lambda,
          * dU/dL), then its statistical weight = exp(kernel * beta) cannot be represented by a double
          * value.
          */
         private double[] kernelValues;
         /**
          * Most recent lambda values for each Walker.
          */
         private double lastReceivedLambda;
         /**
          * Most recent dU/dL value for each walker.
          */
         private double lastReceiveddUdL;
         private final boolean spreadBias;
         /**
          * The HistogramData object contains the parameters for the Histogram that are saved to disk.
          */
         final HistogramData hd;
         /**
          * The LambdaData object contains the parameters for the Lambda particle that are saved to disk.
          */
         final LambdaData ld;
         /**
          * The Stochastic integrator used to propagate the lambda particle.
          */
         private final Stochastic stochastic;
         /**
          * The state of the lambda particle.
          */
         private final SystemState lambdaState;
         /**
          * Number of hills when the last free energy estimate was made.
          */
         private int currentNumberOfHills = 0;
         /**
          * The current free energy difference.
          */
         private double currentFreeEnergyDifference = 0.0;
 
         /**
          * Histogram constructor.
          *
          * @param properties    a CompositeConfiguration used to configure the Histogram.
          * @param histogramData An object containing the values this Histogram will be set to.
          */
         Histogram(CompositeConfiguration properties, HistogramData histogramData, LambdaData lambdaData) {
             hd = histogramData;
             ld = lambdaData;
 
       /*
       double gaussNormalization;
       if (hd.lambdaBiasCutoff == 0 && !hd.wasHistogramRead()) {
         gaussNormalization = ONE_D_NORMALIZATION;
         logger.info(format(" Bias magnitude multiplied by a factor of %.4f "
             + "sqrt(2*pi) to match 1D Gaussian volume to 2D Gaussian volume.", gaussNormalization));
       } else {
         gaussNormalization = TWO_D_NORMALIZATION;
       }
       biasMag = settings.getBiasMag() * gaussNormalization;
        */
 
             deltaT = hd.temperingFactor * R * dynamicsOptions.getTemperature();
 
             // Allocate space to compute the <dU/dL>
             ensembleAveragedUdL = new double[hd.getLambdaBins()];
 
             // Allocate space to regularize kernel values.
             kernelValues = new double[hd.getDUDLBins()];
 
             String propString = properties.getString("ost-integrationType", "SIMPSONS");
             IntegrationType testType;
             try {
                 testType = IntegrationType.valueOf(propString.toUpperCase());
             } catch (Exception ex) {
                 logger.warning(format(" Invalid argument %s to ost-integrationType; resetting to SIMPSONS", propString));
                 testType = SIMPSONS;
             }
             integrationType = testType;
 
             spreadBias = properties.getBoolean("ost-spread-bias", false);
 
       /*
        Set up the multi-walker communication variables for Parallel Java
        communication between nodes.
       */
             world = Comm.world();
             int numProc = world.size();
             rank = world.rank();
             if (hd.asynchronous) {
                 // Use asynchronous communication.
                 sendAsynchronous = new SendAsynchronous(this);
                 sendAsynchronous.start();
                 sendSynchronous = null;
             } else {
                 Histogram[] histograms = new Histogram[numProc];
                 int[] rankToHistogramMap = new int[numProc];
                 for (int i = 0; i < numProc; i++) {
                     histograms[i] = this;
                     rankToHistogramMap[i] = 0;
                 }
                 sendSynchronous = new SendSynchronous(histograms, rankToHistogramMap);
                 sendAsynchronous = null;
             }
 
             lastReceivedLambda = ld.lambda;
             if (hd.discreteLambda) {
                 lastReceivedLambda = mapLambda(lastReceivedLambda);
                 logger.info(format(" Discrete lambda: initializing lambda to nearest bin %.5f", lastReceivedLambda));
                 ld.lambda = lastReceivedLambda;
                 ld.theta = asin(sqrt(lastReceivedLambda));
                 lambdaInterface.setLambda(lastReceivedLambda);
                 lambdaState = null;
                 stochastic = null;
             } else {
                 // Configure the Stochastic integrator.
                 lambdaState = new SystemState(1);
                 stochastic = new Stochastic(lambdaParticleOptions.getLambdaFriction(), lambdaState);
                 stochastic.setTemperature(dynamicsOptions.getTemperature());
                 stochastic.setTimeStep(dynamicsOptions.getDtPsec());
                 double[] mass = lambdaState.getMass();
                 double[] thetaPosition = lambdaState.x();
                 double[] thetaVelocity = lambdaState.v();
                 double[] thetaAccel = lambdaState.a();
                 mass[0] = lambdaParticleOptions.getLambdaMass();
                 thetaPosition[0] = ld.theta;
                 thetaVelocity[0] = ld.thetaVelocity;
                 thetaAccel[0] = ld.thetaAcceleration;
             }
 
             lastReceiveddUdL = getdEdL();
 
             if (hd.wasHistogramRead()) {
                 updateFreeEnergyDifference(true, false);
             }
         }
 
         public void disableResetStatistics() {
             hd.resetHistogram = false;
         }
 
         /**
          * evaluateTotalBias.
          *
          * @param bias If false, return the negative of the Total OST bias.
          * @return A StringBuffer The total OST Bias.
          */
         public StringBuffer evaluateTotalOSTBias(boolean bias) {
             StringBuffer sb = new StringBuffer();
             double[] chainRule = new double[2];
             for (int dUdLBin = 0; dUdLBin < hd.getDUDLBins(); dUdLBin++) {
                 double currentdUdL = dUdLforBin(dUdLBin);
                 sb.append(format(" %16.8f", currentdUdL));
                 for (int lambdaBin = 0; lambdaBin < hd.getLambdaBins(); lambdaBin++) {
                     double currentLambda = lambdaForIndex(lambdaBin);
                     double bias1D = energyAndGradient1D(currentLambda, false);
                     double bias2D = energyAndGradient2D(currentLambda, currentdUdL, chainRule, hd.biasMag);
                     double totalBias = bias1D + bias2D;
                     // If the bias flag is false, turn the total bias into the PMF.
                     if (!bias) {
                         totalBias = -totalBias;
                     }
                     sb.append(format(" %16.8f", totalBias));
                 }
                 sb.append("\n");
             }
             return sb;
         }
 
         /**
          * evaluate2DOSTBias.
          *
          * @param bias If false, return the negative of the 2D OST bias.
          * @return A StringBuffer with the 2D OST bias.
          */
         public StringBuffer evaluate2DOSTBias(boolean bias) {
             StringBuffer sb = new StringBuffer();
             double[] chainRule = new double[2];
             for (int dUdLBin = 0; dUdLBin < hd.getDUDLBins(); dUdLBin++) {
                 double currentdUdL = dUdLforBin(dUdLBin);
                 sb.append(format(" %16.8f", currentdUdL));
                 for (int lambdaBin = 0; lambdaBin < hd.getLambdaBins(); lambdaBin++) {
                     double currentLambda = lambdaForIndex(lambdaBin);
                     double bias2D = energyAndGradient2D(currentLambda, currentdUdL, chainRule, hd.biasMag);
                     if (!bias) {
                         bias2D = -bias2D;
                     }
                     sb.append(format(" %16.8f", bias2D));
                 }
                 sb.append("\n");
             }
             return sb;
         }
 
         /**
          * For MPI parallel jobs, returns true if the walkers are independent (i.e. contribute to only
          * their own histogram).
          *
          * @return True if the walkers are independent.
          */
         public boolean getIndependentWalkers() {
             return hd.independentWalkers;
         }
 
         public double getLambdaResetValue() {
             return hd.resetHistogramAtLambda;
         }
 
         /**
          * For MPI parallel jobs, return the rank of this process.
          *
          * @return The rank of this process.
          */
         public int getRank() {
             return rank;
         }
 
         public boolean getResetStatistics() {
             return hd.resetHistogram;
         }
 
         /**
          * Return the SynchronousSend associated with this Histogram, if any.
          *
          * @return The SynchronousSend, if any.
          */
         public Optional<SendSynchronous> getSynchronousSend() {
             return Optional.ofNullable(sendSynchronous);
         }
 
         public void setLastReceiveddUdL(double lastReceiveddUdL) {
             this.lastReceiveddUdL = lastReceiveddUdL;
         }
 
         public double updateFreeEnergyDifference(boolean print, boolean save) {
             if (hd.metaDynamics) {
                 return updateMetaDynamicsFreeEnergyDifference(print, save);
             } else {
                 return updateOSTFreeEnergyDifference(print, save);
             }
         }
 
         /**
          * Update the free energy estimate for Meta Dynamics.
          *
          * @param print Whether to write the histogram to screen.
          * @param save  Whether to write the histogram to disk.
          * @return Free energy (via integration of ensemble-average dU/dL)
          */
         public double updateMetaDynamicsFreeEnergyDifference(boolean print, boolean save) {
             synchronized (this) {
                 if (!print && currentNumberOfHills == hd.counts) {
                     return currentFreeEnergyDifference;
                 }
                 double freeEnergyDifferenceMeta = 0.0;
                 double freeEnergyDifferenceTI = 0.0;
 
                 double minFL = Double.MAX_VALUE;
 
                 int lambdaBins = hd.getLambdaBins();
                 double[] metaFreeEnergy = new double[lambdaBins];
 
                 // Total histogram weight.
                 double totalWeight = 0;
                 StringBuilder stringBuilder = new StringBuilder();
 
                 // Loop over lambda bins, computing <dU/dL> for each bin.
                 for (int lambdaBin = 0; lambdaBin < lambdaBins; lambdaBin++) {
                     int firstdUdLBin = firstdUdLBin(lambdaBin);
                     int lastdUdLBin = lastdUdLBin(lambdaBin);
                     double lambdaCount = 0;
                     // The dUdL range sampled for lambda.
                     double mindUdLForLambda = 0.0;
                     double maxdUdLforLambda = 0.0;
                     double maxBias = 0;
                     if (firstdUdLBin == -1 || lastdUdLBin == -1) {
                         ensembleAveragedUdL[lambdaBin] = 0.0;
                         minFL = 0.0;
                     } else {
                         double ensembleAverageFLambda = 0.0;
                         double partitionFunction = 0.0;
 
                         // Evaluate and regularize all kernel values for this value of lambda.
                         double offset = evaluateKernelForLambda(lambdaBin, firstdUdLBin, lastdUdLBin);
 
                         for (int dUdLBin = firstdUdLBin; dUdLBin <= lastdUdLBin; dUdLBin++) {
                             double kernel = kernelValues[dUdLBin];
 
                             if (kernel - offset > maxBias) {
                                 maxBias = kernel - offset;
                             }
 
                             // The weight is just the kernel value for no 2D bias.
                             partitionFunction += kernel;
 
                             double currentdUdL = dUdLforBin(dUdLBin);
                             ensembleAverageFLambda += currentdUdL * kernel;
                             lambdaCount += getRecursionKernelValue(lambdaBin, dUdLBin);
                         }
                         if (minFL > maxBias) {
                             minFL = maxBias;
                         }
                         ensembleAveragedUdL[lambdaBin] = (partitionFunction == 0) ? 0 : ensembleAverageFLambda / partitionFunction;
                         mindUdLForLambda = hd.dUdLMinimum + firstdUdLBin * hd.dUdLBinWidth;
                         maxdUdLforLambda = hd.dUdLMinimum + (lastdUdLBin + 1) * hd.dUdLBinWidth;
                     }
 
                     double deltaFreeEnergy = ensembleAveragedUdL[lambdaBin] * deltaForLambdaBin(lambdaBin);
                     freeEnergyDifferenceTI += deltaFreeEnergy;
                     totalWeight += lambdaCount;
 
 
                     double lambdaBinWidth = hd.lambdaBinWidth;
                     double llL = lambdaBin * lambdaBinWidth - hd.lambdaBinWidth_2;
                     double ulL = llL + lambdaBinWidth;
                     if (llL < 0.0) {
                         llL = 0.0;
                     }
                     if (ulL > 1.0) {
                         ulL = 1.0;
                     }
 
                     double midLambda = llL;
                     if (!hd.discreteLambda) {
                         midLambda = (llL + ulL) / 2.0;
                     }
                     double bias1D = energyAndGradient1D(midLambda, false);
                     metaFreeEnergy[lambdaBin] = energyAndGradientMeta(midLambda, false);
                     freeEnergyDifferenceMeta = -(metaFreeEnergy[lambdaBin] - metaFreeEnergy[0]);
 
                     if (print || save) {
                         stringBuilder.append(format(" %6.2e %7.5f %7.1f %7.1f %9.2f %9.2f %9.2f %9.2f %9.2f\n", lambdaCount,
                                 midLambda, mindUdLForLambda, maxdUdLforLambda, ensembleAveragedUdL[lambdaBin],
                                 bias1D, freeEnergyDifferenceTI, metaFreeEnergy[lambdaBin], freeEnergyDifferenceMeta));
                     }
                 }
 
                 double temperingOffset = hd.getTemperingOffset();
                 double temperEnergy = (minFL > temperingOffset) ? temperingOffset - minFL : 0;
                 temperingWeight = exp(temperEnergy / deltaT);
 
                 if (print) {
                     logger.info("  Weight   Lambda      dU/dL Bins   <dU/dL>      g(L) dG_OST(L)  Meta(L) dG_Meta(L)");
                     logger.info(stringBuilder.toString());
                     logger.info(" Histogram Evaluation");
                     logger.info(format("  Free Energy Difference: %12.4f kcal/mol", freeEnergyDifferenceMeta));
                     logger.info(format("  Number of Hills:        %12d", hd.counts));
                     logger.info(format("  Total Bias Added:       %12.4f kcal/mol", totalWeight * hd.biasMag));
                     logger.info(format("  Minimum Bias:           %12.4f kcal/mol", minFL));
                     logger.info(format("  Tempering Percentage:   %12.4f %%\n", temperingWeight * 100));
                 }
 
                 if (save) {
                     String modelFilename = molecularAssembly.getFile().getAbsolutePath();
                     File saveDir = new File(FilenameUtils.getFullPath(modelFilename));
                     String dirName = saveDir + File.separator;
                     String fileName = dirName + "histogram.txt";
                     try {
                         logger.info(" Writing " + fileName);
                         PrintWriter printWriter = new PrintWriter(fileName);
                         printWriter.write(stringBuilder.toString());
                         printWriter.close();
                     } catch (Exception e) {
                         logger.info(format(" Failed to write %s.", fileName));
                     }
                 }
 
                 currentNumberOfHills = hd.counts;
                 currentFreeEnergyDifference = freeEnergyDifferenceMeta;
                 return currentFreeEnergyDifference;
             }
         }
 
         /**
          * Eqs. 7 and 8 from the 2012 Crystal Thermodynamics paper.
          *
          * @param print Whether to write the histogram to screen.
          * @param save  Whether to write the histogram to disk.
          * @return Free energy (via integration of ensemble-average dU/dL)
          */
         public double updateOSTFreeEnergyDifference(boolean print, boolean save) {
             synchronized (this) {
                 if (!print && currentNumberOfHills == hd.counts) {
                     return currentFreeEnergyDifference;
                 }
 
                 double freeEnergyDifferenceOST = 0.0;
                 double minFL = Double.MAX_VALUE;
 
                 // If the bias magnitude is zero, computing <dU/dL> from
                 // counts will not be correct. Assign a temporary non-zero bias magnitude.
                 boolean biasMagZero = hd.biasMag <= 0.0;
 
                 // Total histogram weight.
                 double totalWeight = 0;
                 double beta = 1.0 / (R * dynamicsOptions.getTemperature());
                 StringBuilder stringBuilder = new StringBuilder();
 
                 // Loop over lambda bins, computing <dU/dL> for each bin.
                 for (int lambdaBin = 0; lambdaBin < hd.lambdaBins; lambdaBin++) {
                     int firstdUdLBin = firstdUdLBin(lambdaBin);
                     int lastdUdLBin = lastdUdLBin(lambdaBin);
                     double lambdaCount = 0;
                     // The FL range sampled for lambda bin [iL*dL .. (iL+1)*dL]
                     double mindUdLforLambda = 0.0;
                     double maxdUdLforLambda = 0.0;
                     double maxBias = 0;
                     if (firstdUdLBin == -1 || lastdUdLBin == -1) {
                         ensembleAveragedUdL[lambdaBin] = 0.0;
                         minFL = 0.0;
                     } else {
                         double ensembleAverage = 0.0;
                         double partitionFunction = 0.0;
 
                         // Evaluate and regularize all kernel values for this value of lambda.
                         double offset = evaluateKernelForLambda(lambdaBin, firstdUdLBin, lastdUdLBin);
 
                         for (int dUdLBin = firstdUdLBin; dUdLBin <= lastdUdLBin; dUdLBin++) {
                             double kernel = kernelValues[dUdLBin];
 
                             if (kernel - offset > maxBias) {
                                 maxBias = kernel - offset;
                             }
 
                             double weight;
                             if (biasMagZero) {
                                 // The weight is just the kernel value for no 2D bias.
                                 weight = kernel;
                             } else {
                                 // The Boltzmann weight based on temperature and the 2D bias height.
                                 weight = exp(kernel * beta);
                             }
                             partitionFunction += weight;
 
                             double currentdUdL = dUdLforBin(dUdLBin);
                             ensembleAverage += currentdUdL * weight;
                             lambdaCount += getRecursionKernelValue(lambdaBin, dUdLBin);
                         }
                         if (minFL > maxBias) {
                             minFL = maxBias;
                         }
                         ensembleAveragedUdL[lambdaBin] = (partitionFunction == 0) ? 0 : ensembleAverage / partitionFunction;
                         mindUdLforLambda = hd.dUdLMinimum + firstdUdLBin * hd.dUdLBinWidth;
                         maxdUdLforLambda = hd.dUdLMinimum + (lastdUdLBin + 1) * hd.dUdLBinWidth;
                     }
 
                     double deltaFreeEnergy = ensembleAveragedUdL[lambdaBin] * deltaForLambdaBin(lambdaBin);
                     freeEnergyDifferenceOST += deltaFreeEnergy;
                     totalWeight += lambdaCount;
 
                     if (print || save) {
                         double llL = lambdaBin * hd.lambdaBinWidth - hd.lambdaBinWidth_2;
                         double ulL = llL + hd.lambdaBinWidth;
                         if (llL < 0.0) {
                             llL = 0.0;
                         }
                         if (ulL > 1.0) {
                             ulL = 1.0;
                         }
 
                         double midLambda = llL;
                         if (!hd.discreteLambda) {
                             midLambda = (llL + ulL) / 2.0;
                         }
                         double bias1D = energyAndGradient1D(midLambda, false);
 
                         double bias2D = 0.0;
                         if (!biasMagZero) {
                             bias2D = computeBiasEnergy(midLambda, ensembleAveragedUdL[lambdaBin]) - bias1D;
                         }
 
                         stringBuilder.append(
                                 format(" %6.2e %7.5f %7.1f %7.1f %8.2f %8.2f %8.2f %8.2f %8.2f   %8.2f\n", lambdaCount,
                                         midLambda, mindUdLforLambda, maxdUdLforLambda, ensembleAveragedUdL[lambdaBin],
                                         bias1D, bias2D, bias1D + bias2D, freeEnergyDifferenceOST, bias1D + bias2D + freeEnergyDifferenceOST));
                     }
                 }
 
                 if (!biasMagZero) {
                     double temperingOffset = hd.getTemperingOffset();
                     double temperEnergy = (minFL > temperingOffset) ? temperingOffset - minFL : 0;
                     temperingWeight = exp(temperEnergy / deltaT);
                 }
 
                 freeEnergyDifferenceOST = integrateNumeric(ensembleAveragedUdL, integrationType);
 
                 if (print) {
                     logger.info("  Weight   Lambda      dU/dL Bins  <dU/dL>    g(L)  f(L,<dU/dL>) Bias    dG(L) Bias+dG(L)");
                     logger.info(stringBuilder.toString());
                     logger.info(" Histogram Evaluation");
                     logger.info(format("  Free Energy Difference: %12.4f kcal/mol", freeEnergyDifferenceOST));
                     logger.info(format("  Number of Hills:        %12d", hd.counts));
                     logger.info(format("  Total Bias Added:       %12.4f kcal/mol", totalWeight * hd.biasMag));
                     logger.info(format("  Minimum Bias:           %12.4f kcal/mol", minFL));
                     logger.info(format("  Tempering Percentage:   %12.4f %%\n", temperingWeight * 100));
                 }
 
                 if (save) {
                     String modelFilename = molecularAssembly.getFile().getAbsolutePath();
                     File saveDir = new File(FilenameUtils.getFullPath(modelFilename));
                     String dirName = saveDir + File.separator;
                     String fileName = dirName + "histogram.txt";
                     try {
                         logger.info("\n Writing " + fileName);
                         PrintWriter printWriter = new PrintWriter(fileName);
                         printWriter.write(stringBuilder.toString());
                         printWriter.close();
                     } catch (Exception e) {
                         logger.info(format(" Failed to write %s.", fileName));
                     }
                 }
 
                 currentNumberOfHills = hd.counts;
                 currentFreeEnergyDifference = freeEnergyDifferenceOST;
                 return currentFreeEnergyDifference;
             }
         }
 
         /**
          * For thermodynamic integration, return the integration width for the given Lambda lambdaBin.
          *
          * @param lambdaBin The lambda lambdaBin.
          * @return The integration width.
          */
         private double deltaForLambdaBin(int lambdaBin) {
             if (!hd.discreteLambda && (lambdaBin == 0 || lambdaBin == hd.lambdaBins - 1)) {
                 // The first and last lambda bins are half size for continuous lambda.
                 return hd.lambdaBinWidth_2;
             } else if (hd.discreteLambda && lambdaBin == 0) {
                 // The free energy change to move from L=0 to L=0 is zero.
                 return 0.0;
             }
             // All other cases.
             return hd.lambdaBinWidth;
         }
 
         /**
          * Returns the index of the first dU/dL bin with counts, or -1 if there are no counts for the
          * given lambda bin.
          *
          * @param lambdaBin Lambda bin to consider.
          * @return Index of the first dUdL bin with counts.
          */
         private int firstdUdLBin(int lambdaBin) {
             // Synchronize use of the recursion kernel.
             synchronized (this) {
                 // Find the smallest FL bin that has counts.
                 for (int jFL = 0; jFL < hd.dUdLBins; jFL++) {
                     double count = hd.zHistogram[lambdaBin][jFL];
                     if (count > 0) {
                         return jFL;
                     }
                 }
                 return -1;
             }
         }
 
         /**
          * Returns the index of the last dU/dL bin with counts, or -1 if there are no counts for the
          * given lambda bin.
          *
          * @param lambdaBin Lambda bin to consider.
          * @return Index of the last dUdL bin with counts.
          */
         private int lastdUdLBin(int lambdaBin) {
             // Synchronize use of the recursion kernel.
             synchronized (this) {
                 // Find the largest FL bin that has counts.
                 for (int jFL = hd.dUdLBins - 1; jFL >= 0; jFL--) {
                     double count = hd.zHistogram[lambdaBin][jFL];
                     if (count > 0) {
                         return jFL;
                     }
                 }
                 return -1;
             }
         }
 
         /**
          * Write histogram and lambda restart files.
          * <p>
          * For MPI parallel jobs, only rank 0 writes a histogram restart file (unless the "independentWrite" flag is set).
          */
         void writeRestart() {
             if (rank == 0 || hd.independentWrite) {
                 updateFreeEnergyDifference(true, false);
                 try {
                     hd.writeHistogram();
                     logger.info(format(" Wrote histogram restart to: %s.", hd.getHistogramFileName()));
                 } catch (Exception ex) {
                     String message = format(" Exception writing histogram restart file %s.", hd.getHistogramFileName());
                     logger.log(Level.INFO, Utilities.stackTraceToString(ex));
                     logger.log(Level.SEVERE, message, ex);
                 }
             }
 
             // All ranks write a lambda restart file.
             try {
                 ld.writeLambdaData();
                 logger.info(format(" Wrote lambda restart to:    %s.", ld.getLambdaFileName()));
             } catch (Exception ex) {
                 String message = format(" Exception writing lambda restart file %s.", ld.getLambdaFileName());
                 logger.log(Level.INFO, Utilities.stackTraceToString(ex));
                 logger.log(Level.SEVERE, message, ex);
             }
         }
 
         private double mapLambda(double lambda) {
             if (hd.discreteLambda) {
                 return hd.lambdaLadder[indexForDiscreteLambda(lambda)];
             } else {
                 return lambda;
             }
         }
 
         /**
          * For continuous lambda, the returned value is the lambda bin. For discrete lambda, the returned
          * value is the discrete lambda index.
          *
          * @param lambda a double.
          * @return a int.
          */
         int indexForLambda(double lambda) {
             if (hd.discreteLambda) {
                 return indexForDiscreteLambda(lambda);
             } else {
                 return indexForContinuousLambda(lambda);
             }
         }
 
         int indexForLambda() {
             return indexForLambda(ld.lambda);
         }
 
         private double lambdaForIndex(int bin) {
             if (hd.discreteLambda) {
                 return hd.lambdaLadder[bin];
             } else {
                 return bin * hd.lambdaBinWidth;
             }
         }
 
         private int indexForContinuousLambda(double lambda) {
             int lambdaBin = (int) floor((lambda - hd.minLambda) / hd.lambdaBinWidth);
             if (lambdaBin < 0) {
                 lambdaBin = 0;
             }
             if (lambdaBin >= hd.lambdaBins) {
                 lambdaBin = hd.lambdaBins - 1;
             }
             return lambdaBin;
         }
 
         private int indexForDiscreteLambda(double lambda) {
             assert hd.discreteLambda && hd.lambdaLadder != null && hd.lambdaLadder.length > 0;
 
             int initialGuess = indexForContinuousLambda(lambda);
             double minErr = Double.MAX_VALUE;
             int minErrBin = -1;
             for (int i = -1; i < 2; i++) {
                 int guessBin = i + initialGuess;
                 if (guessBin < 0 || guessBin >= hd.lambdaBins) {
                     continue;
                 }
                 double guessLam = hd.lambdaLadder[guessBin];
                 double guessErr = Math.abs(guessLam - lambda);
                 if (guessErr < minErr) {
                     minErr = guessErr;
                     minErrBin = guessBin;
                 }
             }
 
             assert minErr < 1.0E-6;
             return minErrBin;
         }
 
         /**
          * Find the bin for the supplied dEdLambda.
          * <p>
          * If the supplied dEdL is outside the range of the count matrix, then -1 is returned.
          *
          * @param dUdL a double.
          * @return The dUdL bin.
          */
         int binFordUdL(double dUdL) {
 
             // No counts below mindUdL.
             if (dUdL < hd.dUdLMinimum) {
                 return -1;
             }
 
             // No counts above the maxdUdL.
             if (dUdL > hd.dUdLMaximum) {
                 return -1;
             }
 
             int bin = (int) floor((dUdL - hd.dUdLMinimum) / hd.dUdLBinWidth);
 
             if (bin == hd.dUdLBins) {
                 bin = hd.dUdLBins - 1;
             }
 
             if (bin < 0) {
                 bin = 0;
             }
 
             return bin;
         }
 
         /**
          * Return the value of a recursion kernel bin. Mirror conditions are applied to the lambda bin.
          * If the dUdLBin is outside the range of the count matrix, zero is returned.
          *
          * @param lambdaBin The lambda bin.
          * @param dUdLBin   The dU/dL bin.
          * @return The value of the bin.
          */
         double getRecursionKernelValue(int lambdaBin, int dUdLBin) {
             // Synchronize use of the recursion kernel.
             synchronized (this) {
                 // Apply lambda mirror condition.
                 lambdaBin = lambdaMirror(lambdaBin);
 
                 // For dUdL outside the count matrix the weight is 0.
                 if (dUdLBin < 0 || dUdLBin >= hd.dUdLBins) {
                     return 0.0;
                 }
 
                 return hd.zHistogram[lambdaBin][dUdLBin];
             }
         }
 
         /**
          * Return the value of a recursion kernel bin. Mirror conditions are applied to the lambda bin.
          * If the dUdL is outside the range of the count matrix, zero is returned.
          *
          * @param lambda the lambda value.
          * @param dUdL   the dU/dL value.
          * @return The value of the Histogram.
          */
         double getRecursionKernelValue(double lambda, double dUdL) {
             return getRecursionKernelValue(indexForLambda(lambda), binFordUdL(dUdL));
         }
 
         /**
          * Add to the value of a recursion kernel bin.
          *
          * @param currentLambda The lambda bin.
          * @param currentdUdL   The dU/dL bin.
          * @param value         The value of the bin.
          */
         void addToRecursionKernelValue(double currentLambda, double currentdUdL, double value) {
             synchronized (this) {
                 if (spreadBias) {
                     // Expand the recursion kernel if necessary.
                     int dUdLbiasCutoff = hd.dUdLBiasCutoff;
                     checkRecursionKernelSize(dUdLforBin(binFordUdL(currentdUdL) - dUdLbiasCutoff));
                     checkRecursionKernelSize(dUdLforBin(binFordUdL(currentdUdL) + dUdLbiasCutoff));
 
                     int currentLambdaBin = indexForLambda(currentLambda);
                     int currentdUdLBin = binFordUdL(currentdUdL);
 
                     // Variances are only used when dividing by twice their value.
                     double invLs2 = 0.5 / hd.lambdaVariance;
                     double invFLs2 = 0.5 / hd.dUdLVariance;
 
                     // Compute the normalization.
                     double normalize = 0.0;
                     int lambdaBiasCutoff = hd.lambdaBiasCutoff;
                     for (int iL = -lambdaBiasCutoff; iL <= lambdaBiasCutoff; iL++) {
                         int lambdaBin = currentLambdaBin + iL;
                         double deltaL = currentLambda - lambdaBin * hd.lambdaBinWidth;
                         double deltaL2 = deltaL * deltaL;
                         // Pre-compute the lambda bias.
                         double L2exp = exp(-deltaL2 * invLs2);
                         for (int jFL = -dUdLbiasCutoff; jFL <= dUdLbiasCutoff; jFL++) {
                             int dUdLBin = currentdUdLBin + jFL;
                             double deltaFL = currentdUdL - dUdLforBin(dUdLBin);
                             double deltaFL2 = deltaFL * deltaFL;
                             normalize += L2exp * exp(-deltaFL2 * invFLs2);
                         }
                     }
 
                     for (int iL = -lambdaBiasCutoff; iL <= lambdaBiasCutoff; iL++) {
                         int lambdaBin = currentLambdaBin + iL;
                         double deltaL = currentLambda - lambdaBin * hd.lambdaBinWidth;
                         double deltaL2 = deltaL * deltaL;
                         // Pre-compute the lambda bias.
                         double L2exp = exp(-deltaL2 * invLs2);
                         lambdaBin = lambdaMirror(lambdaBin);
                         for (int jFL = -dUdLbiasCutoff; jFL <= dUdLbiasCutoff; jFL++) {
                             int dUdLBin = currentdUdLBin + jFL;
                             double deltaFL = currentdUdL - dUdLforBin(dUdLBin);
                             double deltaFL2 = deltaFL * deltaFL;
                             double weight = value / normalize * L2exp * exp(-deltaFL2 * invFLs2);
                             hd.zHistogram[lambdaBin][dUdLBin] += weight;
                         }
                     }
                     hd.counts++;
                 } else {
                     // Check the recursion kernel size.
                     checkRecursionKernelSize(currentdUdL);
                     int lambdaBin = indexForLambda(currentLambda);
                     int dUdLBin = binFordUdL(currentdUdL);
                     try {
                         hd.zHistogram[lambdaBin][dUdLBin] += value;
                         hd.counts++;
                     } catch (IndexOutOfBoundsException e) {
                         logger.info(format(" Histogram dimensions %d %d", hd.lambdaBins, hd.dUdLBins));
                         logger.info(format(
                                 " Count skipped in addToRecursionKernelValue due to an index out of bounds exception.\n L=%10.8f (%d), dU/dL=%10.8f (%d) and count=%10.8f",
                                 currentLambda, lambdaBin, currentdUdL, dUdLBin, value));
 
                     }
                 }
             }
         }
 
         /**
          * Allocate memory for the recursion kernel.
          */
         void allocateRecursionKernel() {
             // Synchronize updates of the recursion kernel.
             synchronized (this) {
                 hd.zHistogram = new double[hd.lambdaBins][hd.dUdLBins];
                 kernelValues = new double[hd.dUdLBins];
             }
         }
 
         /**
          * Integrates dUdL over lambda using more sophisticated techniques than midpoint rectangular
          * integration.
          *
          * <p>The ends (from 0 to dL and 1-dL to 1) are integrated with trapezoids for continuous
          * lambda.
          *
          * @param dUdLs dUdL at the midpoint of each bin.
          * @param type  Integration type to use.
          * @return Current delta-G estimate.
          */
         private double integrateNumeric(double[] dUdLs, IntegrationType type) {
             double val;
             if (hd.discreteLambda) {
                 // Integrate between the first bin and the last bin.
                 double[] lams = Integrate1DNumeric.generateXPoints(0.0, 1.0, hd.lambdaBins, false);
                 DataSet dSet = new DoublesDataSet(lams, dUdLs, false);
                 val = Integrate1DNumeric.integrateData(dSet, Integrate1DNumeric.IntegrationSide.LEFT, type);
             } else {
                 // Integrate between the second bin midpoint and the second-to-last bin midpoint.
                 double[] midLams = Integrate1DNumeric.generateXPoints(hd.lambdaBinWidth, 1.0 - hd.lambdaBinWidth,
                         (hd.lambdaBins - 2), false);
                 double[] midVals = Arrays.copyOfRange(dUdLs, 1, (hd.lambdaBins - 1));
                 DataSet dSet = new DoublesDataSet(midLams, midVals, false);
 
                 val = Integrate1DNumeric.integrateData(dSet, Integrate1DNumeric.IntegrationSide.LEFT, type);
 
                 // Everything after this is just adding in the endpoint contributions.
 
                 double dL_4 = hd.lambdaBinWidth_2 * 0.5;
 
                 // Initially, assume dU/dL is exactly 0 at the endpoints. This is sometimes a true
                 // assumption.
                 double val0 = 0;
                 double val1 = 0;
 
                 // If we cannot guarantee that dUdL is exactly 0 at the endpoints, interpolate.
                 if (!lambdaInterface.dEdLZeroAtEnds()) {
                     double recipSlopeLen = 1.0 / (hd.lambdaBinWidth * 0.75);
 
                     double slope = dUdLs[0] - dUdLs[1];
                     slope *= recipSlopeLen;
                     val0 = dUdLs[0] + (slope * dL_4);
 
                     slope = dUdLs[hd.lambdaBins - 1] - dUdLs[hd.lambdaBins - 2];
                     slope *= recipSlopeLen;
                     val1 = dUdLs[hd.lambdaBins - 1] + (slope * dL_4);
                     logger.fine(format(" Inferred dU/dL values at 0 and 1: %10.5g , %10.5g", val0, val1));
                 }
 
                 // Integrate trapezoids from 0 to the second bin midpoint, and from second-to-last bin
                 // midpoint to 1.
                 val += trapezoid(0, dL_4, val0, dUdLs[0]);
                 val += trapezoid(dL_4, hd.lambdaBinWidth, dUdLs[0], dUdLs[1]);
                 val += trapezoid(1.0 - hd.lambdaBinWidth, 1.0 - dL_4, dUdLs[hd.lambdaBins - 2],
                         dUdLs[hd.lambdaBins - 1]);
                 val += trapezoid(1.0 - dL_4, 1.0, dUdLs[hd.lambdaBins - 1], val1);
             }
 
             return val;
         }
 
         /**
          * Integrates a trapezoid.
          *
          * @param x0  First x point
          * @param x1  Second x point
          * @param fx0 First f(x) point
          * @param fx1 Second f(x) point
          * @return The area under a trapezoid.
          */
         private double trapezoid(double x0, double x1, double fx0, double fx1) {
             double val = 0.5 * (fx0 + fx1);
             val *= (x1 - x0);
             return val;
         }
 
         /**
          * Evaluate the kernel across dU/dL values at given value of lambda. The largest kernel value V
          * is used to define an offset (-V), which is added to all to kernel values. Then, the largest
          * kernel value is zero, and will result in a statistical weight of 1.
          *
          * <p>The offset avoids the recursion kernel becoming too large for some combinations of
          * (Lambda, dU/dL). This can result in its statistical weight = exp(kernel * beta) exceeding the
          * maximum representable double value.
          *
          * @param lambda Value of Lambda to evaluate the kernel for.
          * @param llFL   Lower FLambda bin.
          * @param ulFL   Upper FLambda bin.
          * @return the applied offset.
          */
         private double evaluateKernelForLambda(int lambda, int llFL, int ulFL) {
             double maxKernel = Double.MIN_VALUE;
 
             double gaussianBiasMagnitude = hd.biasMag;
             if (hd.biasMag <= 0.0) {
                 gaussianBiasMagnitude = PSEUDO_BIAS_MAGNITUDE;
             }
 
             for (int jFL = llFL; jFL <= ulFL; jFL++) {
                 double value = evaluateKernel(lambda, jFL, gaussianBiasMagnitude);
                 kernelValues[jFL] = value;
                 if (value > maxKernel) {
                     maxKernel = value;
                 }
             }
 
             // Only offset the kernel values for use with the 2D OST bias.
             double offset = 0.0;
             if (hd.biasMag > 0.0 && !hd.metaDynamics) {
                 offset = -maxKernel;
                 for (int jFL = llFL; jFL <= ulFL; jFL++) {
                     kernelValues[jFL] += offset;
                 }
             }
             return offset;
         }
 
         /**
          * Mirror the lambda bin if its less < 0 or greater than the last bin.
          *
          * @param bin Lambda bin to mirror.
          * @return The mirrored lambda bin.
          */
         private int lambdaMirror(int bin) {
             if (bin < 0) {
                 return -bin;
             }
             int lastBin = hd.lambdaBins - 1;
             if (bin > lastBin) {
                 // Number of bins past the last bin.
                 bin -= lastBin;
                 // Return Mirror bin
                 return lastBin - bin;
             }
             // No mirror condition.
             return bin;
         }
 
         /**
          * For continuous lambda, the width of the first and last bins is dLambda_2, so the mirror
          * condition is to double their counts.
          *
          * @param bin Current lambda bin.
          * @return The mirror factor (either 1.0 or 2.0).
          */
         private double mirrorFactor(int bin) {
             if (hd.discreteLambda) {
                 return 1.0;
             }
             if (bin == 0 || bin == hd.lambdaBins - 1) {
                 return 2.0;
             }
             return 1.0;
         }
 
         /**
          * Compute the value of dU/dL for the given Histogram bin.
          *
          * @param dUdLBin The bin index in the dU/dL dimension.
          * @return The value of dU/dL at the center of the bin.
          */
         private double dUdLforBin(int dUdLBin) {
             return hd.dUdLMinimum + dUdLBin * hd.dUdLBinWidth + hd.dUdLBinWidth_2;
         }
 
         /**
          * Evaluate the bias at [currentLambdaBin, cF_lambda]
          */
         private double evaluateKernel(int currentLambdaBin, int currentdUdLBin, double gaussianBiasMagnitude) {
             // Compute the value of L and FL for the center of the current bin.
             double currentLambda = currentLambdaBin * hd.lambdaBinWidth;
             double currentdUdL = dUdLforBin(currentdUdLBin);
 
             // Variances are only used when dividing by twice their value.
             double invLs2 = 0.5 / hd.lambdaVariance;
             double invFLs2 = 0.5 / hd.dUdLVariance;
 
             double sum = 0.0;
             int lambdaBiasCutoff = hd.lambdaBiasCutoff;
             for (int iL = -lambdaBiasCutoff; iL <= lambdaBiasCutoff; iL++) {
                 int lambdaBin = currentLambdaBin + iL;
                 double deltaL = currentLambda - lambdaBin * hd.lambdaBinWidth;
                 double deltaL2 = deltaL * deltaL;
 
                 // Pre-compute the lambda bias.
                 double L2exp = exp(-deltaL2 * invLs2);
 
                 // Mirror condition for Lambda counts.
                 lambdaBin = lambdaMirror(lambdaBin);
                 double mirrorFactor = mirrorFactor(lambdaBin);
 
                 int dUdLbiasCutoff = hd.dUdLBiasCutoff;
                 for (int jFL = -dUdLbiasCutoff; jFL <= dUdLbiasCutoff; jFL++) {
                     int dUdLBin = currentdUdLBin + jFL;
                     double weight = mirrorFactor * getRecursionKernelValue(lambdaBin, dUdLBin);
                     if (weight <= 0.0) {
                         continue;
                     }
                     double deltaFL = currentdUdL - dUdLforBin(dUdLBin);
                     double deltaFL2 = deltaFL * deltaFL;
                     double e = weight * L2exp * exp(-deltaFL2 * invFLs2);
                     sum += e;
                 }
             }
             return gaussianBiasMagnitude * sum;
         }
 
         /**
          * Compute the total Bias energy at (currentLambda, currentdUdL).
          *
          * @param currentLambda The value of lambda.
          * @param currentdUdL   The value of dU/dL.
          * @return The bias energy.
          */
         double computeBiasEnergy(double currentLambda, double currentdUdL) {
             synchronized (this) {
                 int currentLambdaBin = indexForLambda(currentLambda);
                 int currentdUdLBin = binFordUdL(currentdUdL);
 
                 double bias1D;
                 double bias2D = 0.0;
 
                 if (!hd.metaDynamics) {
                     if (hd.biasMag > 0.0) {
                         int lambdaBiasCutoff = hd.lambdaBiasCutoff;
                         for (int iL = -lambdaBiasCutoff; iL <= lambdaBiasCutoff; iL++) {
 
                             int lambdaBin = currentLambdaBin + iL;
                             double deltaL = currentLambda - (lambdaBin * hd.lambdaBinWidth);
                             double deltaL2 = deltaL * deltaL;
                             double expL2 = exp(-deltaL2 / (2.0 * hd.lambdaVariance));
 
                             // Mirror conditions for recursion kernel counts.
                             lambdaBin = lambdaMirror(lambdaBin);
                             double mirrorFactor = mirrorFactor(lambdaBin);
 
                             int dUdLbiasCutoff = hd.dUdLBiasCutoff;
                             for (int iFL = -dUdLbiasCutoff; iFL <= dUdLbiasCutoff; iFL++) {
                                 int dUdLBin = currentdUdLBin + iFL;
 
                                 double weight = mirrorFactor * getRecursionKernelValue(lambdaBin, dUdLBin);
                                 if (weight <= 0.0) {
                                     continue;
                                 }
 
                                 double deltaFL = currentdUdL - dUdLforBin(dUdLBin);
                                 double deltaFL2 = deltaFL * deltaFL;
                                 double bias = weight * hd.biasMag * expL2 * exp(-deltaFL2 / (2.0 * hd.dUdLVariance));
                                 bias2D += bias;
                             }
                         }
                     }
                     // Compute the energy for the recursion worker at F(L) using interpolation.
                     bias1D = energyAndGradient1D(currentLambda, false);
                 } else {
                     bias1D = energyAndGradientMeta(currentLambda, false);
                 }
 
                 // Return the total bias.
                 return bias1D + bias2D;
             }
         }
 
         /**
          * Compute the total Bias energy at (currentLambda, currentdUdL) and its gradient.
          *
          * @param currentdUdLambda The value of dU/dL.
          * @param chainRule        The chain rule terms for the bias.
          * @return The bias energy.
          */
         double energyAndGradient2D(double currentdUdLambda, double[] chainRule) {
             return energyAndGradient2D(ld.lambda, currentdUdLambda, chainRule, hd.biasMag);
         }
 
         /**
          * Compute the total Bias energy at (currentLambda, currentdUdL) and its gradient.
          *
          * @param currentLambda         The value of lambda.
          * @param currentdUdLambda      The value of dU/dL.
          * @param chainRule             The chain rule terms for the bias.
          * @param gaussianBiasMagnitude The magnitude of the Gaussian bias.
          * @return The bias energy.
          */
         double energyAndGradient2D(double currentLambda, double currentdUdLambda, double[] chainRule,
                                    double gaussianBiasMagnitude) {
 
             if (gaussianBiasMagnitude <= 0.0) {
                 chainRule[0] = 0.0;
                 chainRule[1] = 0.0;
                 return 0;
             }
 
             // Calculate recursion kernel G(L, F_L) and its derivatives with respect to L and F_L.
             double gLdEdL = 0.0;
             double dGdLambda = 0.0;
             double dGdFLambda = 0.0;
             int currentLambdaBin = indexForLambda(currentLambda);
             int currentdUdLBin = binFordUdL(currentdUdLambda);
 
             int lambdaBiasCutoff = hd.lambdaBiasCutoff;
             for (int iL = -lambdaBiasCutoff; iL <= lambdaBiasCutoff; iL++) {
                 int lambdaBin = currentLambdaBin + iL;
 
                 // Pass in the bin offset to the weight instance.
                 double deltaL = currentLambda - (lambdaBin * hd.lambdaBinWidth);
                 double deltaL2 = deltaL * deltaL;
                 double expL2 = exp(-deltaL2 / (2.0 * hd.lambdaVariance));
 
                 // Mirror conditions for recursion kernel counts.
                 double mirrorFactor = mirrorFactor(lambdaBin);
                 int dUdLbiasCutoff = hd.dUdLBiasCutoff;
                 for (int iFL = -dUdLbiasCutoff; iFL <= dUdLbiasCutoff; iFL++) {
                     int dUdLBin = currentdUdLBin + iFL;
 
                     double weight = mirrorFactor * getRecursionKernelValue(lambdaBin, dUdLBin);
                     if (weight <= 0.0) {
                         continue;
                     }
 
                     double deltaFL = currentdUdLambda - dUdLforBin(dUdLBin);
                     double deltaFL2 = deltaFL * deltaFL;
                     double bias = weight * gaussianBiasMagnitude * expL2 * exp(-deltaFL2 / (2.0 * hd.dUdLVariance));
                     gLdEdL += bias;
                     dGdLambda -= deltaL / hd.lambdaVariance * bias;
                     dGdFLambda -= deltaFL / hd.dUdLVariance * bias;
                 }
             }
 
             chainRule[0] = dGdLambda;
             chainRule[1] = dGdFLambda;
             return gLdEdL;
         }
 
 
         /**
          * This calculates the 1D OST bias and its derivative with respect to Lambda.
          * <p>
          * See Equation 8 in <a href="http://doi.org/10.1021/ct300035u">
          * The Structure, Thermodynamics, and Solubility of Organic Crystals from Simulation with a
          * Polarizable Force Field</a>.
          * <p>
          *
          * @param gradient If true, compute the gradient.
          * @return a double.
          */
         private double energyAndGradient1D(boolean gradient) {
             return energyAndGradient1D(ld.lambda, gradient);
         }
 
         /**
          * This calculates the 1D OST bias and its derivative with respect to Lambda.
          * <p>
          * See Equation 8 in <a href="http://doi.org/10.1021/ct300035u">
          * The Structure, Thermodynamics, and Solubility of Organic Crystals from Simulation with a
          * Polarizable Force Field</a>.
          * <p>
          *
          * @param currentLambda The current value of lambda.
          * @param gradient      If true, compute the gradient.
          * @return a double.
          */
         private double energyAndGradient1D(double currentLambda, boolean gradient) {
             synchronized (this) {
                 double biasEnergy = 0.0;
                 for (int iL0 = 0; iL0 < hd.lambdaBins - 1; iL0++) {
                     int iL1 = iL0 + 1;
 
                     // Find bin centers and values for interpolation / extrapolation points.
                     double L0 = iL0 * hd.lambdaBinWidth;
                     double L1 = L0 + hd.lambdaBinWidth;
                     double FL0 = ensembleAveragedUdL[iL0];
                     double FL1 = ensembleAveragedUdL[iL1];
                     double deltaFL = FL1 - FL0;
           /*
            If the lambda is less than or equal to the upper limit, this is
            the final interval. Set the upper limit to L, compute the partial
            derivative and break.
           */
                     boolean done = false;
                     if (currentLambda <= L1) {
                         done = true;
                         L1 = currentLambda;
                     }
 
                     // Upper limit - lower limit of the integral of the extrapolation / interpolation.
                     biasEnergy += (FL0 * L1 + deltaFL * L1 * (0.5 * L1 - L0) / hd.lambdaBinWidth);
                     biasEnergy -= (FL0 * L0 + deltaFL * L0 * (-0.5 * L0) / hd.lambdaBinWidth);
                     if (done) {
                         // Compute the gradient d F(L) / dL at L.
                         if (gradient) {
                             dUdLambda -= FL0 + (L1 - L0) * deltaFL / hd.lambdaBinWidth;
                         }
                         break;
                     }
                 }
                 return -biasEnergy;
             }
         }
 
         /**
          * This calculates the 1D Metadynamics bias and its derivative with respect to Lambda.
          *
          * @param gradient If true, compute the gradient.
          * @return a double.
          */
         private double energyAndGradientMeta(boolean gradient) {
             return energyAndGradientMeta(ld.lambda, gradient);
         }
 
         /**
          * This calculates the 1D Metadynamics bias and its derivative with respect to Lambda.
          *
          * @param currentLambda The current value of lambda.
          * @param gradient      If true, compute the gradient.
          * @return a double.
          */
         private double energyAndGradientMeta(double currentLambda, boolean gradient) {
             // Zero out the metadynamics bias energy.
             double biasEnergy = 0.0;
 
             // Current lambda bin.
             int currentBin = indexForLambda(currentLambda);
 
             // Loop over bins within the cutoff.
             int lambdaBiasCutoff = hd.lambdaBiasCutoff;
             for (int iL = -lambdaBiasCutoff; iL <= lambdaBiasCutoff; iL++) {
                 int lambdaBin = currentBin + iL;
                 double deltaL = currentLambda - (lambdaBin * hd.lambdaBinWidth);
                 double deltaL2 = deltaL * deltaL;
                 // Mirror conditions for the lambda bin and count magnitude.
                 lambdaBin = lambdaMirror(lambdaBin);
                 double mirrorFactor = mirrorFactor(lambdaBin);
                 double weight = mirrorFactor * hd.biasMag * countsForLambda(lambdaBin);
                 if (weight > 0) {
                     double e = weight * exp(-deltaL2 / (2.0 * hd.lambdaVariance));
                     biasEnergy += e;
                     // Add the dMeta/dL contribution.
                     if (gradient) {
                         dUdLambda -= deltaL / hd.lambdaVariance * e;
                     }
                 }
             }
 
             return biasEnergy;
         }
 
         /**
          * Marginalize over dU/dL counts for the given lambda bin.
          *
          * @param lambdaBin Lambda bin to marginalize for.
          * @return Total number of counts.
          */
         private double countsForLambda(int lambdaBin) {
             synchronized (this) {
                 double count = 0.0;
                 for (int i = 0; i < hd.dUdLBins; i++) {
                     count += hd.zHistogram[lambdaBin][i];
                 }
                 return count;
             }
         }
 
         /**
          * If necessary, allocate more space.
          */
         void checkRecursionKernelSize(double currentdUdL) {
             // Synchronize updates of the recursion kernel.
             synchronized (this) {
                 if (currentdUdL > hd.dUdLMaximum) {
                     logger.info(format(" Current F_lambda %8.2f > maximum histogram size %8.2f.", currentdUdL, hd.dUdLMaximum));
 
                     double origDeltaG = updateFreeEnergyDifference(false, false);
 
                     int newdUdLBins = hd.dUdLBins;
                     while (hd.dUdLMinimum + newdUdLBins * hd.dUdLBinWidth < currentdUdL) {
                         newdUdLBins += 100;
                     }
                     double[][] newRecursionKernel = new double[hd.lambdaBins][newdUdLBins];
 
                     // We have added bins above the indices of the current counts just copy them into the new array.
                     for (int i = 0; i < hd.lambdaBins; i++) {
                         arraycopy(hd.zHistogram[i], 0, newRecursionKernel[i], 0, hd.dUdLBins);
                     }
 
                     hd.zHistogram = newRecursionKernel;
                     hd.dUdLBins = newdUdLBins;
                     kernelValues = new double[hd.dUdLBins];
                     hd.dUdLMaximum = hd.dUdLMinimum + hd.dUdLBinWidth * hd.dUdLBins;
                     logger.info(format(" New histogram %8.2f to %8.2f with %d bins.\n", hd.dUdLMinimum, hd.dUdLMaximum, hd.dUdLBins));
 
                     double newFreeEnergy = updateFreeEnergyDifference(true, false);
                     assert (origDeltaG == newFreeEnergy);
                 }
                 if (currentdUdL < hd.dUdLMinimum) {
                     logger.info(format(" Current F_lambda %8.2f < minimum histogram size %8.2f.", currentdUdL, hd.dUdLMinimum));
 
                     double origDeltaG = updateFreeEnergyDifference(false, false);
 
                     int offset = 100;
                     while (currentdUdL < hd.dUdLMinimum - offset * hd.dUdLBinWidth) {
                         offset += 100;
                     }
                     int newdUdLBins = hd.dUdLBins + offset;
                     double[][] newRecursionKernel = new double[hd.lambdaBins][newdUdLBins];
 
                     // We have added bins below the current counts,
                     // so their indices must be increased by: offset = newFLBins - FLBins
                     for (int i = 0; i < hd.lambdaBins; i++) {
                         arraycopy(hd.zHistogram[i], 0, newRecursionKernel[i], offset, hd.dUdLBins);
                     }
 
                     hd.zHistogram = newRecursionKernel;
                     hd.dUdLMinimum = hd.dUdLMinimum - offset * hd.dUdLBinWidth;
                     hd.dUdLBins = newdUdLBins;
                     kernelValues = new double[hd.dUdLBins];
                     logger.info(format(" New histogram %8.2f to %8.2f with %d bins.\n", hd.dUdLMinimum, hd.dUdLMaximum, hd.dUdLBins));
 
                     double newFreeEnergy = updateFreeEnergyDifference(true, false);
                     assert (origDeltaG == newFreeEnergy);
                 }
             }
         }
 
         /**
          * Add a Gaussian hill to the Histogram at (lambda, dUdL).
          *
          * @param dUdL The value of dU/dL.
          */
         void addBias(double dUdL) {
             // Communicate adding the bias to all walkers.
             if (hd.asynchronous) {
                 sendAsynchronous.send(ld.lambda, dUdL, temperingWeight);
             } else {
                 sendSynchronous.send(ld.lambda, dUdL, temperingWeight);
             }
         }
 
         /**
          * Pre-force portion of the stochastic integrator.
          */
         private void stochasticPreForce() {
             // Update theta.
             double[] thetaPosition = lambdaState.x();
             thetaPosition[0] = ld.theta;
             stochastic.preForce(OrthogonalSpaceTempering.this);
             ld.theta = thetaPosition[0];
 
             // Maintain theta in the interval -PI to PI.
             if (ld.theta > PI) {
                 ld.theta -= 2.0 * PI;
             } else if (ld.theta <= -PI) {
                 ld.theta += 2.0 * PI;
             }
 
             // Update lambda as sin(theta)^2.
             double sinTheta = sin(ld.theta);
             ld.lambda = sinTheta * sinTheta;
             lambdaInterface.setLambda(ld.lambda);
         }
 
         /**
          * Post-force portion of the stochastic integrator.
          */
         private void stochasticPostForce() {
             double[] thetaPosition = lambdaState.x();
             double[] gradient = lambdaState.gradient();
             // Update theta and dU/dtheta.
             gradient[0] = dUdLambda * sin(2 * ld.theta);
             thetaPosition[0] = ld.theta;
             stochastic.postForce(gradient);
 
             // Load the new theta, velocity and acceleration.
             double[] thetaVelocity = lambdaState.v();
             double[] thetaAcceleration = lambdaState.a();
             ld.theta = thetaPosition[0];
             ld.thetaVelocity = thetaVelocity[0];
             ld.thetaAcceleration = thetaAcceleration[0];
 
             // Maintain theta in the interval -PI to PI.
             if (ld.theta > PI) {
                 ld.theta -= 2.0 * PI;
             } else if (ld.theta <= -PI) {
                 ld.theta += 2.0 * PI;
             }
 
             // Update lambda as sin(theta)^2.
             double sinTheta = sin(ld.theta);
             ld.lambda = sinTheta * sinTheta;
             lambdaInterface.setLambda(ld.lambda);
         }
 
         /**
          * Gets the last lambda value received by this Histogram. This can be out-of-date w.r.t. the
          * OST's current lambda!
          *
          * @return Lambda value of the last bias added to this Histogram.
          */
         double getLastReceivedLambda() {
             return lastReceivedLambda;
         }
 
         public void setLastReceivedLambda(double lastReceivedLambda) {
             this.lastReceivedLambda = lastReceivedLambda;
         }
 
         /**
          * Gets the last dU/dL value received by this Histogram. This can be out-of-date w.r.t. the OST's
          * current dU/dL!
          *
          * @return dU/dL value of the last bias added to this Histogram.
          */
         double getLastReceivedDUDL() {
             return lastReceiveddUdL;
         }
 
         void destroy() {
             if (sendAsynchronous != null && sendAsynchronous.isAlive()) {
                 double[] killMessage = new double[]{Double.NaN, Double.NaN, Double.NaN, Double.NaN};
                 DoubleBuf killBuf = DoubleBuf.buffer(killMessage);
                 try {
                     logger.fine(" Sending the termination message.");
                     world.send(rank, killBuf);
                     logger.fine(" Termination message was sent successfully.");
                     logger.fine(format(" Receive thread alive %b status %s", sendAsynchronous.isAlive(),
                             sendAsynchronous.getState()));
                 } catch (Exception ex) {
                     String message = format(" Asynchronous Multiwalker OST termination signal "
                             + "failed to be sent for process %d.", rank);
                     logger.log(Level.SEVERE, message, ex);
                 }
             } else {
                 logger.fine(
                         " CountReceiveThread was either not initialized, or is not alive. This is the case for the Histogram script.");
             }
         }
 
     }
 
 }