MDP utilities#
- class bindflow.mdp.mdp.StepMDP(step: str | None = None, step_path: PathLike | str | bytes | None = None, **kwargs)[source]#
This subclass will inherit from
bindflow.mdp.mdp.MDPIt is meant to be used in combination with the templates that can be access frombindflow.mdp.templates.TemplatePath. This class define the methodset_new_step. One time initialized, the instance could be used to access other steps on the step_path- Parameters:
MDP (
bindflow.mdp.mdp.MDP) – base MDP class
- __init__(step: str | None = None, step_path: PathLike | str | bytes | None = None, **kwargs)[source]#
Constructor. It is assume a tree directory as:
. ├── emin.mdp ├── npt.mdp ├── npt-norest.mpd ├── nvt.mdp └── prod.mdp
- Parameters:
step (str, optional) – the step, basically the name of the mdp file on templates, by default None
step_path (PathLike, optional) – where to look for the mdp, by default None
Example
In [1]: from bindflow.mdp import mdp In [2]: from bindflow.mdp.templates import TemplatePath In [3]: from pathlib import Path In [4]: my_mdp = mdp.StepMDP(step='00_min', step_path=Path(TemplatePath.ligand.fep)/'coul') In [5]: my_mdp.set_parameters(**{"init-lambda-state": "0", "coul-lambdas": "0 0.5 1"}) In [6]: print(my_mdp) StepMDP({ "define": "-DFLEXIBLE", "integrator": "steep", "nsteps": "10000", "emtol": "100", "emstep": "0.01", "nstcomm": "100", "nstxout": "0 ; save coordinates to .trr every 250 steps", "nstvout": "0 ; don't save velocities to .trr", "nstfout": "0 ; don't save forces to .trr", "nstxout-compressed": "0 ; xtc compressed trajectory output every 500 steps", "compressed-x-precision": "0", "nstlog": "0 ; update log file every 500 steps", "nstenergy": "0 ; save energies every 500 steps", "nstcalcenergy": "100", "cutoff-scheme": "Verlet", "ns-type": "grid", "nstlist": "1", "rlist": "1.0", "constraints": "none", "coulombtype": "PME", "rcoulomb": "1.0", "pme-order": "4", "fourierspacing": "0.10", "ewald-rtol": "1e-6", "vdwtype": "Cut-off", "vdw-modifier": "Potential-shift-Verlet", "verlet-buffer-tolerance": "0.005", "rvdw": "1.0 ; short-range van der Waals cutoff (in nm)", "DispCorr": "EnerPres ; apply long range dispersion corrections for Energy and Pressure", "tcoupl": "no", "pcoupl": "no", "gen-vel": "no", "pbc": "xyz", "free-energy": "yes", "couple-moltype": "LIG", "couple-lambda0": "vdw-q", "couple-lambda1": "vdw", "init-lambda-state": "0", "coul-lambdas": "0 0.5 1", "nstdhdl": "0", "dhdl-print-energy": "total", "calc-lambda-neighbors": "-1", "separate-dhdl-file": "yes", "couple-intramol": "yes" }) In [7]: my_mdp.set_new_step(step='01_nvt') Out[7]: StepMDP({ "define": "-DPOSRES", "integrator": "sd ; stochastic leap-frog integrator", "nsteps": "5000 ; 2 * 5,000 fs = 10 ps", "dt": "0.002 ; 2 fs", "comm-mode": "Linear ; remove center of mass translation", "nstcomm": "50 ; frequency for center of mass motion removal", "nstxout": "0 ; save coordinates to .trr every 100 ps", "nstvout": "0 ; don't save velocities to .trr", "nstfout": "0 ; don't save forces to .trr", "nstxout-compressed": "0 ; xtc compressed trajectory output every 2 ps", "compressed-x-precision": "0 ; precision with which to write to the compressed trajectory file", "nstlog": "0 ; update log file every 2 ps", "nstenergy": "0 ; save energies every 2 ps", "nstcalcenergy": "50 ; calculate energies every 200 fs", "constraint-algorithm": "lincs ; holonomic constraints", "constraints": "all-bonds ; all bonds are constrained (HMR)", "lincs-iter": "1 ; accuracy of LINCS (1 is default)", "lincs-order": "6 ; also related to accuracy (4 is default)", "lincs-warnangle": "30 ; maximum angle that a bond can rotate before LINCS will complain (30 is default)", "continuation": "no ; formerly known as 'unconstrained-start' - useful for exact continuations and reruns", "cutoff-scheme": "Verlet", "ns-type": "grid ; search neighboring grid cells", "nstlist": "10 ; 20 fs (default is 10)", "rlist": "1.2 ; short-range neighborlist cutoff (in nm)", "pbc": "xyz ; 3D PBC", "coulombtype": "PME ; Particle Mesh Ewald for long-range electrostatics", "rcoulomb": "1.0 ; short-range electrostatic cutoff (in nm)", "ewald-geometry": "3d ; Ewald sum is performed in all three dimensions", "pme-order": "4 ; interpolation order for PME (default is 4)", "fourierspacing": "0.10 ; grid spacing for FFT", "ewald-rtol": "1e-6 ; relative strength of the Ewald-shifted direct potential at rcoulomb", "vdwtype": "Cut-off", "vdw-modifier": "Potential-shift-Verlet", "verlet-buffer-tolerance": "0.005", "rvdw": "1.0 ; short-range van der Waals cutoff (in nm)", "DispCorr": "EnerPres ; apply long range dispersion corrections for Energy and Pressure", "tc-grps": "System", "tau-t": "2.0", "ref-t": "298.15", "pcoupl": "no", "gen-vel": "yes ; Velocity generation is on (if gen_vel is 'yes', continuation should be 'no')", "gen-seed": "-1 ; Use random seed", "gen-temp": "298.15", "free-energy": "yes", "couple-moltype": "LIG", "couple-lambda0": "vdw-q", "couple-lambda1": "vdw", "init-lambda-state": "<state>", "coul-lambdas": "<lamRange>", "nstdhdl": "0", "dhdl-print-energy": "total", "calc-lambda-neighbors": "-1", "separate-dhdl-file": "yes", "couple-intramol": "yes" }) In [8]: print(my_mdp.to_string()) define = -DPOSRES integrator = sd ; stochastic leap-frog integrator nsteps = 5000 ; 2 * 5,000 fs = 10 ps dt = 0.002 ; 2 fs comm-mode = Linear ; remove center of mass translation nstcomm = 50 ; frequency for center of mass motion removal nstxout = 0 ; save coordinates to .trr every 100 ps nstvout = 0 ; don't save velocities to .trr nstfout = 0 ; don't save forces to .trr nstxout-compressed = 0 ; xtc compressed trajectory output every 2 ps compressed-x-precision = 0 ; precision with which to write to the compressed trajectory file nstlog = 0 ; update log file every 2 ps nstenergy = 0 ; save energies every 2 ps nstcalcenergy = 50 ; calculate energies every 200 fs constraint-algorithm = lincs ; holonomic constraints constraints = all-bonds ; all bonds are constrained (HMR) lincs-iter = 1 ; accuracy of LINCS (1 is default) lincs-order = 6 ; also related to accuracy (4 is default) lincs-warnangle = 30 ; maximum angle that a bond can rotate before LINCS will complain (30 is default) continuation = no ; formerly known as 'unconstrained-start' - useful for exact continuations and reruns cutoff-scheme = Verlet ns-type = grid ; search neighboring grid cells nstlist = 10 ; 20 fs (default is 10) rlist = 1.2 ; short-range neighborlist cutoff (in nm) pbc = xyz ; 3D PBC coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm) ewald-geometry = 3d ; Ewald sum is performed in all three dimensions pme-order = 4 ; interpolation order for PME (default is 4) fourierspacing = 0.10 ; grid spacing for FFT ewald-rtol = 1e-6 ; relative strength of the Ewald-shifted direct potential at rcoulomb vdwtype = Cut-off vdw-modifier = Potential-shift-Verlet verlet-buffer-tolerance = 0.005 rvdw = 1.0 ; short-range van der Waals cutoff (in nm) DispCorr = EnerPres ; apply long range dispersion corrections for Energy and Pressure tc-grps = System tau-t = 2.0 ref-t = 298.15 pcoupl = no gen-vel = yes ; Velocity generation is on (if gen_vel is 'yes', continuation should be 'no') gen-seed = -1 ; Use random seed gen-temp = 298.15 free-energy = yes couple-moltype = LIG couple-lambda0 = vdw-q couple-lambda1 = vdw init-lambda-state = <state> coul-lambdas = <lamRange> nstdhdl = 0 dhdl-print-energy = total calc-lambda-neighbors = -1 separate-dhdl-file = yes couple-intramol = yes
- bindflow.mdp.mdp.make_fep_dir_structure(sim_dir: PathLike | str | bytes, template_dir: PathLike | str | bytes, lambda_values: List[float], lambda_type: str, sys_type: str, dt_max: float, mdp_extra_kwargs: dict | None = None)[source]#
This function is meant to be used on
ligand_fep_setupandcomplex_fet_setuprules. It will create the structure of the simulation directory:{sim_dir}/simulation/{lambda_type}.{i}/{step}/{step}.mdpWhere:
i: init-lambda-state,
step: the name of the simulation to carry on
- Parameters:
sim_dir (PathLike) – Where the simulation suppose to run
template_dir (PathLike) – This is the directory that storage the mdp templates: bindflow.mdp.templates.TemplatePath.ligand.fep or bindlfow.mdp.templates.TemplatePath.complex.fep
lambda_values (List[float]) – This is a the list of lambda values to be used inside the mdp on the entrance {lambda_type}-lambdas
lambda_type (str) – Must be one of the following strings “vdw”, “coul”, “bonded” (the last is for restraints)
sys_type (str) – Must one of the following strings “ligand” or “complex”. This is used in order to turn on the bonded lambdas for the complex simulations
mdp_extra_kwargs (dict) –
The MDP options for the fep calculations on every step. This dictionary must have the structure:
{ 'vdw':{ 'step1': <mdp options>, 'step2': <mdp options>, ... } 'coul':{ 'step1': <mdp options>, 'step2': <mdp options>, ... 'bonded':{ 'step1': <mdp options>, 'step2': <mdp options>, ... } }
- Raises:
ValueError – In case of an invalid
lambda_typeValueError – In case of an invalid
sys_type