BaseNPTWaterMCSampler

class grandfep.sampler.BaseNPTWaterMCSampler(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', rst_file='md.rst7', dcd_file=None, append=False, set_reporter=True)

Bases: object

Base class for water MC sampling.

This class provides a flexible framework for customizing OpenMM forces so that water molecules can be alchemically inserted (ghost → real) and deleted (real → ghost) at the same time. The last two water molecules in the system will be set as the switching water, which will be used to smoothly turn on/off. 4 global parameters will be added to control the alchemical transformation.

• lambda_vdw_swit6: vdw lambda for the 2nd water molecule (switching water 1)

• lambda_coulomb_swit6: Coulomb lambda for the 2nd water molecule (switching water 1)

• lambda_vdw_swit7: vdw lambda for the last water molecule (switching water 2)

• lambda_coulomb_swit7: Coulomb lambda for the last water molecule (switching water 2)

__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', rst_file='md.rst7', dcd_file=None, append=False, set_reporter=True)

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

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

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:

energy = (
    "U_rest2;"

    # REST2 scaling
    "U_rest2 = U_sterics * k_rest2_sqrt^is_hot;"
    "is_hot = step(3-atom_group1) + step(3-atom_group2);"

    # LJ with softcore + per-particle coupling factors for swit6/swit7
    "U_sterics = 4*epsilon*x*(x-1.0) * coupling1 * coupling2;"

    # ---------------------------
    # 6. softcore activation
    # ---------------------------
    # lambda_alpha is turned on if at least one particle is "possibly dummy"
    #
    "x = 1 / (softcore_alpha*lambda_alpha + (r/sigma)^6);"
    "lambda_alpha = (soft_need1 + soft_need2) / soft_num;"
    "soft_num = max(1, soft_num);" # avoid division by zero
    "soft_num = is_new1 + is_old1 + is_swit6_1 + is_swit7_1 + is_new2 + is_old2 + is_swit6_2 + is_swit7_2;"

    "soft_need1 = is_new1*(1-lambda_sterics_insert) + is_old1*lambda_sterics_delete"
                 " + is_swit6_1*(1-lambda_vdw_swit6) + is_swit7_1*(1-lambda_vdw_swit7);"
    "soft_need2 = is_new2*(1-lambda_sterics_insert) + is_old2*lambda_sterics_delete"
                 " + is_swit6_2*(1-lambda_vdw_swit6) + is_swit7_2*(1-lambda_vdw_swit7);"

    # ---------------------------
    # 5. A/B interpolation
    # ---------------------------
    "epsilon = (1-lambda_sterics)*epsilonA + lambda_sterics*epsilonB;"
    "sigma   = (1-lambda_sterics)*sigmaA   + lambda_sterics*sigmaB;"
    "lambda_sterics = (new_X*lambda_sterics_insert + old_X*lambda_sterics_delete + core_cew*lambda_sterics_core);"

    # ---------------------------
    # 4 coupling for switchable water
    # ---------------------------
    # coupling = 1 for normal particles
    #         = lambda_vdw_swit6 for atom_group==6
    #         = lambda_vdw_swit7 for atom_group==7
    #
    "coupling1 = 1 + is_swit6_1*(lambda_vdw_swit6-1) + is_swit7_1*(lambda_vdw_swit7-1);"
    "coupling2 = 1 + is_swit6_2*(lambda_vdw_swit6-1) + is_swit7_2*(lambda_vdw_swit7-1);"

    # ---------------------------
    # 3. determine interaction types
    # ---------------------------
    "core_cew = delta(old_X + new_X);" # core-core + core-envh + core-envc + core-swit6 + core-swit7
    "new_X    = max(is_new1, is_new2);"
    "old_X    = max(is_old1, is_old2);"


    # ---------------------------
    # 2. determine atom groups
    # ---------------------------
    "is_core1 = delta(0-atom_group1);"
    "is_new1  = delta(1-atom_group1);"
    "is_old1  = delta(2-atom_group1);"
    "is_envh1 = delta(3-atom_group1);"
    "is_envc1 = delta(4-atom_group1);"
    "is_swit6_1 = delta(6-atom_group1);"
    "is_swit7_1 = delta(7-atom_group1);"

    "is_core2 = delta(0-atom_group2);"
    "is_new2  = delta(1-atom_group2);"
    "is_old2  = delta(2-atom_group2);"
    "is_envh2 = delta(3-atom_group2);"
    "is_envc2 = delta(4-atom_group2);"
    "is_swit6_2 = delta(6-atom_group2);"
    "is_swit7_2 = delta(7-atom_group2);"

    # ---------------------------
    # 1. LJ mixing rules
    # ---------------------------
    "epsilonA = sqrt(epsilonA1*epsilonA2);"
    "epsilonB = sqrt(epsilonB1*epsilonB2);"
    "sigmaA   = 0.5*(sigmaA1 + sigmaA2);"
    "sigmaB   = 0.5*(sigmaB1 + sigmaB2);"
)

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

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

dcd_reporter_dict

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

kBT

k_B* 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