NoneqNPTWaterMCSampler

class grandfep.sampler.NoneqNPTWaterMCSampler(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)

Bases: BaseNPTWaterMCSampler

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

Initialize the NPT sampler

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

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

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

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.

report_dcd(state=None)

Append a frame to the DCD trajectory file.

Returns:

None

report_rst(state=None)

Write an Amber rst7 restart file.

Returns:

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.

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.

logger

Logger for the Sampler

rst_reporter_dict

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

temperature: openmm.unit.quantity.Quantity

reference temperature of the system, with unit