utils

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.

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

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)