Running a Molecular Dynamics (MD) simulation of a protein-ligand complex#

⚠️ Make sure you have GROMACS in your path:

Before running this notebook, make sure that you have Gromacs installed and that you have sourced GMXRC to get access to GROMACS,

e.g. source /usr/local/gromacs/bin/GMXRC.

In this notebook we run an MD simulation of benzene bound to T4-lysozyme L99A using Gromacs.image

The Gromacs MD protocol allows the user to run an MD simulation in solvent of e.g. a small molecule, a protein, or a protein-ligand complex. Running an MD simulations can be useful for a variety of things, such as pre-equilibration of a protein-ligand complex prior to running free energy calculations or for getting insights into the general dynamics of the system of interest.

The MD protocol uses different software tools in a series of steps and provides multiple outputs for the user:

Step

Software used

  1. Input handling using gufe

OpenFE, Gufe, RDKit

  1. Parameterization using OpenMMForceFields & OpenFF

OpenFE - OpenMMForceFields - OpenFF

  1. Creating of Gromacs input files (.gro, .top) via Interchange

OpenFE - OpenFF

  1. Minimization

GROMACS

  1. NVT equilibration

GROMACS

  1. NPT MD

GROMACS

1. Defining the ChemicalSystem#

ChemicalSystems are containers which define the various components which exist in a system of interest. Here, we will be passing the SmallMoleculeComponent for benzene, a ProteinComponent generated from a PDB file, and a SolventComponent which will contain the necessary information for OpenMM’s Modeller class to add water and 0.15 M NaCl around the solute.

[1]:
import openfe_gromacs
from openfe_gromacs import ChemicalSystem, ProteinComponent, SmallMoleculeComponent, SolventComponent
from openff.units import unit

# Define the ligand we are interested in
ligand = SmallMoleculeComponent.from_sdf_file('assets/benzene.sdf')

# Define the solvent environment and protein structure
solvent = SolventComponent(ion_concentration=0.15 * unit.molar)
protein = ProteinComponent.from_pdb_file('assets/t4_lysozyme.pdb', name='t4-lysozyme')

# create the ChemicalSystem
system = ChemicalSystem({'ligand': ligand, 'protein': protein, 'solvent': solvent}, name=f"{ligand.name}_{protein.name}")

2. Defining the MD simulation settings#

There are various different parameters which can be set to determine how the MD simulation will take place. To allow for maximum user flexibility, these are defined as a series of settings objects which control the following:

Setting

Description

protocol_repeats

How many repeats of the MD protocol should be run.

simulation_settings_em

Energy minimization simulation settings, e.g. nsteps.

simulation_settings_nvt

NVT MD simulation settings, e.g. nsteps.

simulation_settings_npt

NPT MD simulation settings, e.g. nsteps.

output_settings_em

Energy minimization outputs settings, e.g. file names and write frequencies.

output_settings_nvt

NVT MD output settings, e.g. file names and write frequencies.

output_settings_npt

NPT MD output settings, e.g. file names and write frequencies.

forcefield_settings

Settings controlling how the system is parameterized.

partial_charge_settings

Settings controlling how small molecule partial charges are assigned.

solvation_settings

Settings controlling how the system should be solvated.

thermo_settings

Settings controlling thermodynamic parameters e.g. the temperature and the pressure of the system.

top

Name of the Gromacs topology .top file that should be written (using Interchange).

gro

Name of the Gromacs cordinate .gro file that should be written (using Interchange).

The easiest way to access and change settings is by first importing the default settings, printing them and then changing the settings according to the user’s needs.

[2]:
from openfe_gromacs.protocols.gromacs_md.md_methods import GromacsMDProtocol
from openff.units import unit

settings = GromacsMDProtocol.default_settings()
settings.simulation_settings_em.nsteps = 10 # setting the number of minimization steps to 10
settings.simulation_settings_nvt.nsteps = 10 # setting the number of NVT MD steps to 10
settings.simulation_settings_npt.nsteps = 10 # setting the number of NVT MD steps to 10
settings.solvation_settings.box_shape = 'dodecahedron'
[3]:
settings
{'engine_settings': {'bonded': 'auto',
                     'nb': 'auto',
                     'ntomp': 1,
                     'pme': 'auto',
                     'pmefft': 'auto',
                     'update': 'auto'},
 'forcefield_settings': {'constraints': None,
                         'forcefields': ['amber/ff14SB.xml',
                                         'amber/tip3p_standard.xml',
                                         'amber/tip3p_HFE_multivalent.xml',
                                         'amber/phosaa10.xml'],
                         'hydrogen_mass': 3.0,
                         'nonbonded_cutoff': <Quantity(1.0, 'nanometer')>,
                         'nonbonded_method': 'PME',
                         'rigid_water': True,
                         'small_molecule_forcefield': 'openff-2.1.1'},
 'gro': 'system.gro',
 'integrator_settings': {'barostat_frequency': <Quantity(25, 'timestep')>,
                         'constraint_tolerance': 1e-06,
                         'langevin_collision_rate': <Quantity(1.0, '1 / picosecond')>,
                         'n_restart_attempts': 20,
                         'reassign_velocities': False,
                         'remove_com': False,
                         'timestep': <Quantity(4, 'femtosecond')>},
 'output_settings_em': {'compressed_x_grps': '',
                        'compressed_x_precision': 1000,
                        'cpt_file': 'em.cpt',
                        'edr_file': 'em.edr',
                        'energygrps': '',
                        'forcefield_cache': 'db.json',
                        'gro_file': 'em.gro',
                        'log_file': 'em.log',
                        'mdp_file': 'em.mdp',
                        'nstcalcenergy': 100,
                        'nstenergy': 1000,
                        'nstfout': 0,
                        'nstlog': 1000,
                        'nstvout': 0,
                        'nstxout': 0,
                        'nstxout_compressed': 0,
                        'tpr_file': 'em.tpr',
                        'trr_file': 'em.trr',
                        'xtc_file': 'em.xtc'},
 'output_settings_npt': {'compressed_x_grps': '',
                         'compressed_x_precision': 1000,
                         'cpt_file': 'npt.cpt',
                         'edr_file': 'npt.edr',
                         'energygrps': '',
                         'forcefield_cache': 'db.json',
                         'gro_file': 'npt.gro',
                         'log_file': 'npt.log',
                         'mdp_file': 'npt.mdp',
                         'nstcalcenergy': 100,
                         'nstenergy': 1000,
                         'nstfout': 5000,
                         'nstlog': 1000,
                         'nstvout': 5000,
                         'nstxout': 5000,
                         'nstxout_compressed': 5000,
                         'tpr_file': 'npt.tpr',
                         'trr_file': 'npt.trr',
                         'xtc_file': 'npt.xtc'},
 'output_settings_nvt': {'compressed_x_grps': '',
                         'compressed_x_precision': 1000,
                         'cpt_file': 'nvt.cpt',
                         'edr_file': 'nvt.edr',
                         'energygrps': '',
                         'forcefield_cache': 'db.json',
                         'gro_file': 'nvt.gro',
                         'log_file': 'nvt.log',
                         'mdp_file': 'nvt.mdp',
                         'nstcalcenergy': 100,
                         'nstenergy': 1000,
                         'nstfout': 0,
                         'nstlog': 1000,
                         'nstvout': 0,
                         'nstxout': 0,
                         'nstxout_compressed': 0,
                         'tpr_file': 'nvt.tpr',
                         'trr_file': 'nvt.trr',
                         'xtc_file': 'nvt.xtc'},
 'partial_charge_settings': {'nagl_model': None,
                             'number_of_conformers': None,
                             'off_toolkit_backend': 'ambertools',
                             'partial_charge_method': 'am1bcc'},
 'protocol_repeats': 1,
 'simulation_settings_em': {'DispCorr': 'EnerPres',
                            'constraint_algorithm': 'lincs',
                            'constraints': 'h-bonds',
                            'coulombtype': 'PME',
                            'cutoff_scheme': 'verlet',
                            'ewald_rtol': 1e-05,
                            'integrator': 'steep',
                            'lincs_iter': 1,
                            'lincs_order': 12,
                            'mass_repartition_factor': 1,
                            'nsteps': 10,
                            'nstlist': 10,
                            'pme_order': 4,
                            'rcoulomb': <Quantity(1.2, 'nanometer')>,
                            'rlist': <Quantity(1, 'nanometer')>,
                            'rvdw': <Quantity(1.0, 'nanometer')>,
                            'shake_tol': 0.0001,
                            'vdw_modifier': 'Potential-shift',
                            'vdwtype': 'Cut-off'},
 'simulation_settings_npt': {'DispCorr': 'EnerPres',
                             'constraint_algorithm': 'lincs',
                             'constraints': 'h-bonds',
                             'coulombtype': 'PME',
                             'cutoff_scheme': 'verlet',
                             'dt': <Quantity(0.002, 'picosecond')>,
                             'ewald_rtol': 1e-05,
                             'gen_seed': -1,
                             'gen_temp': <Quantity(298.15, 'kelvin')>,
                             'gen_vel': 'no',
                             'integrator': 'sd',
                             'ld_seed': -1,
                             'lincs_iter': 1,
                             'lincs_order': 12,
                             'mass_repartition_factor': 1,
                             'nsteps': 10,
                             'nstlist': 10,
                             'pcoupl': 'C-rescale',
                             'pme_order': 4,
                             'rcoulomb': <Quantity(1.2, 'nanometer')>,
                             'ref_p': <Quantity(1.01325, 'bar')>,
                             'ref_t': <Quantity(298.15, 'kelvin')>,
                             'refcoord_scaling': 'no',
                             'rlist': <Quantity(1, 'nanometer')>,
                             'rvdw': <Quantity(1.0, 'nanometer')>,
                             'shake_tol': 0.0001,
                             'tcoupl': 'no',
                             'vdw_modifier': 'Potential-shift',
                             'vdwtype': 'Cut-off'},
 'simulation_settings_nvt': {'DispCorr': 'EnerPres',
                             'constraint_algorithm': 'lincs',
                             'constraints': 'h-bonds',
                             'coulombtype': 'PME',
                             'cutoff_scheme': 'verlet',
                             'dt': <Quantity(0.002, 'picosecond')>,
                             'ewald_rtol': 1e-05,
                             'gen_seed': -1,
                             'gen_temp': <Quantity(298.15, 'kelvin')>,
                             'gen_vel': 'yes',
                             'integrator': 'sd',
                             'ld_seed': -1,
                             'lincs_iter': 1,
                             'lincs_order': 12,
                             'mass_repartition_factor': 1,
                             'nsteps': 10,
                             'nstlist': 10,
                             'pcoupl': 'no',
                             'pme_order': 4,
                             'rcoulomb': <Quantity(1.2, 'nanometer')>,
                             'ref_p': <Quantity(1.01325, 'bar')>,
                             'ref_t': <Quantity(298.15, 'kelvin')>,
                             'refcoord_scaling': 'no',
                             'rlist': <Quantity(1, 'nanometer')>,
                             'rvdw': <Quantity(1.0, 'nanometer')>,
                             'shake_tol': 0.0001,
                             'tcoupl': 'no',
                             'vdw_modifier': 'Potential-shift',
                             'vdwtype': 'Cut-off'},
 'solvation_settings': {'box_shape': 'dodecahedron',
                        'box_size': None,
                        'box_vectors': None,
                        'number_of_solvent_molecules': None,
                        'solvent_model': 'tip3p',
                        'solvent_padding': <Quantity(1.2, 'nanometer')>},
 'thermo_settings': {'ph': None,
                     'pressure': <Quantity(0.986923267, 'standard_atmosphere')>,
                     'redox_potential': None,
                     'temperature': <Quantity(298.15, 'kelvin')>},
 'top': 'system.top'}

3. Creating a Protocol#

The actual simulation is performed by a `Protocol <https://docs.openfree.energy/en/stable/guide/models/execution.html#protocols-and-the-execution-model>`__.

With the Settings inspected and adjusted, we can provide these to the Protocol. Here, the Gromacs-based MD Protocol is named GromacsMDProtocol.

[4]:
# Creating the Protocol
from openfe_gromacs.protocols.gromacs_md.md_methods import GromacsMDProtocol
protocol = GromacsMDProtocol(settings=settings)

4. Running the MD simulation#

Here we will show you how to run the MD simulation using our Python API.

The MD simulation can be run by executing the ProtocolDAG. The ProtocolDAG is created using the protocol.create() method and requires as input the ChemicalSystem.

Note: we use the shared_basedir and scratch_basedir argument of execute_DAG in order to set the directory where the simulation files are written to.

[5]:
import gufe
import pathlib

# Creating the Protocol
protocol = GromacsMDProtocol(settings=settings)
# Creating the Protocol DAG
dag = protocol.create(stateA=system, stateB=system, mapping=None)
workdir = pathlib.Path('')
# Running the MD simulations
dagres = gufe.protocols.execute_DAG(
    dag,
    shared_basedir=workdir,
    scratch_basedir=workdir,
    keep_shared=True, # set this to True to save the outputs
    n_retries=3
)
/home/ialibay/github/openfe-gromacs/openfe_gromacs/protocols/gromacs_md/md_methods.py:797: UserWarning: Environment variable INTERCHANGE_EXPERIMENTAL=1 is set for using the interchange functionality 'from_openmm' which is not well tested yet.
  warnings.warn(war)
                :-) GROMACS - gmx grompp, 2024.3-conda_forge (-:

Executable:   /home/ialibay/software/mambaforge/install/envs/openfe_gromacs/bin.AVX2_256/gmx
Data prefix:  /home/ialibay/software/mambaforge/install/envs/openfe_gromacs
Working dir:  /home/ialibay/github/openfe-gromacs/examples/md
Command line:
  gmx grompp -f shared_GromacsMDSetupUnit-9bca6ec36249416695cebd72db5a0d47_attempt_0/em.mdp -c shared_GromacsMDSetupUnit-9bca6ec36249416695cebd72db5a0d47_attempt_0/system.gro -p shared_GromacsMDSetupUnit-9bca6ec36249416695cebd72db5a0d47_attempt_0/system.top -o shared_GromacsMDRunUnit-552d5bfdff8041bc88ad20e8653d514e_attempt_0/em.tpr

Generating 1-4 interactions: fudge = 0.5
Number of degrees of freedom in T-Coupling group rest is 53927.00
The integrator does not provide a ensemble temperature, there is no system ensemble temperature

GROMACS reminds you: "If humanity has fled shivering from the starry spaces, it has become minutely at home in the interstices of the speck that it inhabits for an instant" (George H. Mead)

                :-) GROMACS - gmx mdrun, 2024.3-conda_forge (-:

Executable:   /home/ialibay/software/mambaforge/install/envs/openfe_gromacs/bin.AVX2_256/gmx
Data prefix:  /home/ialibay/software/mambaforge/install/envs/openfe_gromacs
Working dir:  /home/ialibay/github/openfe-gromacs/examples/md/shared_GromacsMDRunUnit-552d5bfdff8041bc88ad20e8653d514e_attempt_0
Command line:
  gmx mdrun -s em.tpr -cpo em.cpt -o em.trr -x em.xtc -c em.gro -e em.edr -g em.log -ntmpi 1 -ntomp 1 -pme auto -pmefft auto -bonded auto -nb auto -update auto

The current CPU can measure timings more accurately than the code in
gmx mdrun was configured to use. This might affect your simulation
speed as accurate timings are needed for load-balancing.
Please consider rebuilding gmx mdrun with the GMX_USE_RDTSCP=ON CMake option.
Reading file em.tpr, VERSION 2024.3-conda_forge (single precision)

Update groups can not be used for this system because there are three or more consecutively coupled constraints

Using 1 MPI thread
Using 1 OpenMP thread

Setting the LD random seed to -337379401

Generated 3457135 of the 3457135 non-bonded parameter combinations

Generated 3457135 of the 3457135 1-4 parameter combinations

Excluding 3 bonded neighbours molecule type 'MOL0'

turning H bonds into constraints...

Excluding 3 bonded neighbours molecule type 'MOL1'

turning H bonds into constraints...

Excluding 3 bonded neighbours molecule type 'MOL2'

turning H bonds into constraints...

Excluding 3 bonded neighbours molecule type 'MOL7874'

turning H bonds into constraints...

Excluding 3 bonded neighbours molecule type 'MOL7895'

turning H bonds into constraints...
Analysing residue names:
There are:   154    Protein residues
There are:     1      Other residues
There are:  7872      Water residues
There are:    50        Ion residues
Analysing Protein...
Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups...

The largest distance between excluded atoms is 0.418 nm between atom 98 and 105
Calculating fourier grid dimensions for X Y Z
Using a fourier grid of 64x64x64, spacing 0.115 0.115 0.115

Estimate for the relative computational load of the PME mesh part: 0.20

This run will generate roughly 2 Mb of data

NOTE: Thread affinity was not set.

Steepest Descents:
   Tolerance (Fmax)   =  1.00000e+01
   Number of steps    =           10

Energy minimization reached the maximum number of steps before the forces
reached the requested precision Fmax < 10.

writing lowest energy coordinates.

Steepest Descents did not converge to Fmax < 10 in 11 steps.
Potential Energy  = -3.1639372e+05
Maximum force     =  1.4821552e+04 on atom 2551
Norm of force     =  4.0669904e+02

GROMACS reminds you: "If humanity has fled shivering from the starry spaces, it has become minutely at home in the interstices of the speck that it inhabits for an instant" (George H. Mead)

                :-) GROMACS - gmx grompp, 2024.3-conda_forge (-:

Executable:   /home/ialibay/software/mambaforge/install/envs/openfe_gromacs/bin.AVX2_256/gmx
Data prefix:  /home/ialibay/software/mambaforge/install/envs/openfe_gromacs
Working dir:  /home/ialibay/github/openfe-gromacs/examples/md
Command line:
  gmx grompp -f shared_GromacsMDSetupUnit-9bca6ec36249416695cebd72db5a0d47_attempt_0/nvt.mdp -c shared_GromacsMDRunUnit-552d5bfdff8041bc88ad20e8653d514e_attempt_0/em.gro -p shared_GromacsMDSetupUnit-9bca6ec36249416695cebd72db5a0d47_attempt_0/system.top -o shared_GromacsMDRunUnit-552d5bfdff8041bc88ad20e8653d514e_attempt_0/nvt.tpr

Generating 1-4 interactions: fudge = 0.5
Number of degrees of freedom in T-Coupling group System is 53927.00

GROMACS reminds you: "Remember, being healthy is basically dying as slowly as possible." (Ricky Gervais)

                :-) GROMACS - gmx mdrun, 2024.3-conda_forge (-:

Executable:   /home/ialibay/software/mambaforge/install/envs/openfe_gromacs/bin.AVX2_256/gmx
Data prefix:  /home/ialibay/software/mambaforge/install/envs/openfe_gromacs
Working dir:  /home/ialibay/github/openfe-gromacs/examples/md/shared_GromacsMDRunUnit-552d5bfdff8041bc88ad20e8653d514e_attempt_0
Command line:
  gmx mdrun -s nvt.tpr -cpo nvt.cpt -o nvt.trr -x nvt.xtc -c nvt.gro -e nvt.edr -g nvt.log -ntmpi 1 -ntomp 1 -pme auto -pmefft auto -bonded auto -nb auto -update auto

Setting the LD random seed to -155189377

Generated 3457135 of the 3457135 non-bonded parameter combinations

Generated 3457135 of the 3457135 1-4 parameter combinations

Excluding 3 bonded neighbours molecule type 'MOL0'

turning H bonds into constraints...

Excluding 3 bonded neighbours molecule type 'MOL1'

turning H bonds into constraints...

Excluding 3 bonded neighbours molecule type 'MOL2'

turning H bonds into constraints...

Excluding 3 bonded neighbours molecule type 'MOL7874'

turning H bonds into constraints...

Excluding 3 bonded neighbours molecule type 'MOL7895'

turning H bonds into constraints...

Setting gen_seed to 1870101402

Velocities were taken from a Maxwell distribution at 298.15 K
Analysing residue names:
There are:   154    Protein residues
There are:     1      Other residues
There are:  7872      Water residues
There are:    50        Ion residues
Analysing Protein...
Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups...

The largest distance between excluded atoms is 0.416 nm between atom 98 and 105

Determining Verlet buffer for a tolerance of 0.005 kJ/mol/ps at 298.15 K

Calculated rlist for 1x1 atom pair-list as 1.232 nm, buffer size 0.032 nm

Set rlist, assuming 4x4 atom pair-list, to 1.200 nm, buffer size 0.000 nm

Note that mdrun will redetermine rlist based on the actual pair-list setup
Calculating fourier grid dimensions for X Y Z
Using a fourier grid of 64x64x64, spacing 0.115 0.115 0.115

Estimate for the relative computational load of the PME mesh part: 0.22

This run will generate roughly 2 Mb of data
The current CPU can measure timings more accurately than the code in
gmx mdrun was configured to use. This might affect your simulation
speed as accurate timings are needed for load-balancing.
Please consider rebuilding gmx mdrun with the GMX_USE_RDTSCP=ON CMake option.
Reading file nvt.tpr, VERSION 2024.3-conda_forge (single precision)
Changing nstlist from 10 to 80, rlist from 1.2 to 1.314


Update groups can not be used for this system because there are three or more consecutively coupled constraints

Using 1 MPI thread
Using 1 OpenMP thread


NOTE: Thread affinity was not set.
starting mdrun 'FOO'
10 steps,      0.0 ps.

Writing final coordinates.

NOTE: 1 % of the run time was spent in domain decomposition,
      12 % of the run time was spent in pair search,
      you might want to increase nstlist (this has no effect on accuracy)

               Core t (s)   Wall t (s)        (%)
       Time:        0.261        0.261      100.0
                 (ns/day)    (hour/ns)
Performance:        7.285        3.294

GROMACS reminds you: "Clickety Clickety Click" (System Manager From Hell)

                :-) GROMACS - gmx grompp, 2024.3-conda_forge (-:

Executable:   /home/ialibay/software/mambaforge/install/envs/openfe_gromacs/bin.AVX2_256/gmx
Data prefix:  /home/ialibay/software/mambaforge/install/envs/openfe_gromacs
Working dir:  /home/ialibay/github/openfe-gromacs/examples/md
Command line:
  gmx grompp -f shared_GromacsMDSetupUnit-9bca6ec36249416695cebd72db5a0d47_attempt_0/npt.mdp -c shared_GromacsMDRunUnit-552d5bfdff8041bc88ad20e8653d514e_attempt_0/nvt.gro -p shared_GromacsMDSetupUnit-9bca6ec36249416695cebd72db5a0d47_attempt_0/system.top -o shared_GromacsMDRunUnit-552d5bfdff8041bc88ad20e8653d514e_attempt_0/npt.tpr

Generating 1-4 interactions: fudge = 0.5
Number of degrees of freedom in T-Coupling group System is 53927.00

GROMACS reminds you: "I managed to get two hours of work done before work" (E. Lindahl)

                :-) GROMACS - gmx mdrun, 2024.3-conda_forge (-:

Executable:   /home/ialibay/software/mambaforge/install/envs/openfe_gromacs/bin.AVX2_256/gmx
Data prefix:  /home/ialibay/software/mambaforge/install/envs/openfe_gromacs
Working dir:  /home/ialibay/github/openfe-gromacs/examples/md/shared_GromacsMDRunUnit-552d5bfdff8041bc88ad20e8653d514e_attempt_0
Command line:
  gmx mdrun -s npt.tpr -cpo npt.cpt -o npt.trr -x npt.xtc -c npt.gro -e npt.edr -g npt.log -ntmpi 1 -ntomp 1 -pme auto -pmefft auto -bonded auto -nb auto -update auto

The current CPU can measure timings more accurately than the code in
gmx mdrun was configured to use. This might affect your simulation
speed as accurate timings are needed for load-balancing.
Please consider rebuilding gmx mdrun with the GMX_USE_RDTSCP=ON CMake option.
Reading file npt.tpr, VERSION 2024.3-conda_forge (single precision)
Setting the LD random seed to -1862836486

Generated 3457135 of the 3457135 non-bonded parameter combinations

Generated 3457135 of the 3457135 1-4 parameter combinations

Excluding 3 bonded neighbours molecule type 'MOL0'

turning H bonds into constraints...

Excluding 3 bonded neighbours molecule type 'MOL1'

turning H bonds into constraints...

Excluding 3 bonded neighbours molecule type 'MOL2'

turning H bonds into constraints...

Excluding 3 bonded neighbours molecule type 'MOL7874'

turning H bonds into constraints...

Excluding 3 bonded neighbours molecule type 'MOL7895'

turning H bonds into constraints...
Analysing residue names:
There are:   154    Protein residues
There are:     1      Other residues
There are:  7872      Water residues
There are:    50        Ion residues
Analysing Protein...
Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups...

The largest distance between excluded atoms is 0.421 nm between atom 1612 and 1622

Determining Verlet buffer for a tolerance of 0.005 kJ/mol/ps at 298.15 K

Calculated rlist for 1x1 atom pair-list as 1.232 nm, buffer size 0.032 nm

Set rlist, assuming 4x4 atom pair-list, to 1.200 nm, buffer size 0.000 nm

Note that mdrun will redetermine rlist based on the actual pair-list setup
Calculating fourier grid dimensions for X Y Z
Using a fourier grid of 64x64x64, spacing 0.115 0.115 0.115

Estimate for the relative computational load of the PME mesh part: 0.22

This run will generate roughly 3 Mb of data
Changing nstlist from 10 to 80, rlist from 1.2 to 1.314


Update groups can not be used for this system because there are three or more consecutively coupled constraints

Using 1 MPI thread
Using 1 OpenMP thread


NOTE: Thread affinity was not set.
starting mdrun 'FOO'
10 steps,      0.0 ps.

Writing final coordinates.

NOTE: 1 % of the run time was spent in domain decomposition,
      12 % of the run time was spent in pair search,
      you might want to increase nstlist (this has no effect on accuracy)

               Core t (s)   Wall t (s)        (%)
       Time:        0.254        0.254      100.0
                 (ns/day)    (hour/ns)
Performance:        7.493        3.203

GROMACS reminds you: "I always think there is something foreign about jolly phrases at breakfast." (Mr. Carson in Downtown Abbey)

Following files were created for the setup of the MD run:

[6]:
!ls shared_GromacsMDSetupUnit-*
db.json  em.mdp  npt.mdp  nvt.mdp  system.gro  system.top

Following files were created by the energy minimization, NVT, and NPT MD simulations:

[7]:
!ls shared_GromacsMDRunUnit-8
ls: cannot access 'shared_GromacsMDRunUnit-8': No such file or directory