NoneqNPTWaterMCSamplerMPI

class grandfep.sampler.NoneqNPTWaterMCSamplerMPI(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', position=None, rst_file='md.rst7', dcd_file=None, append=False, active_site=None, init_lambda_state=None, lambda_dict=None)

Bases: _ReplicaExchangeMixin, NoneqNPTWaterMCSampler

NPT waterMC Sampler class with MPI (replica exchange) support. Only Hamiltonian is allowed to be different.

__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', position=None, rst_file='md.rst7', dcd_file=None, append=False, active_site=None, init_lambda_state=None, lambda_dict=None)

Initialize the NPT waterMC sampler with MPI support

Parameters:

  • system (openmm.System) – OpenMM system object, this system should include a barostat

  • topology (app.Topology) – OpenMM topology object

  • temperature (unit.Quantity) – Reference temperature of the system, with unit

  • collision_rate (unit.Quantity) – Collision rate of the system, with unit. e.g., 1 / (1.0 * unit.picoseconds)

  • timestep (unit.Quantity) – Timestep of the simulation, with unit. e.g., 4 * unit.femtoseconds with Hydrogen Mass Repartitioning

  • log (Union[str, Path]) – Log file path for the simulation

  • platform (openmm.Platform) – OpenMM platform to use for the simulation. Default is ‘CUDA’.

  • water_resname (str) – The residue name of water in the system. Default is “HOH”.

  • water_O_name (str) – The atom name of oxygen in the water residue. Default is “O”.

  • position (openmm.unit.quantity.Quantity) – Initial position of the system.

  • rst_file (str) – Restart file path for the simulation. Default is “md.rst7”.

  • dcd_file (str) – DCD file path for the simulation. Default is None, which means no dcd output.

  • append (bool) – If True, append to the existing dcd file. Default is False, which means overwrite the existing dcd file.

  • active_site (dict) – A dictionary defining the active site.

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

customise_force_hybridREST2(system)

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

Groups

core

new

old

envh

envc

swit6

swit7

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

6 swit6

C_alc

C_alc

C_alc

C_alc

C_alc

C_alc

7 swit7

C_alc

C_alc

C_alc

C_alc

C_alc

C_alc

C_alc
  • C_alc

    CustomNonbondedForce for Alchemical atoms

  • NonB

    NonbondedForce for regular atoms

The energy expression for vdw is the following:

 1energy = (
 2    "U_rest2;"
 3
 4    # REST2 scaling
 5    "U_rest2 = U_sterics * k_rest2_sqrt^is_hot;"
 6    "is_hot = step(3-atom_group1) + step(3-atom_group2);"
 7
 8    # LJ with softcore + per-particle coupling factors for swit6/swit7
 9    "U_sterics = 4*epsilon*x*(x-1.0) * coupling1 * coupling2;"
10
11    # ---------------------------
12    # 6. softcore activation
13    # ---------------------------
14    # lambda_alpha is turned on if at least one particle is "possibly dummy"
15    #
16    "x = 1 / (softcore_alpha*lambda_alpha + (r/sigma)^6);"
17    "lambda_alpha = (soft_need1 + soft_need2) / soft_num;"
18    "soft_num = max(1, soft_num);" # avoid division by zero
19    "soft_num = is_new1 + is_old1 + is_swit6_1 + is_swit7_1 + is_new2 + is_old2 + is_swit6_2 + is_swit7_2;"
20
21    "soft_need1 = is_new1*(1-lambda_sterics_insert) + is_old1*lambda_sterics_delete"
22                 " + is_swit6_1*(1-lambda_vdw_swit6) + is_swit7_1*(1-lambda_vdw_swit7);"
23    "soft_need2 = is_new2*(1-lambda_sterics_insert) + is_old2*lambda_sterics_delete"
24                 " + is_swit6_2*(1-lambda_vdw_swit6) + is_swit7_2*(1-lambda_vdw_swit7);"
25
26    # ---------------------------
27    # 5. A/B interpolation
28    # ---------------------------
29    "epsilon = (1-lambda_sterics)*epsilonA + lambda_sterics*epsilonB;"
30    "sigma   = (1-lambda_sterics)*sigmaA   + lambda_sterics*sigmaB;"
31    "lambda_sterics = (new_X*lambda_sterics_insert + old_X*lambda_sterics_delete + core_cew*lambda_sterics_core);"
32
33    # ---------------------------
34    # 4 coupling for switchable water
35    # ---------------------------
36    # coupling = 1 for normal particles
37    #         = lambda_vdw_swit6 for atom_group==6
38    #         = lambda_vdw_swit7 for atom_group==7
39    #
40    "coupling1 = 1 + is_swit6_1*(lambda_vdw_swit6-1) + is_swit7_1*(lambda_vdw_swit7-1);"
41    "coupling2 = 1 + is_swit6_2*(lambda_vdw_swit6-1) + is_swit7_2*(lambda_vdw_swit7-1);"
42
43    # ---------------------------
44    # 3. determine interaction types
45    # ---------------------------
46    "core_cew = delta(old_X + new_X);" # core-core + core-envh + core-envc + core-swit6 + core-swit7
47    "new_X    = max(is_new1, is_new2);"
48    "old_X    = max(is_old1, is_old2);"
49
50
51    # ---------------------------
52    # 2. determine atom groups
53    # ---------------------------
54    "is_core1 = delta(0-atom_group1);"
55    "is_new1  = delta(1-atom_group1);"
56    "is_old1  = delta(2-atom_group1);"
57    "is_envh1 = delta(3-atom_group1);"
58    "is_envc1 = delta(4-atom_group1);"
59    "is_swit6_1 = delta(6-atom_group1);"
60    "is_swit7_1 = delta(7-atom_group1);"
61
62    "is_core2 = delta(0-atom_group2);"
63    "is_new2  = delta(1-atom_group2);"
64    "is_old2  = delta(2-atom_group2);"
65    "is_envh2 = delta(3-atom_group2);"
66    "is_envc2 = delta(4-atom_group2);"
67    "is_swit6_2 = delta(6-atom_group2);"
68    "is_swit7_2 = delta(7-atom_group2);"
69
70    # ---------------------------
71    # 1. LJ mixing rules
72    # ---------------------------
73    "epsilonA = sqrt(epsilonA1*epsilonA2);"
74    "epsilonB = sqrt(epsilonB1*epsilonB2);"
75    "sigmaA   = 0.5*(sigmaA1 + sigmaA2);"
76    "sigmaB   = 0.5*(sigmaB1 + sigmaB2);"
77)
Parameters:

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

Return type:

None

load_rst(rst_input)

Load positions/boxVector/velocities from a restart file

Parameters:

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

Return type:

None

move_in(l_vdw_list, l_chg_list, n_prop, write_pdb=None)

Perform a MC move to insert a switching7 water into the active site, and delete a switching6 water from the bulk.

Parameters:

  • l_vdw_list (list) – A list of value for lambda_gc_vdw, 0 to 1.

  • l_chg_list (list) – A list of value for lambda_gc_coulomb, 0 to 1.

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

  • write_pdb (str) – Folder to write pdb for debugging. Default is None.

Returns:

  • accept (bool) – Whether the MC move is accepted.

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

  • 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 each element

move_out(l_vdw_list, l_chg_list, n_prop, write_pdb=None)

Perform a MC move to insert a switching7 water into the bulk, and delete a switching6 water from the active site.

Parameters:

  • l_vdw_list (list) – A list of value for lambda_gc_vdw, 0 to 1.

  • l_chg_list (list) – A list of value for lambda_gc_coulomb, 0 to 1.

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

  • write_pdb (str) – Folder to write pdb for debugging. Default is None.

Returns:

  • accept (bool) – Whether the MC move is accepted.

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

  • 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 each element

random_pos_vel_water()

Generate random position and velocity for a water molecule from the cached water molecules.

rank_0_print_log(msg)

Print log message only on rank 0, and log on all ranks.

Parameters:

msg – The message to be printed.

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.

  • exchange (bool) – Whether this replica has exchanged its state with another replica.

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_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_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_mc_count(in_out, work, accept, acc_prob)

Update the MC move count

Parameters:

  • in_out (int) – 0 for in_move, 1 for out_move

  • work (unit.Quantity) – The work done during the MC move, with unit

  • accept (bool) – Whether the MC move is accepted

  • acc_prob (float) – The acceptance probability of the MC move

Return type:

None

work_process(n_prop, l_vdw_list, l_chg_list)

Perform the work process for delete swtich6 and insert switch7.

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.

comm

MPI communicator

dcd_reporter_dict

A dictionary of all the dcd reporter. Call the reporter inside to write the dcd trajectory file.

kBT

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.

rank

Rank of this process

re_step

Number of replica exchanges performed

rst_reporter_dict

A dictionary of all the rst7 reporter. Call the reporter inside to write the rst7 restart file.

size

Size of the communicator (number of processes)

temperature: openmm.unit.quantity.Quantity

reference temperature of the system, with unit