1 2
import logging
2 2
import os
3 2
import subprocess as sp
4 2
from collections import OrderedDict
5 2
from sys import platform
6

7 2
logger = logging.getLogger(__name__)
8

9

10 2
class Simulation(object):
11
    """
12
    A wrapper that can be used to set AMBER simulation parameters.
13
    """
14

15 2
    @property
16
    def path(self):
17
        """os.PathLike: The path for the creation and execution of files for this protocol."""
18 2
        return self._path
19

20 2
    @path.setter
21
    def path(self, value):
22 2
        self._path = value
23

24 2
    @property
25
    def executable(self):
26
        """The AMBER executable that will be used.
27

28
        .. note ::
29
            This could be made safer by making an ``ENUM`` of ``sander``, ``pmemd`` and ``pmemd.cuda``.
30
        """
31 2
        return self._executable
32

33 2
    @executable.setter
34
    def executable(self, value):
35 2
        self._executable = value
36

37 2
    @property
38
    def gpu_devices(self):
39
        """A wrapper around the environmental variable ``CUDA_VISIBLE_DEVICES``."""
40 2
        return self._gpu_devices
41

42 2
    @gpu_devices.setter
43
    def gpu_devices(self, value):
44 0
        self._gpu_devices = value
45

46 2
    @property
47
    def phase(self):
48
        """str: Which phase of the calculation to simulate.
49

50
        .. note ::
51
            This could probably be combined with the ``window`` attribute for simplicity.
52
        """
53 0
        return self._phase
54

55 2
    @phase.setter
56
    def phase(self, value):
57 0
        self._phase = value
58

59 2
    @property
60
    def window(self):
61
        """str: Which window of the calculation to simulate."""
62 0
        return self._window
63

64 2
    @window.setter
65
    def window(self, value):
66 0
        self._window = value
67

68 2
    @property
69
    def topology(self):
70
        """os.PathLike: The topology file to simulate."""
71 2
        return self._topology
72

73 2
    @topology.setter
74
    def topology(self, value):
75 2
        self._topology = value
76

77 2
    @property
78
    def restraint_file(self):
79
        """os.PathLike: The file containing NMR-style restraints for AMBER.
80

81
        .. note ::
82
            When running AMBER simulations, you can only use either an AMBER NMR-style
83
            restraints or a Plumed-style restraints and not both.
84
        """
85 2
        return self._restraint_file
86

87 2
    @restraint_file.setter
88
    def restraint_file(self, value):
89 2
        self._restraint_file = value
90

91 2
    @property
92 2
    def plumed_file(self) -> str:
93
        """os.PathLike: The name of the Plumed-style restraints file for AMBER.
94

95
        .. note ::
96
            When running AMBER simulations, you can only use either an AMBER NMR-style
97
            restraints or a Plumed-style restraints and not both.
98
        """
99 2
        return self._plumed_file
100

101 2
    @plumed_file.setter
102 2
    def plumed_file(self, value: str):
103 0
        self._plumed_file = value
104

105 2
    @property
106 2
    def plumed_kernel_library(self) -> str:
107
        """os.PathLike: The file path of the Plumed kernel library if it is different from the default. The default
108
        name is `libplumedKernel` with its extension determined automatically depending on the operating system and is
109
        located in `os.environ['CONDA_PREFIX']/lib`. This variable is useful for users who did not install or want to
110
        use Plumed from the conda repository."""
111 2
        return self._plumed_kernel_library
112

113 2
    @plumed_kernel_library.setter
114 2
    def plumed_kernel_library(self, value: str):
115 2
        self._plumed_kernel_library = value
116

117 2
    @property
118 2
    def converged(self) -> bool:
119
        """bool: Whether the simulation is converged.
120

121
        .. warning ::
122
            This is just a placeholder and is not implemented yet.
123
        """
124 0
        return self._converged
125

126 2
    @converged.setter
127 2
    def converged(self, value: bool):
128 2
        self._converged = value
129

130 2
    @property
131
    def prefix(self):
132
        """
133
        The prefix for file names generated from this simulation.
134
        """
135 0
        return self._prefix
136

137 2
    @prefix.setter
138
    def prefix(self, new_prefix):
139 2
        self._prefix = new_prefix
140 2
        self.input = new_prefix + ".in"
141 2
        self.inpcrd = new_prefix + ".inpcrd"
142 2
        self.ref = new_prefix + ".inpcrd"
143 2
        self.output = new_prefix + ".out"
144 2
        self.restart = new_prefix + ".rst7"
145 2
        self.mdinfo = new_prefix + ".mdinfo"
146 2
        self.mdcrd = new_prefix + ".nc"
147 2
        self.mden = new_prefix + ".mden"
148

149 2
    @property
150
    def cntrl(self):
151
        """
152
        An ordered dictionary of simulation "control" namelist parameters. The main puprose of this class attribute is
153
        to make it easy to override certain simulation parameters, such as positional restraints on dummy atoms, and
154
        the inclusion of exclusion of the NMR-style APR restraints.
155

156
        As of AMBER18, these are described, in part, in chapter 18 on ``pmemd``.
157

158
        .. note ::
159
            I can't recall why we wanted an ``OrderedDict`` here.
160

161
        .. note ::
162
            **This is fragile** and this could be hardened by making these ``ENUM``s and doing much more type-checking.
163

164
            These must be valid AMBER keywords or else the simulation will crash. The default keywords that are set when
165
            this class is initialized are partially overridden by the specific "config" functions detailed below.
166

167
        The default dictionary keys and values are as follows:
168

169
            - ``imin``       : 0
170
            - ``ntx``        : 1
171
            - ``irest``      : 0
172
            - ``maxcyc``     : 0
173
            - ``ncyc``       : 0
174
            - ``dt``         : 0.002
175
            - ``nstlim``     : 5000
176
            - ``ntpr``       : 500
177
            - ``ntwe``       : 500
178
            - ``ntwr``       : 5000
179
            - ``ntwx``       : 500
180
            - ``ntxo``       : 1
181
            - ``ioutfm``     : 1
182
            - ``ntf``        : 2
183
            - ``ntc``        : 2
184
            - ``cut``        : 8
185
            - ``igb``        : 0
186
            - ``tempi``      : 298.15
187
            - ``tempo``      : 298.15
188
            - ``ntt``        : 3
189
            - ``gamma_ln``   : 1.0
190
            - ``ig``         : -1
191
            - ``ntp``        : 1
192
            - ``barostat``   : 2
193
            - ``ntr``        : ``None``
194
            - ``restraint_wt``  : ``None``
195
            - ``restraintmask`` : ``None``
196
            - ``nmropt``     : 1
197
            - ``pencut``     : -1
198

199
        """
200 2
        return self._cntrl
201

202 2
    @cntrl.setter
203
    def cntrl(self, value):
204 0
        self._cntrl = value
205

206 2
    @property
207
    def ewald(self):
208
        """Additional Ewald summation settings."""
209 2
        return self._ewald
210

211 2
    @ewald.setter
212
    def ewald(self, value):
213 0
        self._ewald = value
214

215 2
    @property
216
    def wt(self):
217
        """ "Weight" assigned to various simulation parameters. Written to the "&wt" line in the input file."""
218 2
        return self._wt
219

220 2
    @wt.setter
221
    def wt(self, value):
222 0
        self._wt = value
223

224 2
    @property
225
    def group(self):
226
        """Group specification for restraints applied to multiple atoms.
227

228
        .. note::
229
            I don't recall ever having to this option and this seems distinct from applying restraints to a collection
230
            of atoms.
231
        """
232 2
        return self._group
233

234 2
    @group.setter
235
    def group(self, value):
236 0
        self._group = value
237

238 2
    def __init__(self):
239

240
        # I/O
241 2
        self._path = "."
242 2
        self._executable = "sander"
243 2
        self._gpu_devices = None
244 2
        self._phase = None
245 2
        self._window = None
246 2
        self._topology = "prmtop"
247 2
        self._restraint_file = None
248 2
        self._plumed_file = None
249 2
        self._plumed_kernel_library = None
250 2
        self.title = "PBC MD Simulation"
251 2
        self.converged = False
252

253
        # File names
254 2
        self._prefix = "md"
255 2
        self.input = self._prefix + ".in"
256 2
        self.inpcrd = self._prefix + ".inpcrd"
257 2
        self.ref = self._prefix + ".inpcrd"
258 2
        self.output = self._prefix + ".out"
259 2
        self.restart = self._prefix + ".rst7"
260 2
        self.mdinfo = self._prefix + ".mdinfo"
261 2
        self.mdcrd = self._prefix + ".nc"
262 2
        self.mden = self._prefix + ".mden"
263

264
        # Input file cntrl settings (Default = NTP)
265 2
        self._cntrl = OrderedDict()
266 2
        self._cntrl["imin"] = 0
267 2
        self._cntrl["ntx"] = 1
268 2
        self._cntrl["irest"] = 0
269 2
        self._cntrl["maxcyc"] = 0
270 2
        self._cntrl["ncyc"] = 0
271 2
        self._cntrl["dt"] = 0.002
272 2
        self._cntrl["nstlim"] = 5000
273 2
        self._cntrl["ntpr"] = 500
274 2
        self._cntrl["ntwe"] = 500
275 2
        self._cntrl["ntwr"] = 5000
276 2
        self._cntrl["ntwx"] = 500
277 2
        self._cntrl["ntxo"] = 1
278 2
        self._cntrl["ioutfm"] = 1
279 2
        self._cntrl["ntf"] = 2
280 2
        self._cntrl["ntc"] = 2
281 2
        self._cntrl["cut"] = 8.0
282 2
        self._cntrl["igb"] = 0
283 2
        self._cntrl["tempi"] = 298.15
284 2
        self._cntrl["temp0"] = 298.15
285 2
        self._cntrl["ntt"] = 3
286 2
        self._cntrl["gamma_ln"] = 1.0
287 2
        self._cntrl["ig"] = -1
288 2
        self._cntrl["ntp"] = 1
289 2
        self._cntrl["barostat"] = 2
290 2
        self._cntrl["ntr"] = None
291 2
        self._cntrl["restraint_wt"] = None
292 2
        self._cntrl["restraintmask"] = None
293 2
        self._cntrl["nmropt"] = 1
294 2
        self._cntrl["pencut"] = -1
295

296
        # Other input file sections
297 2
        self._ewald = None
298 2
        self._wt = None  # or []
299 2
        self._group = None  # or []
300

301 2
    def _config_min(self):
302
        """
303
        Configure input settings for minimization (without periodic boundary conditions).
304

305
        """
306 2
        self.cntrl["imin"] = 1
307 2
        self.cntrl["ntx"] = 1
308 2
        self.cntrl["irest"] = 0
309 2
        self.cntrl["maxcyc"] = 5000
310 2
        self.cntrl["ncyc"] = 1000
311 2
        self.cntrl["dt"] = 0.0
312 2
        self.cntrl["nstlim"] = 0
313 2
        self.cntrl["ntpr"] = 100
314 2
        self.cntrl["ntwr"] = 5000
315 2
        self.cntrl["ntwx"] = 0
316 2
        self.cntrl["ntwe"] = 0
317 2
        self.cntrl["ntxo"] = 1
318 2
        self.cntrl["ntf"] = 1
319 2
        self.cntrl["ntc"] = 1
320 2
        self.cntrl["ntt"] = 0
321 2
        self.cntrl["gamma_ln"] = 0.0
322 2
        self.cntrl["ig"] = 0
323 2
        self.cntrl["ntp"] = 0
324 2
        self.cntrl["barostat"] = 0
325 2
        self.mdcrd = None
326 2
        self.mden = None
327

328 2
    def config_pbc_min(self):
329
        """
330
        Configure input settings for minimization in periodic boundary conditions.
331
        """
332 0
        self._config_min()
333 0
        self.title = "PBC Minimization"
334 0
        self.cntrl["cut"] = 8.0
335 0
        self.cntrl["igb"] = 0
336

337 2
    def config_gb_min(self):
338
        """
339
        Configure input settings for minimization in continuum solvent.
340
        """
341

342 2
        self._config_min()
343 2
        self.title = "GB Minimization"
344 2
        self.cntrl["cut"] = 999.0
345 2
        self.cntrl["igb"] = 1
346

347 2
    def _config_md(self):
348
        """
349
        Configure input settings for MD.
350
        """
351 2
        self.cntrl["imin"] = 0
352 2
        self.cntrl["ntx"] = 1
353 2
        self.cntrl["irest"] = 0
354 2
        self.cntrl["maxcyc"] = 0
355 2
        self.cntrl["ncyc"] = 0
356 2
        self.cntrl["dt"] = 0.002
357 2
        self.cntrl["nstlim"] = 5000
358 2
        self.cntrl["ntpr"] = 500
359 2
        self.cntrl["ntwe"] = 500
360 2
        self.cntrl["ntwr"] = 5000
361 2
        self.cntrl["ntwx"] = 500
362 2
        self.cntrl["ntxo"] = 1
363 2
        self.cntrl["ioutfm"] = 1
364 2
        self.cntrl["ntf"] = 2
365 2
        self.cntrl["ntc"] = 2
366 2
        self.cntrl["ntt"] = 3
367 2
        self.cntrl["gamma_ln"] = 1.0
368 2
        self.cntrl["ig"] = -1
369

370 2
    def config_gb_md(self):
371
        """
372
        Configure input settings for MD in default GB.
373
        """
374

375 2
        self._config_md()
376 2
        self.title = "GB MD Simulation"
377 2
        self.cntrl["cut"] = 999.0
378 2
        self.cntrl["igb"] = 1
379 2
        self.cntrl["ntp"] = 0
380 2
        self.cntrl["barostat"] = 0
381

382 2
    def config_pbc_md(self):
383
        """
384
        Configure input settings for default NTP.
385
        """
386

387 0
        self._config_md()
388 0
        self.title = "PBC MD Simulation"
389 0
        self.cntrl["cut"] = 8.0
390 0
        self.cntrl["igb"] = 0
391 0
        self.cntrl["iwrap"] = 1
392 0
        self.cntrl["ntp"] = 1
393 0
        self.cntrl["barostat"] = 2
394

395 2
    def _write_dict_to_mdin(self, f, dictionary):
396
        """
397
        Write dictionary to file, following AMBER format.
398

399
        Parameters
400
        ----------
401
        f : os.PathLike
402
            File where the dictionary should be written
403
        dictionary : dict
404
            Dictionary of values
405

406
        """
407

408 2
        for key, val in dictionary.items():
409 2
            if val is not None:
410 2
                f.write("  {:15s} {:s},\n".format(key + " =", str(val)))
411

412
        # Write PLUMED option if preferred over Amber NMR restraints
413 2
        if self.plumed_file:
414 0
            f.write("  {:15s} {:s},\n".format("plumed = ", str(1)))
415 0
            f.write("  {:15s} {:s},\n".format("plumedfile =", f"'{self.plumed_file}'"))
416

417 2
        f.write(" /\n")
418

419 2
    def _amber_write_input_file(self):
420
        """
421
        Write the input file specification to file.
422

423
        """
424 2
        logger.debug("Writing {}".format(self.input))
425 2
        with open(os.path.join(self.path, self.input), "w") as f:
426 2
            f.write("{}\n".format(self.title))
427 2
            f.write(" &cntrl\n")
428 2
            self._write_dict_to_mdin(f, self.cntrl)
429

430 2
            if self.ewald is not None:
431 0
                f.write(" &ewald\n")
432 0
                self._write_dict_to_mdin(f, self.ewald)
433

434 2
            if self.cntrl["nmropt"] == 1:
435 2
                if self.wt is not None:
436 0
                    for line in self.wt:
437 0
                        f.write(" " + line + "\n")
438 2
                f.write(" &wt type = 'END', /\n")
439

440
                # Specify Amber NMR file
441 2
                if self.restraint_file is not None:
442 0
                    f.write("DISANG = {}\n".format(self.restraint_file))
443 0
                    f.write("LISTOUT = POUT\n\n")
444

445 2
            if self.group is not None:
446 0
                f.write("{:s}".format(self.group))
447

448 2
    def run(self, soft_minimize=False, overwrite=False, fail_ok=False):
449
        """
450
        Run minimization.
451

452
        Parameters
453
        ----------
454
        soft_minimize: bool
455
            Whether to slowly turn on non-bonded interactions so that restraints get enforced first.
456
        overwrite: bool
457
            Whether to overwrite simulation files.
458
        fail_ok: bool
459
            Whether a failing simulation should stop execution of ``pAPRika``.
460
        """
461

462 2
        if overwrite or not self.has_timings():
463

464
            # These settings hardcoded at the moment ... possibly expose for
465
            # editing in the future
466 2
            if soft_minimize:
467
                # Set a burn in value that is 25% of the way between ncyc and
468
                # maxcyc
469 0
                ncyc = self.cntrl["ncyc"]
470 0
                maxcyc = self.cntrl["maxcyc"]
471 0
                burn_in = int(float(ncyc) + 0.20 * (float(maxcyc) - float(ncyc)))
472
                # If the burn_in value is nuts, then just set it to zero
473 0
                if burn_in < 0 or burn_in >= maxcyc:
474 0
                    burn_in = 0
475
                # Set an end_soft value that is 75% of way between ncyc and
476
                # maxcyc
477 0
                end_soft = int(float(ncyc) + 0.60 * (float(maxcyc) - float(ncyc)))
478 0
                self.wt = [
479
                    "&wt type = 'NB', istep1=0, istep2={:.0f}, value1 = 0.0, value2=0.0, IINC=50, /".format(
480
                        burn_in
481
                    ),
482
                    "&wt type = 'NB', istep1={:.0f}, istep2={:.0f}, value1 = 0.0, value2=1.0, IINC=50, /".format(
483
                        burn_in, end_soft
484
                    ),
485
                ]
486

487
            # Check restraints file
488 2
            if self.restraint_file and self.plumed_file:
489 0
                raise Exception(
490
                    "Cannot use both NMR-style and Plumed-style restraints at the same time."
491
                )
492

493
            # _amber_write_input_file(self.path+'/'+self.input, self.min, title='GB Minimization.')
494 2
            self._amber_write_input_file()
495

496 2
            if self.cntrl["imin"] == 1:
497 2
                logger.info("Running Minimization at {}".format(self.path))
498
            else:
499 2
                logger.info("Running MD at {}".format(self.path))
500

501
            # Create executable list for subprocess
502 2
            exec_list = self.executable.split() + ["-O", "-p", self.topology]
503 2
            if self.ref is not None:
504 2
                exec_list += ["-ref", self.ref]
505 2
            exec_list += [
506
                "-c",
507
                self.inpcrd,
508
                "-i",
509
                self.input,
510
                "-o",
511
                self.output,
512
                "-r",
513
                self.restart,
514
            ]
515 2
            if self.mdcrd is not None:
516 2
                exec_list += ["-x", self.mdcrd]
517 2
            if self.mdinfo is not None:
518 2
                exec_list += ["-inf", self.mdinfo]
519 2
            if self.mden is not None:
520 2
                exec_list += ["-e", self.mden]
521

522 2
            logger.debug("Exec line: " + " ".join(exec_list))
523

524
            # Add CUDA_VISIBLE_DEVICES variable to env dict
525 2
            if self.gpu_devices is not None:
526 0
                os.environ = dict(
527
                    os.environ, CUDA_VISIBLE_DEVICES=str(self.gpu_devices)
528
                )
529

530
            # Set the Plumed kernel library path
531 2
            if self.plumed_kernel_library is None:
532
                # using "startswith" as recommended from https://docs.python.org/3/library/sys.html#sys.platform
533 2
                if platform.startswith("linux"):
534 1
                    self.plumed_kernel_library = os.path.join(
535
                        os.environ["CONDA_PREFIX"], "lib", "libplumedKernel.so"
536
                    )
537 1
                elif platform.startswith("darwin"):
538 1
                    self.plumed_kernel_library = os.path.join(
539
                        os.environ["CONDA_PREFIX"], "lib", "libplumedKernel.dylib"
540
                    )
541
                else:
542 0
                    logger.warning(
543
                        f"Plumed is not supported with the '{platform}' platform."
544
                    )
545

546
            # Add Plumed kernel library to env dict
547 2
            if os.path.isfile(self.plumed_kernel_library):
548 2
                os.environ = dict(os.environ, PLUMED_KERNEL=self.plumed_kernel_library)
549
            else:
550 0
                logger.warning(
551
                    f"Plumed kernel library {self.plumed_kernel_library} does not exist, "
552
                    "PLUMED_KERNEL environment variable is not set."
553
                )
554

555
            # Execute
556 2
            amber_output = sp.Popen(
557
                exec_list,
558
                cwd=self.path,
559
                stdout=sp.PIPE,
560
                stderr=sp.PIPE,
561
                env=os.environ,
562
            )
563

564 2
            amber_stdout = amber_output.stdout.read().splitlines()
565 2
            amber_stderr = amber_output.stderr.read().splitlines()
566

567
            # Report any stdout/stderr which are output from execution
568 2
            if amber_stdout:
569 0
                logger.info("STDOUT received from AMBER execution")
570 0
                for line in amber_stdout:
571 0
                    logger.info(line)
572

573 2
            if amber_stderr:
574 0
                logger.info("STDERR received from AMBER execution")
575 0
                for line in amber_stderr:
576 0
                    logger.info(line)
577

578
            # Check completion status
579 2
            if self.cntrl["imin"] == 1 and self.has_timings():
580 2
                logger.info("Minimization completed...")
581 2
            elif self.has_timings():
582 2
                logger.info("MD completed ...")
583
            else:
584 0
                logger.info(
585
                    "Simulation did not complete when executing the following ...."
586
                )
587 0
                logger.info(" ".join(exec_list))
588 0
                if not fail_ok:
589 0
                    raise Exception(
590
                        "Exiting due to failed simulation! Check logging info."
591
                    )
592

593
        else:
594 0
            logger.info(
595
                "Completed output detected ... Skipping. Use: run(overwrite=True) to overwrite"
596
            )
597

598 2
    def has_timings(self, alternate_file=None):
599
        """
600
        Check for the string "TIMINGS" in ``self.output`` file. If "TIMINGS" is found, then the simulation completed
601
        successfully, as of AMBER18.
602

603
        Parameters
604
        ----------
605
        alternate_file : os.PathLike
606
            If present, check for "TIMINGS" in this file rather than ``self.output``. Default: None
607

608
        Returns
609
        -------
610
        timings : bool
611
            True if "TIMINGS" is found in file. False, otherwise.
612

613
        """
614

615
        # Assume not completed
616 2
        timings = False
617

618 2
        if alternate_file:
619 0
            output_file = alternate_file
620
        else:
621 2
            output_file = os.path.join(self.path, self.output)
622

623 2
        if os.path.isfile(output_file):
624 2
            with open(output_file, "r") as f:
625 2
                strings = f.read()
626 2
                if " TIMINGS" in strings:
627 2
                    timings = True
628 2
        if timings:
629 2
            logger.debug("{} has TIMINGS".format(output_file))
630
        else:
631 2
            logger.debug("{} does not have TIMINGS".format(output_file))
632

633 2
        return timings

Read our documentation on viewing source code .

Loading