gromacs_py package

Subpackages

Submodules

gromacs_py.free_ener module

class gromacs_py.free_ener.FreeEner(mol_name, out_folder, unit='kcal')

Bases: object

Free Energy encapsulation class.

This class can be used to launch and analyze free energy calculations using Free Energy Perturbation.

Parameters

  • out_folder (str) – output folder

  • lambda_coul (list of float) – lambda points for Coulomb

  • lambda_vdw (list of float) – lambda points for Lennard Jones

  • lambda_restr (list of float) – lambda points for restraints

  • xvg_file_list (list of string) – List of free energy xvg files

  • lambda_sys_list (list of string) – List of lambda GmxSys

  • temp (float) – Temperature

  • smile (str) – Ligand SMILE

add_intermol_restr_index(rec_index_list, lig_index_list, ref_coor, k=41.84, temp=300)

Compute and add the intermolecular restraints based on the ref_coor distance and angles.

Give three atoms for each receptor and ligand index list: \(R_0\), \(R_1\), \(R_2\) and \(L_0\), \(L_1\), \(L_2\) Will define:

  • 1 bond:

    • \(R_0 - L_0\)

  • 2 angles:

    • \(R_0 - L_0 - L_1\)

    • \(R_1 - R_0 - L_0\)

  • 2 dihedral angles:

    • \(R_0 - L_0 - L_1 - L_2\)

    • \(R_2 - R_1 - R_0 - L_0\)

Parameters

  • rec_index_list (list) – List of the receptor atom index

  • lig_index_list (list) – List of the ligand atom index

  • ref_coor (str) – Reference coordinates file

  • k (float) – intermolecular force constant, (default= \(41.84 \; kcal \, mol^{-1} \, nm^{-2}\))

  • temp (float) – Temperature defalult=300 K

Object requirement(s):

  • self.gmxsys.coor_file

  • self.gmxsys.top_file

Object field(s) changed:

  • self.gmxsys.top_file

align_ref_traj(rec_group='Protein')
compute_add_intermol_from_traj(ref_coor=None, rec_group='Protein', k=41.84, cutoff_prot=6.0)

Compute intermolecular restraint from the last GmxSys trajectory. Get a coor object Get distance and angles

compute_convergence_alchemlyb(dt=10)
compute_convergence_gbar(dt=10)
static compute_lambda_point(gmx_sys, lambda_iter, mol_name, out_folder, free_ener_option, pbar, mbar, em_steps, nvt_steps, npt_steps, prod_steps, maxwarn=1, monitor_tool={'file_check_ext': 'log', 'function': <function progress_bar>})

Run the different steps of a single lambda point.

Parameters

  • gmx_sys (GmxSys) – Gmx System to start from

  • lambda_iter (int) – lambda point

  • mol_name (str) – Molecule residue name

  • out_folder (str) – path of the output folder

  • free_ener_option (dict) – Mdp options

  • pbar (tqmd bar) – progress bar object

  • mbar (bool) – MBAR flag

  • em_steps (int) – number of minimisation steps

  • nvt_steps (int) – number of NVT equilibration steps

  • npt_steps (int) – number of NPT equilibration steps

  • prod_steps (int) – number of production steps

  • maxwarn (int, default=0) – Maximum number of warnings when using gmx grompp

  • monitor – option to monitor a simulation, if not none monitor should contains two values: function the function to be ran while simulation is running and input parameters for the function

property conv_fac

Conversion factor from kT to self.unit

equilibrate_complex(em_steps=10000, HA_time=0.25, CA_time=0.5, CA_LOW_time=1.0, dt=0.002, dt_HA=0.001, temp=300, receptor_grp='Protein', short_steps=50000)

Equilibrate a receptor-ligand complex system.

Parameters

  • em_steps (int) – Energy minimisation step number, default=10000

  • HA_time (float, default=0.25) – Heavy atoms restraints equilibration time (ns)

  • CA_time (float, default=0.5) – Alpha carbon atoms restraints equilibration time (ns)

  • CA_LOW_time (float, default=1.0) – Low alpha carbon atoms restraints equilibration time (ns)

  • dt (float, default=0.002 ps) – Integration time step (ps)

  • dt_HA (float, default=0.001 ps) – Integration time step (ps)

  • temp (float, default=300.0 K) – Temperature of equilibration (K)

  • receptor_grp (str, default='Protein') – Receptor group (for temperature coupling)

  • short_steps (int, default=50000) – Short equilibration steps

equilibrate_solvent_box(em_steps=10000, dt=0.002, prod_time=10.0, short_steps=50000, temp=300.0)

Equilibrate a solvent (water, octanol) box.

Parameters

  • em_steps (int, default=10000) – Energy minimisation step number

  • dt (float, default=0.002 ps) – Integration time step (ps)

  • prod_time (float, default=10.0 ns) – Equilibration time (ns)

  • short_steps (int, default=50000) – Short equilibration steps

  • temp (float, default=300.0) – Temperature of equilibration

extend_lambda_prod(prod_time)

Extend free energy production.

Parameters

prod_time (float) – production time (ns)

static get_bar(xvg_file_list, bar_xvg='bar.xvg', barint_xvg='barint.xvg', hist_xvg='histogram.xvg', begin_time=0, end_time=- 1, check_file_out=True, keep_ener_file=False)

Get energy of a system using gmx bar.

I don’t know how to compute std like in gmx bar.

static get_conv_fac(unit, temp)

Conversion factor from kT to self.unit

get_free_ener(begin_time=0, end_time=- 1, unit=None)

Show free energy calculation output

NEED TO FIX STD !!

get_ligand_atoms(ref_coor)
get_protein_atoms_from_res(resid, rec_group='Protein')
get_protein_atoms_from_rmsf(ref_coor, lig_atom_list, rec_group='Protein', cutoff_max=6.0)
get_water_restr(temp=300)

Compute ligand restaint energy in water using Boresh et al. equation:

\(\Delta G_{restr} = RT \ln \left( \dfrac{8 \pi^2 V^0} {r^2_0 \sin{\theta_{a 0}} \sin{\theta_{b 0}}} \dfrac{ \sqrt{k_r k_{\theta_a} k_{\theta_b} k_{\tau_\alpha} k_{\tau_\beta} k_{\tau_\gamma}}} {2 \pi KT^3} \right)\)

octanol_box_from_SMILE(smile, method_3d='rdkit', iter_num=5000, box_dist=1.3)

Create an octanol box coordinates and toplogie with a molecule defined by its SMILE.

Parameters

  • smile (str) – Molecule’s SMILE

  • method_3d (str, default='rdkit') – Method to compute 3D coordinates

  • iter_num (int, default=5000) – Iteration step number for 3D coordinate computation

  • box_dist (float, default=1.3 nm) – Ditance to egde box (nm)

Note

Default box dist 1.3 nm was taken as minimal distance for CH4 molecule, to allow domain decomposition with gmx mdrun.

plot_convergence(graph_out=None, dt=10)
plot_convergence_graph(graph_out=None)
plot_intermol_restr(graph_out=None)
prepare_complex_pdb(pdb_in, smile, ff='amber99sb-ildn')

Prepare topologie from a pdb file, create box around and solvate it.

Parameters

  • pdb_in (str) – Input PDB file

  • smile (str) – Ligand’s SMILE

  • ff (str, default='amber99sb-ildn') – Forcefield

run(lambda_coul_list, lambda_vdw_list, lambda_restr_list=[], mbar=False, dir_name='free_ener_run', em_steps=5000, nvt_time=10, npt_time=10, prod_time=100, dt=0.002, temp=300.0, temp_groups='Protein non-Protein', maxwarn=1, monitor_tool={'file_check_ext': 'log', 'function': <function progress_bar>})

Compute free energy to uncouple a molecule to a system.

Parameters

  • dir_name (str, default='free_ener_run') – path of the output folder

  • mol_name (str) – Name of the molecule

  • lambda_coul_list (list) – List lambda points for Coulomb

  • lambda_vdw_list (list) – List lambda points for Lennard Jones

  • lambda_bond_list (list, default=[]) – List lambda points for restraints

  • mbar (bool, default=False) – MBAR flag

  • em_steps (int, default=5000) – number of minimisation steps

  • nvt_time (int, default=10 ps) – Time (ps) of NVT equilibration

  • npt_time (int, default=10 ps) – Time (ps) of NPT equilibration

  • prod_time (float, default=100 ps) – Time (ps) of production run

  • dt (float, default=0.002) – integration time step

  • name (str, default=None) – name of the simulation to run

  • temp (float, default=300.0) – Temperature K

  • temp_groups (str, default='Protein non-Protein') – Group(s) for temperature coupling

  • maxwarn (int, default=0) – Maximum number of warnings when using gmx grompp

  • monitor – option to monitor a simulation, if not none monitor should contains two values: function the function to be ran while simulation is running and input parameters for the function

Object requirement(s):

  • self.gmxsys.coor_file

  • self.gmxsys.top_file

  • self.gmxsys.nt

  • self.gmxsys.ntmpi

  • self.gmxsys.gpu_id

Object field(s) changed:

  • self.lambda_sys_list

  • self.lambda_coul

  • self.lambda_vdw

  • self.lambda_restr

  • self.prod_time

  • self.xvg_file_list

show_intermol_restr()

Show traj with atom implied in intermolecular restraints using nglview library.

static symmetry_correction(smile, temp=300)

Compute symmetry correction \(\Delta_{sym} = −k T ln(\sigma)\)

return value in kcal/mol

> FreeEner.symmetry_correction('c1ccccc1')
-1.4814...
property unit_graph

Return unit as math latex for matplotlib label purpose.

property unit_name
water_box_from_SMILE(smile, method_3d='rdkit', iter_num=5000, box_dist=1.1)

Create a water box coordinates and toplogie with a molecule defined by its SMILE.

Parameters

  • smile (str) – Molecule’s SMILE

  • method_3d (str, default='rdkit') – Method to compute 3D coordinates

  • iter_num (int, default=5000) – Iteration step number for 3D coordinate computation

  • box_dist (float, default=1.1 nm) – Ditance to egde box (nm)

Note

Default box dist 1.1 nm was taken as minimal distance for CH4 molecule, to allow domain decomposition with gmx mdrun.

gromacs_py.free_ener.show_log()

To use only with Doctest !!! Redirect logger output to sys.stdout

Module contents