Cannot write RST7 file after GCMC — coordinates exceed Amber format width

[MartinFuentetaja]

Description

When running OpenBPMD with the --gcmc flag, the pipeline crashes with an OSError when attempting to save the post-GCMC system as an Amber RST7 file via sire.

The error occurs at:

sr.save(mols, os.path.join(out_dir, f"sd_{lig}_final.rst7"))

I have checked that mols = sr.load(input.rst7, input.prmtop) returns Periodic Boundary Box. I wish someone could give me insight why this is happening.

Root Cause

The Amber RST7 format stores coordinates as fixed-width fields of 12 characters with 7 decimal places (AmberFormat( 6 x float[width = 12, precision = 7] )). This format can only represent values in the range (-1000, 10000) Å — values ≤ −1000 Å require a 13-character field and cannot be written.

After GCMC sampling, a subset of water molecules are displaced outside the periodic box, resulting in wrapped coordinates with values such as -1153.06, -1249.21, etc. These out-of-box coordinates trigger the format overflow.

Example errors from the traceback:

Could not write the float at index 7340, value '-1153.06' into the specified format
AmberFormat( 6 x float[width = 12, precision = 7] )

Traceback

OSError: SireError::io_error: Cannot write the files as the following errors occurred:
Failed to write '/scratch/.../sd_JI1_final.rst7' using the parser for fileformat 'RST7'.
Errors converting the system to Amber Rst7 format...
Could not write the float at index 7340, value '-1153.06' into AmberFormat( 6 x float[width = 12, precision = 7] ).
[... 29 further instances of the same error ...]

[lohedges]

Can you check whether it is coordinates, or velocities that are causing the overflow. If the latter, then you can use save_velocities=False with sire.save, or use a binary format.

Note that I am not familiar with OpenBPMD. In using loch directly, or via our own software stack, I have never experienced wrapping issues. It might be worth contacting the OpenBPMD developer if this is only happening with their interface.

[MartinFuentetaja]

Hello!

Thank you for answering so fast. I have tried your suggestion but It did not work. I am able to save .pdb, .gro, .prmtop and .dcd files. After saving them, I am using mdtraj to remove PBC (similar to trjconv in GROMACS), and the resulting files are the ones I use for production (BPMD).

I will leave the code If you are interested in. I have checked that the resulting files from NVT simulation are correct.

#!/usr/bin/env python
descriptn = \
    """
    OpenBPMD - an open source implementation of Binding Pose Metadynamics
    (BPMD) with OpenMM. Replicates the protocol as described by
    Clark et al. 2016 (DOI: 10.1021/acs.jctc.6b00201).

    Runs ten 10 ns metadynamics simulations that biases the RMSD of the ligand.

    The stability of the ligand is calculated using the ligand RMSD (PoseScore)
    and the persistence of the original noncovalent interactions between the
    protein and the ligand (ContactScore). Stable poses have a low RMSD and
    a high fraction of the native contacts preserved until the end of the
    simulation.

    A composite score is calculated using the following formula:
    CompScore = PoseScore - 5 * ContactScore
    """

# OpenMM
from openmm import *
from openmm.app import *
from openmm.unit import *
from openmm.app.metadynamics import *

# The rest
import grand
from loch import GCMCSampler
import sire as sr
import argparse
import numpy as np
import MDAnalysis as mda
from MDAnalysis.analysis import rms, contacts
import mdtraj as md
import pandas as pd
import parmed as pmd
import glob
import os
from loguru import logger
from pathlib import Path

__author__ = "Dominykas Lukauskis"
__version__ = "1.0.0"
__email__ = "dominykas.lukauskis.19@ucl.ac.uk"

logger.add(sys.stderr, format="{time} {level} {message}", filter="my_module", level="INFO")

def main(args):
    """Main entry point of the app. Takes in argparse.Namespace object as
    a function argument. Carries out a sequence of steps required to obtain a
    stability score for a given ligand pose in the provided structure file.

    1. Load the structure and parameter files.
    2. If absent, create an output folder.
    3. Minimization up to ener. tolerance of 10 kJ/mol.
    4. 500 ps equilibration in NVT ensemble with position
       restraints on solute heavy atoms with the force 
       constant of 5 kcal/mol/A^2
    5. Run NREPs (default=10) of binding pose metadynamics simulations,
       writing trajectory files and a time-resolved BPM scores for each
       repeat.
    6. Collect results from the OpenBPMD simulations and
       write a final score for a given protein-ligand
       structure.

    Parameters
    ----------
    args.structure : str, default='solvated.rst7'
        Name of the structure file, either Amber or Gromacs format.
    args.parameters : str, default='solvated.prmtop'
        Name of the parameter or topology file, either Amber or Gromacs
        format.
    args.output : str, default='.'
        Path to and the name of the output directory.
    args.lig_resname : str, default='LIG'
        Residue name of the ligand in the structure/parameter file.
    args.nreps : int, default=10
        Number of repeat OpenBPMD simulations to run in series.
    args.hill_height : float, default=0.3
        Size of the metadynamical hill, in kcal/mol.
    """
    if args.structure.endswith('.gro'):
        coords = GromacsGroFile(args.structure)
        box_vectors = coords.getPeriodicBoxVectors()
        parm = GromacsTopFile(args.parameters, periodicBoxVectors=box_vectors)
    else:
        coords = AmberInpcrdFile(args.structure)
        parm = AmberPrmtopFile(args.parameters)

    if not os.path.isdir(f'{args.output}'):
        os.mkdir(f'{args.output}')

    # Minimize
    min_file_name = 'minimized_system.pdb'
    min_pdb = os.path.join(args.output,min_file_name)

    if not os.path.isfile(os.path.join(args.output,min_file_name)):
        logger.info("Minimizing...")
        minimize(parm, coords.positions, args.output, min_file_name)
    else:
        logger.info(f"Minimizing already done -> {min_pdb}. SKIPING.")


    # Equilibrate NVT
    eq_file_name = 'equil_system_NVT.pdb'
    eq_pdb = os.path.join(args.output,eq_file_name)

    if not os.path.isfile(eq_pdb):
        logger.info("Equilibrating (NVT)...")
        equilibrate(min_pdb, parm, args.output, eq_file_name, args.parameters, npt = False)
    else:
        logger.info(f"Equilibrating (NVT) already done -> {eq_pdb}. SKIPING.")

    eq_prmtop = str(Path(eq_pdb).with_suffix(".prmtop"))
    eq_rst7 = str(Path(eq_pdb).with_suffix(".rst7"))
    cent_eq_pdb = os.path.join(args.output,'centred_'+eq_file_name)
    cent_eq_rst7 = str(Path(cent_eq_pdb).with_suffix('.rst7'))

    if os.path.isfile(eq_pdb) and not os.path.isfile(cent_eq_pdb):
	# mdtraj can't use GMX TOP, so we have to specify the GRO file instead
        if args.structure.endswith('.gro'):
            mdtraj_top = args.structure
        else:
            mdtraj_top = args.parameters
        mdu = md.load(eq_pdb, top=eq_prmtop)#eq_rst7eq_prmtop
        mdu.image_molecules(inplace=True)
        
        mdu.save_pdb(cent_eq_pdb)
        mdu.save_amberrst7(cent_eq_rst7)

    # Equilibrate NPT
    if args.npt:
        
        eq_file_name = 'equil_system_NPT.pdb'
        eq_pdb = os.path.join(args.output,eq_file_name)

        if not os.path.isfile(eq_pdb):
            logger.info("Equilibrating (NPT)...")
            equilibrate(cent_eq_pdb, parm, args.output, eq_file_name, args.parameters, npt = True)
        else:
            logger.info(f"Equilibrating (NPT) already done -> {eq_pdb}. SKIPING.")
        
        eq_prmtop = str(Path(eq_pdb).with_suffix(".prmtop"))
        eq_rst7 = str(Path(eq_pdb).with_suffix(".rst7"))
        cent_eq_pdb = os.path.join(args.output,'centred_'+eq_file_name)
        cent_eq_rst7 = str(Path(cent_eq_pdb).with_suffix('.rst7'))

        if os.path.isfile(eq_pdb) and not os.path.isfile(cent_eq_pdb):
        # mdtraj can't use GMX TOP, so we have to specify the GRO file instead
           
            mdu = md.load(eq_rst7, top=eq_prmtop)
            mdu.image_molecules(inplace=True)
            
            mdu.save_pdb(cent_eq_pdb)
            mdu.save_amberrst7(cent_eq_rst7)

    
    # Grand-Canonical Monte Carlo (GCMC)
    if args.gcmc:
        gcmc_file_name = f"sd_{args.lig_resname}_final"
        gcmc_prmtop, gcmc_rst7, gcmc_pdb = path_to_files(gcmc_file_name, args.output)

        if not os.path.isfile(gcmc_pdb):
            logger.info("GCMC...")
            GCMC(cent_eq_rst7, eq_prmtop, args.output, args.lig_resname)
        else:
            logger.info(f"GCMC already done -> {gcmc_pdb}. SKIPING.")

        cent_gcmc_pdb = os.path.join(args.output,'centred_' + gcmc_file_name + '.pdb')
        
        if os.path.isfile(gcmc_pdb) and not os.path.isfile(cent_gcmc_pdb):
            
            mdu = md.load(gcmc_pdb, top=gcmc_prmtop)
            mdu.image_molecules(make_whole=True, inplace=True)
            
            mdu.save_pdb(cent_gcmc_pdb)

            logger.info("GCMC MDAnalysis completed!")

        # Renaming variables
        cent_eq_pdb = cent_gcmc_pdb
        eq_prmtop = gcmc_prmtop

    # Run NREPS number of production simulations
    for idx in range(0, args.nreps):
        rep_dir = os.path.join(args.output,f'rep_{idx}')
        if not os.path.isdir(rep_dir):
            os.mkdir(rep_dir)

        if os.path.isfile(os.path.join(rep_dir,'bpm_results.csv')):
            logger.info(f"BPMD rep_{idx} already completed! SKIPING!")
            continue
        
        logger.info(f"Initializing BPMD rep_{idx}...")

        produce(
            args.output, 
            idx, 
            args.lig_resname, 
            cent_eq_pdb,
            AmberPrmtopFile(eq_prmtop), #if (args.npt or args.gcmc) else parm,
            eq_prmtop, #if (args.npt or args.gcmc) else args.structure,
            args.hill_height
        )
        
        logger.info(f"BPMD rep_{idx} completed!")
              
        trj_name = os.path.join(rep_dir,'trj.dcd')
                
        PoseScoreArr = get_pose_score(cent_eq_pdb, trj_name, args.lig_resname)

        ContactScoreArr = get_contact_score(cent_eq_pdb, trj_name, args.lig_resname)

        # Calculate the CompScore at every frame
        CompScoreArr = np.zeros(99)
        for index in range(ContactScoreArr.shape[0]):
            ContactScore, PoseScore = ContactScoreArr[index], PoseScoreArr[index]
            CompScore = PoseScore - 5 * ContactScore
            CompScoreArr[index] = CompScore

        Scores = np.stack((CompScoreArr, PoseScoreArr, ContactScoreArr), axis=-1)

        # Save a DataFrame to CSV
        df = pd.DataFrame(Scores, columns=['CompScore', 'PoseScore',
                                           'ContactScore'])
        df.to_csv(os.path.join(rep_dir,'bpm_results.csv'), index=False)
                
    collect_results(args.output, args.output)

    return None
    

def get_contact_score(structure_file, trajectory_file, lig_resname):
    """A function the gets the ContactScore from an OpenBPMD trajectory.

    Parameters
    ----------
    structure_file : str
        The name of the centred equilibrated system PDB file that 
        was used to start the OpenBPMD simulation.
    trajectory_file : str
        The name of the OpenBPMD trajectory file.
    lig_resname : str
        Residue name of the ligand that was biased.

    Returns
    -------
    contact_scores : np.array 
        ContactScore for every frame of the trajectory.
    """
    u = mda.Universe(structure_file, trajectory_file)

    sel_donor = f"resname {lig_resname} and not name *H*"
    sel_acceptor = f"protein and not name H* and \
                     around 5 resname {lig_resname}"

    # reference groups (first frame of the trajectory, but you could also use
    # a separate PDB, eg crystal structure)
    a_donors = u.select_atoms(sel_donor)
    a_acceptors = u.select_atoms(sel_acceptor)

    cont_analysis = contacts.Contacts(u, select=(sel_donor, sel_acceptor),
                                      refgroup=(a_donors, a_acceptors),
                                      radius=3.5)

    cont_analysis.run()
    # print number of average contacts in the first ns
    # NOTE - hard coded number of frames (100 per traj)
    frame_idx_first_ns = int(len(cont_analysis.timeseries)/10)
    first_ns_mean = np.mean(cont_analysis.timeseries[1:frame_idx_first_ns, 1])
    if first_ns_mean == 0:
        normed_contacts = cont_analysis.timeseries[1:, 1]
    else:
        normed_contacts = cont_analysis.timeseries[1:, 1]/first_ns_mean
    contact_scores = np.where(normed_contacts > 1, 1, normed_contacts)

    return contact_scores


def get_pose_score(structure_file, trajectory_file, lig_resname):
    """A function the gets the PoseScore (ligand RMSD) from an OpenBPMD
    trajectory.

    Parameters
    ----------
    'structure_file : str
        The name of the centred equilibrated system
        PDB file that was used to start the OpenBPMD simulation.
    trajectory_file : str
        The name of the OpenBPMD trajectory file.
    lig_resname : str
        Residue name of the ligand that was biased.

    Returns
    -------
    pose_scores : np.array 
        PoseScore for every frame of the trajectory.
    """
    # Load a MDA universe with the trajectory
    u = mda.Universe(structure_file, trajectory_file)
    # Align each frame using the backbone as reference
    # Calculate the RMSD of ligand heavy atoms
    r = rms.RMSD(u, select='backbone',
                 groupselections=[f'resname {lig_resname} and not name H*'],
                 ref_frame=0).run()
    # Get the PoseScores as np.array
    pose_scores = r.rmsd[1:, -1]

    return pose_scores


def minimize(parm, input_positions, out_dir, min_file_name):
    """An energy minimization function down with an energy tolerance
    of 10 kJ/mol.

    Parameters
    ----------
    parm : Parmed or OpenMM parameter file object
        Used to create the OpenMM System object.
    input_positions : OpenMM Quantity
        3D coordinates of the equilibrated system.
    out_dir : str
        Directory to write the outputs.
    min_file_name : str
        Name of the minimized PDB file to write.
    """
    system = parm.createSystem(
        nonbondedMethod=PME,
        nonbondedCutoff=1*nanometers,
        constraints=HBonds,
    )

    # Define platform properties
    platform = Platform.getPlatformByName('CUDA')
    properties = {'CudaPrecision': 'mixed'}

    # Set up the simulation parameters
    # Langevin integrator at 300 K w/ 1 ps^-1 friction coefficient
    # and a 2-fs timestep
    # NOTE - no dynamics performed, but required for setting up
    # the OpenMM system.
    integrator = LangevinIntegrator(300*kelvin, 1/picosecond,
                                    0.002*picoseconds)
    simulation = Simulation(parm.topology, system, integrator, platform,
                            properties)
    simulation.context.setPositions(input_positions)

    # Minimize the system - no predefined number of steps
    simulation.minimizeEnergy()

    # Write out the minimized system to use w/ MDAnalysis
    positions = simulation.context.getState(getPositions=True).getPositions()
    out_file = os.path.join(out_dir,min_file_name)
    PDBFile.writeFile(simulation.topology, positions,
                      open(out_file, 'w'))

    return None


def equilibrate(min_pdb, parm, out_dir, eq_file_name, parm_file : str = None, npt : bool = False):
    """A function that does a 500 ps NVT equilibration with position
    restraints, with a 5 kcal/mol/A**2 harmonic constant on solute heavy
    atoms, using a 2 fs timestep.

    Parameters
    ----------
    min_pdb : str
        Name of the minimized PDB file.
    parm : Parmed or OpenMM parameter file object
        Used to create the OpenMM System object.
    out_dir : str
        Directory to write the outputs to.
    eq_file_name : str
        Name of the equilibrated PDB file to write.
    """
    RESTRAINT_K      = 5 * kilocalories_per_mole / angstrom**2
    EQUIL_STEPS      = 250_000   # 500 ps @ 2 fs/step
    NONSOLVENT       = {'WAT', 'SOL', 'HOH', 'CL', 'NA'}
    PROTOCOL         = "NPT" if npt else "NVT"
    PROTOCOL         = "NPT" if npt else "NVT"

    system = parm.createSystem(
        nonbondedMethod=PME,
        nonbondedCutoff=1*nanometers,
        constraints=HBonds,
    )

    # Add the harmonic restraints on the positions
    # of specified atoms
    restraint = CustomExternalForce('k*periodicdistance(x, y, z, x0, y0, z0)^2')
    restraint.addGlobalParameter('k', RESTRAINT_K)
    restraint.addPerParticleParameter('x0')
    restraint.addPerParticleParameter('y0')
    restraint.addPerParticleParameter('z0')

    system.addForce(restraint)
    
    pdb = PDBFile(min_pdb)
    input_positions = pdb.getPositions()

    for atom in pdb.topology.atoms():
        if atom.residue.name not in NONSOLVENT and atom.element.symbol != 'H':
            restraint.addParticle(atom.index, pdb.positions[atom.index])

    if npt:
        system.addForce(MonteCarloBarostat(1.0 * atmospheres, 300 * kelvin, 25))

    integrator = LangevinIntegrator(300*kelvin, 1/picosecond,
                                    0.002*picoseconds)
    
    platform = Platform.getPlatformByName('CUDA')
    properties = {'CudaPrecision': 'mixed'}
    
    sim = Simulation(parm.topology, system, integrator,
                     platform, properties)
    sim.context.setPositions(input_positions)
    sim.step(EQUIL_STEPS)  # run 500 ps of equilibration
    
    all_positions = sim.context.getState(
        getPositions=True, 
        enforcePeriodicBox=True
    ).getPositions()
    
    # write the positions of atoms
    out_file = os.path.join(out_dir,eq_file_name)
    PDBFile.writeFile(sim.topology, all_positions,
                      open(out_file, 'w'))

    logger.info(f"{PROTOCOL} PDB file saved -> {out_file}")

    
    # ── 6. Write RST7 with ParmEd ───────────────────────────────────────────
    struct = pmd.load_file(parm_file, hasbox = True)
    
    for i, atom in enumerate(struct.atoms):
        pos = all_positions[i].value_in_unit(angstrom)
        atom.xx = pos.x  # nm → Å
        atom.xy = pos.y
        atom.xz = pos.z
    
    box_vectors = sim.context.getState().getPeriodicBoxVectors()
    a, b, c = [v.value_in_unit(angstrom) for v in box_vectors]
    struct.box = [
        np.linalg.norm(a),
        np.linalg.norm(b),
        np.linalg.norm(c),
        np.degrees(np.arccos(np.dot(b, c)/(np.linalg.norm(b)*np.linalg.norm(c)))),
        np.degrees(np.arccos(np.dot(a, c)/(np.linalg.norm(a)*np.linalg.norm(c)))),
        np.degrees(np.arccos(np.dot(a, b)/(np.linalg.norm(a)*np.linalg.norm(b))))
    ]
        
    eq_file_name = eq_file_name.split(".")[0]
    out_rst = os.path.join(out_dir, eq_file_name + '.rst7')
    out_prm = os.path.join(out_dir, eq_file_name + '.prmtop')  # solo hace falta una vez

    struct.save(out_rst, overwrite=True)
    logger.info(f"{PROTOCOL} RST7 file saved -> {out_rst}")
    
    struct.save(out_prm, overwrite=True)
    logger.info(f"{PROTOCOL} prmtop file saved -> {out_prm}")

    return None


def GCMC(eq_rst7, eq_prmtop, out_dir, ligand_name):
    
    # Store the ligand index.
    lig = ligand_name
    
    # Load the scytalone dehydratase system.
    mols = sr.load(eq_prmtop, eq_rst7, show_warnings=True)
    reference = f"resname {lig}"

    ref_atoms = []

    for i, atom in enumerate(mols[lig].atoms()):
        if atom.element().name() == "Carbon":
            ref_atoms.append({
                "name": atom.name().value(),
                "resname": atom.residue().name().value(),
                "resid": str(atom.residue().number().value())
            })

    logger.info(f"GCMC Reference atoms: {ref_atoms}")

    # Create a GCMC sampler.
    sampler = GCMCSampler(
        mols,
        reference=reference,
        batch_size=1000,
        num_attempts=10000,
        cutoff_type="pme",
        cutoff="10 A",
        radius="6.0 A", #4.2 A
        temperature="300 K",
        num_ghost_waters=100,
        bulk_sampling_probability=0.1,
        ghost_file=os.path.join(out_dir,f"ghosts_{lig}.txt"),
        log_file=os.path.join(out_dir,f"gcmc_{lig}.txt"),
        log_level="info",
        platform="auto",
        overwrite=True,
    )

    # Create a dynamics object using the modified GCMC system.
    # This contains a number of ghost waters that can be used
    # for insertion moves.
    d = sampler.system().dynamics(
        cutoff_type="pme",
        cutoff="10 A",
        temperature="300 K",
        integrator="langevin_middle",
        pressure=None,
        constraint="h_bonds",
        timestep="2 fs",
        platform="auto",
    )
    d.randomise_velocities()

    # Bind the GCMC sampler to the dynamics object. This ensures that the GCMC
    # water state can be reset following a crash.
    sampler.bind_dynamics(d)

    # Delete any existing waters from the GCMC region.
    sampler.delete_waters(d.context())

    # Store the frame frequency.
    frame_frequency = 5

    # Store the radius.
    radius = sampler._radius.value()

    # Perform initial GCMC equilibration on the initial structure.
    logger.info("Equilibrating the system with GCMC moves...")
    for i in range(100):
        sampler.move(d.context())

    # Run 1ns dynamics with GCMC moves every 1ps.
    logger.info("Running 1ns of dynamics with GCMC moves...")
    for i in range(1000):
        # Run 1ps of dynamics.
        d.run("1ps", energy_frequency="5ps", frame_frequency="5ps")

        # Perform a GCMC move.
        moves = sampler.move(d.context())

        # If we hit the frame frequency, then save the current ghost residue indices.
        if i > 0 and (i + 1) % frame_frequency == 0:
            sampler.write_ghost_residues()

        if i % 10 == 0:
            logger.info(
                f"Cycle {i}, N = {sampler.num_waters()}, "
                f"insertions = {sampler.num_insertions()}, "
                f"deletions = {sampler.num_deletions()}"
            )
            logger.info(
                f"Current potential energy: {d.current_potential_energy().value():.3f} kcal/mol"
            )

    # Save the trajectory.
    mols = d.commit()
    sr.save(mols, os.path.join(out_dir,f"sd_{lig}_final.pdb"))
    sr.save(mols, os.path.join(out_dir,f"sd_{lig}_final.prmtop"))
    sr.save(mols, os.path.join(out_dir,f"sd_{lig}_final.rst7"))
    sr.save(mols.trajectory(), os.path.join(out_dir,f"sd_{lig}.dcd"))

    logger.info(f"Insertions: {sampler.num_insertions()}")
    logger.info(f"Deletions: {sampler.num_deletions()}")
    logger.info(f"Move acceptance probability: {sampler.move_acceptance_probability():.4f}")
    logger.info(f"Attempt acceptance probability: {sampler.attempt_acceptance_probability():.4f}")
    logger.info(f"GCMC Analysis completed! -> {os.path.join(out_dir,f"sd_{lig}_final.pdb")}")


def produce(out_dir, idx, lig_resname, eq_pdb, parm, parm_file,
            set_hill_height):
    """An OpenBPMD production simulation function. Ligand RMSD is biased with
    metadynamics. The integrator uses a 4 fs time step and
    runs for 10 ns, writing a frame every 100 ps.

    Writes a 'trj.dcd', 'COLVAR.npy', 'bias_*.npy' and 'sim_log.csv' files
    during the metadynamics simulation in the '{out_dir}/rep_{idx}' directory.
    After the simulation is done, it analyses the trajectories and writes a
    'bpm_results.csv' file with time-resolved PoseScore and ContactScore.

    Parameters
    ----------
    out_dir : str
        Directory where your equilibration PDBs and 'rep_*' dirs are at.
    idx : int
        Current replica index.
    lig_resname : str
        Residue name of the ligand.
    eq_pdb : str
        Name of the PDB for equilibrated system.
    parm : Parmed or OpenMM parameter file object
        Used to create the OpenMM System object.
    parm_file : str
        The name of the parameter or topology file of the system.
    set_hill_height : float
        Metadynamic hill height, in kcal/mol.

    """
    # First, assign the replica directory to which we'll write the files
    write_dir = os.path.join(out_dir,f'rep_{idx}')
    # Get the anchor atoms by ...
    universe = mda.Universe(eq_pdb,
                            format='XPDB', in_memory=True)
    # ... finding the protein's COM ...
    prot_com = universe.select_atoms('protein').center_of_mass()
    x, y, z = prot_com[0], prot_com[1], prot_com[2]
    # ... and taking the heavy backbone atoms within 10A of the COM
    sel_str = f'point {x} {y} {z} 10 and backbone and not name H*'
    anchor_atoms = universe.select_atoms(sel_str)
    if len(anchor_atoms) == 0:
        raise ValueError('No Calpha atoms found within 10 ang of the center of mass of the protein. \
                Check your input files.')
    anchor_atom_idx = anchor_atoms.indices.tolist()

    # Get indices of ligand heavy atoms
    lig = universe.select_atoms(f'resname {lig_resname} and not name H*')
    if len(lig) == 0:
        raise ValueError(f"Ligand with resname '{lig_resname}' not found.")

    lig_ha_idx = lig.indices.tolist()

    # Set up the system to run metadynamics
    system = parm.createSystem(
        nonbondedMethod=PME,
        nonbondedCutoff=1*nanometers,
        constraints=HBonds,
        hydrogenMass=4*amu
    )
    # get the atom positions for the system from the equilibrated
    # system
    input_positions = PDBFile(eq_pdb).getPositions()

    # Add an 'empty' flat-bottom restraint to fix the issue with PBC.
    # Without one, RMSDForce object fails to account for PBC.
    k = 0*kilojoules_per_mole  # NOTE - 0 kJ/mol constant
    upper_wall = 10.00*nanometer
    fb_eq = '(k/2)*max(distance(g1,g2) - upper_wall, 0)^2'
    upper_wall_rest = CustomCentroidBondForce(2, fb_eq)
    upper_wall_rest.addGroup(lig_ha_idx)
    upper_wall_rest.addGroup(anchor_atom_idx)
    upper_wall_rest.addBond([0, 1])
    upper_wall_rest.addGlobalParameter('k', k)
    upper_wall_rest.addGlobalParameter('upper_wall', upper_wall)
    upper_wall_rest.setUsesPeriodicBoundaryConditions(True)
    system.addForce(upper_wall_rest)

    alignment_indices = lig_ha_idx + anchor_atom_idx

    rmsd = RMSDForce(input_positions, alignment_indices)
    # Set up the typical metadynamics parameters
    grid_min, grid_max = 0.0, 1.0  # nm
    hill_height = set_hill_height*kilocalories_per_mole
    hill_width = 0.002  # nm, also known as sigma

    grid_width = hill_width / 5
    # 'grid' here refers to the number of grid points
    grid = int(abs(grid_min - grid_max) / grid_width)

    rmsd_cv = BiasVariable(rmsd, grid_min, grid_max, hill_width,
                           False, gridWidth=grid)

    # define the metadynamics object
    # deposit bias every 1 ps, BF = 4, write bias every ns
    meta = Metadynamics(system, [rmsd_cv], 300.0*kelvin, 4.0, hill_height,
                        250, biasDir=write_dir,
                        saveFrequency=250000)

    # Set up and run metadynamics
    integrator = LangevinIntegrator(300*kelvin, 1.0/picosecond,
                                    0.004*picoseconds)
    platform = Platform.getPlatformByName('CUDA')
    properties = {'CudaPrecision': 'mixed'}

    simulation = Simulation(parm.topology, system, integrator, platform,
                            properties)
    simulation.context.setPositions(input_positions)

    trj_name = os.path.join(write_dir,'trj.dcd')

    sim_time = 10  # ns
    steps = 250000 * sim_time

    simulation.reporters.append(DCDReporter(trj_name, 25000))  # every 100 ps
    simulation.reporters.append(StateDataReporter(
                                os.path.join(write_dir,'sim_log.csv'), 250000,
                                step=True, temperature=True, progress=True,
                                remainingTime=True, speed=True,
                                totalSteps=steps, separator=','))  # every 1 ns

    colvar_array = np.array([meta.getCollectiveVariables(simulation)])
    for i in range(0, int(steps), 500):
        if i % 25000 == 0:
            # log the stored COLVAR every 100ps
            np.save(os.path.join(write_dir,'COLVAR.npy'), colvar_array)
        meta.step(simulation, 500)
        current_cvs = meta.getCollectiveVariables(simulation)
        # record the CVs every 2 ps
        colvar_array = np.append(colvar_array, [current_cvs], axis=0)
    np.save(os.path.join(write_dir,'COLVAR.npy'), colvar_array)

    # center everything using MDTraj, to fix any PBC imaging issues
    # mdtraj can't use GMX TOP, so we have to specify the GRO file instead
    
    mdu = md.load(trj_name, top=parm_file)
    mdu.image_molecules(inplace=True)
    mdu.save(trj_name)

    return None


def collect_results(in_dir, out_dir):
    """A function that collects the time-resolved BPM results,
    takes the scores from last 2 ns of the simulation, averages them
    and writes that average as the final score for a given pose.

    Writes a 'results.csv' file in 'out_dir' directory.
    
    Parameters
    ----------
    in_dir : str
        Directory with 'rep_*' directories.
    out_dir : str
        Directory where the 'results.csv' file will be written
    """
    compList = []
    contactList = []
    poseList = []
    # find how many repeats have been run
    glob_str = os.path.join(in_dir,'rep_*')
    nreps = len(glob.glob(glob_str))
    for idx in range(0, nreps):
        f = os.path.join(in_dir,f'rep_{idx}','bpm_results.csv')
        df = pd.read_csv(f)
        # Since we only want last 2 ns, get the index of
        # the last 20% of the data points
        last_2ns_idx = round(len(df['CompScore'].values)/5)  # round up
        compList.append(df['CompScore'].values[-last_2ns_idx:])
        contactList.append(df['ContactScore'].values[-last_2ns_idx:])
        poseList.append(df['PoseScore'].values[-last_2ns_idx:])

    # Get the means of the last 2 ns
    meanCompScore = np.mean(compList)
    meanPoseScore = np.mean(poseList)
    meanContact = np.mean(contactList)
    # Get the standard deviation of the final 2 ns
    meanCompScore_std = np.std(compList)
    meanPoseScore_std = np.std(poseList)
    meanContact_std = np.std(contactList)
    # Format it the Pandas way
    d = {'CompScore': [meanCompScore], 'CompScoreSD': [meanCompScore_std],
         'PoseScore': [meanPoseScore], 'PoseScoreSD': [meanPoseScore_std],
         'ContactScore': [meanContact], 'ContactScoreSD': [meanContact_std]}

    results_df = pd.DataFrame(data=d)
    results_df = results_df.round(3)
    results_df.to_csv(os.path.join(out_dir,'results.csv'), index=False)


def path_to_files(name, outdir):

    parm_prmtop = os.path.join(outdir, f"{name}.prmtop") 
    structure_rst7 = os.path.join(outdir, f"{name}.rst7") 
    structure_pdb = os.path.join(outdir, f"{name}.pdb") 

    return parm_prmtop, structure_rst7, structure_pdb


if __name__ == "__main__":
    """ This is executed when run from the command line """
    # Parse the CLI arguments
    parser = argparse.ArgumentParser(
        formatter_class=argparse.RawDescriptionHelpFormatter,
        description=descriptn)

    parser.add_argument("-s", "--structure", type=str, default='system.gro',
                        help='input structure file name (default: %(default)s)')
    parser.add_argument("-p", "--parameters", type=str, default='system.top',
                        help='input topology file name (default: %(default)s)')
    parser.add_argument("-o", "--output", type=str, default='.',
                        help='output location (default: %(default)s)')
    parser.add_argument("-lig_resname", type=str, default='MOL',
                        help='the name of the ligand (default: %(default)s)')
    parser.add_argument("-nreps", type=int, default=10,
                        help="number of OpenBPMD repeats (default: %(default)i)")
    parser.add_argument("-hill_height", type=float, default=0.05,
                        help="the hill height in kcal/mol (default: %(default)f)")
    parser.add_argument("--npt", action="store_true", 
                        help="Enable NPT")
    parser.add_argument("--no-npt", dest="npt", action="store_false", 
                        help="Disable NPT")
    parser.add_argument("--gcmc", action="store_true", 
                        help="Enable GCMC")
    parser.add_argument("--no-gcmc", dest="gcmc", action="store_false", 
                        help="Disable GCMC")
    args = parser.parse_args()
    main(args)

[lohedges]

No problem. I am away for a few weeks and will take a look when I'm back. If you are able to provide your inputs and the command-line used, that would be great.

[MartinFuentetaja]

Of course! The .gro and .top files are inside .zip folder. -> Inputs.zip

Regarding the command-line used:

• NVT only:

python MY_openbpmd.py \
    -s INPUT.gro \
    -p topol.top \
    -o ./results \
    -lig_resname "JI1" \
    --no-npt \
    --no-gcmc

• NVT + NPT:

python MY_openbpmd.py \
    -s INPUT.gro \
    -p topol.top \
    -o ./results \
    -lig_resname "JI1" \
    --npt \
    --no-gcmc

• NVT + GCMC:

python MY_openbpmd.py \
    -s INPUT.gro \
    -p topol.top \
    -o ./results \
    -lig_resname "JI1" \
    --no-npt \
    --gcmc

• NVT + NPT + GCMC:

python MY_openbpmd.py \
    -s INPUT.gro \
    -p topol.top \
    -o ./results \
    -lig_resname "JI1" \
    --npt \
    --gcmc

The script fails when --gcmc option is activated. My guess is that it is related with how I am using parmed, mdtraj and mdanalysis.

[lohedges]

Thanks. I am away until the 15th but will try to reproduce and debug the issue when I'm back.

[MartinFuentetaja]

Okey! No problem. I have not mentioned, but my mamba enviroment is the following one:

enviroment.zip

[lohedges]

I've run your code and see the same error. However, if I add save_velocities=False to the sr.save() call for the RST7 file, then it writes just fine. However, the simulation then crashes at the next stage, i.e.:

2026-06-15 12:34:39.557 | INFO     | __main__:GCMC:645 - Insertions: 146
2026-06-15 12:34:39.557 | INFO     | __main__:GCMC:646 - Deletions: 136
2026-06-15 12:34:39.557 | INFO     | __main__:GCMC:647 - Move acceptance probability: 1.4100
2026-06-15 12:34:39.557 | INFO     | __main__:GCMC:650 - Attempt acceptance probability: 0.0044
2026-06-15 12:34:39.557 | INFO     | __main__:GCMC:653 - GCMC Analysis completed! -> ./results/sd_JI1_final.pdb
2026-06-15 12:34:40.476 | INFO     | __main__:main:181 - GCMC MDAnalysis completed!
2026-06-15 12:34:40.476 | INFO     | __main__:main:197 - Initializing BPMD rep_0...
INFO:MDAnalysis.topology.PDBParser:Setting segids from chainIDs because no segids found in the PDB file.
/home/lester/Code/openbiosim/sire/.pixi/envs/dev/lib/python3.12/site-packages/MDAnalysis/core/universe.py:939: UserWarning: Reader has no dt information, set to 1.0 ps
  dt=self.trajectory.ts.dt * step,
INFO:MDAnalysis.core.universe:The attribute(s) types have already been read from the topology file. The guesser will only guess empty values for this attribute, if any exists. To overwrite it by completely guessed values, you can pass the attribute to the force_guess parameter instead of the to_guess one
INFO:MDAnalysis.guesser.base:There is no empty types values. Guesser did not guess any new values for types attribute
INFO:MDAnalysis.core.universe:attribute masses has been guessed successfully.
╭─────────────────────────────── Traceback (most recent call last) ────────────────────────────────╮
│ /home/lester/Code/openbiosim/loch/position_overflow/script.py:928 in <module>                    │
│                                                                                                  │
│   925 │   │   "--no-gcmc", dest="gcmc", action="store_false", help="Disable GCMC"                │
│   926 │   )                                                                                      │
│   927 │   args = parser.parse_args()                                                             │
│ ❱ 928 │   main(args)                                                                             │
│   929                                                                                            │
│                                                                                                  │
│ /home/lester/Code/openbiosim/loch/position_overflow/script.py:199 in main                        │
│                                                                                                  │
│   196 │   │                                                                                      │
│   197 │   │   logger.info(f"Initializing BPMD rep_{idx}...")                                     │
│   198 │   │                                                                                      │
│ ❱ 199 │   │   produce(                                                                           │
│   200 │   │   │   args.output,                                                                   │
│   201 │   │   │   idx,                                                                           │
│   202 │   │   │   args.lig_resname,                                                              │
│                                                                                                  │
│ /home/lester/Code/openbiosim/loch/position_overflow/script.py:795 in produce                     │
│                                                                                                  │
│   792 │   │   if i % 25000 == 0:                                                                 │
│   793 │   │   │   # log the stored COLVAR every 100ps                                            │
│   794 │   │   │   np.save(os.path.join(write_dir, "COLVAR.npy"), colvar_array)                   │
│ ❱ 795 │   │   meta.step(simulation, 500)                                                         │
│   796 │   │   current_cvs = meta.getCollectiveVariables(simulation)                              │
│   797 │   │   # record the CVs every 2 ps                                                        │
│   798 │   │   colvar_array = np.append(colvar_array, [current_cvs], axis=0)                      │
│                                                                                                  │
│ /home/lester/Code/openbiosim/sire/.pixi/envs/dev/lib/python3.12/site-packages/openmm/app/metadyn │
│ amics.py:170 in step                                                                             │
│                                                                                                  │
│   167 │   │   │   │   nextSteps = min(nextSteps, self.frequency)                                 │
│   168 │   │   │   else:                                                                          │
│   169 │   │   │   │   nextSteps = min(nextSteps, self.frequency - simulation.currentStep % sel   │
│ ❱ 170 │   │   │   simulation.step(nextSteps)                                                     │
│   171 │   │   │   if simulation.currentStep % self.frequency == 0:                               │
│   172 │   │   │   │   position = self._force.getCollectiveVariableValues(simulation.context)     │
│   173 │   │   │   │   energy = simulation.context.getState(energy=True, groups={forceGroup}).g   │
│                                                                                                  │
│ /home/lester/Code/openbiosim/sire/.pixi/envs/dev/lib/python3.12/site-packages/openmm/app/simulat │
│ ion.py:149 in step                                                                               │
│                                                                                                  │
│   146 │   │   if int(steps) != steps:                                                            │
│   147 │   │   │   raise TypeError(f'Expected an integer for steps, got {type(steps)}')           │
│   148 │   │                                                                                      │
│ ❱ 149 │   │   self._simulate(endStep=self.currentStep+steps)                                     │
│   150 │                                                                                          │
│   151 │   def runForClockTime(self, time, checkpointFile=None, stateFile=None, checkpointInter   │
│   152 │   │   """Advance the simulation by integrating time steps until a fixed amount of cloc   │
│                                                                                                  │
│ /home/lester/Code/openbiosim/sire/.pixi/envs/dev/lib/python3.12/site-packages/openmm/app/simulat │
│ ion.py:228 in _simulate                                                                          │
│                                                                                                  │
│   225 │   │   │   │   │   anyReport = True                                                       │
│   226 │   │   │   stepsToGo = nextSteps                                                          │
│   227 │   │   │   while stepsToGo > 10:                                                          │
│ ❱ 228 │   │   │   │   self.integrator.step(10) # Only take 10 steps at a time, to give Python    │
│   229 │   │   │   │   stepsToGo -= 10                                                            │
│   230 │   │   │   │   if endTime is not None and datetime.now() >= endTime:                      │
│   231 │   │   │   │   │   return                                                                 │
│                                                                                                  │
│ /home/lester/Code/openbiosim/sire/.pixi/envs/dev/lib/python3.12/site-packages/openmm/openmm.py:1 │
│ 6460 in step                                                                                     │
│                                                                                                  │
│   16457 │   │   steps : int                                                                      │
│   16458 │   │   │   the number of time steps to take                                             │
│   16459 │   │   """                                                                              │
│ ❱ 16460 │   │   return _openmm.LangevinMiddleIntegrator_step(self, steps)                        │
│   16461 │                                                                                        │
│   16462 │   def __init__(self, *args):                                                           │
│   16463 │   │   r"""                                                                             │
╰──────────────────────────────────────────────────────────────────────────────────────────────────╯
OpenMMException: NaN encountered during RMSD force calculation

[MartinFuentetaja]

Hello @lohedges !

Firstly, thank you for taking your time. Yes, I have encountered the same error. Instead of saving the .rst7 file, I saved the .pdb, but it raised the same error. Interestingly, If you unselect gcmc option (--no-gcmc), it works correctly. Therefore, I guess that something is happening at GCMC step...

[lohedges]

This is because the GCMC prmtop that is written doesn't know about the ghost waters (these are only in the OpenMM context, not the Sire topology). When you start the production simulation, the topology mismatches with the water coordinates in the PDB file, so you will likely be turning on full interactions for many water molecules that overlap with other atoms, hence the blow up. What you'll need to do is loop over the water molecules and update their force field parameters in Sire so that they match the correct state, i.e. real or ghost. I'll add a convenience function to the code to do this for you.

[lohedges]

I've added a fix in #26. Once merged, you can use an updated version of your script.

[MartinFuentetaja]

THANK YOU!!!! It works perfectly now! 👍

[lohedges]

No problem. Thanks for reporting, and for confirming the fix.