NoneqGrandCanonicalMonteCarloSamplerMPI

class grandfep.sampler.NoneqGrandCanonicalMonteCarloSamplerMPI(system, topology, temperature, collision_rate, timestep, log, platform=<openmm.openmm.Platform; proxy of <Swig Object of type 'OpenMM::Platform *'> >, water_resname='HOH', water_O_name='O', position=None, chemical_potential=None, standard_volume=None, sphere_radius=Quantity(value=10.0, unit=angstrom), reference_atoms=None, rst_file='md.rst7', dcd_file=None, append_dcd=False, jsonl_file='md.jsonl', init_lambda_state=None, lambda_dict=None)

Bases: _ReplicaExchangeMixin, NoneqGrandCanonicalMonteCarloSampler

Nonequilibrium Grand Canonical Monte Carlo (Noneq-GCMC) sampler.

In this class, GCMC is achieved by performing insertion/deletion in a nonequilibrium candidate Monte Carlo manner. The work value of the insertion/deletion will be used to evaluate the acceptance ratio. Insertion/deletion can either be performed in the whole box or in a sub-volume (active site) of the box. In an equilibrium sampling, when the water is free to move in/out of the sub-volume, I recommend alternating between GCMC and box.

Parameters:

  • system (openmm.openmm.System) – The OpenMM System object. Must include CustomNonbondedForce and NonbondedForce with appropriate per-particle parameters and global parameter definitions. Call system.setDefaultPeriodicBoxVectors() before passing it to this class, the newly created context will use this box vectors. Only a rectangular box is supported.

  • topology (openmm.app.topology.Topology) – The OpenMM Topology object. Must contain water molecules with the specified residue and atom names. Must have the correct boxVector. Only a rectangular box is supported.

  • temperature (openmm.unit.quantity.Quantity) – The reference temperature for the system, with proper units (e.g., kelvin).

  • collision_rate (openmm.unit.quantity.Quantity) – The collision rate (friction) for the Langevin integrator, with time units.

  • timestep (openmm.unit.quantity.Quantity) – The timestep for the integrator, with time units (e.g., femtoseconds).

  • log (typing.Union[str, pathlib.Path]) – Path to the log file. This file will be opened in appended mode.

  • platform (openmm.openmm.Platform) – The OpenMM computational platform to use. Default is CUDA.

  • water_resname (str) – The residue name of water in the topology. Default is ‘HOH’.

  • water_O_name (str) – The atom name of oxygen in water. Default is ‘O’.

  • position (openmm.unit.quantity.Quantity) – Initial position of the system. Need to be provided for box Vectors. Default is None.

  • chemical_potential – Chemical potential of the system, with units. Default is None.

  • standard_volume – Standard volume of a water molecule in the reservoir. with units. Default is None.

  • sphere_radius (openmm.unit.quantity.Quantity) – Radius of the GCMC sphere. Default is 10.0 * unit.angstroms.

  • reference_atoms (list) – A list of atom indices in the topology that will be set as the center of the GCMC sphere. Default is None.

  • rst_file (str) – File name for the restart file.

  • dcd_file (str) – File name for the DCD trajectory file. Default is None, no DCD reporter.

  • append_dcd (bool) – Whether to append to the DCD file. Default is True.

  • jsonl_file (str) – File name for the JSONL file. Default is “md.jsonl”. This file saves the ghost_list. rst7 file needs this jsonl to restart the ghost list, and dcd file needs this jsonl to remove the ghost water.

  • init_lambda_state (int) – Lambda state index for this replica, counting from 0

  • lambda_dict (dict) – A dictionary of mapping from global parameters to their values in all the sampling states.

static get_particle_parameter_index_cust_nb_force(custom_nonbonded_force)

Get the indices of is_real and is_switching parameters in the perParticleParameters of self.custom_nonbonded_force.

Parameters:

custom_nonbonded_force (openmm.openmm.CustomNonbondedForce) – The CustomNonbondedForce to be checked.

Returns:

A tuple (is_real_index, is_switching_index) indicating the indices of the corresponding parameters in self.custom_nonbonded_force.getPerParticleParameterName(i).

Return type:

tuple of int

check_ghost_list()

Loop over all water particles in the system to validate that self.ghost_list correctly reflects the current ghost and switching water configuration.

  • Ghost water:

    1. is_real = 0.0 and is_switching = 0.0 in every CustomNonbondedForce in self.custom_nonbonded_force_list

    2. Charge = 0.0 in self.nonbonded_force

  • Real and not switching water:

    1. is_real = 1.0 in every CustomNonbondedForce in self.custom_nonbonded_force_list

    2. Charge = Proper_water_charge in self.nonbonded_force

  • Switching water:

    1. is_real = 1.0 in every CustomNonbondedForce in self.custom_nonbonded_force_list

    2. Charge = 0.0 in self.nonbonded_force

The switching water should not be present in self.ghost_list.

Raises:

ValueError – If any condition for ghost, real, or switching water is violated.

Return type:

None

check_switching()

Loop over all water particles in the system to make sure the one water given by self.switching_water is the only one that is switching.

  • Switching water:

    1. is_real = 1.0 and is_switching = 1.0 in self.custom_nonbonded_force

    2. Charge = 0.0 in self.nonbonded_force (ParticleParameters)

    3. chargeScale = proper_water_charge in ParticleParameterOffsets

  • Non-switching waters:

    1. is_switching = 0.0 in self.custom_nonbonded_force

    2. No ParticleParameterOffsets in self.nonbonded_force

Raises:

ValueError – If any switching or non-switching water fails to meet the expected criteria.

Return type:

None

customise_force_amber(system)

In Amber, NonbondedForce handles both electrostatics and vdW. This function will remove vdW from NonbondedForce and create a list of CustomNonbondedForce (in this case, only 1 in the list) to handle vdW, so that the interaction can be switched off for certain water.

Parameters:

system (openmm.openmm.System) – openmm.System The system to be converted.

Return type:

None

Returns:

None

customise_force_charmm(system)

In Charmm, NonbondedForce handles electrostatics, and CustomNonbondedForce handles vdW. For vdW, this function will add perParticle parameters ‘is_real’, ‘is_switching’, global parameter lambda_gc_vdw to the CustomNonbondedForce. For Coulomb, this function will add ParticleParameterOffset to the switching water and lambda_gc_coulomb to the NonbondedForce.

The CustomNonbondedForce should have the following energy expression:

‘(a/r6)^2-b/r6; r6=r^6;a=acoef(type1, type2);b=bcoef(type1, type2)’ or ‘acoef(type1, type2)/r^12 - bcoef(type1, type2)/r^6;’

Parameters:

system (openmm.openmm.System) – openmm.System The system to be converted.

Return type:

None

Returns:

None

customise_force_hybrid(system)

If the system is Hybrid, this function will add perParticleParameters is_real and is_switching to the custom_nonbonded_force (openmm.openmm.CustomNonbondedForce) for vdw.

Groups

old

core

new

fix

wat

swit

old

C1

core

C1

C1

new

None

C1

C1

fix

C1

C1

C1

C4

wat

C1

C1

C1

C3

C3

swit

C1

C1

C1

C2

C2

C2
Parameters:

system (openmm.openmm.System) – The system to be converted.

Return type:

None

get_ghost_list(check_system=False)

Get a copy of the current ghost water list.

If check_system is True, self.check_ghost_list will be called to validate the consistency of self.ghost_list with self.custom_nonbonded_force and self.nonbonded_force. If the validation fails, a ValueError will be raised.

Parameters:

check_system (bool) – If True, self.check_ghost_list() will be called to perform consistency checks in the system forces to ensure that the ghost list is correct. Default is False.

Returns:

A copy of self.ghost_list, which contains the residue indices of ghost waters.

Return type:

list

get_particle_parameter_offset_dict()

Retrieve all ParticleParameterOffset entries from self.nonbonded_force.

Returns:

A dictionary mapping atom indices to their corresponding ParticleParameterOffset data.

The key is the atom index. Each value is a list with the following structure: [ param_offset_index, global_parameter_name, atom_index, chargeScale, sigmaScale, epsilonScale ]

Return type:

dict

get_sphere_center(positions)

Average the reference atoms positions

Returns:

sphere_center – The center of the GCMC sphere, with unit

Return type:

unit.Quantity

get_water_state(positions)

Check whether the water molecule is inside the box (1) or not (0).

Parameters:

positions (unit.Quantity) – The current positions of the system.

Return type:

tuple[dict, numpy.array]

Returns:

  • water_state_dict (dict) – A dictionary label whether the water molecule is inside the box (1) or not (0).

  • dist_all_o (np.array) – All the distance bewteen water oxygen and the center of the sphere.

load_rst(rst_input)

Load positions/velocities from a restart file. boxVector will be checked, if not match, raise an error.

Parameters:

rst_input (Union[str, Path]) – Amber inpcrd file

Return type:

None

move_deletion_GCMC(l_vdw_list, l_chg_list, n_prop)

Perform a non-equilibrium deletion, according the given lambda schedule

Parameters:

  • l_vdw_list (list) – A list of lambda_gc_vdw value that defines the path of insertion. It should start with 1.0 and end with 0.0.

  • l_chg_list (list) – A list of lambda_gc_coulomb value that defines the path of insertion. It should start with 1.0 and end with 0.0.

  • n_prop (int) – Number of propagation steps (equilibrium MD) between each lambda switching.

Return type:

tuple[bool, float, openmm.unit.quantity.Quantity, list, int, bool]

Returns:

  • accept (bool) – Whether the insertion is accepted.

  • acc_prob (float) – The acceptance probability of the insertion.

  • protocol_work (unit.Quantity) – The work done during the insertion process, with unit.

  • protocol_work_list (list) – A list of work values during each perturbation step, in kcal/mol. no unit in this list.

  • n_water (int) – The number of water molecules in the system after the insertion.

  • switching_water_inside (bool) – Whether the swithching water is inside the GCMC sphere in the end point of non-eq process

move_deletion_box(l_vdw_list, l_chg_list, n_prop)

Perform a non-equilibrium deletion, according the given lambda schedule

Parameters:

  • l_vdw_list (list) – A list of lambda_gc_vdw value that defines the path of insertion. It should start with 1.0 and end with 0.0.

  • l_chg_list (list) – A list of lambda_gc_coulomb value that defines the path of insertion. It should start with 1.0 and end with 0.0.

  • n_prop (int) – Number of propagation steps (equilibrium MD) between each lambda switching.

Return type:

tuple[bool, float, openmm.unit.quantity.Quantity, list, int]

Returns:

  • accept (bool) – Whether the insertion is accepted.

  • acc_prob (float) – The acceptance probability of the insertion.

  • protocol_work (unit.Quantity) – The work done during the insertion process, with unit.

  • protocol_work_list (list) – A list of work values during each perturbation step, in kcal/mol. no unit in this list.

  • n_water (int) – The number of water molecules in the system after the insertion.

move_insert_delete(l_vdw_list, l_chg_list, n_prop, box)

Randomly do insertion or deletion in 1:1 probability

Parameters:

  • l_vdw_list (list) – A list of lambda_gc_vdw value that defines the path of insertion. It should start with 0.0 and end with 1.0.

  • l_chg_list (list) – A list of lambda_gc_coulomb value that defines the path of insertion. It should start with 0.0 and end with 1.0.

  • n_prop (int) – Number of propagation steps (equilibrium MD) between each lambda switching.

  • box (bool) – GC in box or GCMC sphere

Return type:

tuple[bool, float, openmm.unit.quantity.Quantity, list, int, bool]

Returns:

  • accept (bool) – Whether the insertion is accepted.

  • acc_prob (float) – The acceptance probability of the insertion.

  • protocol_work (unit.Quantity) – The work done during the insertion process, with unit.

  • protocol_work_list (list) – A list of work values during each perturbation step, in kcal/mol. no unit in this list.

  • n_water (int) – The number of water molecules in the system after the insertion.

  • switching_water_inside (bool) – Whether the swithching water is inside reasonable region in the end point of non-eq process

move_insertion_GCMC(l_vdw_list, l_chg_list, n_prop)

Perform a non-equilibrium insertion, according the given lambda schedule

Parameters:

  • l_vdw_list (list) – A list of lambda_gc_vdw value that defines the path of insertion. It should start with 0.0 and end with 1.0.

  • l_chg_list (list) – A list of lambda_gc_coulomb value that defines the path of insertion. It should start with 0.0 and end with 1.0.

  • n_prop (int) – Number of propagation steps (equilibrium MD) between each lambda switching.

Return type:

tuple[bool, float, openmm.unit.quantity.Quantity, list, int, bool]

Returns:

  • accept (bool) – Whether the insertion is accepted.

  • acc_prob (float) – The acceptance probability of the insertion.

  • protocol_work (unit.Quantity) – The work done during the insertion process, with unit.

  • protocol_work_list (list) – A list of work values during each perturbation step, in kcal/mol. no unit in this list.

  • n_water (int) – The number of water molecules in the system after the insertion.

  • switching_water_inside (bool) – Whether the swithching water is inside the GCMC sphere in the end point of non-eq process

move_insertion_box(l_vdw_list, l_chg_list, n_prop)

Perform a non-equilibrium insertion, according the given lambda schedule

Parameters:

  • l_vdw_list (list) – A list of lambda_gc_vdw value that defines the path of insertion. It should start with 0.0 and end with 1.0.

  • l_chg_list (list) – A list of lambda_gc_coulomb value that defines the path of insertion. It should start with 0.0 and end with 1.0.

  • n_prop (int) – Number of propagation steps (equilibrium MD) between each lambda switching.

Return type:

tuple[bool, float, openmm.unit.quantity.Quantity, list, int]

Returns:

  • accept (bool) – Whether the insertion is accepted.

  • acc_prob (float) – The acceptance probability of the insertion.

  • protocol_work (unit.Quantity) – The work done during the insertion process, with unit.

  • protocol_work_list (list) – A list of work values during each perturbation step, in kcal/mol. no unit in this list.

  • n_water (int) – The number of water molecules in the system after the insertion.

random_place_water(state, res_index, sphere_center=None)

Shift the coordinate+orientation of a water molecule to a random place. If the center is None, a random position in the box will be selected. If the center is not None, a randomed position in the GCMC sphere will be selected.

Parameters:

  • state (openmm.State) – The current state of the simulation context. The positions will be read from this state.

  • res_index (int) – The residue index of the water molecule to be shifted.

  • sphere_center (unit.Quantity) – The center of the GCMC sphere. If None, a random position in the box will be selected. Default is None.

Return type:

tuple[openmm.unit.quantity.Quantity, openmm.unit.quantity.Quantity]

Returns:

  • positions (unit.Quantity) – The new positions with the water molecule shifted.

  • velocities (unit.Quantity) – The new velocities with the water molecule shifted.

replica_exchange(calc_neighbor_only=True)

Perform one neighbor swap replica exchange. In odd RE steps, attempt exchange between 0-1, 2-3, … In even RE steps, attempt exchange between 1-2, 3-4, … If RE is accepted, update the position/boxVector/velocity

Parameters:

calc_neighbor_only (bool) – If True, only the nearest neighbor will be calculated.

Returns:

  • reduced_energy_matrix (np.array) – A numpy array with the size of (n_lambda_states, n_simulated_states)

  • re_decision (dict) – The exchange decision between all the pairs. e.g., {(0,1):(True, ratio_0_1), (2,3):(True, ratio_2_3)}. The key is the MPI rank pairs, and the value is the decision and the acceptance ratio.

replica_exchange_global_param(calc_neighbor_only=True)

Perform one neighbor swap replica exchange. In odd RE steps, attempt exchange between 0-1, 2-3, … In even RE steps, attempt exchange between 1-2, 3-4, … If RE is accepted, update the position/boxVector/velocity

Parameters:

calc_neighbor_only (bool) – If True, only the nearest neighbor will be calculated.

Returns:

  • reduced_energy_matrix (np.array) – A numpy array with the size of (n_lambda_states, n_simulated_states)

  • re_decision (dict) – The exchange decision between all the pairs. e.g., {(0,1):(True, ratio_0_1), (2,3):(True, ratio_2_3)}. The key is the MPI rank pairs, and the value is the decision and the acceptance ratio.

report_dcd(state=None)

gather the coordinates from all MPI rank, only rank 0 writes the file. Report both dcd and rst7 files.

Parameters:

state (openmm.openmm.State) – State of the simulation. If None, it will get the current state from the simulation context.

report_dcd_rank0(state=None)

gather the coordinates from all MPI rank, only rank 0 writes the file.

Parameters:

state (openmm.openmm.State) – State of the simulation. If None, it will get the current state from the simulation context.

report_jsonl(dcd)

Append a line to the jsonl file. An example line is like this. GC_count is the current number of GC step. ghost_list is the ghost water residue 0-index. dcd tells whether there is a dcd frame saved with this line of recording

{"GC_count": 50, "ghost_list": [11, 12, 13], "dcd": 0}
Parameters:

dcd – Whether there is a dcd writing in this step for appending a line to jsonl. 1: yes, 0: no

Return type:

None

report_rst(state=None)

gather the coordinates from all MPI rank, only rank 0 write the file.

Parameters:

state (openmm.openmm.State) – XXX

report_rst_rank0(state=None)

gather the coordinates from all MPI rank, only rank 0 write the file.

Parameters:

state (openmm.openmm.State) – XXX

set_ghost_from_jsonl(jsonl_file)

Read the last line of a jsonl file and set the current ghost_list,gc_count.

Parameters:

jsonl_file

Return type:

None

set_ghost_list(ghost_list, check_system=True)

Update the water residues in ghost_list to ghost.

If check_system is True, self.check_ghost_list will be called to validate the consistency of self.ghost_list with self.custom_nonbonded_force_list and self.nonbonded_force. If the validation fails, a ValueError will be raised.

Parameters:

  • ghost_list (list) – A list of residue indices (integers) that should be marked as ghost waters.

  • check_system (bool, optional) – If True, perform validation to ensure the internal force parameters are consistent with the updated ghost_list. Default is True.

Return type:

None

set_lambda_dict(init_lambda_state, lambda_dict)
Set internal attributes :

  • init_lambda_state: (int) Lambda state index for this replica, counting from 0

  • lambda_dict: (dict) Global parameter values for all the sampling states

  • lambda_states_list: (list) Index of the Lambda state which is simulated. All the init_lambda_state in this MPI run

  • n_lambda_states: (int) Number of Lambda state to be sampled

You can have 10 states defined in the lambda_dict, and only simulate 4 replicas. In this case, each value in the lambda_dict will have a length of 10, and n_lambda_states=10, and the lambda_states_list will be a list of 4 integers, continuously increasing between 0 and 9, and init_lambda_state should be one of the values in lambda_states_list.

Parameters:

  • init_lambda_state (int) – The lambda state for this replica, counting from 0

  • lambda_dict (dict) – A dictionary of mapping from global parameters to their values in all the sampling states. For example, {"lambda_gc_vdw": [0.0, 0.5, 1.0, 1.0, 1.0], "lambda_gc_coulomb": [0.0, 0.0, 1.0, 0.5, 1.0]}

Return type:

None

set_lambda_state(i)

Set the lambda state index for this replica.

Parameters:

i (int) – The lambda state index for this replica, counting from 0

Return type:

None

set_re_step(re_step)
Parameters:

re_step (int) – Number of replica exchanges to perform

Return type:

None

set_re_step_from_log(log_input)

Read the previous log file and determine the RE step.

Parameters:

log_input (Union[str, Path]) – previous log file where the RE step counting should continue.

Return type:

None

update_gc_count(insert_delete, work, box_GCMC, n_water, accept, acc_prob)

Update the internal gc_count dictionary that tracks GCMC insertion/deletion statistics.

Parameters:

  • insert_delete (int) – Indicator for the type of move: 0 for insertion, 1 for deletion.

  • work (unit.Quantity) – Work value (e.g., from nonequilibrium switching), with units. The value in kcal/mol will be recorded.

  • box_GCMC (int) – Region type indicator: 0 for the whole simulation box, 1 for a GCMC sub-volume.

  • n_water (int) – Number of water molecules after this move in the whole box or in the GCMC sub-volume, depending on box_GCMC.

  • accept (bool) – Indicator of whether the move was accepted

  • acc_prob (float) – Acceptance probability.

Return type:

None

work_process(n_prop, l_vdw_list, l_chg_list)

Perform the work process for insertion/deletion.

Parameters:

  • n_prop (int) – Number of propagation steps (equilibrium MD) between each lambda switching.

  • l_vdw_list (list) – A list of value for lambda_gc_vdw.

  • l_chg_list (list) – A list of value for lambda_gc_coulomb

Return type:

tuple[bool, openmm.unit.quantity.Quantity, list]

Returns:

  • explosion (bool) – Whether the non-equilibrium process get nan.

  • protocol_work (unit.Quantity) – The work done during the insertion process.

  • protocol_work_list (list) – A list of work done during each perturbation step, in kcal/mol. no unit in this list.

Adam_GCMC: openmm.unit.quantity.Quantity

Adam value, considering the selected sphere and water molecules in the sphere.

Adam_box: openmm.unit.quantity.Quantity

Adam value, considering the whole simulation box and all water molecules.

box_vectors

The box vectors of the system. The box has to be rectangular.

chemical_potential: openmm.unit.quantity.Quantity

Chemical potential of the GC particle.

comm

MPI communicator

compound_integrator: openmm.openmm.CompoundIntegrator

2 integrators in this attribute. The 1st one is for Canonical simulation, the 2nd one is for non-equilibirum insertion/deletion.

custom_nonbonded_force_list: list[openmm.openmm.CustomNonbondedForce]
ghost_list

A list of residue indices of water that are set to ghost. Should only be modified with set_ghost_list(), and get by get_ghost_list().

kBT: openmm.unit.quantity.Quantity

kB* T, with unit.

lambda_dict: dict

A dictionary of mapping from global parameters to their values in all the sampling states.

lambda_state_index: int

Lambda state index for this replica, counting from 0

lambda_states_list: list

Index of the Lambda state which is simulated. All the init_lambda_state in this MPI run can be fetched with self.lambda_states_list[rank]

logger

Logger for the Sampler

n_lambda_states: int

Number of Lambda state to be sampled. It should be equal to the length of the values in lambda_dict.

nonbonded_force: openmm.openmm.NonbondedForce

This force handles Coulomb. The switching water has ParticleParameterOffset with global parameter lambda_gc_coulomb to control the switching water.

num_of_points_water: int

The number of points in the water model. 3 for TIP3P, 4 for OPC.

rank

Rank of this process

re_step

Number of replica exchanges performed

reference_atoms

A list of atom indices that will be set as the center of the GCMC sphere.

simulation: openmm.app.simulation.Simulation

Simulation ties together Topology, System, Integrator, and Context in this sampler.

size

Size of the communicator (number of processes)

standard_volume: openmm.unit.quantity.Quantity

Standard volume of a water molecule in the reservoir.

switching_water: int

The residue index of the switching water. The switching water will be set as the last water during initialization. It should not be changed during the simulation, as the ParticleParameterOffset can not be updated in NonbondedForce.

system: openmm.openmm.System

The OpenMM System object.

system_type: str

The type of the system. Can be Amber, Charmm or Hybrid depending on the system and energy expression in the given CustomNonbondedForce.

temperature: openmm.unit.quantity.Quantity

The reference temperature for the system, with proper units (e.g., kelvin).

topology: openmm.app.topology.Topology

The OpenMM Topology object. All the res_name, atom_index, atom_name, etc. are in this topology.

wat_params: dict

A dictionary to track the nonbonded parameter of water. The keys are charge, sigma, epsilon, The values are a list of parameters with unit.

water_res_2_O: dict

A dictionary of residue index to the atom index of water oxygen.

water_res_2_atom: dict

A dictionary of residue index to a list of atom indices of water.