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

maxh¶

The maximum run time. Unit in hour

Type:: float

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.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, water_resname='HOH')¶

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.
water_resname (str) – The residue name of water molecules to be excluded from the restraints. Default is “HOH”.

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

Navigation

Related Topics

utils¶

Example:¶

Parameter¶