utils

grandfep.utils

class grandfep.utils.ActiveSite(center_index)

Bases: object

Class to define an active site in the system. Only works for cubic box now.

get_atom_states(atom_indices, box_vectors, positions)

random_position(box_vectors, positions, in_active_site=True)

rap_box(box_vectors, positions)

class grandfep.utils.ActiveSiteCube(center_index, box_abc)

Bases: ActiveSite

Class to define a spherical active site in the system. only works for cubic box now.

get_atom_states(atom_indices, box_vectors, positions)

class grandfep.utils.ActiveSiteSphere(center_index, radius)

Bases: ActiveSite

Class to define a spherical active site in the system.

get_atom_states(atom_indices, box_vectors, positions)

get_volume(state)

random_position(box_vectors, positions, in_active_site=True)

class grandfep.utils.ActiveSiteSphereRelative(center_index, radius, box_vectors)

Bases: ActiveSite

Class to define a spherical active site in the system.

get_atom_states(atom_indices, box_vectors, positions)

get_volume(state)

random_position(box_vectors, positions, in_active_site=True)

class grandfep.utils.FreeEAnalysis(file_list, keyword, separator, drop_equil=True, begin=0)

Bases: object

Class to analyze free energy calculations using BAR/MBAR.

static print_res_all(res_all)

static read_energy(log_file, keyword='Reduced Energy U_i(x):', separator=',', begin=0)

Parse the energy from a file

Parameters:

• log_file (typing.Union[str, pathlib.Path]) – File to be parsed

• keyword (str) – The keyword to locate the line. The energy should immediately follow this keyword.

• separator (str) – The separator to use between energy in the file

• begin (int) – The first frame to be analyzed. The first frame is 0.

Return type:

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

Returns:

• e_array (np.array) – The energy array.

• temperature (unit.Quantity) – The temperature in Kelvin. If not found, None is returned.

bar_U_all()

mbar_U_all()

print_uncorrelate()

Print the number of uncorrelated sample in each file

set_temperature(temperature)

Set the temperature for the analysis.

Parameters:

temperature (unit.Quantity) – Temperature in Kelvin, with unit

sub_sample(drop_equil=True)

class grandfep.utils.dcd_reporter(file, reportInterval, append=False, enforcePeriodicBox=None)

Bases: DCDReporter

A custom DCD reporter that can report from coordinate np.array.

Parameters:

• file (string) – The file to write to

• reportInterval (int) – The interval (in time steps) at which to write frames

• append (bool=False) – If True, open an existing DCD file to append to. If False, create a new file.

• enforcePeriodicBox (bool) – Specifies whether particle positions should be translated so the center of every molecule lies in the same periodic box. If None (the default), it will automatically decide whether to translate molecules based on whether the system being simulated uses periodic boundary conditions.

report_positions(simulation, periodicBoxVectors, positions_nm)

Write one frame to the DCD file.

Parameters:

• simulation (openmm.app.simulation.Simulation) – The Simulation to generate a report for

• periodicBoxVectors (typing.Union[openmm.vec3.Vec3, numpy.array]) – You can prepare this np.array with state.getPeriodicBoxVectors()

• positions_nm (typing.Union[openmm.unit.quantity.Quantity, numpy.ndarray]) – You can prepare this np.array with state.getPositions(asNumpy=True).value_in_unit(unit.nanometers). It should either be unit.Quantity or value in nm.

class grandfep.utils.md_params_yml(yml_file=None)

Bases: object

Class to manage MD parameters with default values and YAML overrides. This class reads whatever in the yml file, and it does not guarantee the correct usage of the parameters anywhere else. The built-in parameters in self._unit_map will be automatically added with units. If you want a unit to be attached, you can add it to self._unit_map before loading the yml file.

Example:

>>> from grandfep import utils
>>> mdp = utils.md_params_yml("test/Water_Chemical_Potential/OPC/multidir/0/md.yml")

integrator

Name of the integrator. Can be “BAOABIntegrator”, “LangevinMiddleIntegrator”, etc.

Type:

str

dt

Time step. Unit in ps

Type:

unit.Quantity

nsteps

Number of steps.

Type:

int

nst_dcd

Number of steps per dcd trajectory output.

Type:

int

nst_csv

Number of steps per csv energy file output.

Type:

int

ncycle_dcd

Number of cycles per dcd trajectory output.

Type:

int

ncycle_csv

Number of cycles per csv energy file output.

Type:

int

tau_t

Temperature coupling time constant. Unit in ps

Type:

unit.Quantity

ref_t

Reference temperature. Unit in K

Type:

unit.Quantity

gen_vel

Generate velocities.

Type:

bool

gen_temp

Temperature for random velocity generation. Unit in K

Type:

unit.Quantity

restraint

Whether to apply restraints.

Type:

bool

restraint_fc

Restraint force constant. Unit in kJ/mol/nm^2

Type:

unit.Quantity

pcoupltype

Specifies the kind of pressure coupling used. MonteCarloBarostat, MonteCarloMembraneBarostat

Type:

str

ref_p

Reference pressure. Unit in bar

Type:

unit.Quantity

nstpcouple

Pressure coupling frequency.

Type:

int

surface_tension

Surface tension. Unit in bar*nm

Type:

unit.Quantity

ex_potential

Excess potential in GC. Unit in kcal/mol

Type:

unit.Quantity

standard_volume

Standard volume in GC. Unit in nm^3

Type:

unit.Quantity

n_propagation

Number of propagation steps bewteen each lambda switching.

Type:

int

init_lambda_state

The lambda state index to simulate

Type:

int

calc_neighbor_only

Whether to calculate the energy of the nearest neighbor only when performing replica exchange.

Type:

bool

md_gc_re_protocol

MD, Grand Canonical, Replica Exchange protocol. The default is [("MD", 200),("GC", 1),("MD", 200),("RE", 1),("MD", 200),("RE", 1),("MD", 200),("RE", 1)]

Type:

list

system_setting

Nonbonded parameters for the OpenMM createSystem. The default is {"nonbondedMethod": app.PME, "nonbondedCutoff": 1.0 * unit.nanometer, "constraints": app.HBonds}. openmm has to be imported in the way the the value string can be evaluated.

>>> from openmm import app, openmm, unit

Type:

dict

sphere_radius

The radius of the GCMC sphere. Unit in nm

Type:

unit.Quantity

lambda_gc_vdw

This is the ghobal parameter for controlling the vdw on the switching water. You can give a list of lambda values for the path of GC insertion.

Type:

list

lambda_gc_coulomb

This is the ghobal parameter for controlling the Coulomb on the switching water. You can give a list of lambda values for the path of GC insertion.

Type:

list

lambda_angles

Lambda

Type:

list

lambda_bonds

Lambda

Type:

list

lambda_sterics_core

Lambda

Type:

list

lambda_electrostatics_core

Lambda

Type:

list

lambda_sterics_delete

Lambda

Type:

list

lambda_electrostatics_delete

Lambda

Type:

list

lambda_sterics_insert

Lambda

Type:

list

lambda_electrostatics_insert

Lambda

Type:

list

lambda_torsions

Lambda

Type:

list

get_lambda_dict()

Returns:

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

Return type:

lambda_dict

get_system_setting()

Evaluate system_setting, and return them in a dictionary

class grandfep.utils.rst7_reporter(file, reportInterval, write_multiple=False, netcdf=False, write_velocities=True)

Bases: RestartReporter

A reporter to handle writing Amber rst7 restart files.

Parameters:

• file (str) – Name of the file to write the restart to.

• reportInterval (int) – Number of steps between writing restart files

• write_multiple (bool=False) – Either write a separate restart each time (appending the step number in the format .# to the file name given above) if True, or overwrite the same file each time if False

• netcdf (bool=False) – Use the Amber NetCDF restart file format

• write_velocities (bool=True) – Write velocities to the restart file. You can turn this off for passing in, for instance, a minimized structure.

report_positions_velocities(sim, state, periodicBoxVectors, positions_A, velocities)

Generate a report.

Parameters:

• sim (openmm.app.simulation.Simulation) – The Simulation to generate a report for

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

• periodicBoxVectors (typing.Union[openmm.vec3.Vec3, numpy.array]) – You can prepare this np.array with state.getPeriodicBoxVectors()

• positions_A (typing.Union[openmm.unit.quantity.Quantity, numpy.ndarray]) – You can prepare this np.array with state.getPositions(asNumpy=True).value_in_unit(unit.angstrom). It should either be unit.Quantity or value in Angstrom.

• velocities (typing.Union[openmm.unit.quantity.Quantity, numpy.ndarray]) – You can prepare this np.array with state.getVelocities(asNumpy=True).value_in_unit(unit.angstrom/unit.picosecond). It should either be unit.Quantity or value in Angstrom/ps.

grandfep.utils.check_system_type(system, no_barostat=True)

Check if system can be converted to a GCMC system.

Parameters:

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

Return type:

str

Returns:

(str) The type of the system. Can be Amber, Charmm or Hybrid Amber should have NonbondedForce without CustomNonbondedForce Charmm should have NonbondedForce and CustomNonbondedForce. The EnergyFunction that is allowed in CustomNonbondedForce is (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;

Hybrid should have

CustomBondForce For perturbed bonds HarmonicBondForce CustomAngleForce For perturbed angles HarmonicAngleForce CustomTorsionForce For perturbed dihedrals PeriodicTorsionForce NonbondedForce For perturbed Coulomb CustomNonbondedForce For perturbed Lennard-Jones CustomBondForce_exceptions For perturbed exceptions, mostly 1-4 nonbonded

Hybrid_REST2 should have

CustomBondForce For perturbed bonds HarmonicBondForce CustomAngleForce For perturbed angles HarmonicAngleForce CustomTorsionForce For perturbed dihedrals PeriodicTorsionForce NonbondedForce For perturbed Coulomb CustomNonbondedForce For perturbed Lennard-Jones CustomBondForce_exceptions_1D For perturbed exceptions, mostly 1-4 nonbonded

grandfep.utils.check_water_points(topology, resname)

Check if the water model is 3-point or 4-point in topology. :return: int

grandfep.utils.copy_exclusion_c2c(c1, c2)

Copy exclusion from one CustomNonbondedForce to another.

Parameters:

• c1 (openmm.openmm.CustomNonbondedForce) – Source CustomNonbondedForce from which to copy exclusion.

• c2 (openmm.openmm.CustomNonbondedForce) – Destination CustomNonbondedForce to update in-place.

Returns:

The function updates c2 in-place.

Return type:

None

grandfep.utils.copy_nonbonded_setting_c2c(c1, c2)

Copy nonbonded settings from one CustomNonbondedForce to another.

Parameters:

• c1 (openmm.openmm.CustomNonbondedForce) – Source CustomNonbondedForce from which to copy settings. Copied settings include: nonbonded method, cutoff distance, switching function usage, switching distance, and long-range dispersion correction usage.

• c2 (openmm.openmm.CustomNonbondedForce) – Destination CustomNonbondedForce to update in-place.

Returns:

The function updates c2 in-place.

Return type:

None

grandfep.utils.copy_nonbonded_setting_n2c(nonbonded_force, custom_nonbonded_force)

Copy nonbonded settings from a NonbondedForce to a CustomNonbondedForce.

Parameters:

• nonbonded_force (openmm.openmm.NonbondedForce) – Source NonbondedForce from which to copy settings (nonbonded method, cutoff distance, switching function and switching distance, dispersion correction).

• custom_nonbonded_force (openmm.openmm.CustomNonbondedForce) – Destination CustomNonbondedForce to update with the copied settings.

Returns:

The function updates custom_nonbonded_force in-place.

Return type:

None

grandfep.utils.find_all_water(topology, resname, water_O_name)

Identify water residues and their oxygen atom indices in the topology.

Parameters:

• topology (openmm.app.Topology) – OpenMM topology to inspect.

• resname (str) – Residue name used to identify water molecules in the topology.

• water_O_name (str) – Atom name used to identify the oxygen atom within a water residue.

Returns:

• water_res_2_atom (dict) – Mapping from water residue index (int) to a list of atom indices (list of int) that belong to that water residue.

• water_res_2_O (dict) – Mapping from water residue index (int) to the atom index (int) of the oxygen atom within that residue.

Raises:

ValueError – If a water residue matching resname does not contain an atom named water_O_name, or if no residues matching resname are present in the topology.

grandfep.utils.find_mapping(map_list, resA, resB)

Check if residues A and B match with a mapping in the mapping list.

Parameter

map_list :

A list of mapping. For example

[{'res_nameA': 'MOL',
  'res_nameB': 'MOL',
  'index_map': {1: 1, 0: 2, 2: 0}
 }]

resA :

residue from state A

resB :

residue from state B

rtype:

tuple[bool, dict]

returns:

• match_flag – resA and resB can be matched

• index_map – A dictionay mapping the atoms from state A to state B with their 0-index

grandfep.utils.find_reference_atom_indices(topology, ref_atoms_list)

Find atom indices in the topology that match the given reference atom definitions.

Parameters:

• topology (openmm.app.topology.Topology) – OpenMM topology object

• ref_atoms_list (list) –

  A list of dictionaries specifying reference atoms. Each dictionary can contain any combination of the following keys:

  • chain_index: int

    Index of the chain in the topology.

  • res_name: str

    Residue name (e.g., “HOH”).

  • res_id: str

    In openmm topology, res_id is string.

  • res_index: int

    0-based index of the residue in the topology.

  • atom_name: str

    Atom name (e.g., “O”, “H1”).

Returns:

A list of integer atom indices matching the provided atom specifications.

Return type:

list

Examples

>>> from grandfep import utils
>>> top = "test/KcsA_5VKE_SF/step1_pdbreader_12WAT.psf"
>>> topology, _ = utils.load_top(top)
>>> ref_atoms_list = [{"res_id":"71", "res_name":"GLU", "atom_name":"O"}]
>>> utils.find_reference_atom_indices(topology, ref_atoms_list)
[21, 164, 307, 450]

grandfep.utils.get_water_parameters(topology, system, water_resname='HOH')

Get charge, vdw of water parameters from the system.

Parameters:

• topology (openmm.app.topology.Topology) – The OpenMM topology object. Water residues will be identified by water_resname.

• system (openmm.openmm.System) – The OpenMM system object containing force field parameters. The water parameters will be extracted from the NonbondedForce in this system.

• water_resname (str) – The residue name used to identify water molecules in the topology. Default is “HOH”.

Returns:

water_params – A dictionary containing the water parameters with keys: - ‘charge’: charge of all the atoms in a water molecule. - ‘sigma’: sigma of all the atoms in a water molecule. - ‘epsilon’: epsilon of all the atoms in a water molecule.

Return type:

dict

grandfep.utils.load_rst(context, rst_input)

Load coordinates and box vectors from an AMBER restart file (.rst7) into an OpenMM context.

Parameters:

• context (openmm.Context) – The OpenMM context to load the coordinates and box vectors into.

• rst_input (str or Path) – Path to the AMBER restart file (.rst7).

Return type:

None

Examples

from grandfep import utils
from openmm import app, openmm, unit

# Load system and create context
pdb = app.PDBFile("input.pdb")
forcefield = app.ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')
system = forcefield.createSystem(pdb.topology, nonbondedMethod=app.PME,
                                 nonbondedCutoff=1 * unit.nanometer, constraints=app.HBonds)
integrator = openmm.LangevinMiddleIntegrator(300*unit.kelvin, 1/unit.picosecond, 0.004*unit.picoseconds)
simulation = app.Simulation(pdb.topology, system, integrator)
simulation.context.setPositions(pdb.positions)

# Load restart file into context
utils.load_rst(simulation.context, "restart.rst7")

grandfep.utils.load_sys(sys_file)

Load a serialized OpenMM system from a xml or xml.gz file.

Parameters:

sys_file (typing.Union[str, pathlib.Path]) – Path to the serialized OpenMM system file (.xml or .xml.gz).

Returns:

The deserialized OpenMM System object.

Return type:

openmm.System

Examples

from grandfep import utils
system = utils.load_sys("system.xml.gz")

grandfep.utils.load_top(top_file)

Load a topology file in PSF (CHARMM) or PRMTOP/PARM7 (AMBER) format.

After loading a top file, you should call .setPeriodicBoxVectors() on the topology object to define the box. Without setting the box, trajectory files may lack box information.

Parameters:

top_file (typing.Union[str, pathlib.Path]) – Path to the topology file (either .psf, .prmtop, or .parm7).

Return type:

tuple[openmm.app.topology.Topology, typing.Union[openmm.app.charmmpsffile.CharmmPsfFile, openmm.app.amberprmtopfile.AmberPrmtopFile]]

Returns:

• topology (openmm.app.Topology) – The OpenMM Topology object for use in simulation setup.

• top_object (openmm.app.CharmmPsfFile or openmm.app.AmberPrmtopFile) – The loaded OpenMM file object used to construct the topology.

Raises:

ValueError – If the file format is not supported. Only the extensions .psf, .prmtop, and .parm7 are supported.

grandfep.utils.period_from_k_mu(k, mu)

Calculate the period T from spring constant k and reduced mass mu. Please make sure k and mu have the proper units.

Parameters:

• k (spring constant (energy/distance^2), e.g. kJ/mol/nm^2)

• mu (reduced mass (mass), e.g. dalton)

Return type:

T

grandfep.utils.prepare_atom_map(topologyA, topologyB, map_list)

With given topology A and B, prepare the atom mapping.

Parameters:

• topologyA (openmm.app.topology.Topology) – topology for state A

• topologyB (openmm.app.topology.Topology) – topology for state B

• map_list (list) –

  A list of mapping. For example

  [{'res_nameA': 'MOL', 'res_nameB': 'MOL',
    'index_map': {1: 1, 0: 2, 2: 0, 25: 26, 22: 29, 23: 28, 24: 27}
  

Return type:

tuple[dict, dict]

Returns:

• old_to_new_all – A dictionay for all the atoms that should map from old (state A) to new (state B)

• old_to_new_core – A dictionay for Alchemical (core) atoms that should map from old (state A) to new (state B)

grandfep.utils.prepare_restraints_force(topology, positions, fc, solvent_resname=None)

Prepare a force to add position restraints to heavy atoms. 1/2 * k * (x - x0)^2

Parameters:

• topology (openmm.app.Topology) – The OpenMM topology object.

• positions (openmm.unit.Quantity) – The reference positions of the atoms in the system.

• fc (openmm.unit.Quantity) – The force constant for the position restraints. Unit in kJ/mol/nm^2.

• solvent_resname (list) – The residue name of water molecules to be excluded from the restraints. Default is “HOH”, “Na+”, “K+”, “Cl-“.

Returns:

• posres – An OpenMM CustomExternalForce object for position restraints.

• res_atom_count (int) – The number of atoms in the system that are restrained.

• res_residue_list (list) – A list of residue names that are restrained.

grandfep.utils.random_rotation_matrix()

Generate a random rotation matrix using Shoemake method.

Returns:

A 3x3 rotation matrix.

Return type:

np.ndarray

Examples

import numpy as np
import matplotlib.pyplot as plt
from grandfep import utils
def gen_random_vec():
    axis = np.random.normal(0, 1, 3)
    axis /= np.linalg.norm(axis)
    return axis
res_new = []
res_ref = []
vec_init = gen_random_vec() # regardless of the initial vector, the rotated vector should have uniform distribution on x,y,z
for i in range(100000):
    rot_matrix = utils.random_rotation_matrix()
    res_new.append(np.dot(rot_matrix, vec_init))
    res_ref.append(gen_random_vec())
res_new = np.array(res_new)
res_ref = np.array(res_ref)
fig, axes = plt.subplots(2, 3, dpi=300, figsize = (9,6))
for res, ax_list in zip([res_new, res_ref], axes):
    for i, ax in enumerate(ax_list):
        ax.hist(res[:, i], orientation='horizontal', density=True)

grandfep.utils.random_rotation_matrix_protoms()

Copied from https://github.com/essex-lab/grand/blob/master/grand/utils.py. Possibly be wrong.

Returns:

A 3x3 rotation matrix.

Return type:

np.ndarray

grandfep.utils.reduced_mass(m1, m2)

(m1*m2)/(m1+m2). Reduced mass for frequency calculation.

Parameters:

• m1 (mass of particle 1)

• m2 (mass of particle 2)

Return type:

reduced_mass

grandfep.utils.seconds_to_hms(seconds)

Convert seconds to hours, minutes, and seconds.

Parameters:

seconds (float) – Time in seconds.

Return type:

tuple[int, int, float]

Returns:

• hours

• minutes

• seconds

grandfep.utils.HybridTopologyFactoryREST2

class grandfep.utils.HybridTopologyFactoryREST2(old_system, old_positions, old_topology, new_system, new_positions, new_topology, old_to_new_atom_map, old_to_new_core_atom_map, old_rest2_atom_indices=None, use_dispersion_correction=False, softcore_alpha=0.5, scale_dihe=None)

Bases: object

This class can add REST2 when creating a hybrid topology.

Atoms in the resulting hybrid system are treated as being from one of four possible types:

unique_old_atomThese atoms are not mapped and only present in the old

system. Their interactions will be on for lambda=0, off for lambda=1

unique_new_atomThese atoms are not mapped and only present in the new

system. Their interactions will be off for lambda=0, on for lambda=1

core_atomThese atoms are mapped between the two end states, and are

part of a residue that is changing alchemically. Their interactions will be those corresponding to the old system at lambda=0, and those corresponding to the new system at lambda=1

environment_atomThese atoms are mapped between the two end states, and

are not part of a residue undergoing an alchemical change. Their interactions are always on and are alchemically unmodified.

All the alchemically changed atoms (core, unique_old, unique_new) will be considered “hot” in REST2. Other atoms in the environment can also be added by old_rest2_region_atom_indices

Properties

hybrid_systemopenmm.System

The hybrid system for simulation

new_to_hybrid_atom_mapdict of intint

The mapping of new system atoms to hybrid atoms

old_to_hybrid_atom_mapdict of intint

The mapping of old system atoms to hybrid atoms

hybrid_positions[n, 3] np.ndarray

The positions of the hybrid system

hybrid_topologymdtraj.Topology

The topology of the hybrid system

omm_hybrid_topologyopenmm.app.Topology

The OpenMM topology object corresponding to the hybrid system

Here are the forces in the hybrid system:

• CustomBondForce

  handles bonds involving core_atoms (these are interpolated). REST2 does not change bonds

• HarmonicBondForce

  handles bonds involving environment_atoms, unique_old_atoms and unique_new_atoms (these are never scaled)

• CustomAngleForce

  handles angles involving core_atoms (these are interpolated). REST2 does not change angles

• HarmonicAngleForce

  handles angles involving environment_atoms, unique_old_atoms and unique_new_atoms (these are never scaled)

• CustomTorsionForce

  handles torsions involving core_atoms, rest2_atoms, unique_old_atoms and unique_new_atoms (these are interpolated). We are planning to scale d-d-d-d, r-d-d-d and r-r-d-d to 0 in the future

• PeriodicTorsionForce

  handles torsions involving environment_atoms (these are never scaled)

• NonbondedForce

  handles all electrostatic interactions, non-rest2 environment-environment steric interactions

• CustomNonbondedForce

  handle all non environment-environment sterics

• CustomBondForce_exceptions

  handles all electrostatics and sterics exceptions involving unique old/new atoms when interpolate_14s is True, otherwise the electrostatics/sterics exception is in the NonbondedForce

where interactions refers to any pair of atoms that is not 1-2, 1-3, 1-4

vdW table +———-+——–+——–+——–+——–+——–+ | Groups | core | new | old | envh | envc | +==========+========+========+========+========+========+ | 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 | C_alc | | +———-+——–+——–+——–+——–+——–+ | 4 envc | C_alc | C_alc | C_alc | C_alc | NonB | +———-+——–+——–+——–+——–+——–+

__init__(old_system, old_positions, old_topology, new_system, new_positions, new_topology, old_to_new_atom_map, old_to_new_core_atom_map, old_rest2_atom_indices=None, use_dispersion_correction=False, softcore_alpha=0.5, scale_dihe=None)

Parameters:

• old_system – The OpenMM System object of the old state (state A).

• old_positions – The positions of the old state (state A).

• old_topology – The OpenMM Topology object of the old state (state A).

• new_system – The OpenMM System object of the new state (state B).

• new_positions – The positions of the new state (state B).

• new_topology – The OpenMM Topology object of the new state (state B).

• old_to_new_atom_map – All atoms that should map from the old state (A) to new state (B).

• old_to_new_core_atom_map – Alchemical Atoms that should map from A to B

• old_rest2_atom_indices – The indices of the atoms in the old system that should be considered as hot in REST2.

• use_dispersion_correction (bool, default False) – Whether to use the long range correction in the custom sterics force. This can be very expensive for derivative.

• softcore_alpha (float, default None) – “alpha” parameter of softcore sterics, default 0.5.

• scale_dihe (dict) – A dictionary defining the dummy atom scaling factors for dihedrals with different periodicities. Keys are periodicities (1, 2, 3, etc.), and values are the corresponding scaling factors. default is scale_dihe={ 1: 0.1, “i1”: 0.1, 2: 1.0, “i2”: 1.0, 3: 0.1, “i3”: 0.1, 4: 0.1, “i4”: 0.1, 5: 0.1, “i5”: 0.1,}.

_handle_nonbonded()

Handle the nonbonded interactions defined in the new and old systems. k_rest2 = T_cold / T_hot

Coulomb

1. Old atoms (atoms that only show in state A).

   chg = chgA * (1-λ) * k_rest2_sqrt One offsets are needed:

2. New atoms (atoms that only show in state B).

   chg = chgB * λ * k_rest2_sqrt One offset is needed:

3. Core atoms (atoms that are present in both states A and B).

   chg = (chgA * (1-λ) + chgB * λ) * k_rest2_sqrt = chgA * (1-λ) * k_rest2_sqrt + chgB * λ * k_rest2_sqrt Two offsets are needed:

4. EnvH atoms (atoms that are not alchemically modified, but included in hot region).

   chg = chg * k_rest2_sqrt One offset is needed:

5. remove all 1-4 involving hot atoms

X	param	math	chargeScale
Old	lam_ele_del_x_k_rest2_sqrt	(1-λ) * k_rest2_sqrt	chgA
New	lam_ele_ins_x_k_rest2_sqrt	λ * k_rest2_sqrt	chgB
Core	lam_ele_coreA_x_k_rest2_sqrt	(1-λ) * k_rest2_sqrt	chgA
Core	lam_ele_coreB_x_k_rest2_sqrt	λ * k_rest2_sqrt	chgB
EnvH	k_rest2_sqrt	k_rest2_sqrt	chg

vdW

• NonbondedForce

  env-env vdW

• CustomNonbondedForce

  non env-env vdW. U_rest2 = U_sterics*k_rest2_sqrt^(is_hot1+is_hot2)