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

8 26
import numpy as np
9 26
import parmed as pmd
10

11 26
from paprika.utils import get_dict_without_keys
12

13 26
from .simulation import Simulation
14

15 26
logger = logging.getLogger(__name__)
16

17

18 26
class NAMD(Simulation, abc.ABC):
19
    """
20
    A wrapper that can be used to set NAMD simulation parameters.
21

22
    .. todo ::
23
        Add wrapper support for Drude polarizable FF, steered-MD and alchemical calculations
24

25
    Below is an example of the configuration file (``namd.conf``) generated by the wrapper. The class property
26
    associated
27
    with defining the configuration variables is shown in brackets.
28

29
    .. code ::
30

31
        ## NPT MD Simulation [self.title]
32
        # AMBER files [self.topology, self.coordinates]
33
        amber                          yes
34
        parmfile                       cb6-but-dum-sol.prmtop
35
        ambercoor                      cb6-but-dum-sol.rst7
36
        readexclusions                 yes
37
        scnb                           2.0
38
        # Previous state [self.checkpoint]
39
        bincoordinates                 equilibration.restart.coor
40
        binvelocities                  equilibration.restart.vel
41
        extendedSystem                 equilibration.restart.xsc
42
        # Output control [self.control]
43
        outputname                     equilibration
44
        restartFreq                    5000
45
        dcdFreq                        500
46
        xstFreq                        500
47
        outputEnergies                 500
48
        outputPressure                 500
49
        timestep                       2.0
50
        # Nonbonded options [self.nb_method]
51
        exclude                        scaled1-4
52
        1-4scaling                     0.833333
53
        cutoff                         9.0
54
        LJCorrection                   yes
55
        nonbondedFreq                  1
56
        fullElectFrequency             2
57
        switching                      off
58
        wrapAll                        on
59
        margin                         5
60
        PME                            yes
61
        PMEGridSpacing                 1.0
62
        stepspercycle                  10
63
        # constraints [self.constraints]
64
        rigidBonds                     all
65
        rigidIterations                100
66
        rigidTolerance                 1.0e-8
67
        usesettle                      on
68
        # Temperature coupling [self.thermostat]
69
        langevin                       on
70
        langevinDamping                1.0
71
        langevinTemp                   298.15
72
        langevinHydrogen               off
73
        # Pressure coupling [self.barostat]
74
        useGroupPressure               yes
75
        useFlexibleCell                no
76
        useConstantArea                no
77
        langevinPiston                 on
78
        langevinPistonTarget           1.01325
79
        langevinPistonPeriod           200.0
80
        langevinPistonDecay            20.0
81
        langevinPistonTemp             298.15
82
        # Collective variables [self.colvars_file]
83
        colvars                        on
84
        colvarsConfig                  colvars.tcl
85
        # Molecular dynamics [self.control]
86
        run                            1500000
87
    """
88

89 26
    class Thermostat(Enum):
90
        """
91
        An enumeration of the different themostat implemented in NAMD.
92
        """
93

94 26
        Off = "no thermostat"
95 26
        Langevin = "Langevin dynamics"
96 26
        tCouple = "Temperature Coupling"
97 26
        tRescale = "Temperature Rescaling"
98 26
        vRescale = "Velocity Rescaling"
99 26
        tReassign = "Temperature Reassignment"
100 26
        LoweAnderson = "Lowe-Andersen dynamics"
101

102 26
    class Barostat(Enum):
103
        """
104
        An enumeration of the different barostat implemented in NAMD.
105
        """
106

107 26
        Off = "no barostat"
108 26
        Berendsen = "Berendsen pressure bath"
109 26
        Langevin = "Nosé-Hoover Langevin piston"
110

111 26
    @property
112 26
    def colvars_file(self) -> os.PathLike:
113
        """os.PathLike: The name of the `Colvars` restraints file."""
114 26
        return self._colvars_file
115

116 26
    @colvars_file.setter
117 26
    def colvars_file(self, value: os.PathLike):
118 0
        self._colvars_file = value
119

120 26
    @property
121 26
    def tcl_forces(self) -> os.PathLike:
122
        """
123
        os.PathLike: The name of a `TclForces` scripts for applying a user defined force.
124
        https://www.ks.uiuc.edu/Research/namd/2.14/ug/node50.html
125

126
        .. note ::
127

128
            ``TclForces`` can by used to apply forces and/or on-the-fly analysis. An example use if this feature
129
            could be to generate a funnel potential. https://github.com/jeff231li/funnel_potential
130
        """
131 26
        return self._tcl_forces
132

133 26
    @tcl_forces.setter
134 26
    def tcl_forces(self, value: os.PathLike):
135 0
        self._tcl_forces = value
136

137 26
    @property
138 26
    def control(self) -> dict:
139
        """dict: Dictionary for the output control of the MD simulation (timestep, frequency of energy, trajectory
140
        etc)."""
141 26
        return self._control
142

143 26
    @control.setter
144 26
    def control(self, value: dict):
145 0
        self._control = value
146

147 26
    @property
148 26
    def nb_method(self) -> dict:
149
        """dict:Dictionary for the non-bonded method options (cutoffs and methods)."""
150 26
        return self._nb_method
151

152 26
    @nb_method.setter
153 26
    def nb_method(self, value: dict):
154 0
        self._nb_method = value
155

156 26
    @property
157 26
    def constraints(self) -> dict:
158
        """dict: Dictionary for the bond constraint options (rigid bonds and settle/shake)."""
159 26
        return self._constraints
160

161 26
    @constraints.setter
162
    def constraints(self, value):
163 0
        self._constraints = value
164

165 26
    @property
166 26
    def implicit(self) -> dict:
167
        """dict: Dictionary for implicit solvent options (Generalized Born Implicit Solvent model)."""
168 26
        return self._implicit
169

170 26
    @implicit.setter
171 26
    def implicit(self, value: dict):
172 0
        self._implicit = value
173

174 26
    @property
175 26
    def charmm_parameters(self) -> list:
176
        """list: A list CHARMM parameter file(s) if running simulations with the CHARMM force field."""
177 0
        return self._charmm_parameters
178

179 26
    @charmm_parameters.setter
180 26
    def charmm_parameters(self, value: list):
181 0
        self._charmm_parameters = value
182

183 26
    @property
184 26
    def prefix(self) -> str:
185
        """
186
        str: The prefix for file names generated from this simulation.
187
        """
188 26
        return self._prefix
189

190 26
    @prefix.setter
191 26
    def prefix(self, new_prefix: str):
192 26
        self._prefix = new_prefix
193 26
        self.input = new_prefix + ".conf"
194 26
        self.logfile = new_prefix + ".log"
195 26
        self.control["outputname"] = new_prefix
196

197 26
    @property
198 26
    def checkpoint(self) -> str:
199
        """
200
        str: Specify the prefix file name for the previous simulation (as checkpoint for continuing a simulation).
201
        This is needed to get the coordinates ``checkpoint.coor``, velocities ``checkpoint.vel`` and box information
202
        ``checkpoint.xsc``. Unlike AMBER and GROMACS, NAMD requires the restart files specified inside the
203
        configuration file instead of the command-line argument.
204
        """
205 26
        return self._checkpoint
206

207 26
    @checkpoint.setter
208 26
    def checkpoint(self, value: str):
209 26
        self._checkpoint = value
210

211 26
    @property
212
    def cell_basis_vectors(self):
213
        """
214
        dict: Dictionary for the cell basis vectors (``cellBasisVector1``, ``cellBasisVector2``, ``cellBasisVector3``
215
        and ``cellOrigin``). If not defined, the basis vectors and the origin can be calculated from the positions
216
        of the system from the coordinates files.
217
        """
218 26
        return self._cell_basis_vectors
219

220 26
    @cell_basis_vectors.setter
221
    def cell_basis_vectors(self, value):
222 0
        self._cell_basis_vectors = value
223

224 26
    @property
225 26
    def custom_run_commands(self) -> list:
226
        """
227
        list: A list of NAMD-supported Tcl commands if a custom run is preferred. The code snippet below shows an
228
        example use of this command for running a simulation with a `slow-heating` protocol.
229

230
        .. code-block :: python
231

232
            >>> simulation.custom_run_commands = [
233
            >>>     "for {set temp 50} {$temp <= 300} {incr temp 50} {",
234
            >>>     "langevinTemp $temp",
235
            >>>     "reinitvels $temp",
236
            >>>     "run 50000",
237
            >>>     "}",
238
            >>> ]
239
        """
240 26
        return self._custom_run_commands
241

242 26
    @custom_run_commands.setter
243 26
    def custom_run_commands(self, value: list):
244 0
        self._custom_run_commands = value
245

246 26
    @property
247 26
    def qmForces(self) -> dict:
248
        """dict: Dictionary for setting QM/MM simulation parameters in NAMD
249
        https://www.ks.uiuc.edu/Research/namd/2.14/ug/node84.html.
250

251
        .. note ::
252

253
            ``qmForces`` in NAMD can also be used to run simulations with Neural Network Potentials (NNP)
254
            https://github.com/RowleyGroup/NNP-MM.
255

256
        """
257 26
        return self._qmForces
258

259 26
    @qmForces.setter
260 26
    def qmForces(self, value: dict):
261 0
        self._qmForces = value
262

263 26
    def __init__(self):
264

265 26
        super().__init__()
266

267
        # I/O
268 26
        self._colvars_file = None
269

270
        # File names
271 26
        self.input = self._prefix + ".conf"
272 26
        self.logfile = self._prefix + ".log"
273 26
        self._checkpoint = None
274 26
        self._custom_run_commands = None
275 26
        self._tcl_forces = None
276 26
        self._charmm_parameters = []
277

278 26
        self._qmForces = {}
279

280 26
        self._constraints = OrderedDict()
281 26
        self._constraints["rigidBonds"] = "all"
282 26
        self._constraints["rigidTolerance"] = 1.0e-8
283 26
        self._constraints["rigidIterations"] = 100
284 26
        self._constraints["usesettle"] = "on"
285

286 26
        self._nb_method = OrderedDict()
287 26
        self._nb_method["exclude"] = "scaled1-4"
288 26
        self._nb_method["1-4scaling"] = 0.833333
289 26
        self._nb_method["cutoff"] = 9.0
290 26
        self._nb_method["LJCorrection"] = "yes"
291 26
        self._nb_method["nonbondedFreq"] = 1
292 26
        self._nb_method["fullElectFrequency"] = 2
293 26
        self._nb_method["switching"] = "off"
294 26
        self._nb_method["wrapAll"] = "on"
295 26
        self._nb_method["margin"] = 5.0
296

297
        # Placeholder dict
298 26
        self._control = {}
299 26
        self._implicit = OrderedDict()
300 26
        self._cell_basis_vectors = OrderedDict()
301

302 26
    def _config_min(self):
303
        """
304
        Configure input settings for minimization.
305
        """
306 0
        self.control["minimize"] = 5000
307

308 0
        if ".prmtop" in self.topology:
309 0
            self.control["readexclusions"] = "yes"
310 0
            self.control["scnb"] = 2.0
311

312 0
        elif ".psf" in self.topology:
313 0
            self.nb_method["1-4scaling"] = 1.0
314 0
            self.nb_method["cutoff"] = 12.0
315 0
            self.nb_method["switching"] = "on"
316 0
            self.nb_method["switchdist"] = 10.0
317 0
            self.nb_method["exclude"] = "scaled1-4"
318 0
            self.nb_method["LJCorrection"] = "off"
319

320 0
        self.control["restartFreq"] = 5000
321 0
        self.control["dcdFreq"] = 500
322 0
        self.control["xstFreq"] = 500
323 0
        self.control["outputEnergies"] = 500
324 0
        self.control["outputPressure"] = 500
325

326 0
        if "namd3" in self.executable:
327 0
            self.nb_method["CUDASOAintegrate"] = "off"
328

329 26
    def _config_md(self, thermostat):
330
        """
331
        Configure input settings for MD.
332

333
        Parameters
334
        ----------
335
        thermostat: :class:`NAMD.Thermostat`
336
            Option to choose one of five thermostat implemented in NAMD: (1) `Langevin`, (2) `tCouple`, (3) `tRescale`,
337
            (4) `tReassign`, (5) `Lowe-Anderson`.
338

339
        """
340 26
        self.control["restartFreq"] = 5000
341 26
        self.control["dcdFreq"] = 500
342 26
        self.control["xstFreq"] = 500
343 26
        self.control["outputEnergies"] = 500
344 26
        self.control["outputPressure"] = 500
345 26
        self.control["run"] = 5000
346 26
        self.control["timestep"] = 2.0
347

348 26
        if ".prmtop" in self.topology:
349 26
            self.control["readexclusions"] = "yes"
350 26
            self.control["scnb"] = 2.0
351

352 0
        elif ".psf" in self.topology:
353 0
            self.nb_method["1-4scaling"] = 1.0
354 0
            self.nb_method["cutoff"] = 12.0
355 0
            self.nb_method["switching"] = "on"
356 0
            self.nb_method["switchdist"] = 10.0
357 0
            self.nb_method["exclude"] = "scaled1-4"
358 0
            self.nb_method["LJCorrection"] = "no"
359

360 26
        if "namd3" in self.executable:
361 0
            self.nb_method["CUDASOAintegrate"] = "on"
362
        else:
363 26
            self.nb_method["margin"] = 1.0
364 26
            self.nb_method["stepspercycle"] = 10
365

366
        # Select Thermotat for the simulation
367 26
        if thermostat == self.Thermostat.Langevin:
368 26
            self.thermostat["langevin"] = "on"
369 26
            self.thermostat["langevinDamping"] = 1.0
370 26
            self.thermostat["langevinTemp"] = self.temperature
371 26
            self.thermostat["langevinHydrogen"] = "off"
372

373 26
        elif thermostat == self.Thermostat.tCouple:
374 0
            self.thermostat["tCouple"] = "on"
375 0
            self.thermostat["tCoupleTemp"] = self.temperature
376

377 26
        elif thermostat == self.Thermostat.tRescale:
378 0
            self.thermostat["rescaleFreq"] = 50
379 0
            self.thermostat["rescaleTemp"] = self.temperature
380

381 26
        elif thermostat == self.Thermostat.vRescale:
382 0
            self.thermostat["stochRescale"] = "on"
383 0
            self.thermostat["stochRescaleTemp"] = self.temperature
384 0
            self.thermostat["stochRescalePeriod"] = 1.0
385

386 26
        elif thermostat == self.Thermostat.tReassign:
387 0
            self.thermostat["reassignFreq"] = 50
388 0
            self.thermostat["reassignTemp"] = self.temperature
389

390 26
        elif thermostat == self.Thermostat.LoweAnderson:
391 26
            self.thermostat["loweAndersen"] = "on"
392 26
            self.thermostat["loweAndersenCutoff"] = 2.7
393 26
            self.thermostat["loweAndersenTemp"] = self.temperature
394 26
            self.thermostat["loweAndersenRate"] = 50
395

396 26
    def config_vac_min(self):
397
        """
398
        Configure a reasonable input setting for a MD run in vacuum. `Users can override the parameters set by this
399
        method.`
400
        """
401 0
        self.title = "Vacuum solvent Minimization"
402

403 0
        self._config_min()
404

405 0
        self.nb_method["cutoff"] = 999.0
406 0
        self.nb_method["LJCorrection"] = "no"
407 0
        self.nb_method["nonbondedFreq"] = 1
408 0
        self.nb_method["fullElectFrequency"] = 2
409 0
        self.nb_method["wrapAll"] = "off"
410 0
        self.nb_method["staticAtomAssignment"] = "on"
411

412 0
        if ".psf" in self.topology:
413 0
            self.nb_method["switching"] = "off"
414 0
            del self.nb_method["switchdist"]
415

416 26
    def config_vac_md(self, thermostat=Thermostat.Langevin):
417
        """
418
        Configure a reasonable input settings for MD in vacuum. `Users can override the parameters set by this method.`
419

420
        Parameters
421
        ----------
422
        thermostat: :class:`NAMD.Thermostat`, default=Thermostat.Langevin
423
            Option to choose one of six thermostat implemented in NAMD: **(1)** `Langevin`, **(2)** `tCouple`,
424
            **(3)** `tRescale`, **(4)** `vRescale`, **(5)** `tReassign`, **(6)** `LoweAnderson`.
425
        """
426 26
        self.title = "Vacuum MD Simulation"
427

428 26
        self._config_md(thermostat)
429

430 26
        self.control["timestep"] = 1.0
431

432 26
        self.implicit["GBIS"] = "off"
433

434 26
        self.nb_method["PME"] = "off"
435 26
        self.nb_method["cutoff"] = 999.0
436 26
        self.nb_method["LJCorrection"] = "no"
437 26
        self.nb_method["nonbondedFreq"] = 1
438 26
        self.nb_method["fullElectFrequency"] = 2
439 26
        self.nb_method["wrapAll"] = "off"
440 26
        self.nb_method["stepspercycle"] = 10
441 26
        self.nb_method["staticAtomAssignment"] = "on"
442

443 26
        if ".psf" in self.topology:
444 0
            self.nb_method["switching"] = "off"
445 0
            del self.nb_method["switchdist"]
446

447 26
    def config_gb_min(self):
448
        """
449
        Configure a reasonable input settings for minimization with implicit solvent. `Users can override the parameters
450
        set by this method.`
451
        """
452 0
        self.title = "Implicit solvent Minimization"
453

454 0
        self._config_min()
455

456 0
        self.implicit["GBIS"] = "on"
457 0
        self.implicit["solventDielectric"] = 78.5
458 0
        self.implicit["ionConcentration"] = 0.0
459

460 0
        self.nb_method["cutoff"] = 999.0
461 0
        self.nb_method["LJCorrection"] = "no"
462 0
        self.nb_method["nonbondedFreq"] = 2
463 0
        self.nb_method["fullElectFrequency"] = 4
464 0
        self.nb_method["wrapAll"] = "off"
465

466 0
        if ".psf" in self.topology:
467 0
            self.nb_method["switching"] = "off"
468 0
            del self.nb_method["switchdist"]
469

470 26
    def config_gb_md(self, thermostat=Thermostat.Langevin):
471
        """
472
        Configure a reasonable input settings for MD with implicit solvent. `Users can override the parameters set by
473
        this method.`
474

475
        Parameters
476
        ----------
477
        thermostat: :class:`NAMD.Thermostat`, default=Thermostat.Langevin
478
            Option to choose one of six thermostat implemented in NAMD: **(1)** `Langevin`, **(2)** `tCouple`,
479
            **(3)** `tRescale`, **(4)** `vRescale`, **(5)** `tReassign`, **(6)** `LoweAnderson`.
480
        """
481 26
        self.title = "Implicit solvent MD Simulation"
482

483 26
        self._config_md(thermostat)
484

485 26
        self.implicit["GBIS"] = "on"
486 26
        self.implicit["solventDielectric"] = 78.5
487 26
        self.implicit["ionConcentration"] = 0.0
488

489 26
        self.control["timestep"] = 1.0
490

491 26
        self.nb_method["PME"] = "off"
492 26
        self.nb_method["cutoff"] = 999.0
493 26
        self.nb_method["LJCorrection"] = "no"
494 26
        self.nb_method["nonbondedFreq"] = 2
495 26
        self.nb_method["fullElectFrequency"] = 4
496 26
        self.nb_method["wrapAll"] = "off"
497 26
        self.nb_method["stepspercycle"] = 20
498

499 26
        if ".psf" in self.topology:
500 0
            self.nb_method["switching"] = "off"
501 0
            del self.nb_method["switchdist"]
502

503 26
    def config_pbc_min(self, calculate_cell_vectors=True):
504
        """
505
        Configure a reasonable input setting for an energy minimization run with periodic boundary conditions. `Users
506
        can override the parameters set by this method.`
507

508
        Parameters
509
        ----------
510
        calculate_cell_vectors: bool, default=True
511
            Calculate cell basis vectors based on ``self.coordinates``. This is usually needed at the start of a
512
            simulation. When continuing from a previous a simulation NAMD reads in the ``*.xsc`` file to get the basis
513
            vectors.
514
        """
515 0
        self.title = "PBC Minimization"
516

517 0
        self._config_min()
518

519 0
        self.nb_method["PME"] = "yes"
520 0
        self.nb_method["PMEGridSpacing"] = 1.0
521

522 0
        if calculate_cell_vectors:
523 0
            self._get_cell_basis_vectors()
524

525 26
    def config_pbc_md(
526
        self,
527
        ensemble=Simulation.Ensemble.NPT,
528
        thermostat=Thermostat.Langevin,
529
        barostat=Barostat.Langevin,
530
        calculate_cell_vectors=False,
531
    ):
532
        """
533
        Configure a reasonable input setting for a MD run with periodic boundary conditions. `Users can override the
534
        parameters set by this method.`
535

536
        Parameters
537
        ----------
538
        ensemble: :class:`Simulation.Ensemble`, default=Ensemble.NPT
539
            Configure a MD simulation with NVE, NVT or NPT thermodynamic ensemble.
540
        thermostat: :class:`NAMD.Thermostat`, default=Thermostat.Langevin
541
            Option to choose one of six thermostat implemented in NAMD: **(1)** `Langevin`, **(2)** `tCouple`,
542
            **(3)** `tRescale`, **(4)** `vRescale`, **(5)** `tReassign`, **(6)** `LoweAnderson`.
543
        barostat: :class:`NAMD.Barostat`, default=Barostat.Langevin
544
            Option to choose one of two barostat implemented in NAMD: **(1)** `Langevin` piston or **(2)** `Berendsen`.
545
        calculate_cell_vectors: bool, default=False
546
            Calculate cell basis vectors based on ``self.coordinates``. This is usually needed at the start of a
547
            simulation. When continuing from a previous a simulation NAMD reads in the ``*.xsc`` file to get the cell
548
            basis vectors.
549
        """
550 26
        self.title = f"{ensemble.value} MD Simulation"
551

552 26
        if ensemble == self.Ensemble.NVE:
553 0
            self._config_md(thermostat=self.Thermostat.Off)
554 0
            self.control["timestep"] = 1.0
555
        else:
556 26
            self._config_md(thermostat)
557

558 26
        self.nb_method["PME"] = "yes"
559 26
        self.nb_method["PMEGridSpacing"] = 1.0
560

561
        # Select Barostat for the simulation
562 26
        if ensemble == self.Ensemble.NPT:
563 26
            self.barostat["useGroupPressure"] = "yes"
564 26
            self.barostat["useFlexibleCell"] = "no"
565 26
            self.barostat["useConstantArea"] = "no"
566

567 26
            if barostat == self.Barostat.Langevin:
568 0
                self.barostat["langevinPiston"] = "on"
569 0
                self.barostat["langevinPistonTarget"] = self.pressure
570 0
                self.barostat["langevinPistonPeriod"] = 200.0
571 0
                self.barostat["langevinPistonDecay"] = 100.0
572 0
                self.barostat["langevinPistonTemp"] = self.temperature
573

574 26
            elif barostat == self.Barostat.Berendsen:
575 26
                self.barostat["BerendsenPressure"] = "on"
576 26
                self.barostat["BerendsenPressureTarget"] = self.pressure
577 26
                self.barostat["BerendsenPressureCompressibility"] = 4.57e-5
578 26
                self.barostat["BerendsenPressureRelaxationTime"] = 100
579 26
                self.barostat["BerendsenPressureFreq"] = 10
580

581 26
        if calculate_cell_vectors:
582 0
            self._get_cell_basis_vectors()
583

584 26
    def _get_cell_basis_vectors(self):
585
        """
586
        Function to calculate the PBC cell basis vectors (needed when running a simulation for the first time).
587
        """
588 0
        structure = pmd.load_file(
589
            os.path.join(self.path, self.topology),
590
            os.path.join(self.path, self.coordinates),
591
            structure=True,
592
        )
593 0
        coordinates = structure.coordinates
594 0
        masses = np.ones(len(coordinates))
595 0
        center = pmd.geometry.center_of_mass(coordinates, masses)
596

597 0
        self.cell_basis_vectors["cellOrigin"] = list(center)
598

599
        # Check if box info exists
600 0
        if structure.box_vectors:
601 0
            import simtk.unit as unit
602

603 0
            self.cell_basis_vectors["cellBasisVector1"] = structure.box_vectors[
604
                0
605
            ].value_in_unit(unit.angstroms)
606 0
            self.cell_basis_vectors["cellBasisVector2"] = structure.box_vectors[
607
                1
608
            ].value_in_unit(unit.angstroms)
609 0
            self.cell_basis_vectors["cellBasisVector3"] = structure.box_vectors[
610
                2
611
            ].value_in_unit(unit.angstroms)
612

613
        else:
614 0
            self.cell_basis_vectors["cellBasisVector1"] = [
615
                structure.coordinates[:, 0].max() - structure.coordinates[:, 0].min(),
616
                0.0,
617
                0.0,
618
            ]
619 0
            self.cell_basis_vectors["cellBasisVector2"] = [
620
                0.0,
621
                structure.coordinates[:, 1].max() - structure.coordinates[:, 1].min(),
622
                0.0,
623
            ]
624 0
            self.cell_basis_vectors["cellBasisVector3"] = [
625
                0.0,
626
                0.0,
627
                structure.coordinates[:, 2].max() - structure.coordinates[:, 2].min(),
628
            ]
629

630 26
    @staticmethod
631
    def _write_dict_to_conf(f, dictionary):
632
        """
633
        Write dictionary to file, following NAMD format.
634

635
        Parameters
636
        ----------
637
        f : TextIO
638
            File where the dictionary should be written.
639
        dictionary : dict
640
            Dictionary of values.
641

642
        """
643 26
        for key, val in dictionary.items():
644 26
            if val is not None and not isinstance(val, list):
645 26
                f.write("{:35s} {:s}\n".format(key, str(val)))
646 0
            elif isinstance(val, list):
647 0
                f.write("{:35s} ".format(key))
648 0
                for i in val:
649 0
                    f.write("{:s} ".format(str(i)))
650 0
                f.write("\n")
651

652 26
    def _write_input_file(self):
653
        """
654
        Write the MD configuration to file.
655
        """
656

657 26
        logger.debug("Writing input file: {}".format(self.input))
658

659 26
        with open(os.path.join(self.path, self.input), "w") as conf:
660 26
            conf.write("## {:s}\n".format(self.title))
661

662
            # Topology and Coordinates files
663
            # AMBER PRMTOP/RST7
664 26
            if ".prmtop" in self.topology:
665 26
                conf.write("# AMBER files\n")
666 26
                conf.write("{:35s} {:s}\n".format("amber", "yes"))
667 26
                conf.write("{:35s} {:s}\n".format("parmfile", self.topology))
668

669 26
                if any([ext in self.coordinates for ext in [".rst7", ".inpcrd"]]):
670 26
                    conf.write("{:35s} {:s}\n".format("ambercoor", self.coordinates))
671 0
                elif ".pdb" in self.coordinates:
672 0
                    conf.write("{:35s} {:s}\n".format("coordinates", self.coordinates))
673
                else:
674 0
                    raise FileExistsError(
675
                        f"Coordinates file {self.coordinates} not does not exist."
676
                    )
677

678 26
                conf.write(
679
                    "{:35s} {:s}\n".format(
680
                        "readexclusions", self.control["readexclusions"]
681
                    )
682
                )
683 26
                conf.write("{:35s} {:s}\n".format("scnb", str(self.control["scnb"])))
684

685
            # GROMACS TOP/GRO
686 0
            elif ".top" in self.topology:
687 0
                conf.write("# GROMACS files\n")
688 0
                conf.write("{:35s} {:s}\n".format("gromacs", "on"))
689 0
                conf.write("{:35s} {:s}\n".format("grotopfile", self.topology))
690

691 0
                if ".gro" in self.coordinates:
692 0
                    conf.write("{:35s} {:s}\n".format("grocoorfile", self.coordinates))
693 0
                elif ".pdb" in self.coordinates:
694 0
                    conf.write("{:35s} {:s}\n".format("coordinates", self.coordinates))
695
                else:
696 0
                    raise FileExistsError(
697
                        f"Coordinates file {self.coordinates} not does not exist."
698
                    )
699

700
            # CHARMM PSF/PDB
701 0
            elif ".psf" in self.topology:
702 0
                conf.write("# CHARMM files\n")
703 0
                conf.write("{:35s} {:s}\n".format("paraTypeXplor", "on"))
704 0
                conf.write("{:35s} {:s}\n".format("structure", self.topology))
705

706 0
                if ".pdb" in self.coordinates:
707 0
                    conf.write("{:35s} {:s}\n".format("coordinates", self.coordinates))
708
                else:
709 0
                    raise FileExistsError(
710
                        f"Coordinates file {self.coordinates} not does not exist."
711
                    )
712

713 0
                if not self.charmm_parameters:
714 0
                    raise FileExistsError(
715
                        "CHARMM parameter file(s) not specified, please specify the parameters "
716
                        "with the variable 'namd.charmm_parameters'."
717
                    )
718

719 0
                conf.write("{:35s} {:s}\n".format("paraTypeCharmm ", "on"))
720 0
                for parameter in self.charmm_parameters:
721 0
                    conf.write("{:35s} {:s}\n".format("parameters", parameter))
722

723
            # Load in previous files (.coor, .vel and/or .xsc)
724 26
            if self.checkpoint is not None:
725 26
                conf.write("# Previous state\n")
726

727
                # Coordinates
728 26
                bincoordinates = self.checkpoint + ".restart.coor"
729 26
                if not os.path.isfile(os.path.join(self.path, bincoordinates)):
730 26
                    bincoordinates = self.checkpoint + ".coor"
731

732 26
                conf.write("{:35s} {:s}\n".format("bincoordinates", bincoordinates))
733

734
                # Velocities
735 26
                binvelocities = self.checkpoint + ".restart.vel"
736 26
                if not os.path.isfile(os.path.join(self.path, binvelocities)):
737 26
                    binvelocities = self.checkpoint + ".vel"
738

739 26
                conf.write("{:35s} {:s}\n".format("binvelocities", binvelocities))
740

741
                # PBC box information
742 26
                if not self.implicit:
743 26
                    extendedsystem = self.checkpoint + ".restart.xsc"
744 26
                    if not os.path.isfile(os.path.join(self.path, extendedsystem)):
745 26
                        extendedsystem = self.checkpoint + ".xsc"
746

747 26
                    conf.write("{:35s} {:s}\n".format("extendedSystem", extendedsystem))
748

749
            else:
750 0
                conf.write("{:35s} {:s}\n".format("temperature", str(self.temperature)))
751

752 26
            conf.write("# Output control\n")
753 26
            self._write_dict_to_conf(
754
                conf,
755
                get_dict_without_keys(
756
                    self.control,
757
                    "minimize",
758
                    "run",
759
                    "readexclusions",
760
                    "scnb",
761
                ),
762
            )
763

764 26
            if self.implicit:
765 26
                conf.write("# Generalized Born Implicit Solvent\n")
766 26
                self._write_dict_to_conf(conf, self.implicit)
767

768 26
            if self.cell_basis_vectors:
769 0
                conf.write("# Cell basis vectors\n")
770 0
                self._write_dict_to_conf(conf, self.cell_basis_vectors)
771

772 26
            conf.write("# Nonbonded options\n")
773 26
            self._write_dict_to_conf(conf, self.nb_method)
774

775 26
            conf.write("# constraints\n")
776 26
            self._write_dict_to_conf(conf, self.constraints)
777

778 26
            if self.thermostat:
779 26
                conf.write("# Temperature coupling\n")
780 26
                if self.custom_run_commands:
781 0
                    if any(
782
                        ["langevinTemp" in line for line in self.custom_run_commands]
783
                    ):
784 0
                        self._write_dict_to_conf(
785
                            conf, get_dict_without_keys(self.thermostat, "langevinTemp")
786
                        )
787
                    else:
788 0
                        self._write_dict_to_conf(conf, self.thermostat)
789
                else:
790 26
                    self._write_dict_to_conf(conf, self.thermostat)
791

792 26
            if self.barostat:
793 26
                conf.write("# Pressure coupling\n")
794 26
                self._write_dict_to_conf(conf, self.barostat)
795

796 26
            if self.colvars_file:
797 0
                conf.write("# Collective variables\n")
798 0
                conf.write("{:35s} {:s}\n".format("colvars", "on"))
799 0
                conf.write("{:35s} {:s}\n".format("colvarsConfig", self.colvars_file))
800

801 26
            if self.plumed_file:
802 0
                conf.write("# Plumed Restraints\n")
803 0
                conf.write("{:35s} {:s}\n".format("plumed", "on"))
804 0
                conf.write("{:35s} {:s}\n".format("plumedfile", self.plumed_file))
805

806 26
            if self.tcl_forces:
807 0
                conf.write("# TclForces\n")
808 0
                conf.write("{:35s} {:s}\n".format("tclForces", "on"))
809 0
                conf.write("{:35s} {:s}\n".format("tclForcesScript", self.tcl_forces))
810

811 26
            if self.qmForces:
812 0
                conf.write("# QM Forces\n")
813 0
                self._write_dict_to_conf(conf, self.qmForces)
814

815 26
            if not self.custom_run_commands:
816 26
                if "minimize" in self.control:
817 0
                    conf.write("# Minimization\n")
818 0
                    conf.write(
819
                        "{:35s} {:s}\n".format(
820
                            "minimize", str(self.control["minimize"])
821
                        )
822
                    )
823 26
                elif "run" in self.control:
824 26
                    conf.write("# Molecular dynamics\n")
825 26
                    conf.write("{:35s} {:s}\n".format("run", str(self.control["run"])))
826
            else:
827 0
                conf.write("# Custom Run\n")
828 0
                for line in self.custom_run_commands:
829 0
                    conf.write(line + "\n")
830

831 26
    def run(self, overwrite=True, fail_ok=False):
832
        """
833
        Method to run Molecular Dynamics simulation with NAMD.
834

835
        Parameters
836
        ----------
837
        overwrite: bool, optional, default=False
838
            Whether to overwrite simulation files.
839
        fail_ok: bool, optional, default=False
840
            Whether a failing simulation should stop execution of ``pAPRika``.
841

842
        """
843

844 0
        if overwrite or not self.check_complete():
845 0
            if "Minimization" in self.title:
846 0
                logger.info("Running Minimization at {}".format(self.path))
847 0
            elif "MD Simulation" in self.title:
848 0
                if self.thermostat and self.barostat:
849 0
                    logger.info("Running NPT MD at {}".format(self.path))
850 0
                elif not self.barostat:
851 0
                    logger.info("Running NVT MD at {}".format(self.path))
852
                else:
853 0
                    logger.info("Running NVE MD at {}".format(self.path))
854

855
            # Check restraints file
856 0
            if self.colvars_file and self.plumed_file:
857 0
                raise Exception(
858
                    "Cannot use both Colvars-style and Plumed-style restraints at the same time."
859
                )
860

861
            # Set Plumed kernel library to path
862 0
            if self.plumed_file:
863 0
                self._set_plumed_kernel()
864

865
            # Write input file
866 0
            self._write_input_file()
867

868
            # Create executable list
869 0
            exec_list = self.executable.split()
870

871 0
            if "+p" not in self.executable:
872 0
                exec_list += ["+p", str(self.n_threads)]
873

874 0
            if "+devices" not in self.executable and self.gpu_devices is not None:
875 0
                exec_list += ["+devices", str(self.gpu_devices)]
876

877 0
            exec_list += [self.input]
878

879 0
            logger.debug("Exec line: " + " ".join(exec_list))
880

881
            # Execute
882 0
            namd_output = sp.Popen(
883
                exec_list,
884
                cwd=self.path,
885
                stdout=open(os.path.join(self.path, self.logfile), "w"),
886
                stderr=sp.PIPE,
887
                env=os.environ,
888
            )
889 0
            namd_stderr = namd_output.stderr.read().splitlines()
890

891 0
            if namd_stderr and any(
892
                ["ERROR:" in line.decode("utf-8").strip() for line in namd_stderr]
893
            ):
894 0
                logger.info("STDERR received from NAMD execution")
895 0
                for line in namd_stderr:
896 0
                    logger.error(line)
897

898
            # Check completion status
899 0
            if "minimize" in self.control and self.check_complete():
900 0
                logger.info("Minimization completed...")
901 0
            elif self.check_complete():
902 0
                logger.info("MD completed ...")
903
            else:
904 0
                logger.info(
905
                    "Simulation did not complete when executing the following ...."
906
                )
907 0
                logger.info(" ".join(exec_list))
908 0
                if not fail_ok:
909 0
                    raise Exception(
910
                        "Exiting due to failed simulation! Check logging info."
911
                    )
912
        else:
913 0
            logger.info(
914
                "Completed output detected ... Skipping. Use: run(overwrite=True) to overwrite"
915
            )
916

917 26
    def check_complete(self, alternate_file=None):
918
        """
919
        Check for the string "WallClock" in ``self.logfile`` file. If "WallClock" is found, then
920
        the simulation completed.
921

922
        Parameters
923
        ----------
924
        alternate_file : os.PathLike, optional, default=None
925
            If present, check for "WallClock" in this file rather than ``self.logfile``.
926
            Default: None
927

928
        Returns
929
        -------
930
        complete : bool
931
            True if "WallClock" is found in file. False, otherwise.
932
        """
933
        # Assume not completed
934 0
        complete = False
935

936 0
        if alternate_file:
937 0
            output_file = alternate_file
938
        else:
939 0
            output_file = os.path.join(self.path, self.logfile)
940

941 0
        if os.path.isfile(output_file):
942 0
            with open(output_file, "r") as f:
943 0
                strings = f.read()
944 0
                if "WallClock:" in strings:
945 0
                    complete = True
946

947 0
        if complete:
948 0
            logger.debug("{} has TIMINGS".format(output_file))
949
        else:
950 0
            logger.debug("{} does not have TIMINGS".format(output_file))
951

952 0
        return complete

Read our documentation on viewing source code .

Loading