BaseGrandCanonicalMonteCarloSampler

class grandfep.sampler.BaseGrandCanonicalMonteCarloSampler(system, topology, temperature, collision_rate, timestep, log, integrator_str='BAOABIntegrator', platform=<openmm.openmm.Platform; proxy of <Swig Object of type 'OpenMM::Platform *'> >, water_resname='HOH', water_O_name='O', create_simulation=True, optimization='O3', n_split_water='log')

Bases: object

Base class for Grand Canonical Monte Carlo (GCMC) sampling.

This class provides a flexible framework for customizing OpenMM forces so that water molecules can be alchemically inserted (ghost → real) or deleted (real → ghost). Each water is assigned two properties:

  • is_real:

    • 1.0 : real water

    • 0.0 : ghost water.

  • is_switching:

    • 1.0 for the switching water (used for NEQ insertion/deletion)

    • 0.0 for all other waters.

The last water in the system is selected as the switching water (is_real=1, is_switching=1). During an insertion or deletion attempt, this water is swapped (coordinate and velocity) with a ghost or real water, enabling NEQ perturbation moves.

vdW interactions are handled via self.custom_nonbonded_force (openmm.CustomNonbondedForce), where per-particle parameters is_real and is_switching interact with the global parameter lambda_gc_vdw.

Electrostatic interactions are handled by self.nonbonded_force (openmm.NonbondedForce). Ghost waters are given zero charge, and switching waters use ParticleParameterOffset with a global parameter lambda_gc_coulomb.

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

__init__(system, topology, temperature, collision_rate, timestep, log, integrator_str='BAOABIntegrator', platform=<openmm.openmm.Platform; proxy of <Swig Object of type 'OpenMM::Platform *'> >, water_resname='HOH', water_O_name='O', create_simulation=True, optimization='O3', n_split_water='log')
Parameters:

  • system (openmm.openmm.System) – The OpenMM System object. Must include CustomNonbondedForce and NonbondedForce with appropriate per-particle parameters and global parameter definitions.

  • topology (openmm.app.topology.Topology) – The OpenMM Topology object. Must contain water molecules with the specified residue and atom names.

  • 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 append mode.

  • integrator_str (str) – “BAOABIntegrator” or “LangevinMiddleIntegrator”

  • 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’.

  • create_simulation (bool) – Whether to create a system inside this class. When you only want to customize the system using this class, you can set this to False, in order to avoid unnecessary memory usage.

  • optimization (typing.Literal['O3', 'O1']) – The optimization level for this system. Options are ‘O1’ and ‘O3’. Default is ‘O3’. In ‘O3’, we will try to hardcode water-water vdw if only Oxygen has vdw parameters.

  • n_split_water (typing.Union[int, str]) – Number of split water for Hybrid_REST2 system. If ‘log’, will use log10(N_water).

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_dict

    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_dict

    2. Charge = Proper_water_charge in self.nonbonded_force

  • Switching water:

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

    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

customise_force_hybridREST2(system, optimization, n_split='log')

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

core

new

old

envh

envc

wat

swit

0 core

C_alc

1 new

C_alc

C_alc

2 old

C_alc

None

C_alc

3 envh

C_alc

C_alc

C_alc

NonB

4 envc

C_alc

C_alc

C_alc

NonB

NonB

5 wat

C_alc

C_alc

C_alc

C_wat1

C_wat2

C_wat3

6 swit

C_alc

C_alc

C_alc

C_alc

C_alc

C_alc

C_alc
  • C_alc

    CustomNonbondedForce for Alchemical atoms

  • C_wat1

    CustomNonbondedForce for

  • C_wat2

    CustomNonbondedForce for

  • C_wat3

    CustomNonbondedForce for

The most basic energy expression for vdw is the following. All the other forces are simplified from this one.:

 1energy  = (
 2    "U_rest2;"
 3
 4    # 9. REST2
 5    "U_rest2 = U_sterics * k_rest2_sqrt^is_hot;"
 6    "is_hot = step(3-atom_group1) + step(3-atom_group2);"
 7
 8    # 8. vdw with 1. switching water 2. water real/dummy
 9    "U_sterics = 4*epsilon*x*(x-1.0) * lambda_gc * is_real1 * is_real2;"
10
11    # 7. introduce softcore when new/old/switching atoms are involved
12    "x = 1 / ((softcore_alpha*lambda_alpha + (r/sigma)^6));"
13    "lambda_alpha = (new_X*(1-lambda_sterics_insert) + old_X*lambda_sterics_delete + swit_X*lambda_gc_vdw) / (swit_X + 1);"
14
15    # 6. Interpolating between states A and B
16    "epsilon = (1-lambda_sterics)*epsilonA + lambda_sterics*epsilonB;"
17    "sigma = (1-lambda_sterics)*sigmaA + lambda_sterics*sigmaB;"
18
19    # 5. select lambda
20    "lambda_sterics = (new_X*lambda_sterics_insert + old_X*lambda_sterics_delete + core_cew*lambda_sterics_core);" # only one lambda_XXX will take effect
21    "lambda_gc = ((1 - swit_X) + lambda_gc_vdw * swit_X);" # if swit_X=0 -> lambda_gc=1, if swit_X=1 -> lambda_gc=lambda_gc_vdw
22
23    # 4. determine interaction types
24    "swit_X   = max(is_swit1, is_swit2);"
25    "core_cew = delta(old_X + new_X);"    # core-core + core-envh + core-envc + core-wat + core-swit
26    "new_X    = max(is_new1, is_new2);"
27    "old_X    = max(is_old1, is_old2);"
28
29    # 3. determine atom groups
30    ## 3.1 check atom1
31    "is_core1 = delta(0-atom_group1);"
32    "is_new1  = delta(1-atom_group1);"
33    "is_old1  = delta(2-atom_group1);"
34    "is_envh1 = delta(3-atom_group1);"
35    "is_envc1 = delta(4-atom_group1);"
36    "is_wat1  = delta(5-atom_group1);"
37    "is_swit1 = delta(6-atom_group1);"
38    ## 3.2 check atom2
39    "is_core2 = delta(0-atom_group2);"
40    "is_new2  = delta(1-atom_group2);"
41    "is_old2  = delta(2-atom_group2);"
42    "is_envh2 = delta(3-atom_group2);"
43    "is_envc2 = delta(4-atom_group2);"
44    "is_wat2  = delta(5-atom_group2);"
45    "is_swit2 = delta(6-atom_group2);"
46
47    # 1. LJ mixing rules
48    "epsilonA = sqrt(epsilonA1*epsilonA2);"
49    "epsilonB = sqrt(epsilonB1*epsilonB2);"
50    "sigmaA = 0.5*(sigmaA1 + sigmaA2);"
51    "sigmaB = 0.5*(sigmaB1 + sigmaB2);"
52    )
Parameters:

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

  • optimization (typing.Literal['O3', 'O1']) – Level of optimization. It can be O3 or O1. When O3 is selected, We will try to hard code sigma and epsilon into the energy expression for wat-wat interaction if there is only one water atom has LJ. This can be applied to TIP3P, OPC. O3 will make no difference if charmm TIP3P is used, because H also has LJ parameters.

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

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

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.

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.

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.

logger

Logger for the Sampler

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.

simulation: openmm.app.simulation.Simulation

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

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.