CI update -> use miniconda v2
jeff231li
| 1 | 38 |
import abc |
| 2 | 38 |
import glob |
| 3 | 38 |
import logging |
| 4 | 38 |
import os |
| 5 | 38 |
import subprocess as sp |
| 6 | 38 |
from collections import OrderedDict |
| 7 | 38 |
from enum import Enum |
| 8 |
|
|
| 9 | 38 |
from paprika.utils import get_dict_without_keys |
| 10 |
|
|
| 11 | 38 |
from .simulation import Simulation |
| 12 |
|
|
| 13 | 38 |
logger = logging.getLogger(__name__) |
| 14 |
|
|
| 15 |
|
|
| 16 | 38 |
class GROMACS(Simulation, abc.ABC): |
| 17 |
"""
|
|
| 18 |
A wrapper that can be used to set GROMACS simulation parameters.
|
|
| 19 |
|
|
| 20 |
.. todo ::
|
|
| 21 |
possibly modify this module to use the official python wrapper of GROMACS.
|
|
| 22 |
|
|
| 23 |
Below is an example of the configuration file (``gromacs.mdp``) generated by the wrapper. The class property
|
|
| 24 |
associated with defining the configuration variables is shown in brackets.
|
|
| 25 |
|
|
| 26 |
.. code ::
|
|
| 27 |
|
|
| 28 |
title = NPT MD Simulation ; [self.title]
|
|
| 29 |
; Run control [self.control]
|
|
| 30 |
nsteps = 1500000
|
|
| 31 |
nstxout = 500
|
|
| 32 |
nstlog = 500
|
|
| 33 |
nstenergy = 500
|
|
| 34 |
nstcalcenergy = 500
|
|
| 35 |
dt = 0.002
|
|
| 36 |
integrator = md
|
|
| 37 |
; Nonbonded options [self.nb_method]
|
|
| 38 |
cutoff-scheme = Verlet
|
|
| 39 |
ns_type = grid
|
|
| 40 |
nstlist = 10
|
|
| 41 |
rlist = 0.9
|
|
| 42 |
rcoulomb = 0.9
|
|
| 43 |
rvdw = 0.9
|
|
| 44 |
coulombtype = PME
|
|
| 45 |
pme_order = 4
|
|
| 46 |
fourierspacing = 0.16
|
|
| 47 |
vdwtype = Cut-off
|
|
| 48 |
DispCorr = EnerPres
|
|
| 49 |
pbc = xyz
|
|
| 50 |
; Bond constraints [self.constraints]
|
|
| 51 |
constraint-algorithm = lincs
|
|
| 52 |
constraints = h-bonds
|
|
| 53 |
lincs_iter = 1
|
|
| 54 |
lincs_order = 4
|
|
| 55 |
; Temperature coupling [self.thermostat]
|
|
| 56 |
tcoupl = v-rescale
|
|
| 57 |
tc-grps = System
|
|
| 58 |
ref_t = 298.15
|
|
| 59 |
tau_t = 0.1
|
|
| 60 |
gen_vel = no
|
|
| 61 |
; Pressure coupling [self.barostat]
|
|
| 62 |
pcoupl = Berendsen
|
|
| 63 |
pcoupltype = isotropic
|
|
| 64 |
tau_p = 2.0
|
|
| 65 |
ref_p = 1.01325
|
|
| 66 |
compressibility = 4.5e-05
|
|
| 67 |
"""
|
|
| 68 |
|
|
| 69 | 38 |
class Thermostat(Enum): |
| 70 |
"""
|
|
| 71 |
An enumeration of the different themostat implemented in GROMACS.
|
|
| 72 |
"""
|
|
| 73 |
|
|
| 74 | 38 |
Off = "no" |
| 75 | 38 |
Berendsen = "berendsen" |
| 76 | 38 |
NoseHoover = "nose-hoover" |
| 77 | 38 |
Andersen1 = "andersen" |
| 78 | 38 |
Andersen2 = "andersen-massive" |
| 79 | 38 |
VelocityRescaling = "v-rescale" |
| 80 |
|
|
| 81 | 38 |
class Barostat(Enum): |
| 82 |
"""
|
|
| 83 |
An enumeration of the different barostat implemented in GROMACS.
|
|
| 84 |
"""
|
|
| 85 |
|
|
| 86 | 38 |
Off = "no" |
| 87 | 38 |
Berendsen = "Berendsen" |
| 88 | 38 |
ParrinelloRahman = "Parrinello-Rahman" |
| 89 | 38 |
MMTK = "MTTK" |
| 90 |
|
|
| 91 | 38 |
class Integrator(Enum): |
| 92 |
"""
|
|
| 93 |
An enumeration of the different integrators implemented in GROMACS.
|
|
| 94 |
"""
|
|
| 95 |
|
|
| 96 | 38 |
LeapFrog = "md" |
| 97 | 38 |
VelocityVerlet = "md-vv" |
| 98 | 38 |
VelocityVerletAveK = "md-vv-avek" |
| 99 | 38 |
LangevinDynamics = "sd" |
| 100 | 38 |
BrownianDynamics = "bd" |
| 101 |
|
|
| 102 | 38 |
class Optimizer(Enum): |
| 103 |
"""
|
|
| 104 |
An enumeration of the different minimization algorithm implemented in GROMACS.
|
|
| 105 |
"""
|
|
| 106 |
|
|
| 107 | 38 |
SteepestDescent = "steep" |
| 108 | 38 |
ConjugateGradient = "cg" |
| 109 | 38 |
Broyden = "l-bfgs" |
| 110 |
|
|
| 111 | 38 |
class BoxScaling(Enum): |
| 112 |
"""
|
|
| 113 |
An enumeration of the different PBC scaling options when running constant pressure simulations in GROMACS.
|
|
| 114 |
"""
|
|
| 115 |
|
|
| 116 | 38 |
Isotropic = "isotropic" |
| 117 | 38 |
Semiisotropic = "semiisotropic" |
| 118 | 38 |
Anisotropic = "anisotropic" |
| 119 | 38 |
SurfaceTension = "surface-tension" |
| 120 |
|
|
| 121 | 38 |
class Constraints(Enum): |
| 122 |
"""
|
|
| 123 |
An enumeration of the different bond constraint options in GROMACS.
|
|
| 124 |
"""
|
|
| 125 |
|
|
| 126 | 38 |
Off = "none" |
| 127 | 38 |
HBonds = "h-bonds" |
| 128 | 38 |
AllBonds = "all-bonds" |
| 129 | 38 |
HAngles = "h-angles" |
| 130 | 38 |
AllAngles = "all-angles" |
| 131 |
|
|
| 132 | 38 |
@property
|
| 133 | 38 |
def index_file(self) -> str: |
| 134 |
"""os.PathLike: GROMACS index file that specifies ``groups`` in the system. This is optional in a GROMACS
|
|
| 135 |
simulation."""
|
|
| 136 |
return self._index_file |
|
| 137 |
|
|
| 138 | 38 |
@index_file.setter
|
| 139 | 38 |
def index_file(self, value: str): |
| 140 |
self._index_file = value |
|
| 141 |
|
|
| 142 | 38 |
@property
|
| 143 | 38 |
def checkpoint(self) -> str: |
| 144 |
"""os.PathLike: Checkpoint file (extension is ``.cpt``) for starting a simulation from a previous state."""
|
|
| 145 | 38 |
return self._checkpoint |
| 146 |
|
|
| 147 | 38 |
@checkpoint.setter
|
| 148 | 38 |
def checkpoint(self, value: str): |
| 149 |
self._checkpoint = value |
|
| 150 |
|
|
| 151 | 38 |
@property
|
| 152 |
def control(self): |
|
| 153 |
"""dict: Dictionary for the output control of the MD simulation (frequency of energy, trajectory etc)."""
|
|
| 154 | 38 |
return self._control |
| 155 |
|
|
| 156 | 38 |
@control.setter
|
| 157 |
def control(self, value): |
|
| 158 |
self._control = value |
|
| 159 |
|
|
| 160 | 38 |
@property
|
| 161 |
def nb_method(self): |
|
| 162 |
"""dict: Dictionary for the non-bonded method options (cutoffs and methods)."""
|
|
| 163 | 38 |
return self._nb_method |
| 164 |
|
|
| 165 | 38 |
@nb_method.setter
|
| 166 |
def nb_method(self, value): |
|
| 167 |
self._nb_method = value |
|
| 168 |
|
|
| 169 | 38 |
@property
|
| 170 |
def constraints(self): |
|
| 171 |
"""dict: Dictionary for the bond constraint options (LINCS or SHAKE)."""
|
|
| 172 | 38 |
return self._constraints |
| 173 |
|
|
| 174 | 38 |
@constraints.setter
|
| 175 |
def constraints(self, value): |
|
| 176 |
self._constraints = value |
|
| 177 |
|
|
| 178 | 38 |
@property
|
| 179 | 38 |
def tc_groups(self) -> list: |
| 180 |
"""
|
|
| 181 |
list: List of groups to apply thermostat "separately" based on the groups defined in the ``index_file``.
|
|
| 182 |
Below is an example of applying the thermostat for different groups separately in a GROMACS input file
|
|
| 183 |
|
|
| 184 |
.. code ::
|
|
| 185 |
|
|
| 186 |
tcoupl = v-rescale
|
|
| 187 |
tc-grps = HOST GUEST HOH
|
|
| 188 |
tau-t = 0.1 0.1 0.1
|
|
| 189 |
ref-t = 300 300 300
|
|
| 190 |
"""
|
|
| 191 | 38 |
return self._tc_groups |
| 192 |
|
|
| 193 | 38 |
@tc_groups.setter
|
| 194 | 38 |
def tc_groups(self, value: list): |
| 195 | 38 |
self._tc_groups = value |
| 196 |
|
|
| 197 | 38 |
@property
|
| 198 |
def prefix(self): |
|
| 199 |
"""str: The prefix for file names generated from this simulation."""
|
|
| 200 | 38 |
return self._prefix |
| 201 |
|
|
| 202 | 38 |
@prefix.setter
|
| 203 |
def prefix(self, new_prefix): |
|
| 204 | 38 |
self._prefix = new_prefix |
| 205 | 38 |
self.input = new_prefix + ".mdp" |
| 206 | 38 |
self.output = new_prefix + ".mdout" |
| 207 | 38 |
self.logfile = new_prefix + ".log" |
| 208 | 38 |
self.tpr = new_prefix + ".tpr" |
| 209 |
|
|
| 210 | 38 |
@property
|
| 211 | 38 |
def custom_mdrun_command(self) -> str: |
| 212 |
"""Custom commands for ``mdrun``. The default commands parsed to ``mdrun`` if all the variables are defined is
|
|
| 213 |
|
|
| 214 |
.. code::
|
|
| 215 |
|
|
| 216 |
gmx mdrun -deffnm ``prefix`` -nt ``n_threads`` -gpu_id ``gpu_devices`` -plumed ``plumed.dat``
|
|
| 217 |
|
|
| 218 |
This is useful depending on how GROMACS was compiled, e.g. if GROMACS is compiled with the MPI library the
|
|
| 219 |
you will need to use the command below:
|
|
| 220 |
|
|
| 221 |
.. code::
|
|
| 222 |
|
|
| 223 |
mpirun -np 6 gmx_mpi mdrun -deffnm ``prefix`` -ntomp 1 -gpu_id 0 -plumed ``plumed.dat``
|
|
| 224 |
"""
|
|
| 225 |
return self._custom_mdrun_command |
|
| 226 |
|
|
| 227 | 38 |
@custom_mdrun_command.setter
|
| 228 | 38 |
def custom_mdrun_command(self, value: str): |
| 229 |
self._custom_mdrun_command = value |
|
| 230 |
|
|
| 231 | 38 |
@property
|
| 232 | 38 |
def grompp_maxwarn(self) -> int: |
| 233 |
"""int: Maximum number of warnings for GROMPP to ignore. default=1."""
|
|
| 234 |
return self._grompp_maxwarn |
|
| 235 |
|
|
| 236 | 38 |
@grompp_maxwarn.setter
|
| 237 | 38 |
def grompp_maxwarn(self, value: int): |
| 238 |
self._grompp_maxwarn = value |
|
| 239 |
|
|
| 240 | 38 |
def __init__(self): |
| 241 |
|
|
| 242 | 38 |
super().__init__() |
| 243 |
|
|
| 244 |
# I/O
|
|
| 245 | 38 |
self._index_file = None |
| 246 | 38 |
self._custom_mdrun_command = None |
| 247 | 38 |
self._tc_groups = None |
| 248 | 38 |
self._grompp_maxwarn = 1 |
| 249 |
|
|
| 250 |
# File names
|
|
| 251 | 38 |
self.input = self._prefix + ".mdp" |
| 252 | 38 |
self.output = self._prefix + ".mdout" |
| 253 | 38 |
self._checkpoint = None |
| 254 | 38 |
self.logfile = self._prefix + ".log" |
| 255 | 38 |
self.tpr = self._prefix + ".tpr" |
| 256 |
|
|
| 257 |
# Input file
|
|
| 258 | 38 |
self._control = OrderedDict() |
| 259 | 38 |
self._control["nsteps"] = 5000 |
| 260 | 38 |
self._control["nstxout"] = 500 |
| 261 | 38 |
self._control["nstlog"] = 500 |
| 262 | 38 |
self._control["nstenergy"] = 500 |
| 263 | 38 |
self._control["nstcalcenergy"] = 500 |
| 264 |
|
|
| 265 | 38 |
self._constraints = OrderedDict() |
| 266 | 38 |
self._constraints["constraint-algorithm"] = "lincs" |
| 267 | 38 |
self._constraints["constraints"] = self.Constraints.HBonds.value |
| 268 | 38 |
self._constraints["lincs_iter"] = 1 |
| 269 | 38 |
self._constraints["lincs_order"] = 4 |
| 270 |
|
|
| 271 | 38 |
self._nb_method = OrderedDict() |
| 272 | 38 |
self._nb_method["cutoff-scheme"] = "Verlet" |
| 273 | 38 |
self._nb_method["ns-type"] = "grid" |
| 274 | 38 |
self._nb_method["nstlist"] = 10 |
| 275 | 38 |
self._nb_method["rlist"] = 0.9 |
| 276 | 38 |
self._nb_method["rcoulomb"] = 0.9 |
| 277 | 38 |
self._nb_method["rvdw"] = 0.9 |
| 278 | 38 |
self._nb_method["coulombtype"] = "PME" |
| 279 | 38 |
self._nb_method["pme_order"] = 4 |
| 280 | 38 |
self._nb_method["fourierspacing"] = 0.16 |
| 281 | 38 |
self._nb_method["vdwtype"] = "Cut-off" |
| 282 | 38 |
self._nb_method["DispCorr"] = "EnerPres" |
| 283 | 38 |
self._nb_method["pbc"] = "xyz" |
| 284 |
|
|
| 285 | 38 |
def _config_min(self, optimizer): |
| 286 |
"""
|
|
| 287 |
Configure input settings for a minimization run.
|
|
| 288 |
|
|
| 289 |
Parameters
|
|
| 290 |
----------
|
|
| 291 |
optimizer: :class:`GROMACS.Optimizer`, default=Optimizer.SteepestDescent
|
|
| 292 |
Algorithm for energy minimization, keyword in the parenthesis are the options for the input file.
|
|
| 293 |
**(1)** `SteepestDescent` (``steep``), **(2)** `ConjugateGradient` (``cg``), and **(3)** `Broyden`
|
|
| 294 |
(``l-bfgs``).
|
|
| 295 |
"""
|
|
| 296 | 38 |
self.constraints["continuation"] = "no" |
| 297 | 38 |
self.control["integrator"] = optimizer.value |
| 298 | 38 |
self.control["emtol"] = 10.0 |
| 299 | 38 |
self.control["emstep"] = 0.01 |
| 300 | 38 |
self.control["nsteps"] = 5000 |
| 301 |
|
|
| 302 | 38 |
def _config_md(self, integrator, thermostat): |
| 303 |
"""
|
|
| 304 |
Configure input setting for a MD.
|
|
| 305 |
|
|
| 306 |
Parameters
|
|
| 307 |
----------
|
|
| 308 |
integrator: :class:`GROMACS.Integrator`, default=Integrator.LeapFrog
|
|
| 309 |
Option to choose the integrator for the MD simulations, keywords in the parenthesis are the options for the
|
|
| 310 |
input file. **(1)** `LeapFrog` (``md``), **(2)** `VelocityVerlet` (``md-vv``),
|
|
| 311 |
**(3)** `VelocityVerletAveK` (``md-vv-avek``), **(4)** `LangevinDynamics` (``sd``), and **(5)**
|
|
| 312 |
`Brownian Dynamics` (``bd``).
|
|
| 313 |
integrator: :class:`GROMACS.Integrator`, default=Integrator.LeapFrog
|
|
| 314 |
Option to choose the integrator for the MD simulations, keywords in the parenthesis are the options for the
|
|
| 315 |
input file. **(1)** `LeapFrog` (``md``), **(2)** `VelocityVerlet` (``md-vv``),
|
|
| 316 |
**(3)** `VelocityVerletAveK` (``md-vv-avek``), **(4)** `LangevinDynamics` (``sd``), and **(5)**
|
|
| 317 |
`Brownian Dynamics` (``bd``).
|
|
| 318 |
"""
|
|
| 319 | 38 |
self.control["dt"] = 0.002 |
| 320 | 38 |
self.control["integrator"] = integrator.value |
| 321 | 38 |
self.constraints["continuation"] = "yes" |
| 322 | 38 |
self.thermostat["tc-grps"] = "System" |
| 323 | 38 |
self.thermostat["ref_t"] = self.temperature |
| 324 |
|
|
| 325 | 38 |
if ( |
| 326 |
integrator != self.Integrator.LangevinDynamics |
|
| 327 |
and integrator != self.Integrator.BrownianDynamics |
|
| 328 |
):
|
|
| 329 | 38 |
self.thermostat["tcoupl"] = thermostat.value |
| 330 | 38 |
self.thermostat["tau_t"] = 1.0 |
| 331 |
else: |
|
| 332 | 38 |
self.thermostat["tau_t"] = 0.1 |
| 333 |
|
|
| 334 | 38 |
def config_vac_min(self, optimizer=Optimizer.SteepestDescent): |
| 335 |
"""
|
|
| 336 |
Configure a reasonable input setting for a MD run in vacuum. `Users can override the parameters set by this
|
|
| 337 |
method.`
|
|
| 338 |
|
|
| 339 |
.. note ::
|
|
| 340 |
Newer versions of GMX no longer support a "True" vacuum simulation so we have to do this by creating a
|
|
| 341 |
"pseudo-PBC" environment. Make sure the coordinates ``.gro`` file has an expanded box, which you can do
|
|
| 342 |
using ``gmx editconf``. See the discussion on
|
|
| 343 |
https://gromacs.bioexcel.eu/t/minimization-in-vacuum-without-pbc/110/2.
|
|
| 344 |
|
|
| 345 |
Parameters
|
|
| 346 |
----------
|
|
| 347 |
optimizer: :class:`GROMACS.Optimizer`, default=Optimizer.SteepestDescent
|
|
| 348 |
Algorithm for energy minimization, keyword in the parenthesis are the options for the input file.
|
|
| 349 |
**(1)** `SteepestDescent` (``steep``), **(2)** `ConjugateGradient` (``cg``), and **(3)** `Broyden`
|
|
| 350 |
(``l-bfgs``).
|
|
| 351 |
"""
|
|
| 352 |
|
|
| 353 |
self.title = "Vacuum Minimization" |
|
| 354 |
|
|
| 355 |
self._config_min(optimizer) |
|
| 356 |
|
|
| 357 |
self.nb_method["pbc"] = "xyz" |
|
| 358 |
self.nb_method["ns_type"] = "grid" |
|
| 359 |
self.nb_method["nstlist"] = 10 |
|
| 360 |
self.nb_method["rlist"] = 333.3 |
|
| 361 |
self.nb_method["coulombtype"] = "Cut-off" |
|
| 362 |
self.nb_method["rcoulomb"] = 333.3 |
|
| 363 |
self.nb_method["vdwtype"] = "Cut-off" |
|
| 364 |
self.nb_method["rvdw"] = 333.3 |
|
| 365 |
self.nb_method["DispCorr"] = "no" |
|
| 366 |
|
|
| 367 | 38 |
def config_vac_md( |
| 368 |
self, integrator=Integrator.LeapFrog, thermostat=Thermostat.VelocityRescaling |
|
| 369 |
):
|
|
| 370 |
"""
|
|
| 371 |
Configure a reasonable input setting for a MD run in vacuum. `Users can override the parameters set by this
|
|
| 372 |
method.`
|
|
| 373 |
|
|
| 374 |
.. note ::
|
|
| 375 |
Newer versions of GMX no longer support a "True" vacuum simulation so we have to do this by creating a
|
|
| 376 |
"pseudo-PBC" environment. Make sure the coordinates ``.gro`` file has an expanded box, which you set
|
|
| 377 |
using ``gmx editconf``. See the discussion on
|
|
| 378 |
https://gromacs.bioexcel.eu/t/minimization-in-vacuum-without-pbc/110/2.
|
|
| 379 |
|
|
| 380 |
Parameters
|
|
| 381 |
----------
|
|
| 382 |
integrator: :class:`GROMACS.Integrator`, default=Integrator.LeapFrog
|
|
| 383 |
Option to choose the integrator for the MD simulations, keywords in the parenthesis are the options for the
|
|
| 384 |
input file. **(1)** `LeapFrog` (``md``), **(2)** `VelocityVerlet` (``md-vv``),
|
|
| 385 |
**(3)** `VelocityVerletAveK` (``md-vv-avek``), **(4)** `LangevinDynamics` (``sd``), and **(5)**
|
|
| 386 |
`Brownian Dynamics` (``bd``).
|
|
| 387 |
thermostat: :class:`GROMACS.Thermostat`, default=Thermostat.VelocityRescaling
|
|
| 388 |
Option to choose one of five thermostat implemented in GROMACS, keywords in the parenthesis are the options
|
|
| 389 |
for the input file. **(1)** `Off` (``no``), **(2)** `Berendsen` (``berendsen``), **(3)** `NoseHoover`
|
|
| 390 |
(``nose-hoover``), **(4)** `Andersen1` (``andersen``), **(5)** `Andersen2` (``andersen-massive``),
|
|
| 391 |
and **(6)** `VelocityRescaling` (``v-rescale``).
|
|
| 392 |
"""
|
|
| 393 | 38 |
self.title = "Vacuum MD Simulation" |
| 394 |
|
|
| 395 | 38 |
self._config_md(integrator, thermostat) |
| 396 |
|
|
| 397 | 38 |
if self.checkpoint is None: |
| 398 | 38 |
self.constraints["continuation"] = "no" |
| 399 |
else: |
|
| 400 |
self.constraints["continuation"] = "yes" |
|
| 401 |
|
|
| 402 | 38 |
self.nb_method["pbc"] = "xyz" |
| 403 | 38 |
self.nb_method["ns_type"] = "grid" |
| 404 | 38 |
self.nb_method["nstlist"] = 10 |
| 405 | 38 |
self.nb_method["rlist"] = 333.3 |
| 406 | 38 |
self.nb_method["coulombtype"] = "Cut-off" |
| 407 | 38 |
self.nb_method["rcoulomb"] = 333.3 |
| 408 | 38 |
self.nb_method["vdwtype"] = "Cut-off" |
| 409 | 38 |
self.nb_method["rvdw"] = 333.3 |
| 410 | 38 |
self.nb_method["DispCorr"] = "no" |
| 411 |
|
|
| 412 | 38 |
def config_pbc_min(self, optimizer=Optimizer.SteepestDescent): |
| 413 |
"""
|
|
| 414 |
Configure a reasonable input setting for an energy minimization run with periodic boundary conditions. `Users
|
|
| 415 |
can override the parameters set by this method.`
|
|
| 416 |
|
|
| 417 |
Parameters
|
|
| 418 |
----------
|
|
| 419 |
optimizer: :class:`GROMACS.Optimizer`, default=Optimizer.SteepestDescent
|
|
| 420 |
Algorithm for energy minimization, keywords in the parenthesis are the options for the input file.
|
|
| 421 |
**(1)** `SteepestDescent` (``steep``), **(2)** `ConjugateGradient` (``cg``), and **(3)** `Broyden`
|
|
| 422 |
(``l-bfgs``).
|
|
| 423 |
"""
|
|
| 424 |
|
|
| 425 | 38 |
self.title = "PBC Minimization" |
| 426 |
|
|
| 427 | 38 |
self._config_min(optimizer) |
| 428 |
|
|
| 429 | 38 |
self.nb_method["nstlist"] = 10 |
| 430 |
|
|
| 431 | 38 |
def config_pbc_md( |
| 432 |
self, |
|
| 433 |
ensemble=Simulation.Ensemble.NPT, |
|
| 434 |
integrator=Integrator.LeapFrog, |
|
| 435 |
thermostat=Thermostat.VelocityRescaling, |
|
| 436 |
barostat=Barostat.Berendsen, |
|
| 437 |
):
|
|
| 438 |
"""
|
|
| 439 |
Configure a reasonable input setting for a MD run with periodic boundary conditions. `Users can override the
|
|
| 440 |
parameters set by this method.`
|
|
| 441 |
|
|
| 442 |
Parameters
|
|
| 443 |
----------
|
|
| 444 |
ensemble: :class:`Simulation.Ensemble`, default=Ensemble.NPT
|
|
| 445 |
Configure a MD simulation with NVE, NVT or NPT thermodynamic ensemble.
|
|
| 446 |
integrator: :class:`GROMACS.Integrator`, default=Integrator.LeapFrog
|
|
| 447 |
Option to choose the integrator for the MD simulations, keywords in the parenthesis are the options for the
|
|
| 448 |
input file. **(1)** `LeapFrog` (``md``), **(2)** `VelocityVerlet` (``md-vv``),
|
|
| 449 |
**(3)** `VelocityVerletAveK` (``md-vv-avek``), **(4)** `LangevinDynamics` (``sd``), and **(5)**
|
|
| 450 |
`Brownian Dynamics` (``bd``).
|
|
| 451 |
thermostat: :class:`GROMACS.Thermostat`, default=Thermostat.VelocityRescaling
|
|
| 452 |
Option to choose one of five thermostat implemented in GROMACS, keywords in the parenthesis are the options
|
|
| 453 |
for the input file. **(1)** `Off` (``no``), **(2)** `Berendsen` (``berendsen``), **(3)** `NoseHoover`
|
|
| 454 |
(``nose-hoover``), **(4)** `Andersen1` (``andersen``), **(5)** `Andersen2` (``andersen-massive``),
|
|
| 455 |
and **(6)** `VelocityRescaling` (``v-rescale``).
|
|
| 456 |
barostat: :class:`GROMACS.Barostat`, default=Barostat.Berendsen
|
|
| 457 |
Option to choose one of three barostat implemented in GROMACS, keywords in the parenthesis are the options
|
|
| 458 |
for the input file. **(1)** `Off` (``no``), **(2)** `Berendsen` (``berendsen``), **(3)** `ParrinelloRahman`
|
|
| 459 |
(``Parrinello-Rahman``), and **(4)** `MMTK` (``MTTK``).
|
|
| 460 |
"""
|
|
| 461 | 38 |
self.title = f"{ensemble.value} MD Simulation" |
| 462 |
|
|
| 463 | 38 |
self._config_md(integrator, thermostat) |
| 464 |
|
|
| 465 | 38 |
if self.checkpoint is None: |
| 466 | 38 |
self.constraints["continuation"] = "no" |
| 467 |
else: |
|
| 468 |
self.constraints["continuation"] = "yes" |
|
| 469 |
|
|
| 470 | 38 |
if ensemble == self.Ensemble.NVE: |
| 471 |
self.thermostat["tcoupl"] = self.Thermostat.Off.value |
|
| 472 |
self.barostat["pcoupl"] = self.Barostat.Off.value |
|
| 473 |
del self.thermostat["tc-grps"] |
|
| 474 |
del self.thermostat["ref_t"] |
|
| 475 |
del self.thermostat["tau_t"] |
|
| 476 |
|
|
| 477 | 38 |
elif ensemble == self.Ensemble.NVT: |
| 478 |
self.thermostat["gen_vel"] = "yes" |
|
| 479 |
self.thermostat["gen_temp"] = self.temperature |
|
| 480 |
self.thermostat["gen_seed"] = -1 |
|
| 481 |
self.barostat["pcoupl"] = self.Barostat.Off.value |
|
| 482 |
|
|
| 483 | 38 |
elif ensemble == self.Ensemble.NPT: |
| 484 | 38 |
self.thermostat["gen_vel"] = "no" |
| 485 | 38 |
self.barostat["pcoupl"] = barostat.value |
| 486 | 38 |
if barostat.value != self.Barostat.Off: |
| 487 | 38 |
self.barostat["pcoupltype"] = self.BoxScaling.Isotropic.value |
| 488 | 38 |
self.barostat["tau_p"] = 2.0 |
| 489 | 38 |
self.barostat["ref_p"] = self.pressure |
| 490 | 38 |
self.barostat["compressibility"] = 4.5e-5 |
| 491 |
|
|
| 492 | 38 |
@staticmethod
|
| 493 |
def _write_dict_to_mdp(f, dictionary): |
|
| 494 |
"""
|
|
| 495 |
Write dictionary to file, following GROMACS format.
|
|
| 496 |
|
|
| 497 |
Parameters
|
|
| 498 |
----------
|
|
| 499 |
f : TextIO
|
|
| 500 |
File where the dictionary should be written.
|
|
| 501 |
dictionary : dict
|
|
| 502 |
Dictionary of values.
|
|
| 503 |
"""
|
|
| 504 | 38 |
for key, val in dictionary.items(): |
| 505 | 38 |
if val is not None and not isinstance(val, list): |
| 506 | 38 |
f.write("{:25s} {:s}\n".format(key, "= " + str(val))) |
| 507 | 38 |
elif isinstance(val, list): |
| 508 | 38 |
f.write("{:25s} {:s}".format(key, "= ")) |
| 509 | 38 |
for i in val: |
| 510 | 38 |
f.write("{:s} ".format(str(i))) |
| 511 | 38 |
f.write("\n") |
| 512 |
|
|
| 513 | 38 |
def _write_input_file(self): |
| 514 |
"""
|
|
| 515 |
Write the input file specification to file.
|
|
| 516 |
"""
|
|
| 517 | 38 |
logger.debug("Writing {}".format(self.input)) |
| 518 | 38 |
with open(os.path.join(self.path, self.input), "w") as mdp: |
| 519 | 38 |
mdp.write("{:25s} {:s}\n".format("title", "= " + self.title)) |
| 520 |
|
|
| 521 | 38 |
mdp.write("; Run control\n") |
| 522 | 38 |
self._write_dict_to_mdp(mdp, self.control) |
| 523 |
|
|
| 524 | 38 |
mdp.write("; Nonbonded options\n") |
| 525 | 38 |
self._write_dict_to_mdp(mdp, self.nb_method) |
| 526 |
|
|
| 527 | 38 |
mdp.write("; Bond constraints\n") |
| 528 | 38 |
if self.constraints["constraint-algorithm"].lower() == "shake": |
| 529 |
self._write_dict_to_mdp( |
|
| 530 |
mdp, |
|
| 531 |
get_dict_without_keys( |
|
| 532 |
self.constraints, "lincs_iter", "lincs_order" |
|
| 533 |
),
|
|
| 534 |
)
|
|
| 535 |
else: |
|
| 536 | 38 |
self._write_dict_to_mdp(mdp, self.constraints) |
| 537 |
|
|
| 538 | 38 |
if self.thermostat: |
| 539 | 38 |
mdp.write("; Temperature coupling\n") |
| 540 |
|
|
| 541 |
# Check if users specify different temperature groups
|
|
| 542 | 38 |
if self.tc_groups: |
| 543 | 38 |
tau_t = self.thermostat["tau_t"] |
| 544 | 38 |
self.thermostat["tc-grps"] = self.tc_groups |
| 545 | 38 |
self.thermostat["tau_t"] = [tau_t] * len(self.tc_groups) |
| 546 | 38 |
self.thermostat["ref_t"] = [self.temperature] * len(self.tc_groups) |
| 547 |
|
|
| 548 | 38 |
self._write_dict_to_mdp(mdp, self.thermostat) |
| 549 |
|
|
| 550 | 38 |
if self.barostat: |
| 551 | 38 |
mdp.write("; Pressure coupling\n") |
| 552 | 38 |
self._write_dict_to_mdp(mdp, self.barostat) |
| 553 |
|
|
| 554 | 38 |
def run(self, run_grompp=True, overwrite=False, fail_ok=False): |
| 555 |
"""
|
|
| 556 |
Method to run Molecular Dynamics simulation with GROMACS.
|
|
| 557 |
|
|
| 558 |
Parameters
|
|
| 559 |
----------
|
|
| 560 |
run_grompp: bool, optional, default=True
|
|
| 561 |
Run GROMPP to generate ``.tpr`` file before running MDRUN
|
|
| 562 |
overwrite: bool, optional, default=False
|
|
| 563 |
Whether to overwrite simulation files.
|
|
| 564 |
fail_ok: bool, optional, default=False
|
|
| 565 |
Whether a failing simulation should stop execution of ``pAPRika``.
|
|
| 566 |
"""
|
|
| 567 |
|
|
| 568 |
if overwrite or not self.check_complete(): |
|
| 569 |
# Check the type of simulation: Minimization, NVT or NPT
|
|
| 570 |
if self.control["integrator"] in [ |
|
| 571 |
self.Optimizer.SteepestDescent.value, |
|
| 572 |
self.Optimizer.ConjugateGradient.value, |
|
| 573 |
self.Optimizer.Broyden.value, |
|
| 574 |
]:
|
|
| 575 |
logger.info("Running Minimization at {}".format(self.path)) |
|
| 576 |
elif self.control["integrator"] in [ |
|
| 577 |
self.Integrator.LeapFrog.value, |
|
| 578 |
self.Integrator.VelocityVerlet.value, |
|
| 579 |
self.Integrator.VelocityVerletAveK.value, |
|
| 580 |
self.Integrator.LangevinDynamics.value, |
|
| 581 |
self.Integrator.BrownianDynamics.value, |
|
| 582 |
]:
|
|
| 583 |
if self.thermostat and self.barostat: |
|
| 584 |
logger.info("Running NPT MD at {}".format(self.path)) |
|
| 585 |
elif not self.barostat: |
|
| 586 |
logger.info("Running NVT MD at {}".format(self.path)) |
|
| 587 |
else: |
|
| 588 |
logger.info("Running NVE MD at {}".format(self.path)) |
|
| 589 |
|
|
| 590 |
# Set Plumed kernel library to path
|
|
| 591 |
self._set_plumed_kernel() |
|
| 592 |
|
|
| 593 |
# create executable list for GROMPP
|
|
| 594 |
# gmx grompp -f npt.mdp -c coordinates.gro -p topology.top -t checkpoint.cpt -o npt.tpr -n index.ndx
|
|
| 595 |
if run_grompp: |
|
| 596 |
|
|
| 597 |
# Clean previously generated files
|
|
| 598 |
for file in glob.glob(os.path.join(self.path, f"{self.prefix}*")): |
|
| 599 |
os.remove(file) |
|
| 600 |
|
|
| 601 |
# Write MDF input file
|
|
| 602 |
self._write_input_file() |
|
| 603 |
|
|
| 604 |
# GROMPP list
|
|
| 605 |
grompp_list = [self.executable, "grompp"] |
|
| 606 |
|
|
| 607 |
grompp_list += [ |
|
| 608 |
"-f", |
|
| 609 |
self.input, |
|
| 610 |
"-p", |
|
| 611 |
self.topology, |
|
| 612 |
"-c", |
|
| 613 |
self.coordinates, |
|
| 614 |
"-o", |
|
| 615 |
self.tpr, |
|
| 616 |
"-po", |
|
| 617 |
self.output, |
|
| 618 |
"-maxwarn", |
|
| 619 |
str(self.grompp_maxwarn), |
|
| 620 |
]
|
|
| 621 |
if self.checkpoint: |
|
| 622 |
grompp_list += ["-t", self.checkpoint] |
|
| 623 |
|
|
| 624 |
if self.index_file: |
|
| 625 |
grompp_list += ["-n", self.index_file] |
|
| 626 |
|
|
| 627 |
# Run GROMPP
|
|
| 628 |
grompp_output = sp.Popen( |
|
| 629 |
grompp_list, |
|
| 630 |
cwd=self.path, |
|
| 631 |
stdout=sp.PIPE, |
|
| 632 |
stderr=sp.PIPE, |
|
| 633 |
env=os.environ, |
|
| 634 |
)
|
|
| 635 |
grompp_stdout = grompp_output.stdout.read().splitlines() |
|
| 636 |
grompp_stderr = grompp_output.stderr.read().splitlines() |
|
| 637 |
|
|
| 638 |
# Report any stdout/stderr which are output from execution
|
|
| 639 |
if grompp_stdout: |
|
| 640 |
logger.info("STDOUT received from GROMACS execution") |
|
| 641 |
for line in grompp_stdout: |
|
| 642 |
logger.info(line) |
|
| 643 |
|
|
| 644 |
# Not sure how to do this more efficiently/elegantly, "subprocess" seems to treat everything
|
|
| 645 |
# Gromacs spits out from "grompp" as an error.
|
|
| 646 |
if grompp_stderr and any( |
|
| 647 |
["Error" in line.decode("utf-8").strip() for line in grompp_stderr] |
|
| 648 |
):
|
|
| 649 |
logger.info("STDERR received from GROMACS execution") |
|
| 650 |
for line in grompp_stderr: |
|
| 651 |
logger.error(line) |
|
| 652 |
|
|
| 653 |
# create executable list for MDRUN
|
|
| 654 |
# gmx_mpi mdrun -v -deffnm npt -nt 6 -gpu_id 0 -plumed plumed.dat
|
|
| 655 |
mdrun_list = [] |
|
| 656 |
|
|
| 657 |
# Add any user specified command
|
|
| 658 |
if self.custom_mdrun_command is not None: |
|
| 659 |
if self.executable not in self.custom_mdrun_command: |
|
| 660 |
mdrun_list += [self.executable] |
|
| 661 |
|
|
| 662 |
if "mdrun" not in self.custom_mdrun_command: |
|
| 663 |
mdrun_list += ["mdrun"] |
|
| 664 |
|
|
| 665 |
mdrun_list += self.custom_mdrun_command.split() |
|
| 666 |
|
|
| 667 |
# Output prefix
|
|
| 668 |
if "-deffnm" not in self.custom_mdrun_command: |
|
| 669 |
mdrun_list += ["-deffnm", self.prefix] |
|
| 670 |
|
|
| 671 |
# Add number of threads if not already specified in custom
|
|
| 672 |
if not any( |
|
| 673 |
[
|
|
| 674 |
cpu in self.custom_mdrun_command |
|
| 675 |
for cpu in ["-nt", "-ntomp", "-ntmpi", "-ntomp_pme"] |
|
| 676 |
]
|
|
| 677 |
):
|
|
| 678 |
mdrun_list += [ |
|
| 679 |
"-ntomp" if "mpi" in self.executable else "-nt", |
|
| 680 |
str(self.n_threads), |
|
| 681 |
]
|
|
| 682 |
|
|
| 683 |
# Add gpu id if not already specified in custom
|
|
| 684 |
if ( |
|
| 685 |
self.gpu_devices is not None |
|
| 686 |
and "-gpu_id" not in self.custom_mdrun_command |
|
| 687 |
):
|
|
| 688 |
mdrun_list += ["-gpu_id", str(self.gpu_devices)] |
|
| 689 |
|
|
| 690 |
# Add plumed file if not already specified in custom
|
|
| 691 |
if self.plumed_file and "-plumed" not in self.custom_mdrun_command: |
|
| 692 |
mdrun_list += ["-plumed", self.plumed_file] |
|
| 693 |
|
|
| 694 |
else: |
|
| 695 |
mdrun_list += [self.executable, "mdrun", "-deffnm", self.prefix] |
|
| 696 |
|
|
| 697 |
# Add number of threads
|
|
| 698 |
mdrun_list += [ |
|
| 699 |
"-ntomp" if "mpi" in self.executable else "-nt", |
|
| 700 |
str(self.n_threads), |
|
| 701 |
]
|
|
| 702 |
|
|
| 703 |
# Add gpu id
|
|
| 704 |
if self.gpu_devices is not None: |
|
| 705 |
mdrun_list += ["-gpu_id", str(self.gpu_devices)] |
|
| 706 |
|
|
| 707 |
# Add plumed file
|
|
| 708 |
if self.plumed_file is not None: |
|
| 709 |
mdrun_list += ["-plumed", self.plumed_file] |
|
| 710 |
|
|
| 711 |
# Run MDRUN
|
|
| 712 |
mdrun_output = sp.Popen( |
|
| 713 |
mdrun_list, |
|
| 714 |
cwd=self.path, |
|
| 715 |
stdout=sp.PIPE, |
|
| 716 |
stderr=sp.PIPE, |
|
| 717 |
env=os.environ, |
|
| 718 |
)
|
|
| 719 |
mdrun_out = mdrun_output.stdout.read().splitlines() |
|
| 720 |
mdrun_err = mdrun_output.stderr.read().splitlines() |
|
| 721 |
|
|
| 722 |
# Report any stdout/stderr which are output from execution
|
|
| 723 |
if mdrun_out: |
|
| 724 |
logger.info("STDOUT received from MDRUN execution") |
|
| 725 |
for line in mdrun_out: |
|
| 726 |
logger.info(line) |
|
| 727 |
|
|
| 728 |
# Same reasoning as before for "grompp".
|
|
| 729 |
if mdrun_err and any( |
|
| 730 |
["Error" in line.decode("utf-8").strip() for line in mdrun_err] |
|
| 731 |
):
|
|
| 732 |
logger.info("STDERR received from MDRUN execution") |
|
| 733 |
for line in mdrun_err: |
|
| 734 |
logger.error(line) |
|
| 735 |
|
|
| 736 |
# Check completion status
|
|
| 737 |
if ( |
|
| 738 |
self.control["integrator"] |
|
| 739 |
in [ |
|
| 740 |
self.Optimizer.SteepestDescent.value, |
|
| 741 |
self.Optimizer.ConjugateGradient.value, |
|
| 742 |
self.Optimizer.Broyden.value, |
|
| 743 |
]
|
|
| 744 |
and self.check_complete() |
|
| 745 |
):
|
|
| 746 |
logger.info("Minimization completed...") |
|
| 747 |
elif self.check_complete(): |
|
| 748 |
logger.info("Simulation completed...") |
|
| 749 |
else: |
|
| 750 |
logger.info( |
|
| 751 |
"Simulation did not complete when executing the following ...."
|
|
| 752 |
)
|
|
| 753 |
logger.info(" ".join(mdrun_list)) |
|
| 754 |
if not fail_ok: |
|
| 755 |
raise Exception( |
|
| 756 |
"Exiting due to failed simulation! Check logging info."
|
|
| 757 |
)
|
|
| 758 |
|
|
| 759 |
else: |
|
| 760 |
logger.info( |
|
| 761 |
"Completed output detected ... Skipping. Use: run(overwrite=True) to overwrite"
|
|
| 762 |
)
|
|
| 763 |
|
|
| 764 | 38 |
def check_complete(self, alternate_file=None): |
| 765 |
"""
|
|
| 766 |
Check for the string "step N" in ``self.output`` file. If "step N" is found, then
|
|
| 767 |
the simulation completed.
|
|
| 768 |
|
|
| 769 |
Parameters
|
|
| 770 |
----------
|
|
| 771 |
alternate_file : os.PathLike, optional, default=None
|
|
| 772 |
If present, check for "step N" in this file rather than ``self.output``.
|
|
| 773 |
Default: None
|
|
| 774 |
|
|
| 775 |
Returns
|
|
| 776 |
-------
|
|
| 777 |
complete : bool
|
|
| 778 |
True if "step N" is found in file. False, otherwise.
|
|
| 779 |
"""
|
|
| 780 |
# Assume not completed
|
|
| 781 |
complete = False |
|
| 782 |
|
|
| 783 |
if alternate_file: |
|
| 784 |
output_file = alternate_file |
|
| 785 |
else: |
|
| 786 |
output_file = os.path.join(self.path, self.logfile) |
|
| 787 |
|
|
| 788 |
if os.path.isfile(output_file): |
|
| 789 |
with open(output_file, "r") as f: |
|
| 790 |
strings = f.read() |
|
| 791 |
if ( |
|
| 792 |
f" step {self.control['nsteps']} " in strings |
|
| 793 |
or "Finished mdrun" in strings |
|
| 794 |
):
|
|
| 795 |
complete = True |
|
| 796 |
|
|
| 797 |
if complete: |
|
| 798 |
logger.debug("{} has TIMINGS".format(output_file)) |
|
| 799 |
else: |
|
| 800 |
logger.debug("{} does not have TIMINGS".format(output_file)) |
|
| 801 |
|
|
| 802 |
return complete |
Read our documentation on viewing source code .