essex-lab / grand Goto Github PK

View Code? Open in Web Editor NEW

53.0 5.0 19.0 1.88 MB

A Python module for carrying out GCMC insertions and deletions of water molecules in OpenMM.

Home Page: https://grand.readthedocs.io

License: MIT License

Python 100.00%

gcmc water openmm monte-carlo molecular-dynamics water-molecules python

grand's People

Contributors

Stargazers

Watchers

Forkers

chebuu msoroush layeqa pk-organics biocheming bcossins biochunan qyuan7 anthonialcaraz byun-jinyoung janikkokot hong-rui dgazgalis vickie02736 warrenheww kpark97

grand's Issues

Expand to other fragments?

Hi, can this module easily expand to other fragments?
I want to use GCMC for tasks like SILCS.

Thank you!

Particle coordinate is NaN when using GCNCMC with OpenMM 8.0.0

Hi!

I was running pytests with OpenMM 8.0.0 with a few tweaks (chaged simtk imports to openmm). Tests for NonequilibriumGCMC samplers fail with openmm.OpenMMException: Particle coordinate is NaN. The earliest OpenMM version I was able to get it working with was OpenMM 7.7.0

I've looked into code and OpenMM changelogs and could not figure out the reason behind such behavior. Even a single step of NonequilibriumLangevinIntegrator leads to explosions at both insertion and deletion moves.

Are you planning to provide support for newer OpenMM versions? Folks at @cresset-group have a fork with various enhancements and OpenMM 8.0.0 compatibility, but I was not able to get it working either.

Issue setting positions from AMBER rst7 file

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)

No registered platform called "CUDA"

Hello, I encountered the following error after installation while running the example.
It looks like related to CUDA, any chance you have encoutered wittth this before?

Traceback (most recent call last):
  File "equil-uvt1.py", line 60, in <module>
    platform = Platform.getPlatformByName('CUDA')
  File "/opt/conda/envs/grand/lib/python3.7/site-packages/simtk/openmm/openmm.py", line 17865, in getPlatformByName
    return _openmm.Platform_getPlatformByName(name)
Exception: There is no registered Platform called "CUDA"

GCMC simulation explode with "Particle coordinate is nan" or "Energy is NaN"

Hi developers,

Recently we have tried the grand gcmc for sampling waters, and I installed the latest grand from the master branch and the OpenMM installed from conda with version 7.7.0.
I followed the example scripts in https://grand.readthedocs.io/en/latest/examples.scytalone.html but it does not work for my system.
The system is a urokinase with docked ligand. I've followed the tutorials to prepare the input files as following:

10302_inputs.tar.gz

The ligand pdb, prepi and prmtop files were generated by antechamber using GAFF2 from AmberTools 22, and the whole sovated pdb were generated by tleap.

The simulation scripts were like:

os.chdir('/root/Projects/GCMC_Test/gcmc_sampling_1')

ligand_prmtop_file_name = '/root/Projects/GCMC_Test/gcmc_sampling_1/10302.prmtop'
ligand_prepi_file_name = '/root/Projects/GCMC_Test/gcmc_sampling_1/10302.prepi'
ligand_pdb_file_name = '/root/Projects/GCMC_Test/gcmc_sampling_1/10302.pdb'
ligand_xml_file_name = '/root/Projects/GCMC_Test/gcmc_sampling_1/10302.xml'
solvated_pdb_file_name = '/root/Projects/GCMC_Test/gcmc_sampling_1/system_solvated_10302.pdb'
solvated_pdb_conect_file_name = '/root/Projects/GCMC_Test/gcmc_sampling_1/system_solvated_10302_conect.pdb'

gcmc_ghost_pdb_file_name = '/root/Projects/GCMC_Test/gcmc_sampling_1/10302_gcmc_extra_waters.pdb'
gcmc_ghost_txt_file_name = '/root/Projects/GCMC_Test/gcmc_sampling_1/10302_gcmc_extra_waters.txt'
gcmc_log_file_name = '/root/Projects/GCMC_Test/gcmc_sampling_1/10302_gcmc.log'
gcmc_dcd_file_name = '/root/Projects/GCMC_Test/gcmc_sampling_1/10302_gcmc.dcd'
gcmc_rst_file_name = '/root/Projects/GCMC_Test/gcmc_sampling_1/10302_gcmc.rst7'
gcmc_out_file_name = '/root/Projects/GCMC_Test/gcmc_sampling_1/10302_gcmc.out'

# Write CONECT lines to PDB
grand.utils.write_conect(solvated_pdb_file_name, 'MOL', ligand_prepi_file_name, solvated_pdb_conect_file_name)

# Write ligand XML
grand.utils.create_ligand_xml(ligand_prmtop_file_name, ligand_prepi_file_name, 'MOL', ligand_xml_file_name)

system_pdb = PDBFile(solvated_pdb_conect_file_name)
system_pdb.topology, system_pdb.positions, ghosts = grand.utils.add_ghosts(system_pdb.topology, system_pdb.positions, n=10, pdb=gcmc_ghost_pdb_file_name)
system_ff = ForceField('amber14-all.xml', 'amber14/tip3p.xml', ligand_xml_file_name)
system = system_ff.createSystem(system_pdb.topology,
                                nonbondedMethod=PME,
                                nonbondedCutoff=12.0 * angstrom,
                                switchDistance=10.0 * angstrom,
                                constraints=HBonds,
                                rigidWater=True,
                                removeCMMotion=True,
                                hydrogenMass=None,
                                ewaldErrorTolerance=0.0005)

gcmc_reference_atom_list = []
for atom in system_pdb.topology.atoms():
    if atom.residue.name == 'MOL':
        atom_info_dict = {'name': atom.name, 'resname': atom.residue.name, 'resid': atom.residue.id}
        gcmc_reference_atom_list.append(atom_info_dict)

gcmc_mover = grand.samplers.StandardGCMCSphereSampler(system=system,
                                                      topology=system_pdb.topology,
                                                      temperature=298.15 * kelvin,
                                                      referenceAtoms=gcmc_reference_atom_list,
                                                      sphereRadius=4.2 * angstroms,
                                                      ghostFile=gcmc_ghost_txt_file_name,
                                                      log=gcmc_log_file_name,
                                                      dcd=gcmc_dcd_file_name,
                                                      rst=gcmc_rst_file_name,
                                                      overwrite=False)

integrator = LangevinIntegrator(298.15 * kelvin, 1.0 / picosecond, 0.002 * picoseconds)
platform = Platform.getPlatformByName('CUDA')
platform_property = {'Precision': 'mixed', 'DisablePmeStream': 'true', 'DeviceIndex': '0'}
simulation = Simulation(system_pdb.topology, system, integrator, platform, platform_property)
simulation.context.setPositions(system_pdb.positions)

simulation.context.setPeriodicBoxVectors(*system_pdb.topology.getPeriodicBoxVectors())
simulation.context.setVelocitiesToTemperature(298.15)

simulation. Step(250000)

# Prepare the GCMC sphere
gcmc_mover.initialise(simulation.context, ghosts)
gcmc_mover.deleteWatersInGCMCSphere()

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))



simulation. Reporters = []
simulation.context.setTime(0)
simulation.currentStep = 0

state_data_reporter = StateDataReporter(gcmc_out_file_name,
                                        100,
                                        step=True,
                                        time=True,
                                        potentialEnergy=True,
                                        temperature=True,
                                        kineticEnergy=True,
                                        totalEnergy=True,
                                        speed=True,
                                        separator='\t',
                                        progress=True,
                                        remainingTime=True,
                                        totalSteps=25000000)

simulation.reporters.append(state_data_reporter)

num_gcmc_moves = 100
total_n_steps = 25000
gcmc_n_steps = 1000
num_gcmc_round = int(total_n_steps / gcmc_n_steps)

gcmc_mover.reset()
for i in range(num_gcmc_round):
    # Carry out 100 GCMC moves per 2 ps of MD
    simulation. Step(gcmc_n_steps)
    gcmc_mover.move(simulation. Context, num_gcmc_moves)
    gcmc_mover.report(simulation)

The simulation explodes even when it's just equilibrating the system using GCMC. However, if I remove all gcmc steps, the simulation runs smoothly without errors. I've tried to do energy minimization before the gcmc moves, but it still does not work.

I've inspect the error, and found that the system explode during the InsertonMoves. Normally the simulation would run for like 1 or 2 minutes and then crash.

It's very helpful for be to hear your response and see how to debug this. Thanks!!

Paper scripts give error

Hi,
Many thanks for making your great GCMC code available.
I have installed with the following:

bash ./Miniconda-py37_4.9.2-Linux-x86_64.sh -b -p /opt/MD-software/miniconda3-grand
conda install --yes -c omnia -c anaconda -c conda-forge -c essexlab grand
conda install --yes -c omnia/label/cuda92 -c conda-forge openmm

Then downloaded and tried the paper BPTI scripts from bpti gcmc-md.
When I do:

cd gcmc-md/equil/
python equil-uvt1.py

The simulation runs fine and carried out 110k std GCMC moves but it dies before it can write out the final pdb - bpti-uvt1.pdb. Has the latest code changes in relation to this command?

pdb.topology, pdb.positions = grand.utils.remove_ghosts(pdb.topology, positions,
                                                        ghosts=ghost_resids,
                                                        pdb='bpti-uvt1.pdb')

Protein-ligand GCMC simulations

Hi,

I'm trying to use grand to equilibrate a protein-ligand system, but ran into some issues with getting grand to understand the ligand parameters. I've followed the scytalone example, but I didn't have the prepi files, so I used MOL2 instead. I get the XML file but when I try to define the system, I get the following error:

ValueError: No template found for residue 299 (UNK).  The set of atoms matches UNK-0, but the bonds are different.

I'm attaching all the files and the NB with the code.

GCMC water box starting from zero density

Dear developers,

I've been playing with the water example script, testing to see if GCMC can properly solvate a box initialized at close to zero density (a single non-GCMC water was present with state=2). Unfortunately, I was never able maintain any GCMC waters:

Equilibration:
765/10000 equilibration GCMC moves accepted. N = 0
Production:
5000 move(s) completed (408 accepted (8.1600 %)). Current N = 0. Average N = 0.090

Fundamentally, the way that the Adams formulation of GCMC may be flawed in the extreme case. Consider the Adams formulation for $P_{insert}$ and $P_{delete}$

$P_{insert}=\min \left(1, \dfrac{1}{N+1} \exp(B)\exp(-\beta \Delta U) \right)$
$P_{delete}=\min \left(1, N \exp(-B)\exp(-\beta \Delta U) \right)$

A typical adams parameter is around $-3$ (as corroborated by the various papers). Starting from a fully empty box $N=0$, after some effort, we succeed in inserting a single GCMC water so that $N=1$. However, the next $P_{delete}$ move will then always succeed, since it evaluates to $P_{delete}=\min \left(1, 1 \exp(3)\exp(0) \right) = 1$ since $\Delta U=0$ as we're moving only a single water from ghost to real (but now in an empty box).

Questions:

Is this the intended behavior?
Are there critical points by which we may suddenly evaporate or fill the box?
Are there reasonable mitigations to this problem?

ValueError and AttributeError when running test

Hello, I encountered with the following errors while running python setup.py test after installation.

Could you help me with these? It looks like there are some value packing issues.

I installed it with conda and python version 3.8

============================================================================================= test session starts =============================================================================================
platform linux -- Python 3.8.12, pytest-7.0.1, pluggy-1.0.0
rootdir: /root/software/grand
collected 37 items                                                                                                                                                                                            

grand/tests/test_potential.py ...                                                                                                                                                                       [  8%]
grand/tests/test_samplers.py ...FF.F.F.FFFF.F..FFF                                                                                                                                                      [ 64%]
grand/tests/test_utils.py .............                                                                                                                                                                 [100%]

================================================================================================== FAILURES ===================================================================================================
________________________________________________________________________________ TestGCMCSphereSampler.test_deleteRandomWater _________________________________________________________________________________

self = <grand.tests.test_samplers.TestGCMCSphereSampler testMethod=test_deleteRandomWater>

    def test_deleteRandomWater(self):
        """
        Make sure the GCMCSphereSampler.deleteRandomWater() method works correctly
        """
        # Insert a random water
>       gcmc_id, wat_id, atom_ids = gcmc_sphere_sampler.deleteRandomWater()
E       ValueError: not enough values to unpack (expected 3, got 2)

grand/tests/test_samplers.py:509: ValueError
--------------------------------------------------------------------------------------------- Captured log setup ----------------------------------------------------------------------------------------------
INFO     grand.samplers:samplers.py:82 kT = 0.5961612775922495 kcal/mol
INFO     grand.samplers:samplers.py:153 BaseGrandCanonicalMonteCarloSampler object initialised
INFO     grand.samplers:samplers.py:586 GCMC sphere is based on reference atom IDs: [146, 676]
INFO     grand.samplers:samplers.py:596 GCMC sphere radius is 4 A
INFO     grand.samplers:samplers.py:607 Simulating at an Adams (B) value of -8.036693247651913
INFO     grand.samplers:samplers.py:609 GCMCSphereSampler object initialised
_____________________________________________________________________________ TestGCMCSphereSampler.test_deleteWatersInGCMCSphere _____________________________________________________________________________

self = <grand.tests.test_samplers.TestGCMCSphereSampler testMethod=test_deleteWatersInGCMCSphere>

    def test_deleteWatersInGCMCSphere(self):
        """
        Make sure the GCMCSphereSampler.deleteWatersInGCMCSphere() method works correctly
        """
        # Now delete the waters in the sphere
        gcmc_sphere_sampler.deleteWatersInGCMCSphere()
>       new_ghosts = [gcmc_sphere_sampler.gcmc_resids[id] for id in np.where(gcmc_sphere_sampler.gcmc_status == 0)[0]]
E       AttributeError: 'GCMCSphereSampler' object has no attribute 'gcmc_status'

grand/tests/test_samplers.py:444: AttributeError
________________________________________________________________________________ TestGCMCSphereSampler.test_insertRandomWater _________________________________________________________________________________

self = <grand.tests.test_samplers.TestGCMCSphereSampler testMethod=test_insertRandomWater>

    def test_insertRandomWater(self):
        """
        Make sure the GCMCSphereSampler.insertRandomWater() method works correctly
        """
        # Insert a random water
>       new_positions, gcmc_id, wat_id, atom_ids = gcmc_sphere_sampler.insertRandomWater()
E       ValueError: not enough values to unpack (expected 4, got 3)

grand/tests/test_samplers.py:489: ValueError
_________________________________________________________________________________ TestGCMCSphereSampler.test_updateGCMCSphere _________________________________________________________________________________

self = <grand.tests.test_samplers.TestGCMCSphereSampler testMethod=test_updateGCMCSphere>

    def test_updateGCMCSphere(self):
        """
        Make sure the GCMCSphereSampler.updateGCMCSphere() method works correctly
        """
        # Get initial gcmc_resids and status
>       gcmc_resids = deepcopy(gcmc_sphere_sampler.gcmc_resids)
E       AttributeError: 'GCMCSphereSampler' object has no attribute 'gcmc_resids'

grand/tests/test_samplers.py:458: AttributeError
____________________________________________________________________________ TestNonequilibriumGCMCSphereSampler.test_deletionMove ____________________________________________________________________________

self = <grand.tests.test_samplers.TestNonequilibriumGCMCSphereSampler testMethod=test_deletionMove>

    def test_deletionMove(self):
        """
        Make sure the NonequilibriumGCMCSphereSampler.deletionMove() method works correctly
        """
        # Prep for a move
        # Read in positions
        neq_gcmc_sphere_sampler.context = neq_gcmc_sphere_simulation.context
        state = neq_gcmc_sphere_sampler.context.getState(getPositions=True, enforcePeriodicBox=True, getVelocities=True)
        neq_gcmc_sphere_sampler.positions = deepcopy(state.getPositions(asNumpy=True))
        neq_gcmc_sphere_sampler.velocities = deepcopy(state.getVelocities(asNumpy=True))
    
        # Update GCMC region based on current state
        neq_gcmc_sphere_sampler.updateGCMCSphere(state)
    
        # Set to NCMC integrator
        neq_gcmc_sphere_sampler.compound_integrator.setCurrentIntegrator(1)
    
        # Just run one move to make sure it doesn't crash
>       neq_gcmc_sphere_sampler.deletionMove()

grand/tests/test_samplers.py:656: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
grand/samplers.py:1415: in deletionMove
    self.ncmc_integrator.step(self.n_prop_steps_per_pert)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

self = <openmmtools.integrators.NonequilibriumLangevinIntegrator; proxy of <Swig Object of type 'OpenMM::CustomIntegrator *' at 0x7f9d72456e40> >, steps = 10

    def step(self, steps):
        r"""
        step(self, steps)
        Advance a simulation through time by taking a series of time steps.
    
        Parameters
        ----------
        steps : int
            the number of time steps to take
        """
>       return _openmm.CustomIntegrator_step(self, steps)
E       openmm.OpenMMException: Particle coordinate is nan

/root/anaconda3/envs/py38-grand/lib/python3.8/site-packages/openmm/openmm.py:6011: OpenMMException
--------------------------------------------------------------------------------------------- Captured log setup ----------------------------------------------------------------------------------------------
INFO     grand.samplers:samplers.py:82 kT = 0.5961612775922495 kcal/mol
INFO     grand.samplers:samplers.py:153 BaseGrandCanonicalMonteCarloSampler object initialised
INFO     grand.samplers:samplers.py:586 GCMC sphere is based on reference atom IDs: [146, 676]
INFO     grand.samplers:samplers.py:596 GCMC sphere radius is 4 A
INFO     grand.samplers:samplers.py:607 Simulating at an Adams (B) value of -8.036693247651913
INFO     grand.samplers:samplers.py:609 GCMCSphereSampler object initialised
INFO     grand.samplers:samplers.py:1255 Each NCMC move will be executed over a total of 0.04 ps
INFO     grand.samplers:samplers.py:1274 NonequilibriumGCMCSphereSampler object initialised
___________________________________________________________________________ TestNonequilibriumGCMCSphereSampler.test_insertionMove ____________________________________________________________________________

self = <grand.tests.test_samplers.TestNonequilibriumGCMCSphereSampler testMethod=test_insertionMove>

    def test_insertionMove(self):
        """
        Make sure the NonequilibriumGCMCSphereSampler.insertionMove() method works correctly
        """
        # Prep for a move
        # Read in positions
        neq_gcmc_sphere_sampler.context = neq_gcmc_sphere_simulation.context
        state = neq_gcmc_sphere_sampler.context.getState(getPositions=True, enforcePeriodicBox=True, getVelocities=True)
        neq_gcmc_sphere_sampler.positions = deepcopy(state.getPositions(asNumpy=True))
        neq_gcmc_sphere_sampler.velocities = deepcopy(state.getVelocities(asNumpy=True))
    
        # Update GCMC region based on current state
        neq_gcmc_sphere_sampler.updateGCMCSphere(state)
    
        # Set to NCMC integrator
        neq_gcmc_sphere_sampler.compound_integrator.setCurrentIntegrator(1)
    
        # Just run one move to make sure it doesn't crash
>       neq_gcmc_sphere_sampler.insertionMove()

grand/tests/test_samplers.py:631: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
grand/samplers.py:1332: in insertionMove
    self.ncmc_integrator.step(self.n_prop_steps_per_pert)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

self = <openmmtools.integrators.NonequilibriumLangevinIntegrator; proxy of <Swig Object of type 'OpenMM::CustomIntegrator *' at 0x7f9d72456e40> >, steps = 10

    def step(self, steps):
        r"""
        step(self, steps)
        Advance a simulation through time by taking a series of time steps.
    
        Parameters
        ----------
        steps : int
            the number of time steps to take
        """
>       return _openmm.CustomIntegrator_step(self, steps)
E       openmm.OpenMMException: Particle coordinate is nan

/root/anaconda3/envs/py38-grand/lib/python3.8/site-packages/openmm/openmm.py:6011: OpenMMException
________________________________________________________________________________ TestNonequilibriumGCMCSphereSampler.test_move ________________________________________________________________________________

self = <grand.tests.test_samplers.TestNonequilibriumGCMCSphereSampler testMethod=test_move>

    def test_move(self):
        """
        Make sure the NonequilibriumGCMCSphereSampler.move() method works correctly
        """
        neq_gcmc_sphere_sampler.reset()
    
        # Just run one move, as they are a bit more expensive
>       neq_gcmc_sphere_sampler.move(neq_gcmc_sphere_simulation.context, 1)

grand/tests/test_samplers.py:596: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
grand/samplers.py:1307: in move
    self.deletionMove()
grand/samplers.py:1415: in deletionMove
    self.ncmc_integrator.step(self.n_prop_steps_per_pert)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

self = <openmmtools.integrators.NonequilibriumLangevinIntegrator; proxy of <Swig Object of type 'OpenMM::CustomIntegrator *' at 0x7f9d72456e40> >, steps = 10

    def step(self, steps):
        r"""
        step(self, steps)
        Advance a simulation through time by taking a series of time steps.
    
        Parameters
        ----------
        steps : int
            the number of time steps to take
        """
>       return _openmm.CustomIntegrator_step(self, steps)
E       openmm.OpenMMException: Particle coordinate is nan

/root/anaconda3/envs/py38-grand/lib/python3.8/site-packages/openmm/openmm.py:6011: OpenMMException
---------------------------------------------------------------------------------------------- Captured log call ----------------------------------------------------------------------------------------------
INFO     grand.samplers:samplers.py:1487 Resetting any tracked variables...
________________________________________________________________________________ TestGCMCSystemSampler.test_deleteRandomWater _________________________________________________________________________________

self = <grand.tests.test_samplers.TestGCMCSystemSampler testMethod=test_deleteRandomWater>

    def test_deleteRandomWater(self):
        """
        Make sure the GCMCSystemSampler.deleteRandomWater() method works correctly
        """
        # Insert a random water
>       gcmc_id, wat_id, atom_ids = gcmc_system_sampler.deleteRandomWater()
E       ValueError: not enough values to unpack (expected 3, got 2)

grand/tests/test_samplers.py:734: ValueError
--------------------------------------------------------------------------------------------- Captured log setup ----------------------------------------------------------------------------------------------
INFO     grand.samplers:samplers.py:82 kT = 0.5961612775922495 kcal/mol
INFO     grand.samplers:samplers.py:153 BaseGrandCanonicalMonteCarloSampler object initialised
INFO     grand.samplers:samplers.py:1571 Simulating at an Adams (B) value of -2.561049941969956
INFO     grand.samplers:samplers.py:1573 GCMCSystemSampler object initialised
________________________________________________________________________________ TestGCMCSystemSampler.test_insertRandomWater _________________________________________________________________________________

self = <grand.tests.test_samplers.TestGCMCSystemSampler testMethod=test_insertRandomWater>

    def test_insertRandomWater(self):
        """
        Make sure the GCMCSystemSampler.insertRandomWater() method works correctly
        """
        # Insert a random water
>       new_positions, gcmc_id, wat_id, atom_ids = gcmc_system_sampler.insertRandomWater()
E       ValueError: not enough values to unpack (expected 4, got 3)

grand/tests/test_samplers.py:714: ValueError
____________________________________________________________________________ TestNonequilibriumGCMCSystemSampler.test_deletionMove ____________________________________________________________________________

self = <grand.tests.test_samplers.TestNonequilibriumGCMCSystemSampler testMethod=test_deletionMove>

    def test_deletionMove(self):
        """
        Make sure the NonequilibriumGCMCSystemSampler.deletionMove() method works correctly
        """
        # Prep for a move
        # Read in positions
        neq_gcmc_system_sampler.context = neq_gcmc_system_simulation.context
        state = neq_gcmc_system_sampler.context.getState(getPositions=True, enforcePeriodicBox=True, getVelocities=True)
        neq_gcmc_system_sampler.positions = deepcopy(state.getPositions(asNumpy=True))
        neq_gcmc_system_sampler.velocities = deepcopy(state.getVelocities(asNumpy=True))
    
        # Set to NCMC integrator
        neq_gcmc_system_sampler.compound_integrator.setCurrentIntegrator(1)
    
        # Just run one move to make sure it doesn't crash
>       neq_gcmc_system_sampler.deletionMove()

grand/tests/test_samplers.py:874: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
grand/samplers.py:2053: in deletionMove
    self.ncmc_integrator.step(self.n_prop_steps_per_pert)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

self = <openmmtools.integrators.NonequilibriumLangevinIntegrator; proxy of <Swig Object of type 'OpenMM::CustomIntegrator *' at 0x7f9d643f1f60> >, steps = 1

    def step(self, steps):
        r"""
        step(self, steps)
        Advance a simulation through time by taking a series of time steps.
    
        Parameters
        ----------
        steps : int
            the number of time steps to take
        """
>       return _openmm.CustomIntegrator_step(self, steps)
E       openmm.OpenMMException: Particle coordinate is nan

/root/anaconda3/envs/py38-grand/lib/python3.8/site-packages/openmm/openmm.py:6011: OpenMMException
--------------------------------------------------------------------------------------------- Captured log setup ----------------------------------------------------------------------------------------------
INFO     grand.samplers:samplers.py:82 kT = 0.5961612775922495 kcal/mol
INFO     grand.samplers:samplers.py:153 BaseGrandCanonicalMonteCarloSampler object initialised
INFO     grand.samplers:samplers.py:1571 Simulating at an Adams (B) value of -2.561049941969956
INFO     grand.samplers:samplers.py:1573 GCMCSystemSampler object initialised
INFO     grand.samplers:samplers.py:1919 Each NCMC move will be executed over a total of 0.004 ps
INFO     grand.samplers:samplers.py:1939 NonequilibriumGCMCSystemSampler object initialised
___________________________________________________________________________ TestNonequilibriumGCMCSystemSampler.test_insertionMove ____________________________________________________________________________

self = <grand.tests.test_samplers.TestNonequilibriumGCMCSystemSampler testMethod=test_insertionMove>

    def test_insertionMove(self):
        """
        Make sure the NonequilibriumGCMCSystemSampler.insertionMove() method works correctly
        """
        # Prep for a move
        # Read in positions
        neq_gcmc_system_sampler.context = neq_gcmc_system_simulation.context
        state = neq_gcmc_system_sampler.context.getState(getPositions=True, enforcePeriodicBox=True, getVelocities=True)
        neq_gcmc_system_sampler.positions = deepcopy(state.getPositions(asNumpy=True))
        neq_gcmc_system_sampler.velocities = deepcopy(state.getVelocities(asNumpy=True))
    
        # Set to NCMC integrator
        neq_gcmc_system_sampler.compound_integrator.setCurrentIntegrator(1)
    
        # Just run one move to make sure it doesn't crash
>       neq_gcmc_system_sampler.insertionMove()

grand/tests/test_samplers.py:852: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
grand/samplers.py:1991: in insertionMove
    self.ncmc_integrator.step(self.n_prop_steps_per_pert)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

self = <openmmtools.integrators.NonequilibriumLangevinIntegrator; proxy of <Swig Object of type 'OpenMM::CustomIntegrator *' at 0x7f9d643f1f60> >, steps = 1

    def step(self, steps):
        r"""
        step(self, steps)
        Advance a simulation through time by taking a series of time steps.
    
        Parameters
        ----------
        steps : int
            the number of time steps to take
        """
>       return _openmm.CustomIntegrator_step(self, steps)
E       openmm.OpenMMException: Particle coordinate is nan

/root/anaconda3/envs/py38-grand/lib/python3.8/site-packages/openmm/openmm.py:6011: OpenMMException
________________________________________________________________________________ TestNonequilibriumGCMCSystemSampler.test_move ________________________________________________________________________________

self = <grand.tests.test_samplers.TestNonequilibriumGCMCSystemSampler testMethod=test_move>

    def test_move(self):
        """
        Make sure the NonequilibriumGCMCSystemSampler.move() method works correctly
        """
        neq_gcmc_system_sampler.reset()
    
        # Just run one move, as they are a bit more expensive
>       neq_gcmc_system_sampler.move(neq_gcmc_system_simulation.context, 1)

grand/tests/test_samplers.py:821: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
grand/samplers.py:1969: in move
    self.deletionMove()
grand/samplers.py:2053: in deletionMove
    self.ncmc_integrator.step(self.n_prop_steps_per_pert)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

self = <openmmtools.integrators.NonequilibriumLangevinIntegrator; proxy of <Swig Object of type 'OpenMM::CustomIntegrator *' at 0x7f9d643f1f60> >, steps = 1

    def step(self, steps):
        r"""
        step(self, steps)
        Advance a simulation through time by taking a series of time steps.
    
        Parameters
        ----------
        steps : int
            the number of time steps to take
        """
>       return _openmm.CustomIntegrator_step(self, steps)
E       openmm.OpenMMException: Particle coordinate is nan

/root/anaconda3/envs/py38-grand/lib/python3.8/site-packages/openmm/openmm.py:6011: OpenMMException
---------------------------------------------------------------------------------------------- Captured log call ----------------------------------------------------------------------------------------------
INFO     grand.samplers:samplers.py:2106 Resetting any tracked variables...
============================================================================================== warnings summary ===============================================================================================
grand/tests/test_potential.py::TestPotential::test_calc_mu
  /root/anaconda3/envs/py38-grand/lib/python3.8/site-packages/pymbar/mbar.py:526: PendingDeprecationWarning: the matrix subclass is not the recommended way to represent matrices or deal with linear algebra (see https://docs.scipy.org/doc/numpy/user/numpy-for-matlab-users.html). Please adjust your code to use regular ndarray.
    f_i = np.matrix(self.f_k)

grand/tests/test_potential.py::TestPotential::test_calc_mu
  /root/anaconda3/envs/py38-grand/lib/python3.8/site-packages/pymbar/mbar.py:1717: PendingDeprecationWarning: the matrix subclass is not the recommended way to represent matrices or deal with linear algebra (see https://docs.scipy.org/doc/numpy/user/numpy-for-matlab-users.html). Please adjust your code to use regular ndarray.
    Ndiag = np.matrix(np.diag(N_k), dtype=np.float64)

grand/tests/test_potential.py::TestPotential::test_calc_mu
  /root/anaconda3/envs/py38-grand/lib/python3.8/site-packages/pymbar/mbar.py:1718: PendingDeprecationWarning: the matrix subclass is not the recommended way to represent matrices or deal with linear algebra (see https://docs.scipy.org/doc/numpy/user/numpy-for-matlab-users.html). Please adjust your code to use regular ndarray.
    W = np.matrix(W, dtype=np.float64)

grand/tests/test_potential.py::TestPotential::test_calc_mu
grand/tests/test_potential.py::TestPotential::test_calc_mu
grand/tests/test_potential.py::TestPotential::test_calc_mu
grand/tests/test_potential.py::TestPotential::test_calc_mu
grand/tests/test_potential.py::TestPotential::test_calc_mu
grand/tests/test_potential.py::TestPotential::test_calc_mu
grand/tests/test_potential.py::TestPotential::test_calc_mu
grand/tests/test_potential.py::TestPotential::test_calc_mu
grand/tests/test_potential.py::TestPotential::test_calc_mu
  /root/anaconda3/envs/py38-grand/lib/python3.8/site-packages/numpy/matrixlib/defmatrix.py:69: PendingDeprecationWarning: the matrix subclass is not the recommended way to represent matrices or deal with linear algebra (see https://docs.scipy.org/doc/numpy/user/numpy-for-matlab-users.html). Please adjust your code to use regular ndarray.
    return matrix(data, dtype=dtype, copy=False)

grand/tests/test_potential.py::TestPotential::test_calc_mu
  /root/anaconda3/envs/py38-grand/lib/python3.8/site-packages/pymbar/mbar.py:1728: PendingDeprecationWarning: the matrix subclass is not the recommended way to represent matrices or deal with linear algebra (see https://docs.scipy.org/doc/numpy/user/numpy-for-matlab-users.html). Please adjust your code to use regular ndarray.
    Sigma = np.matrix(np.diag(np.sqrt(S2)))

grand/tests/test_potential.py::TestPotential::test_calc_mu
  /root/anaconda3/envs/py38-grand/lib/python3.8/site-packages/pymbar/mbar.py:1729: PendingDeprecationWarning: the matrix subclass is not the recommended way to represent matrices or deal with linear algebra (see https://docs.scipy.org/doc/numpy/user/numpy-for-matlab-users.html). Please adjust your code to use regular ndarray.
    V = np.matrix(V)

grand/tests/test_potential.py::TestPotential::test_calc_mu
  /root/anaconda3/envs/py38-grand/lib/python3.8/site-packages/pymbar/mbar.py:1581: PendingDeprecationWarning: the matrix subclass is not the recommended way to represent matrices or deal with linear algebra (see https://docs.scipy.org/doc/numpy/user/numpy-for-matlab-users.html). Please adjust your code to use regular ndarray.
    diag = np.matrix(cov.diagonal())

grand/tests/test_samplers.py::TestBaseGrandCanonicalMonteCarloSampler::test_report
  /root/anaconda3/envs/py38-grand/lib/python3.8/site-packages/numpy/core/fromnumeric.py:3440: RuntimeWarning: Mean of empty slice.
    return _methods._mean(a, axis=axis, dtype=dtype,

grand/tests/test_samplers.py::TestBaseGrandCanonicalMonteCarloSampler::test_report
  /root/anaconda3/envs/py38-grand/lib/python3.8/site-packages/numpy/core/_methods.py:189: RuntimeWarning: invalid value encountered in double_scalars
    ret = ret.dtype.type(ret / rcount)

grand/tests/test_utils.py::TestUtils::test_cluster_waters
  /root/software/grand/grand/utils.py:901: UserWarning: Chains are not specified for at least one reference atom, will select the first instance of Tyr10 found
    warnings.warn("Chains are not specified for at least one reference atom, will select the first instance"

grand/tests/test_utils.py::TestUtils::test_cluster_waters
  /root/software/grand/grand/utils.py:901: UserWarning: Chains are not specified for at least one reference atom, will select the first instance of Asn43 found
    warnings.warn("Chains are not specified for at least one reference atom, will select the first instance"

grand/tests/test_utils.py::TestUtils::test_write_sphere_traj
  /root/software/grand/grand/utils.py:797: UserWarning: Chains are not specified for at least one reference atom, will select the first instance of Tyr9 found
    warnings.warn("Chains are not specified for at least one reference atom, will select the first instance"

grand/tests/test_utils.py::TestUtils::test_write_sphere_traj
  /root/software/grand/grand/utils.py:797: UserWarning: Chains are not specified for at least one reference atom, will select the first instance of Tyr10 found
    warnings.warn("Chains are not specified for at least one reference atom, will select the first instance"

grand/tests/test_utils.py::TestUtils::test_write_sphere_traj
  /root/software/grand/grand/utils.py:797: UserWarning: Chains are not specified for at least one reference atom, will select the first instance of Asn43 found
    warnings.warn("Chains are not specified for at least one reference atom, will select the first instance"

-- Docs: https://docs.pytest.org/en/stable/how-to/capture-warnings.html
=========================================================================================== short test summary info ===========================================================================================
FAILED grand/tests/test_samplers.py::TestGCMCSphereSampler::test_deleteRandomWater - ValueError: not enough values to unpack (expected 3, got 2)
FAILED grand/tests/test_samplers.py::TestGCMCSphereSampler::test_deleteWatersInGCMCSphere - AttributeError: 'GCMCSphereSampler' object has no attribute 'gcmc_status'
FAILED grand/tests/test_samplers.py::TestGCMCSphereSampler::test_insertRandomWater - ValueError: not enough values to unpack (expected 4, got 3)
FAILED grand/tests/test_samplers.py::TestGCMCSphereSampler::test_updateGCMCSphere - AttributeError: 'GCMCSphereSampler' object has no attribute 'gcmc_resids'
FAILED grand/tests/test_samplers.py::TestNonequilibriumGCMCSphereSampler::test_deletionMove - openmm.OpenMMException: Particle coordinate is nan
FAILED grand/tests/test_samplers.py::TestNonequilibriumGCMCSphereSampler::test_insertionMove - openmm.OpenMMException: Particle coordinate is nan
FAILED grand/tests/test_samplers.py::TestNonequilibriumGCMCSphereSampler::test_move - openmm.OpenMMException: Particle coordinate is nan
FAILED grand/tests/test_samplers.py::TestGCMCSystemSampler::test_deleteRandomWater - ValueError: not enough values to unpack (expected 3, got 2)
FAILED grand/tests/test_samplers.py::TestGCMCSystemSampler::test_insertRandomWater - ValueError: not enough values to unpack (expected 4, got 3)
FAILED grand/tests/test_samplers.py::TestNonequilibriumGCMCSystemSampler::test_deletionMove - openmm.OpenMMException: Particle coordinate is nan
FAILED grand/tests/test_samplers.py::TestNonequilibriumGCMCSystemSampler::test_insertionMove - openmm.OpenMMException: Particle coordinate is nan
FAILED grand/tests/test_samplers.py::TestNonequilibriumGCMCSystemSampler::test_move - openmm.OpenMMException: Particle coordinate is nan
============================================================================ 12 failed, 25 passed, 22 warnings in 98.14s (0:01:38) ============================================================================

Cannot specify chain information in analysis functions

Hi,

I have a simulated system that has two chains. When I define ref_atoms for cluster_waters function, I cannot define two atoms from both chains. Attached is the pdb file of the system and you can see residues in two chains have the same resid. So without specifying chain information, I don't think I can make the scripts understand my reference atom definition. But in the current version, I cannot specify chainID in these analysis functions.

I notice there is a function called getReferenceAtomIndices to deal with chain/resid/atom_id information. But both write_sphere_traj and cluster_waters do not use this function for reference atom definition. Maybe a quick fix of this could be make both functions use getReferenceAtomIndices so that it is possible to specify chainID as well?

1ec0-ghosts.pdb.zip

Do probabilities of insertion and deletion have to be equal?

Hi,

I noticed in the grand code and some of the papers that the attempt of insertion or deletion was determined with uniform probability, i.e. the probability of attempting insertion or deletion was both ~0.5 for a GC(NC)MC move. My question is that can I change the probability of insertion/deletion for a GC(NC)MC move during simulation, for example favouring water insertion with higher probability (or the other way round), will the difference in probabilities of insertion/deletion be wrong theoretically?

Thanks very much,

Running NVT MD on GCMC equilibrated structure

Hello,

In the paper for grand (J. Chem. Inf. Model. 2020, 60, 10, 4436–4441), comparison were made between GCMC/MD and NVT MD on BPTI, where all steps of GCMC/MD workflow were replaced by NVT MD. I'm wondering would the result be any different if I GCMC/MD for equilibration first, and run NVT MD for the production run only, in order to reduce the computational cost?

A related question would be: if all the GCMC moves in the GCMC/MD were rejected (especially during later stages of GCMC/MD when the ), would such GCMC/MD be of any different compared to an NVT MD?

Thanks very much,

Adams volume calibration

I've been having difficulty re-producing the density data in (https://pubs.acs.org/doi/10.1021/acs.jcim.0c00648) using the water scripts from this repo and I wanted to clarify the settings on the # of waters in states ${0,1,2}$, where $0$ - GCMC Ghost water, $1$ - GCMC Real Water, $2$ - Non-GCMC Real Water. Notably, when converged, I was seeing densities that are consistently higher than that was reported.

When reading up on the Adams parameter definition in the paper(s), I noticed it was defined as:

$B=\beta \mu'_{sol} + \ln \left( \dfrac{V_G}{V^0} \right)$ where $V_G$ is the GCMC volume. Intuitively, the ratio $\dfrac{V_G}{V^0}$ can be interpreted as the average number of GCMC Real waters, $\bar{N}$, also from the Adams paper. This is important, because the acceptance probabilities themselves for insertion and deletion define the $N$ and $\dfrac{1}{N+1}$ prefactors in terms of the number of GCMC Real waters (state=1), also shown in the code here, as opposed to the total number of waters.

However, in the example water scripts, there's usually a large number of non-GCMC waters. If non-GCMC waters are used, shouldn't the calibration of the Adams parameter then be offset by the number of non-GCMC-waters? ie:

$B=\beta \mu'_{sol} + \ln \left( \dfrac{V_G}{V^0} - N_W \right)$ where $N_W$ is the number of non-GCMC waters.

This Integrator is not bound to a context!

Hello,

I'm trying to run GCNCMC but I get the following error:

Traceback (most recent call last):
  File "/biggin/b196/scro4068/GCNCMC/run_gcncmc.py", line 82, in <module>
    gcncmc_mover.move(simulation.context, 1)
  File "/biggin/b196/scro4068/miniconda3/lib/python3.10/site-packages/grand/samplers.py", line 1307, in move
    self.insertionMove()
  File "/biggin/b196/scro4068/miniconda3/lib/python3.10/site-packages/grand/samplers.py", line 1335, in insertionMove
    self.ncmc_integrator.step(self.n_prop_steps_per_pert)
  File "/biggin/b196/scro4068/miniconda3/lib/python3.10/site-packages/openmm/openmm.py", line 18830, in step
    return _openmm.CustomIntegrator_step(self, steps)
openmm.OpenMMException: This Integrator is not bound to a context!

My script is based on one of the examples from the paper:

from simtk.openmm import *
from simtk.openmm.app import *
from simtk.unit import *
from openmmtools.integrators import BAOABIntegrator
from sys import stdout
import numpy as np
import mdtraj
import argparse
import grand
from openmmforcefields.generators import GAFFTemplateGenerator
from openff.toolkit.topology import Molecule
from openmmforcefields.generators import SystemGenerator
from parmed.gromacs import GromacsTopologyFile
import parmed

parser = argparse.ArgumentParser()
parser.add_argument('--pert', type=int, default=99, help='Number of perturbation steps per NCMC move')
parser.add_argument('--prop', type=int, default=50, help='Number of propagation steps between each perturbation')
args = parser.parse_args()

pdb = PDBFile('minim.pdb')
                                

# Add ghost waters
pdb.topology, pdb.positions, ghosts = grand.utils.add_ghosts(pdb.topology, pdb.positions,
                                                             n=15, pdb='prod-ghosts.pdb')

# Create system
forcefield = ForceField('amber14-all.xml','amber14/protein.ff14SB.xml', 'amber14/tip3p.xml', 'NS.xml')



system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME, 
                                 nonbondedCutoff=12*angstrom,
                                 switchDistance=10*angstrom, 
                                 constraints=HBonds)

# Define reference atoms for the GCMC sphere
ref_atoms = [{'name': 'CA', 'resname': 'TRP', 'resid': '156', 'chain': 0},
             {'name': 'CA', 'resname': 'GLU', 'resid': '198', 'chain': 0}]



# Define integrator
integrator = BAOABIntegrator(310*kelvin, 1.0/picosecond, 0.002*picosecond)


# Define platform and precision
platform = Platform.getPlatformByName('CUDA')
platform.setPropertyDefaultValue('Precision', 'mixed')


# Create GCMC sampler object
gcncmc_mover = grand.samplers.NonequilibriumGCMCSphereSampler(system=system, 
                                                              topology=pdb.topology, 
                                                              temperature=310*kelvin, 
                                                              integrator=integrator,
                                                              nPertSteps=args.pert,
                                                              nPropStepsPerPert=args.prop,
                                                              referenceAtoms=ref_atoms, 
                                                              sphereRadius=6*angstroms, 
                                                              ghostFile="prod-ghosts.txt", 
                                                              log="prod.log", 
                                                              rst="prod.rst7", 
                                                              dcd="prod.dcd", 
                                                              overwrite=True)



# Create simulation object and set positions, velocities and box vectors
simulation = Simulation(pdb.topology, system, integrator, platform)
simulation.context.setPositions(pdb.positions)
simulation.context.setVelocitiesToTemperature(310*kelvin)
simulation.context.setPeriodicBoxVectors(*pdb.topology.getPeriodicBoxVectors())

# Initialise the GCMC mover
gcncmc_mover.initialise(simulation.context, ghosts)

# Run GCNCMC/MD production - 4000 iterations of 1 NCMC move and 5 ps MD
for i in range(4000):
    simulation.step(2500)
    gcncmc_mover.move(simulation.context, 1)
    gcncmc_mover.report(simulation)

I have no idea what may be causing this, it is worth noting that my system was prepared in GROMACS and I converted the ligand topology using grand.utils create_ligand_xml and write_conect (as my ligand topology was originally generated using antechamber).

'No ghost water molecules left' error

I'm using the equilibration protocol from the BPTI example, with two UVT and NPT simulations, to equilibrate one of my systems. I'm getting a strange grand error that stochstically pops up during the first UVT sim.

2021-11-24 21:32:37,131 - ERROR: No ghost water molecules left, so insertion moves cannot occur - add more ghost waters

I'm using quite a large sphere (10 angstrom) and lots of ghosts (n = 15) to be on the safe side of getting all the waters right in the binding site. Input and log files are here.

Recommend Projects

React

A declarative, efficient, and flexible JavaScript library for building user interfaces.
Vue.js

🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
Typescript

TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
TensorFlow

An Open Source Machine Learning Framework for Everyone
Django

The Web framework for perfectionists with deadlines.
Laravel

A PHP framework for web artisans
D3

Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

javascript

JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
web

Some thing interesting about web. New door for the world.
server

A server is a program made to process requests and deliver data to clients.
Machine learning

Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Visualization

Some thing interesting about visualization, use data art
Game

Some thing interesting about game, make everyone happy.

Recommend Org

Facebook

We are working to build community through open source technology. NB: members must have two-factor auth.
Microsoft

Open source projects and samples from Microsoft.
Google

Google ❤️ Open Source for everyone.
Alibaba

Alibaba Open Source for everyone
D3

Data-Driven Documents codes.
Tencent

China tencent open source team.