1 26
import abc
2 26
import glob
3 26
import logging
4 26
import os
5 26
import subprocess as sp
6 26
from collections import OrderedDict
7 26
from enum import Enum
8

9 26
from paprika.utils import get_dict_without_keys
10

11 26
from .simulation import Simulation
12

13 26
logger = logging.getLogger(__name__)
14

15

16 26
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 26
    class Thermostat(Enum):
70
        """
71
        An enumeration of the different themostat implemented in GROMACS.
72
        """
73

74 26
        Off = "no"
75 26
        Berendsen = "berendsen"
76 26
        NoseHoover = "nose-hoover"
77 26
        Andersen1 = "andersen"
78 26
        Andersen2 = "andersen-massive"
79 26
        VelocityRescaling = "v-rescale"
80

81 26
    class Barostat(Enum):
82
        """
83
        An enumeration of the different barostat implemented in GROMACS.
84
        """
85

86 26
        Off = "no"
87 26
        Berendsen = "Berendsen"
88 26
        ParrinelloRahman = "Parrinello-Rahman"
89 26
        MMTK = "MTTK"
90

91 26
    class Integrator(Enum):
92
        """
93
        An enumeration of the different integrators implemented in GROMACS.
94
        """
95

96 26
        LeapFrog = "md"
97 26
        VelocityVerlet = "md-vv"
98 26
        VelocityVerletAveK = "md-vv-avek"
99 26
        LangevinDynamics = "sd"
100 26
        BrownianDynamics = "bd"
101

102 26
    class Optimizer(Enum):
103
        """
104
        An enumeration of the different minimization algorithm implemented in GROMACS.
105
        """
106

107 26
        SteepestDescent = "steep"
108 26
        ConjugateGradient = "cg"
109 26
        Broyden = "l-bfgs"
110

111 26
    class BoxScaling(Enum):
112
        """
113
        An enumeration of the different PBC scaling options when running constant pressure simulations in GROMACS.
114
        """
115

116 26
        Isotropic = "isotropic"
117 26
        Semiisotropic = "semiisotropic"
118 26
        Anisotropic = "anisotropic"
119 26
        SurfaceTension = "surface-tension"
120

121 26
    class Constraints(Enum):
122
        """
123
        An enumeration of the different bond constraint options in GROMACS.
124
        """
125

126 26
        Off = "none"
127 26
        HBonds = "h-bonds"
128 26
        AllBonds = "all-bonds"
129 26
        HAngles = "h-angles"
130 26
        AllAngles = "all-angles"
131

132 26
    @property
133 26
    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 0
        return self._index_file
137

138 26
    @index_file.setter
139 26
    def index_file(self, value: str):
140 0
        self._index_file = value
141

142 26
    @property
143 26
    def checkpoint(self) -> str:
144
        """os.PathLike: Checkpoint file (extension is ``.cpt``) for starting a simulation from a previous state."""
145 26
        return self._checkpoint
146

147 26
    @checkpoint.setter
148 26
    def checkpoint(self, value: str):
149 0
        self._checkpoint = value
150

151 26
    @property
152
    def control(self):
153
        """dict: Dictionary for the output control of the MD simulation (frequency of energy, trajectory etc)."""
154 26
        return self._control
155

156 26
    @control.setter
157
    def control(self, value):
158 0
        self._control = value
159

160 26
    @property
161
    def nb_method(self):
162
        """dict: Dictionary for the non-bonded method options (cutoffs and methods)."""
163 26
        return self._nb_method
164

165 26
    @nb_method.setter
166
    def nb_method(self, value):
167 0
        self._nb_method = value
168

169 26
    @property
170
    def constraints(self):
171
        """dict: Dictionary for the bond constraint options (LINCS or SHAKE)."""
172 26
        return self._constraints
173

174 26
    @constraints.setter
175
    def constraints(self, value):
176 0
        self._constraints = value
177

178 26
    @property
179 26
    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 26
        return self._tc_groups
192

193 26
    @tc_groups.setter
194 26
    def tc_groups(self, value: list):
195 26
        self._tc_groups = value
196

197 26
    @property
198
    def prefix(self):
199
        """str: The prefix for file names generated from this simulation."""
200 26
        return self._prefix
201

202 26
    @prefix.setter
203
    def prefix(self, new_prefix):
204 26
        self._prefix = new_prefix
205 26
        self.input = new_prefix + ".mdp"
206 26
        self.output = new_prefix + ".mdout"
207 26
        self.logfile = new_prefix + ".log"
208 26
        self.tpr = new_prefix + ".tpr"
209

210 26
    @property
211 26
    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 0
        return self._custom_mdrun_command
226

227 26
    @custom_mdrun_command.setter
228 26
    def custom_mdrun_command(self, value: str):
229 0
        self._custom_mdrun_command = value
230

231 26
    @property
232 26
    def grompp_maxwarn(self) -> int:
233
        """int: Maximum number of warnings for GROMPP to ignore. default=1."""
234 0
        return self._grompp_maxwarn
235

236 26
    @grompp_maxwarn.setter
237 26
    def grompp_maxwarn(self, value: int):
238 0
        self._grompp_maxwarn = value
239

240 26
    def __init__(self):
241

242 26
        super().__init__()
243

244
        # I/O
245 26
        self._index_file = None
246 26
        self._custom_mdrun_command = None
247 26
        self._tc_groups = None
248 26
        self._grompp_maxwarn = 1
249

250
        # File names
251 26
        self.input = self._prefix + ".mdp"
252 26
        self.output = self._prefix + ".mdout"
253 26
        self._checkpoint = None
254 26
        self.logfile = self._prefix + ".log"
255 26
        self.tpr = self._prefix + ".tpr"
256

257
        # Input file
258 26
        self._control = OrderedDict()
259 26
        self._control["nsteps"] = 5000
260 26
        self._control["nstxout"] = 500
261 26
        self._control["nstlog"] = 500
262 26
        self._control["nstenergy"] = 500
263 26
        self._control["nstcalcenergy"] = 500
264

265 26
        self._constraints = OrderedDict()
266 26
        self._constraints["constraint-algorithm"] = "lincs"
267 26
        self._constraints["constraints"] = self.Constraints.HBonds.value
268 26
        self._constraints["lincs_iter"] = 1
269 26
        self._constraints["lincs_order"] = 4
270

271 26
        self._nb_method = OrderedDict()
272 26
        self._nb_method["cutoff-scheme"] = "Verlet"
273 26
        self._nb_method["ns-type"] = "grid"
274 26
        self._nb_method["nstlist"] = 10
275 26
        self._nb_method["rlist"] = 0.9
276 26
        self._nb_method["rcoulomb"] = 0.9
277 26
        self._nb_method["rvdw"] = 0.9
278 26
        self._nb_method["coulombtype"] = "PME"
279 26
        self._nb_method["pme_order"] = 4
280 26
        self._nb_method["fourierspacing"] = 0.16
281 26
        self._nb_method["vdwtype"] = "Cut-off"
282 26
        self._nb_method["DispCorr"] = "EnerPres"
283 26
        self._nb_method["pbc"] = "xyz"
284

285 26
    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 26
        self.constraints["continuation"] = "no"
297 26
        self.control["integrator"] = optimizer.value
298 26
        self.control["emtol"] = 10.0
299 26
        self.control["emstep"] = 0.01
300 26
        self.control["nsteps"] = 5000
301

302 26
    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 26
        self.control["dt"] = 0.002
320 26
        self.control["integrator"] = integrator.value
321 26
        self.constraints["continuation"] = "yes"
322 26
        self.thermostat["tc-grps"] = "System"
323 26
        self.thermostat["ref_t"] = self.temperature
324

325 26
        if (
326
            integrator != self.Integrator.LangevinDynamics
327
            and integrator != self.Integrator.BrownianDynamics
328
        ):
329 26
            self.thermostat["tcoupl"] = thermostat.value
330 26
            self.thermostat["tau_t"] = 1.0
331
        else:
332 26
            self.thermostat["tau_t"] = 0.1
333

334 26
    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 0
        self.title = "Vacuum Minimization"
354

355 0
        self._config_min(optimizer)
356

357 0
        self.nb_method["pbc"] = "xyz"
358 0
        self.nb_method["ns_type"] = "grid"
359 0
        self.nb_method["nstlist"] = 10
360 0
        self.nb_method["rlist"] = 333.3
361 0
        self.nb_method["coulombtype"] = "Cut-off"
362 0
        self.nb_method["rcoulomb"] = 333.3
363 0
        self.nb_method["vdwtype"] = "Cut-off"
364 0
        self.nb_method["rvdw"] = 333.3
365 0
        self.nb_method["DispCorr"] = "no"
366

367 26
    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 26
        self.title = "Vacuum MD Simulation"
394

395 26
        self._config_md(integrator, thermostat)
396

397 26
        if self.checkpoint is None:
398 26
            self.constraints["continuation"] = "no"
399
        else:
400 0
            self.constraints["continuation"] = "yes"
401

402 26
        self.nb_method["pbc"] = "xyz"
403 26
        self.nb_method["ns_type"] = "grid"
404 26
        self.nb_method["nstlist"] = 10
405 26
        self.nb_method["rlist"] = 333.3
406 26
        self.nb_method["coulombtype"] = "Cut-off"
407 26
        self.nb_method["rcoulomb"] = 333.3
408 26
        self.nb_method["vdwtype"] = "Cut-off"
409 26
        self.nb_method["rvdw"] = 333.3
410 26
        self.nb_method["DispCorr"] = "no"
411

412 26
    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 26
        self.title = "PBC Minimization"
426

427 26
        self._config_min(optimizer)
428

429 26
        self.nb_method["nstlist"] = 10
430

431 26
    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 26
        self.title = f"{ensemble.value} MD Simulation"
462

463 26
        self._config_md(integrator, thermostat)
464

465 26
        if self.checkpoint is None:
466 26
            self.constraints["continuation"] = "no"
467
        else:
468 0
            self.constraints["continuation"] = "yes"
469

470 26
        if ensemble == self.Ensemble.NVE:
471 0
            self.thermostat["tcoupl"] = self.Thermostat.Off.value
472 0
            self.barostat["pcoupl"] = self.Barostat.Off.value
473 0
            del self.thermostat["tc-grps"]
474 0
            del self.thermostat["ref_t"]
475 0
            del self.thermostat["tau_t"]
476

477 26
        elif ensemble == self.Ensemble.NVT:
478 0
            self.thermostat["gen_vel"] = "yes"
479 0
            self.thermostat["gen_temp"] = self.temperature
480 0
            self.thermostat["gen_seed"] = -1
481 0
            self.barostat["pcoupl"] = self.Barostat.Off.value
482

483 26
        elif ensemble == self.Ensemble.NPT:
484 26
            self.thermostat["gen_vel"] = "no"
485 26
            self.barostat["pcoupl"] = barostat.value
486 26
            if barostat.value != self.Barostat.Off:
487 26
                self.barostat["pcoupltype"] = self.BoxScaling.Isotropic.value
488 26
                self.barostat["tau_p"] = 2.0
489 26
                self.barostat["ref_p"] = self.pressure
490 26
                self.barostat["compressibility"] = 4.5e-5
491

492 26
    @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 26
        for key, val in dictionary.items():
505 26
            if val is not None and not isinstance(val, list):
506 26
                f.write("{:25s} {:s}\n".format(key, "= " + str(val)))
507 26
            elif isinstance(val, list):
508 26
                f.write("{:25s} {:s}".format(key, "= "))
509 26
                for i in val:
510 26
                    f.write("{:s} ".format(str(i)))
511 26
                f.write("\n")
512

513 26
    def _write_input_file(self):
514
        """
515
        Write the input file specification to file.
516
        """
517 26
        logger.debug("Writing {}".format(self.input))
518 26
        with open(os.path.join(self.path, self.input), "w") as mdp:
519 26
            mdp.write("{:25s} {:s}\n".format("title", "= " + self.title))
520

521 26
            mdp.write("; Run control\n")
522 26
            self._write_dict_to_mdp(mdp, self.control)
523

524 26
            mdp.write("; Nonbonded options\n")
525 26
            self._write_dict_to_mdp(mdp, self.nb_method)
526

527 26
            mdp.write("; Bond constraints\n")
528 26
            if self.constraints["constraint-algorithm"].lower() == "shake":
529 0
                self._write_dict_to_mdp(
530
                    mdp,
531
                    get_dict_without_keys(
532
                        self.constraints, "lincs_iter", "lincs_order"
533
                    ),
534
                )
535
            else:
536 26
                self._write_dict_to_mdp(mdp, self.constraints)
537

538 26
            if self.thermostat:
539 26
                mdp.write("; Temperature coupling\n")
540

541
                # Check if users specify different temperature groups
542 26
                if self.tc_groups:
543 26
                    tau_t = self.thermostat["tau_t"]
544 26
                    self.thermostat["tc-grps"] = self.tc_groups
545 26
                    self.thermostat["tau_t"] = [tau_t] * len(self.tc_groups)
546 26
                    self.thermostat["ref_t"] = [self.temperature] * len(self.tc_groups)
547

548 26
                self._write_dict_to_mdp(mdp, self.thermostat)
549

550 26
            if self.barostat:
551 26
                mdp.write("; Pressure coupling\n")
552 26
                self._write_dict_to_mdp(mdp, self.barostat)
553

554 26
    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 0
        if overwrite or not self.check_complete():
569
            # Check the type of simulation: Minimization, NVT or NPT
570 0
            if self.control["integrator"] in [
571
                self.Optimizer.SteepestDescent.value,
572
                self.Optimizer.ConjugateGradient.value,
573
                self.Optimizer.Broyden.value,
574
            ]:
575 0
                logger.info("Running Minimization at {}".format(self.path))
576 0
            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 0
                if self.thermostat and self.barostat:
584 0
                    logger.info("Running NPT MD at {}".format(self.path))
585 0
                elif not self.barostat:
586 0
                    logger.info("Running NVT MD at {}".format(self.path))
587
                else:
588 0
                    logger.info("Running NVE MD at {}".format(self.path))
589

590
            # Set Plumed kernel library to path
591 0
            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 0
            if run_grompp:
596

597
                # Clean previously generated files
598 0
                for file in glob.glob(os.path.join(self.path, f"{self.prefix}*")):
599 0
                    os.remove(file)
600

601
                # Write MDF input file
602 0
                self._write_input_file()
603

604
                # GROMPP list
605 0
                grompp_list = [self.executable, "grompp"]
606

607 0
                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 0
                if self.checkpoint:
622 0
                    grompp_list += ["-t", self.checkpoint]
623

624 0
                if self.index_file:
625 0
                    grompp_list += ["-n", self.index_file]
626

627
                # Run GROMPP
628 0
                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 0
                grompp_stdout = grompp_output.stdout.read().splitlines()
636 0
                grompp_stderr = grompp_output.stderr.read().splitlines()
637

638
                # Report any stdout/stderr which are output from execution
639 0
                if grompp_stdout:
640 0
                    logger.info("STDOUT received from GROMACS execution")
641 0
                    for line in grompp_stdout:
642 0
                        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 0
                if grompp_stderr and any(
647
                    ["Error" in line.decode("utf-8").strip() for line in grompp_stderr]
648
                ):
649 0
                    logger.info("STDERR received from GROMACS execution")
650 0
                    for line in grompp_stderr:
651 0
                        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 0
            mdrun_list = []
656

657
            # Add any user specified command
658 0
            if self.custom_mdrun_command is not None:
659 0
                if self.executable not in self.custom_mdrun_command:
660 0
                    mdrun_list += [self.executable]
661

662 0
                if "mdrun" not in self.custom_mdrun_command:
663 0
                    mdrun_list += ["mdrun"]
664

665 0
                mdrun_list += self.custom_mdrun_command.split()
666

667
                # Output prefix
668 0
                if "-deffnm" not in self.custom_mdrun_command:
669 0
                    mdrun_list += ["-deffnm", self.prefix]
670

671
                # Add number of threads if not already specified in custom
672 0
                if not any(
673
                    [
674
                        cpu in self.custom_mdrun_command
675
                        for cpu in ["-nt", "-ntomp", "-ntmpi", "-ntomp_pme"]
676
                    ]
677
                ):
678 0
                    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 0
                if (
685
                    self.gpu_devices is not None
686
                    and "-gpu_id" not in self.custom_mdrun_command
687
                ):
688 0
                    mdrun_list += ["-gpu_id", str(self.gpu_devices)]
689

690
                # Add plumed file if not already specified in custom
691 0
                if self.plumed_file and "-plumed" not in self.custom_mdrun_command:
692 0
                    mdrun_list += ["-plumed", self.plumed_file]
693

694
            else:
695 0
                mdrun_list += [self.executable, "mdrun", "-deffnm", self.prefix]
696

697
                # Add number of threads
698 0
                mdrun_list += [
699
                    "-ntomp" if "mpi" in self.executable else "-nt",
700
                    str(self.n_threads),
701
                ]
702

703
                # Add gpu id
704 0
                if self.gpu_devices is not None:
705 0
                    mdrun_list += ["-gpu_id", str(self.gpu_devices)]
706

707
                # Add plumed file
708 0
                if self.plumed_file is not None:
709 0
                    mdrun_list += ["-plumed", self.plumed_file]
710

711
            # Run MDRUN
712 0
            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 0
            mdrun_out = mdrun_output.stdout.read().splitlines()
720 0
            mdrun_err = mdrun_output.stderr.read().splitlines()
721

722
            # Report any stdout/stderr which are output from execution
723 0
            if mdrun_out:
724 0
                logger.info("STDOUT received from MDRUN execution")
725 0
                for line in mdrun_out:
726 0
                    logger.info(line)
727

728
            # Same reasoning as before for "grompp".
729 0
            if mdrun_err and any(
730
                ["Error" in line.decode("utf-8").strip() for line in mdrun_err]
731
            ):
732 0
                logger.info("STDERR received from MDRUN execution")
733 0
                for line in mdrun_err:
734 0
                    logger.error(line)
735

736
            # Check completion status
737 0
            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 0
                logger.info("Minimization completed...")
747 0
            elif self.check_complete():
748 0
                logger.info("Simulation completed...")
749
            else:
750 0
                logger.info(
751
                    "Simulation did not complete when executing the following ...."
752
                )
753 0
                logger.info(" ".join(mdrun_list))
754 0
                if not fail_ok:
755 0
                    raise Exception(
756
                        "Exiting due to failed simulation! Check logging info."
757
                    )
758

759
        else:
760 0
            logger.info(
761
                "Completed output detected ... Skipping. Use: run(overwrite=True) to overwrite"
762
            )
763

764 26
    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 0
        complete = False
782

783 0
        if alternate_file:
784 0
            output_file = alternate_file
785
        else:
786 0
            output_file = os.path.join(self.path, self.logfile)
787

788 0
        if os.path.isfile(output_file):
789 0
            with open(output_file, "r") as f:
790 0
                strings = f.read()
791 0
                if (
792
                    f" step {self.control['nsteps']} " in strings
793
                    or "Finished mdrun" in strings
794
                ):
795 0
                    complete = True
796

797 0
        if complete:
798 0
            logger.debug("{} has TIMINGS".format(output_file))
799
        else:
800 0
            logger.debug("{} does not have TIMINGS".format(output_file))
801

802 0
        return complete

Read our documentation on viewing source code .

Loading