def _geometry_forward(self, topology_proposal, old_sampler_state):
        """
        Run geometry engine to propose new positions and compute logP

        Parameters
        ----------
        topology_proposal : TopologyProposal
            Contains old/new Topology and System objects and atom mappings.
        old_sampler_state : openmmtools.states.SamplerState
            Configurational properties of the old system atoms.

        Returns
        -------
        new_sampler_state : openmmtools.states.SamplerState
            Configurational properties of new atoms proposed by geometry engine calculation.
        geometry_logp_propose : float
            The log probability of the forward-only proposal
        """
        if self.verbose: print("Geometry engine proposal...")
        # Generate coordinates for new atoms and compute probability ratio of old and new probabilities.
        initial_time = time.time()
        new_positions, geometry_logp_propose = self.geometry_engine.propose(topology_proposal, old_sampler_state.positions, self.sampler.thermodynamic_state.beta)
        if self.verbose: print('proposal took %.3f s' % (time.time() - initial_time))

        if self.geometry_pdbfile is not None:
            print("Writing proposed geometry...")
            from simtk.openmm.app import PDBFile
            PDBFile.writeFile(topology_proposal.new_topology, new_positions, file=self.geometry_pdbfile)
            self.geometry_pdbfile.flush()

        new_sampler_state = SamplerState(new_positions, box_vectors=old_sampler_state.box_vectors)  

        return new_sampler_state, geometry_logp_propose

def add_hydrogens_to_mol(mol):
  """
  Add hydrogens to a molecule object
  TODO (LESWING) see if there are more flags to add here for default
  :param mol: Rdkit Mol
  :return: Rdkit Mol
  """
  molecule_file = None
  try:
    pdbblock = Chem.MolToPDBBlock(mol)
    pdb_stringio = StringIO()
    pdb_stringio.write(pdbblock)
    pdb_stringio.seek(0)
    fixer = PDBFixer(pdbfile=pdb_stringio)
    fixer.addMissingHydrogens(7.4)

    hydrogenated_io = StringIO()
    PDBFile.writeFile(fixer.topology, fixer.positions, hydrogenated_io)
    hydrogenated_io.seek(0)
    return Chem.MolFromPDBBlock(
        hydrogenated_io.read(), sanitize=False, removeHs=False)
  except ValueError as e:
    logging.warning("Unable to add hydrogens", e)
    raise MoleculeLoadException(e)
  finally:
    try:
      os.remove(molecule_file)
    except (OSError, TypeError):
      pass

def write_trajectory_dcd(netcdf_filename, topology, pdb_trajectory_filename, dcd_trajectory_filename):
    """
    Write trajectory.

    Parameters
    ----------
    netcdf_filename : str
        NetCDF filename.
    topology : Topology
        Topology object
    pdb_trajectory_filename : str
        PDB trajectory output filename
    dcd_trajectory_filename : str
        Output trajectory filename.

    """
    ncfile = netCDF4.Dataset(netcdf_filename, 'r')
    [nsamples, nstates] = ncfile.variables['logZ'].shape

    # Write reference.pdb file
    from simtk.openmm.app import PDBFile
    outfile = open(pdb_trajectory_filename, 'w')
    positions = unit.Quantity(ncfile.variables['positions'][0,:,:], unit.angstroms)
    PDBFile.writeFile(topology, positions, file=outfile)
    outfile.close()

    # TODO: Export as DCD trajectory with MDTraj
    from mdtraj.formats import DCDTrajectoryFile
    with DCDTrajectoryFile(dcd_trajectory_filename, 'w') as f:
        f.write(ncfile.variables['positions'][:,:,:])

    def from_pdb(cls, path, forcefield=None, **kwargs):
        """
        Loads topology, positions and, potentially, velocities and vectors,
        from a PDB file

        Parameters
        ----------
        path : str
            Path to PDB file
        forcefields : list of str
            Paths to FFXML and/or FRCMOD forcefields. REQUIRED.

        Returns
        -------
        pdb : SystemHandler
            SystemHandler with topology, positions, and, potentially, velocities and
            box vectors. Forcefields are embedded in the `master` attribute.
        """
        pdb = PDBFile(path)
        box = kwargs.pop('box', pdb.topology.getPeriodicBoxVectors())
        positions = kwargs.pop('positions', pdb.positions)
        velocities = kwargs.pop('velocities', getattr(pdb, 'velocities', None))

        if not forcefield:
            from .md import FORCEFIELDS as forcefield
            print('INFO: Forcefields for PDB not specified. Using default:\n ',
                  ', '.join(forcefield))
        pdb.forcefield = ForceField(*list(process_forcefield(*forcefield)))

        return cls(master=pdb, topology=pdb.topology, positions=positions,
                   velocities=velocities, box=box, path=path, **kwargs)

    def __init__(self, file):
        """Load a prmtop file."""
        top = Topology()
        ## The Topology read from the prmtop file
        self.topology = top

        # Load the prmtop file

        prmtop = amber_file_parser.PrmtopLoader(file)
        self._prmtop = prmtop

        # Add atoms to the topology

        PDBFile._loadNameReplacementTables()
        lastResidue = None
        c = top.addChain()
        for index in range(prmtop.getNumAtoms()):
            resNumber = prmtop.getResidueNumber(index)
            if resNumber != lastResidue:
                lastResidue = resNumber
                resName = prmtop.getResidueLabel(iAtom=index).strip()
                if resName in PDBFile._residueNameReplacements:
                    resName = PDBFile._residueNameReplacements[resName]
                r = top.addResidue(resName, c)
                if resName in PDBFile._atomNameReplacements:
                    atomReplacements = PDBFile._atomNameReplacements[resName]
                else:
                    atomReplacements = {}
            atomName = prmtop.getAtomName(index).strip()
            if atomName in atomReplacements:
                atomName = atomReplacements[atomName]

            # Try to guess the element.

            upper = atomName.upper()
            if upper.startswith('CL'):
                element = elem.chlorine
            elif upper.startswith('NA'):
                element = elem.sodium
            elif upper.startswith('MG'):
                element = elem.magnesium
            else:
                try:
                    element = elem.get_by_symbol(atomName[0])
                except KeyError:
                    element = None
            top.addAtom(atomName, element, r)

        # Add bonds to the topology

        atoms = list(top.atoms())
        for bond in prmtop.getBondsWithH():
            top.addBond(atoms[bond[0]], atoms[bond[1]])
        for bond in prmtop.getBondsNoH():
            top.addBond(atoms[bond[0]], atoms[bond[1]])

        # Set the periodic box size.

        if prmtop.getIfBox():
            top.setUnitCellDimensions(tuple(x.value_in_unit(unit.nanometer) for x in prmtop.getBoxBetaAndDimensions()[1:4])*unit.nanometer)

 def write_pdb(self, path):
     """
     Outputs a PDB file with the current contents of the system
     """
     if self.master is None and self.positions is None:
         raise ValueError('Topology and positions are needed to write output files.')
     with open(path, 'w') as f:
         PDBFile.writeFile(self.topology, self.positions, f)

    def update(self):
        """
        Update the sampler with one step of sampling.
        """
        if not self._initialized:
            self._initialize()

        if self.verbose:
            print("." * 80)
            print("MCMC sampler iteration %d" % self.iteration)

        initial_time = time.time()

        # Reset statistics
        if hasattr(self.integrator, 'setGlobalVariableByName'):
            self.integrator.setGlobalVariableByName('naccept', 0)

        # Take some steps
        self.integrator.step(self.nsteps)

        # Get new sampler state.
        self.sampler_state = SamplerState.createFromContext(self.context)

        # Report statistics
        if hasattr(self.integrator, 'getGlobalVariableByName'):
            naccept = self.integrator.getGlobalVariableByName('naccept')
            fraction_accepted = float(naccept) / float(self.nsteps)
            if self.verbose: print("Accepted %d / %d GHMC steps (%.2f%%)." % (naccept, self.nsteps, fraction_accepted * 100))

        final_time = time.time()
        elapsed_time = final_time - initial_time
        self._timing['sample positions'] = elapsed_time

        if self.verbose:
            final_energy = self.context.getState(getEnergy=True).getPotentialEnergy() * self.thermodynamic_state.beta
            print('Final energy is %12.3f kT' % (final_energy))
            print('elapsed time %8.3f s' % elapsed_time)

        if self.ncfile:
            self.ncfile.variables['positions'][self.iteration,:,:] = self.sampler_state.positions[:,:] / unit.nanometers
            for k in range(3):
                self.ncfile.variables['box_vectors'][self.iteration,k,:] = self.sampler_state.box_vectors[k,:] / unit.nanometers
            self.ncfile.variables['potential'][self.iteration] = self.thermodynamic_state.beta * self.context.getState(getEnergy=True).getPotentialEnergy()
            self.ncfile.variables['sample_positions_time'][self.iteration] = elapsed_time

        # Increment iteration count
        self.iteration += 1

        if self.verbose:
            print("." * 80)
        if self.pdbfile is not None:
            print("Writing frame...")
            from simtk.openmm.app import PDBFile
            PDBFile.writeModel(self.topology, self.sampler_state.positions, self.pdbfile, self.iteration)
            self.pdbfile.flush()

    def report(self, simulation, state):
        """Generate a report.

        Parameters:
         - simulation (Simulation) The Simulation to generate a report for
         - state (State) The current state of the simulation
        """
        if self._nextModel == 0:
            PDBFile.writeHeader(simulation.topology, self._out)
            self._topology = simulation.topology
            self._nextModel += 1
        PDBFile.writeModel(simulation.topology, state.getPositions(), self._out, self._nextModel)
        self._nextModel += 1

    def __init__(self, **kwargs):
        super(AlanineDipeptideExplicitSimulatedTempering, self).__init__(**kwargs)
        self.description = 'Alanine dipeptide in explicit solvent simulated tempering simulation'

        # Create topology, positions, and system.
        from openmmtools.testsystems import AlanineDipeptideExplicit
        testsystem = AlanineDipeptideExplicit(nonbondedMethod=app.CutoffPeriodic)
        self.topology = testsystem.topology
        self.positions = testsystem.positions
        self.system = testsystem.system

        # DEBUG: Write PDB
        from simtk.openmm.app import PDBFile
        outfile = open('initial.pdb', 'w')
        PDBFile.writeFile(self.topology, self.positions, outfile)
        outfile.close()

        # Add a MonteCarloBarostat
        temperature = 270 * unit.kelvin # will be replaced as thermodynamic state is updated
        pressure = 1.0 * unit.atmospheres
        barostat = openmm.MonteCarloBarostat(pressure, temperature)
        self.system.addForce(barostat)

        # Create thermodynamic states.
        Tmin = 270 * unit.kelvin
        Tmax = 600 * unit.kelvin
        ntemps = 256 # number of temperatures
        from sams import ThermodynamicState
        temperatures = unit.Quantity(np.logspace(np.log10(Tmin / unit.kelvin), np.log10(Tmax / unit.kelvin), ntemps), unit.kelvin)
        self.thermodynamic_states = [ ThermodynamicState(system=self.system, temperature=temperature, pressure=pressure) for temperature in temperatures ]

        # Create SAMS samplers
        from sams.samplers import SamplerState, MCMCSampler, ExpandedEnsembleSampler, SAMSSampler
        thermodynamic_state_index = 0 # initial thermodynamic state index
        thermodynamic_state = self.thermodynamic_states[thermodynamic_state_index]
        sampler_state = SamplerState(positions=self.positions)
        self.mcmc_sampler = MCMCSampler(sampler_state=sampler_state, thermodynamic_state=thermodynamic_state, ncfile=self.ncfile)
        #self.mcmc_sampler.pdbfile = open('output.pdb', 'w')
        self.mcmc_sampler.topology = self.topology
        self.mcmc_sampler.nsteps = 500
        self.mcmc_sampler.timestep = 2.0 * unit.femtoseconds
        self.mcmc_sampler.verbose = True
        self.exen_sampler = ExpandedEnsembleSampler(self.mcmc_sampler, self.thermodynamic_states)
        self.exen_sampler.verbose = True
        self.sams_sampler = SAMSSampler(self.exen_sampler)
        self.sams_sampler.verbose = True

def check_hydrogens(molecule, ID):
    # Check that Hydrogens are in structure
    if len(molecule.top.select("name == H")) == 0:
        # If absent, then add Hydrogens using the Amber99sb force-field
        try:
            from simtk.openmm.app import PDBFile, Modeller, ForceField
            pdb = PDBFile(ID + ".pdb")
            modeller = Modeller(pdb.topology, pdb.positions)
            forcefield = ForceField('amber99sb.xml','tip3p.xml')
            modeller.addHydrogens(forcefield)
            PDBFile.writeFile(modeller.topology, modeller.positions, open(ID + ".pdb", 'w'))
            molecule = md.load(ID + ".pdb").remove_solvent()
        except:
            warnings.warn("""PDB topology missing Hydrogens. Either manually add
            or install OpenMM through SIMTK to automatically correct.""")
            pass
    return molecule

def download_pdb(pdbid, file_pathway):
    """

    Args:
        pdbid: 4 letter string specifying the PDB ID of the file yoou want to fix
        file_pathway: a string containing the pathway specifying how you want to organize the PDB files once written

    Returns: nothing, but it does write the PDB file

    ***Note: this function does NOT fix any mistakes with the PDB file

    """

    if not os.path.exists(file_pathway):
        os.makedirs(file_pathway)
    fixer = PDBFixer(pdbid=pdbid)
    PDBFile.writeFile(fixer.topology, fixer.positions, open(os.path.join(file_pathway, '%s.pdb' % pdbid), 'w'))

    def update(self):
        """
        Update the sampler with one step of sampling.
        """
        if self.verbose:
            print("-" * 80)
            print("Expanded Ensemble sampler iteration %8d" % self.iteration)
        self.update_positions()
        self.update_state()
        self.iteration += 1
        if self.verbose:
            print("-" * 80)

        if self.pdbfile is not None:
            print("Writing frame...")
            from simtk.openmm.app import PDBFile
            PDBFile.writeModel(self.topology, self.sampler.sampler_state.positions, self.pdbfile, self.iteration)
            self.pdbfile.flush()

    def report(self, simulation, state):
        """Generate a report.

        Parameters
        ----------
        simulation : Simulation
            The Simulation to generate a report for
        state : State
            The current state of the simulation
        """
        if self._nextModel == 0:
            PDBFile.writeHeader(simulation.topology, self._out)
            self._topology = simulation.topology
            self._nextModel += 1
        PDBFile.writeModel(simulation.topology, state.getPositions(), self._out, self._nextModel)
        self._nextModel += 1
        if hasattr(self._out, 'flush') and callable(self._out.flush):
            self._out.flush()

    def report(self, simulation, _):
        """Generate a report.

        Parameters
        ----------
        simulation : Simulation
            The Simulation to generate a report for
        _ : State
            The current state of the simulation
        """
        state = simulation.context.getState(getPositions=True, enforcePeriodicBox=self._enforcePeriodicBox)
        if self._nextModel == 0:
            PDBFile.writeHeader(simulation.topology, self._out)
            self._topology = simulation.topology
            self._nextModel += 1
        PDBFile.writeModel(simulation.topology, state.getPositions(), self._out, self._nextModel)
        self._nextModel += 1
        if hasattr(self._out, 'flush') and callable(self._out.flush):
            self._out.flush()

def pdb_fix_pdbfixer(pdbid, file_pathway, ph, chains_to_remove):
    """

    Args:
        pdbid: 4 letter string specifying the PDB ID of the file yoou want to fix
        file_pathway: a string containing the pathway specifying how you want to organize the PDB files once written
        ph: the pH at which hydrogens will be determined and added
        chains_to_remove: dictionary containing pdbs with chains to remove
    Returns: nothing, but it does right PDB files

    """
    print(pdbid)

    # Download the topology from rcsb based on pdbod
    fixer = PDBFixer(pdbid=pdbid)

    # Remove chains based on hand curated .csv file
    if pdbid in chains_to_remove['pdbid']:
        chains = chains_to_remove['chain_to_remove'][chain_to_remove['pdbid'].index(pdbid)]
        chains_list = chains.split()
        fixer.removeChains(chainIds=chains_list)

    # Determine the first and last residue resolved in chain 0
    chains = [chain for chain in fixer.topology.chains()]
    resindices = [residue.index for residue in chains[0].residues()]
    resindices = natsorted(resindices)
    first_resindex = resindices[0]
    last_resindex = resindices[-1]

    # Find Missing residues and determine if they are C or N terminal fragments (which will be removed)

    fixer.findMissingResidues()
    if len(fixer.missingResidues) > 0:
        if sorted(fixer.missingResidues.keys())[0][-1] <= first_resindex:
            fixer.missingResidues.pop((sorted(fixer.missingResidues.keys())[0]))

        if sorted(fixer.missingResidues.keys())[-1][-1] >= last_resindex:
            fixer.missingResidues.pop((sorted(fixer.missingResidues.keys())[-1]))

    fixer.findNonstandardResidues()
    fixer.replaceNonstandardResidues()
    fixer.findMissingAtoms()
    fixer.addMissingAtoms()
    fixer.addMissingHydrogens(ph)
    # Write fixed PDB file, with all of the waters and ligands
    PDBFile.writeFile(fixer.topology, fixer.positions, open(os.path.join(file_pathway,
                                                                         '%s_fixed_ph%s.pdb' % (pdbid, ph)), 'w'),
                      keepIds=keepNumbers)

    # Remove the ligand and write a pdb file
    fixer.removeHeterogens(True)
    PDBFile.writeFile(fixer.topology, fixer.positions, open(os.path.join(file_pathway,
                                                                         '%s_fixed_ph%s_apo.pdb' % (pdbid, ph)), 'w'),
                      keepIds=keepNumbers)
    # Remove the waters and write a pdb file
    fixer.removeHeterogens(False)
    PDBFile.writeFile(fixer.topology, fixer.positions, open(os.path.join(file_pathway,
                                                                         '%s_fixed_ph%s_apo_nowater.pdb' % (pdbid, ph)),
                                                            'w'), keepIds=keepNumbers)

def fix_pdb(pdb_id, pdb_file, pdb_group):
    chains_to_retain = get_required_chains(pdb_group)
    chains_to_remove = []

    for chain in PDBParser().get_structure(pdb_id, pdb_file)[0]:
        if chain.get_id() not in chains_to_retain:
            chains_to_remove.append(chain.get_id())

    fixer = PDBFixer(filename=pdb_file)

    fixer.removeChains(chainIds=chains_to_remove)

    fixer.findMissingResidues()
    fixer.findMissingAtoms()
    fixer.addMissingAtoms()
    fixer.removeHeterogens(True)

    # KeepIds flag is critical here, otherwise we loose all information binding
    pdb_file = dirname(pdb_file) + '/' + pdb_id + '.pdb'
    PDBFile.writeFile(fixer.topology, fixer.positions, open(pdb_file, 'w'), keepIds=True)

    return pdb_file

 def __del__(self):
     if self._topology is not None:
         PDBFile.writeFooter(self._topology, self._out)
     self._out.close()

nsteps = 5 # number of timesteps per iteration
niterations = 500 # number of iterations
mctrials = 10 # number of Monte Carlo trials per iteration
nsalt = 0 # current number of salt pairs
tol = 1e-6 # constraint tolerance

# Determine number of molecules
nmolecules = 0
for residue in topology.residues():
    nmolecules += 1
print('system originally has %d water molecules' % nmolecules)

# Open PDB file for writing.
from simtk.openmm.app import PDBFile
pdbfile = open('output.pdb', 'w')
PDBFile.writeHeader(topology, file=pdbfile)
PDBFile.writeModel(topology, positions, file=pdbfile, modelIndex=0)

# Simulate
for iteration in range(niterations):
    print('iteration %5d / %5d' % (iteration, niterations))

    # Create a simulation
    from openmmtools.integrators import VelocityVerletIntegrator
    integrator = VelocityVerletIntegrator(timestep)
    integrator.setConstraintTolerance(tol)
    context = openmm.Context(system, integrator)
    context.setPositions(positions)

    # Propagate dynamics at constant counterion number.
    print('propagating dynamics for %d steps...' % nsteps)

def run(options):
        fixer = PDBFixer(options['pdb'])
        fixer.addMissingHydrogens(7.0)
        fixer.addSolvent(boxSize=Vec3(2.62,2.62,2.62)*nanometers, padding=None, 
                    positiveIon='Na+', negativeIon='Cl-', ionicStrength=0.0*molar)
        PDBFile.writeFile(fixer.topology, fixer.positions, open(options['outfile'], 'w'))

    def __init__(self, file, unitCellDimensions=None, includeDir='/usr/local/gromacs/share/gromacs/top', defines={}):
        """Load a top file.

        Parameters:
         - file (string) the name of the file to load
         - unitCellDimensions (Vec3=None) the dimensions of the crystallographic unit cell
         - includeDir (string=/usr/local/gromacs/share/gromacs/top) a directory in which to look for other files
           included from the top file
         - defines (map={}) preprocessor definitions that should be predefined when parsing the file
         """
        self._includeDirs = (os.path.dirname(file), includeDir)
        self._defines = defines

        # Parse the file.

        self._currentCategory = None
        self._ifStack = []
        self._moleculeTypes = {}
        self._molecules = []
        self._currentMoleculeType = None
        self._atomTypes = {}
        self._bondTypes= {}
        self._angleTypes = {}
        self._dihedralTypes = {}
        self._implicitTypes = {}
        self._pairTypes = {}
        self._cmapTypes = {}
        self._processFile(file)

        # Create the Topology from it.

        top = Topology()
        ## The Topology read from the prmtop file
        self.topology = top
        top.setUnitCellDimensions(unitCellDimensions)
        PDBFile._loadNameReplacementTables()
        for moleculeName, moleculeCount in self._molecules:
            if moleculeName not in self._moleculeTypes:
                raise ValueError("Unknown molecule type: "+moleculeName)
            moleculeType = self._moleculeTypes[moleculeName]

            # Create the specified number of molecules of this type.

            for i in range(moleculeCount):
                atoms = []
                lastResidue = None
                c = top.addChain()
                for index, fields in enumerate(moleculeType.atoms):
                    resNumber = fields[2]
                    if resNumber != lastResidue:
                        lastResidue = resNumber
                        resName = fields[3]
                        if resName in PDBFile._residueNameReplacements:
                            resName = PDBFile._residueNameReplacements[resName]
                        r = top.addResidue(resName, c)
                        if resName in PDBFile._atomNameReplacements:
                            atomReplacements = PDBFile._atomNameReplacements[resName]
                        else:
                            atomReplacements = {}
                    atomName = fields[4]
                    if atomName in atomReplacements:
                        atomName = atomReplacements[atomName]

                    # Try to guess the element.

                    upper = atomName.upper()
                    if upper.startswith('CL'):
                        element = elem.chlorine
                    elif upper.startswith('NA'):
                        element = elem.sodium
                    elif upper.startswith('MG'):
                        element = elem.magnesium
                    else:
                        try:
                            element = elem.get_by_symbol(atomName[0])
                        except KeyError:
                            element = None
                    atoms.append(top.addAtom(atomName, element, r))

                # Add bonds to the topology

                for fields in moleculeType.bonds:
                    top.addBond(atoms[int(fields[0])-1], atoms[int(fields[1])-1])