Issue setting positions from AMBER rst7 file

[althea-hansel]

Hello,

I am trying to run grand on a protein-ligand complex I have previously prepared with openMM. I am trying to import the prmtop and rst7 files, but get the following error when trying to set positions from the rst7:

INFO:grand.samplers:StandardGCMCSphereSampler object initialised
Traceback (most recent call last):
  File "/home/althea/Documents/IMP2_project/e1_allsites/rrm12/input/../../test_grand.py", line 43, in <module>
    simulation.context.setPositions(inpcrd.positions)
  File "/home/althea/programs/miniconda3/envs/openmm/lib/python3.9/site-packages/openmm/openmm.py", line 3681, in setPositions
    return _openmm.Context_setPositions(self, positions)
openmm.OpenMMException: Called setPositions() on a Context with the wrong number of positions

Would you be able to assist me in debugging this issue? I can provide the input prmtop and rst7 if needed.
My script is below:

#!/usr/bin/env python

from openmm.app import *
from openmm import *
from openmm.unit import *
from sys import stdout

from openmmtools.integrators import BAOABIntegrator

import grand

if __name__ == "__main__":
    prmtop = AmberPrmtopFile('e1_TYR154_0_system.prmtop')
    inpcrd = AmberInpcrdFile('e1_TYR154_0_system.rst7')

    prmtop.topology, inpcrd.positions, ghosts = grand.utils.add_ghosts(prmtop.topology, inpcrd.positions, n=10, pdb='sd-ghosts.pdb')

    system = prmtop.createSystem(nonbondedMethod=PME, nonbondedCutoff=1.2*nanometer, constraints=HBonds)

    # Define reference atoms around which the GCMC sphere is based
    ref_atoms = [{'name': 'OH', 'resname': 'TYR', 'resid': '154'},
                #{'name': 'OH', 'resname': 'TYR', 'resid': '43'}
                ]

    # Create GCMC Sampler object
    gcmc_mover = grand.samplers.StandardGCMCSphereSampler(system=system,
                                                        topology=prmtop.topology,
                                                        temperature=300*kelvin,
                                                        referenceAtoms=ref_atoms,
                                                        sphereRadius=4*angstroms,
                                                        log='sd-gcmc.log',
                                                        dcd='sd-raw.dcd',
                                                        rst='sd-gcmc.rst7',
                                                        overwrite=False)

    # BAOAB Langevin integrator (important)
    integrator = BAOABIntegrator(300*kelvin, 1.0/picosecond, 0.002*picoseconds)

    platform = Platform.getPlatformByName('CUDA')
    platform.setPropertyDefaultValue('Precision', 'mixed')

    simulation = Simulation(prmtop.topology, system, integrator, platform)
    simulation.context.setPositions(inpcrd.positions)
    simulation.context.setVelocitiesToTemperature(300*kelvin)
    simulation.context.setPeriodicBoxVectors(*prmtop.topology.getPeriodicBoxVectors())

    # Switch off ghost waters and in sphere
    gcmc_mover.initialise(simulation.context, ghosts)
    gcmc_mover.deleteWatersInGCMCSphere()

    # Equilibrate water distribution - 10k moves over 5 ps
    print("Equilibration...")
    for i in range(50):
        # Carry out 200 moves every 100 fs
        gcmc_mover.move(simulation.context, 200)
        simulation.step(50)
    print("{}/{} equilibration GCMC moves accepted. N = {}".format(gcmc_mover.n_accepted,
                                                                gcmc_mover.n_moves,
                                                                gcmc_mover.N))

    # Add StateDataReporter for production
    simulation.reporters.append(StateDataReporter(stdout,
                                                1000,
                                                step=True,
                                                potentialEnergy=True,
                                                temperature=True,
                                                volume=True))
    # Reset GCMC statistics
    gcmc_mover.reset()

    # Run simulation - 5k moves over 50 ps
    print("\nProduction")
    for i in range(50):
        # Carry out 100 GCMC moves per 1 ps of MD
        simulation.step(500)
        gcmc_mover.move(simulation.context, 100)
        # Write data out
        gcmc_mover.report(simulation)

    #
    # Need to process the trajectory for visualisation
    #

    # Move ghost waters out of the simulation cell
    trj = grand.utils.shift_ghost_waters(ghost_file='gcmc-ghost-wats.txt',
                                        topology='sd-ghosts.pdb',
                                        trajectory='sd-raw.dcd')

    # Recentre the trajectory on a particular residue
    trj = grand.utils.recentre_traj(t=trj, resname='TYR', resid=154)

    # Align the trajectory to the protein
    grand.utils.align_traj(t=trj, output='sd-gcmc.dcd')

    # Write out a trajectory of the GCMC sphere
    grand.utils.write_sphere_traj(radius=4.0,
                                ref_atoms=ref_atoms,
                                topology='sd-ghosts.pdb',
                                trajectory='sd-gcmc.dcd',
                                output='gcmc_sphere.pdb',
                                initial_frame=True)

[olliemelling]

Hi,

Please could you send the .inpcrd and .prmtop files and I'll look into what's causing this error.

Many thanks,

Ollie

[althea-hansel]

grand_system.tgz
This should be them

[olliemelling]

Thanks for sending your input files, the issue is that Amber input files (.prmtop, .inpcrd, .rst7) are difficult to manipulate in OpenMM so when the ghost waters are added to the topology, this information isn't retained when the system is created as OpenMM re-reads the original .prmtop file - meaning the new water molecules aren't in the system, hence the error when you set the positions.

The easiest solution here is to use a PDB file and a .xml file for your ligand and then to load them as in the examples - creating a PDBFile object and a ForceField object.

The ambpdb command within the ambertools package can convert your .prmtop and .inpcrd files into a PDB file with the following command:

ambpdb -p file.prmtop -c file.inpcrd > output.pdb

There are functions within grand.utils to create an .xml file for the ligand and add CONECT lines to your PDB if needed, parmed should also be able to do this as well.

Let me know if you have any further issues.

[althea-hansel]

Thank you for these tips, I will try this. I did also try running the same grand code on a system built from a PDB and RDKit ligand (not via amber input files) and had the same error/issue.

[olliemelling]

Would you be able to send input files and script for the PDB + ligand system that you tried previously?

A protein-ligand system loaded as a PDB and XML file should work fine so would need to see the script to know what's going wrong here.

[noahharrison64]

Are there any other options for providing input to grand, rather than PDBFile (with conect records) + ligand XML?
I've been trying to use .gro and .top input files but have also been receiving the error:
'openmm.OpenMMException: Called setPositions() on a Context with the wrong number of positions'

Is this likely for the same reason? It's a bit of a faff having to try to convert my systems to a XML, a .prepi and a .pdb file. Mostly because .prepi is a pretty old and unsupported format.

Input files incase you want to try to recreate
input.zip