@@ -2,10 +2,12 @@
Loading
2 2
Tests the restraints utilities.
3 3
"""
4 4
5 -
import pytest
6 5
import os
7 6
7 +
import pytest
8 +
8 9
from paprika.restraints.restraints import *
10 +
from paprika.restraints.utils import *
9 11
10 12
logger = logging.getLogger(__name__)
11 13
@@ -17,7 +19,9 @@
Loading
17 19
    rest1.amber_index = True
18 20
    rest1.continuous_apr = False
19 21
    rest1.auto_apr = False
20 -
    rest1.topology = os.path.join(os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb")
22 +
    rest1.topology = os.path.join(
23 +
        os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb"
24 +
    )
21 25
    rest1.mask1 = ":CB6@O,O2,O4,O6,O8,O10"
22 26
    rest1.mask2 = ":BUT@C3"
23 27
    rest1.attach["target"] = 3.0
@@ -73,7 +77,9 @@
Loading
73 77
    rest2.amber_index = True
74 78
    rest2.continuous_apr = False
75 79
    rest2.auto_apr = False
76 -
    rest2.topology = os.path.join(os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb")
80 +
    rest2.topology = os.path.join(
81 +
        os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb"
82 +
    )
77 83
    rest2.mask1 = ":CB6@O,O2,O4,O6,O8,O10"
78 84
    rest2.mask2 = ":BUT@C3"
79 85
    rest2.mask3 = ":BUT@C"
@@ -131,7 +137,9 @@
Loading
131 137
    rest3.amber_index = True
132 138
    rest3.continuous_apr = False
133 139
    rest3.auto_apr = True
134 -
    rest3.topology = os.path.join(os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb")
140 +
    rest3.topology = os.path.join(
141 +
        os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb"
142 +
    )
135 143
    rest3.mask1 = ":CB6@O2"
136 144
    rest3.mask2 = ":CB6@O"
137 145
    rest3.mask3 = ":BUT@C3"
@@ -188,7 +196,9 @@
Loading
188 196
    rest4.amber_index = True
189 197
    rest4.continuous_apr = False
190 198
    rest4.auto_apr = False
191 -
    rest4.topology = os.path.join(os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb")
199 +
    rest4.topology = os.path.join(
200 +
        os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb"
201 +
    )
192 202
    rest4.mask1 = ":CB6@O2"
193 203
    rest4.mask2 = ":CB6@O"
194 204
    rest4.mask3 = ":BUT@C3"
@@ -243,7 +253,9 @@
Loading
243 253
    rest5.amber_index = True
244 254
    rest5.continuous_apr = False
245 255
    rest5.auto_apr = False
246 -
    rest5.topology = os.path.join(os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb")
256 +
    rest5.topology = os.path.join(
257 +
        os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb"
258 +
    )
247 259
    rest5.mask1 = ":CB6@O,O2,O4,O6,O8,O10"
248 260
    rest5.mask2 = ":BUT@C*"
249 261
    rest5.attach["target"] = 0.0
@@ -295,7 +307,9 @@
Loading
295 307
    rest6.amber_index = True
296 308
    rest6.continuous_apr = False
297 309
    rest6.auto_apr = False
298 -
    rest6.topology = os.path.join(os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb")
310 +
    rest6.topology = os.path.join(
311 +
        os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb"
312 +
    )
299 313
    rest6.mask1 = ":CB6@O,O2,O4,O6,O8,O10"
300 314
    rest6.mask2 = ":BUT@C*"
301 315
    rest6.attach["target"] = 0.0
@@ -352,7 +366,9 @@
Loading
352 366
    rest7.amber_index = True
353 367
    rest7.continuous_apr = True
354 368
    rest7.auto_apr = False
355 -
    rest7.topology = os.path.join(os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb")
369 +
    rest7.topology = os.path.join(
370 +
        os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb"
371 +
    )
356 372
    rest7.mask1 = ":1@O,O1,:BUT@H1"
357 373
    rest7.mask2 = ":CB6@N"
358 374
    rest7.attach["target"] = 0.0
@@ -400,7 +416,9 @@
Loading
400 416
    rest8.amber_index = True
401 417
    rest8.continuous_apr = False
402 418
    rest8.auto_apr = False
403 -
    rest8.topology = os.path.join(os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb")
419 +
    rest8.topology = os.path.join(
420 +
        os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb"
421 +
    )
404 422
    rest8.mask1 = ":CB6@O"
405 423
    rest8.mask2 = ":BUT@C3"
406 424
    rest8.attach["target"] = 0.0
@@ -429,7 +447,9 @@
Loading
429 447
    rest9.amber_index = True
430 448
    rest9.continuous_apr = False
431 449
    rest9.auto_apr = False
432 -
    rest9.topology = os.path.join(os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb")
450 +
    rest9.topology = os.path.join(
451 +
        os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb"
452 +
    )
433 453
    rest9.mask1 = ":CB6@O"
434 454
    rest9.mask2 = ":BUT@C3"
435 455
    rest9.pull["fc"] = 3.0
@@ -458,7 +478,9 @@
Loading
458 478
    rest10.amber_index = True
459 479
    rest10.continuous_apr = False
460 480
    rest10.auto_apr = False
461 -
    rest10.topology = os.path.join(os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb")
481 +
    rest10.topology = os.path.join(
482 +
        os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb"
483 +
    )
462 484
    rest10.mask1 = ":CB6@O"
463 485
    rest10.mask2 = ":BUT@C3"
464 486
    rest10.release["target"] = 0.0
@@ -487,4 +509,274 @@
Loading
487 509
488 510
    # Test inconsistent windows:
489 511
    with pytest.raises(Exception) as e_info:
490 -
        window_list = create_window_list([rest1, rest10])

@@ -2,13 +2,14 @@
Loading
2 2
import os
3 3
import subprocess as sp
4 4
from collections import OrderedDict
5 +
from sys import platform
5 6
6 7
logger = logging.getLogger(__name__)
7 8
8 9
9 10
class Simulation(object):
10 11
    """
11 -
    A wrapper that can be used to set AMBER simulation parameters .
12 +
    A wrapper that can be used to set AMBER simulation parameters.
12 13
    """
13 14
14 15
    @property
@@ -34,13 +35,13 @@
Loading
34 35
        self._executable = value
35 36
36 37
    @property
37 -
    def CUDA_VISIBLE_DEVICES(self):
38 +
    def gpu_devices(self):
38 39
        """A wrapper around the environmental variable ``CUDA_VISIBLE_DEVICES``."""
39 -
        return self._CUDA_VISIBLE_DEVICES
40 +
        return self._gpu_devices
40 41
41 -
    @CUDA_VISIBLE_DEVICES.setter
42 -
    def CUDA_VISIBLE_DEVICES(self, value):
43 -
        self._CUDA_VISIBLE_DEVICES = value
42 +
    @gpu_devices.setter
43 +
    def gpu_devices(self, value):
44 +
        self._gpu_devices = value
44 45
45 46
    @property
46 47
    def phase(self):
@@ -75,13 +76,44 @@
Loading
75 76
76 77
    @property
77 78
    def restraint_file(self):
78 -
        """os.PathLike: The file containing NMR-style restraints for AMBER."""
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 +
        """
79 85
        return self._restraint_file
80 86
81 87
    @restraint_file.setter
82 88
    def restraint_file(self, value):
83 89
        self._restraint_file = value
84 90
91 +
    @property
92 +
    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 +
        return self._plumed_file
100 +
101 +
    @plumed_file.setter
102 +
    def plumed_file(self, value: str):
103 +
        self._plumed_file = value
104 +
105 +
    @property
106 +
    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 +
        return self._plumed_kernel_library
112 +
113 +
    @plumed_kernel_library.setter
114 +
    def plumed_kernel_library(self, value: str):
115 +
        self._plumed_kernel_library = value
116 +
85 117
    @property
86 118
    def converged(self) -> bool:
87 119
        """bool: Whether the simulation is converged.
@@ -105,14 +137,14 @@
Loading
105 137
    @prefix.setter
106 138
    def prefix(self, new_prefix):
107 139
        self._prefix = new_prefix
108 -
        self.input = new_prefix+".in"
109 -
        self.inpcrd = new_prefix+".inpcrd"
110 -
        self.ref = new_prefix+".inpcrd"
111 -
        self.output = new_prefix+".out"
112 -
        self.restart = new_prefix+".rst7"
113 -
        self.mdinfo = new_prefix+".mdinfo"
114 -
        self.mdcrd = new_prefix+".nc"
115 -
        self.mden = new_prefix+".mden"
140 +
        self.input = new_prefix + ".in"
141 +
        self.inpcrd = new_prefix + ".inpcrd"
142 +
        self.ref = new_prefix + ".inpcrd"
143 +
        self.output = new_prefix + ".out"
144 +
        self.restart = new_prefix + ".rst7"
145 +
        self.mdinfo = new_prefix + ".mdinfo"
146 +
        self.mdcrd = new_prefix + ".nc"
147 +
        self.mden = new_prefix + ".mden"
116 148
117 149
    @property
118 150
    def cntrl(self):
@@ -208,24 +240,26 @@
Loading
208 240
        # I/O
209 241
        self._path = "."
210 242
        self._executable = "sander"
211 -
        self._CUDA_VISIBLE_DEVICES = None
243 +
        self._gpu_devices = None
212 244
        self._phase = None
213 245
        self._window = None
214 246
        self._topology = "prmtop"
215 -
        self._restraint_file = "restraints.in"
247 +
        self._restraint_file = None
248 +
        self._plumed_file = None
249 +
        self._plumed_kernel_library = None
216 250
        self.title = "PBC MD Simulation"
217 251
        self.converged = False
218 252
219 253
        # File names
220 254
        self._prefix = "md"
221 -
        self.input = self._prefix+".in"
222 -
        self.inpcrd = self._prefix+".inpcrd"
223 -
        self.ref = self._prefix+".inpcrd"
224 -
        self.output = self._prefix+".out"
225 -
        self.restart = self._prefix+".rst7"
226 -
        self.mdinfo = self._prefix+".mdinfo"
227 -
        self.mdcrd = self._prefix+".nc"
228 -
        self.mden = self._prefix+".mden"
255 +
        self.input = self._prefix + ".in"
256 +
        self.inpcrd = self._prefix + ".inpcrd"
257 +
        self.ref = self._prefix + ".inpcrd"
258 +
        self.output = self._prefix + ".out"
259 +
        self.restart = self._prefix + ".rst7"
260 +
        self.mdinfo = self._prefix + ".mdinfo"
261 +
        self.mdcrd = self._prefix + ".nc"
262 +
        self.mden = self._prefix + ".mden"
229 263
230 264
        # Input file cntrl settings (Default = NTP)
231 265
        self._cntrl = OrderedDict()
@@ -373,7 +407,13 @@
Loading
373 407
374 408
        for key, val in dictionary.items():
375 409
            if val is not None:
376 -
                f.write("  {:15s} {:s},\n".format(key+" =", str(val)))
410 +
                f.write("  {:15s} {:s},\n".format(key + " =", str(val)))
411 +
412 +
        # Write PLUMED option if preferred over Amber NMR restraints
413 +
        if self.plumed_file:
414 +
            f.write("  {:15s} {:s},\n".format("plumed = ", str(1)))
415 +
            f.write("  {:15s} {:s},\n".format("plumedfile =", f"'{self.plumed_file}'"))
416 +
377 417
        f.write(" /\n")
378 418
379 419
    def _amber_write_input_file(self):
@@ -394,11 +434,14 @@
Loading
394 434
            if self.cntrl["nmropt"] == 1:
395 435
                if self.wt is not None:
396 436
                    for line in self.wt:
397 -
                        f.write(" "+line+"\n")
437 +
                        f.write(" " + line + "\n")
398 438
                f.write(" &wt type = 'END', /\n")
439 +
440 +
                # Specify Amber NMR file
399 441
                if self.restraint_file is not None:
400 442
                    f.write("DISANG = {}\n".format(self.restraint_file))
401 443
                    f.write("LISTOUT = POUT\n\n")
444 +
402 445
            if self.group is not None:
403 446
                f.write("{:s}".format(self.group))
404 447
@@ -425,13 +468,13 @@
Loading
425 468
                # maxcyc
426 469
                ncyc = self.cntrl["ncyc"]
427 470
                maxcyc = self.cntrl["maxcyc"]
428 -
                burn_in = int(float(ncyc)+0.20 * (float(maxcyc)-float(ncyc)))
471 +
                burn_in = int(float(ncyc) + 0.20 * (float(maxcyc) - float(ncyc)))
429 472
                # If the burn_in value is nuts, then just set it to zero
430 473
                if burn_in < 0 or burn_in >= maxcyc:
431 474
                    burn_in = 0
432 475
                # Set an end_soft value that is 75% of way between ncyc and
433 476
                # maxcyc
434 -
                end_soft = int(float(ncyc)+0.60 * (float(maxcyc)-float(ncyc)))
477 +
                end_soft = int(float(ncyc) + 0.60 * (float(maxcyc) - float(ncyc)))
435 478
                self.wt = [
436 479
                    "&wt type = 'NB', istep1=0, istep2={:.0f}, value1 = 0.0, value2=0.0, IINC=50, /".format(
437 480
                        burn_in
@@ -441,6 +484,12 @@
Loading
441 484
                    ),
442 485
                ]
443 486
487 +
            # Check restraints file
488 +
            if self.restraint_file and self.plumed_file:
489 +
                raise Exception(
490 +
                    "Cannot use both NMR-style and Plumed-style restraints at the same time."
491 +
                )
492 +
444 493
            # _amber_write_input_file(self.path+'/'+self.input, self.min, title='GB Minimization.')
445 494
            self._amber_write_input_file()
446 495
@@ -450,7 +499,7 @@
Loading
450 499
                logger.info("Running MD at {}".format(self.path))
451 500
452 501
            # Create executable list for subprocess
453 -
            exec_list = self.executable.split()+["-O", "-p", self.topology]
502 +
            exec_list = self.executable.split() + ["-O", "-p", self.topology]
454 503
            if self.ref is not None:
455 504
                exec_list += ["-ref", self.ref]
456 505
            exec_list += [
@@ -470,30 +519,60 @@
Loading
470 519
            if self.mden is not None:
471 520
                exec_list += ["-e", self.mden]
472 521
473 -
            logger.debug("Exec line: "+" ".join(exec_list))
522 +
            logger.debug("Exec line: " + " ".join(exec_list))
474 523
475 -
            # Execute
476 -
            if self.CUDA_VISIBLE_DEVICES:
477 -
                amber_output = sp.Popen(
478 -
                    exec_list,
479 -
                    cwd=self.path,
480 -
                    stdout=sp.PIPE,
481 -
                    stderr=sp.PIPE,
482 -
                    env=dict(
483 -
                        os.environ, CUDA_VISIBLE_DEVICES=str(self.CUDA_VISIBLE_DEVICES)
484 -
                    ),
524 +
            # Add CUDA_VISIBLE_DEVICES variable to env dict
525 +
            if self.gpu_devices is not None:
526 +
                os.environ = dict(
527 +
                    os.environ, CUDA_VISIBLE_DEVICES=str(self.gpu_devices)
485 528
                )
529 +
530 +
            # Set the Plumed kernel library path
531 +
            if self.plumed_kernel_library is None:
532 +
                # using "startswith" as recommended from https://docs.python.org/3/library/sys.html#sys.platform
533 +
                if platform.startswith("linux"):
534 +
                    self.plumed_kernel_library = os.path.join(
535 +
                        os.environ["CONDA_PREFIX"], "lib", "libplumedKernel.so"
536 +
                    )
537 +
                elif platform.startswith("darwin"):
538 +
                    self.plumed_kernel_library = os.path.join(
539 +
                        os.environ["CONDA_PREFIX"], "lib", "libplumedKernel.dylib"
540 +
                    )
541 +
                else:
542 +
                    logger.warning(
543 +
                        f"Plumed is not supported with the '{platform}' platform."
544 +
                    )
545 +
546 +
            # Add Plumed kernel library to env dict
547 +
            if os.path.isfile(self.plumed_kernel_library):
548 +
                os.environ = dict(os.environ, PLUMED_KERNEL=self.plumed_kernel_library)
486 549
            else:
487 -
                amber_output = sp.Popen(
488 -
                    exec_list, cwd=self.path, stdout=sp.PIPE, stderr=sp.PIPE
550 +
                logger.warning(
551 +
                    f"Plumed kernel library {self.plumed_kernel_library} does not exist, "
552 +
                    "PLUMED_KERNEL environment variable is not set."
489 553
                )
490 554
491 -
            amber_output = amber_output.stdout.read().splitlines()
555 +
            # Execute
556 +
            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 +
            amber_stdout = amber_output.stdout.read().splitlines()
565 +
            amber_stderr = amber_output.stderr.read().splitlines()
492 566
493 567
            # Report any stdout/stderr which are output from execution
494 -
            if amber_output:
495 -
                logger.info("STDOUT/STDERR received from AMBER execution")
496 -
                for line in amber_output:
568 +
            if amber_stdout:
569 +
                logger.info("STDOUT received from AMBER execution")
570 +
                for line in amber_stdout:
571 +
                    logger.info(line)
572 +
573 +
            if amber_stderr:
574 +
                logger.info("STDERR received from AMBER execution")
575 +
                for line in amber_stderr:
497 576
                    logger.info(line)
498 577
499 578
            # Check completion status

@@ -1,6 +1,6 @@
Loading
1 1
import logging
2 2
3 -
from paprika.restraints.utils import parse_window
3 +
from paprika.restraints.utils import get_restraint_values, parse_window
4 4
5 5
logger = logging.getLogger(__name__)
6 6
@@ -40,13 +40,13 @@
Loading
40 40
    """
41 41
42 42
    window, phase = parse_window(window)
43 -
    if restraint.phase[phase]["force_constants"] is None and restraint.phase[phase]["targets"] is None:
43 +
    if (
44 +
        restraint.phase[phase]["force_constants"] is None
45 +
        and restraint.phase[phase]["targets"] is None
46 +
    ):
44 47
        return ""
45 48
46 -
    if not restraint.index1:
47 -
        iat1 = " "
48 -
        raise Exception("There must be at least two atoms in a restraint.")
49 -
    elif not restraint.group1:
49 +
    if not restraint.group1:
50 50
        iat1 = "{},".format(restraint.index1[0])
51 51
    else:
52 52
        iat1 = "-1,"
@@ -54,10 +54,7 @@
Loading
54 54
        for index in restraint.index1:
55 55
            igr1 += "{},".format(index)
56 56
57 -
    if not restraint.index2:
58 -
        iat2 = " "
59 -
        raise Exception("There must be at least two atoms in a restraint.")
60 -
    elif not restraint.group2:
57 +
    if not restraint.group2:
61 58
        iat2 = "{},".format(restraint.index2[0])
62 59
    else:
63 60
        iat2 = "-1,"
@@ -65,51 +62,26 @@
Loading
65 62
        for index in restraint.index2:
66 63
            igr2 += "{},".format(index)
67 64
68 -
    if not restraint.index3:
69 -
        iat3 = ""
70 -
    elif not restraint.group3:
65 +
    iat3 = ""
66 +
    if restraint.index3 and not restraint.group3:
71 67
        iat3 = "{},".format(restraint.index3[0])
72 -
    else:
68 +
    elif restraint.group3:
73 69
        iat3 = "-1,"
74 70
        igr3 = ""
75 71
        for index in restraint.index3:
76 72
            igr3 += "{},".format(index)
77 73
78 -
    if not restraint.index4:
79 -
        iat4 = ""
80 -
    elif not restraint.group4:
74 +
    iat4 = ""
75 +
    if restraint.index4 and not restraint.group4:
81 76
        iat4 = "{},".format(restraint.index4[0])
82 -
    else:
77 +
    elif restraint.group4:
83 78
        iat4 = "-1,"
84 79
        igr4 = ""
85 80
        for index in restraint.index4:
86 81
            igr4 += "{},".format(index)
87 82
88 -
89 -
    # Set upper/lower bounds depending on whether distance, angle, or torsion
90 -
    lower_bound = 0.0
91 -
    upper_bound = 999.0
92 -
93 -
    if restraint.mask3 and not restraint.mask4:
94 -
        upper_bound = 180.0
95 -
96 -
    if restraint.mask3 and restraint.mask4:
97 -
        lower_bound = restraint.phase[phase]["targets"][window] - 180.0
98 -
        upper_bound = restraint.phase[phase]["targets"][window] + 180.0
99 -
100 -
    amber_restraint_values = {
101 -
        "r1": lower_bound,
102 -
        "r2": restraint.phase[phase]["targets"][window],
103 -
        "r3": restraint.phase[phase]["targets"][window],
104 -
        "r4": upper_bound,
105 -
        "rk2": restraint.phase[phase]["force_constants"][window],
106 -
        "rk3": restraint.phase[phase]["force_constants"][window],
107 -
    }
108 -
109 -
    for key, value in restraint.custom_restraint_values.items():
110 -
        if value is not None:
111 -
            logger.debug("Overriding {} = {}".format(key, value))
112 -
            amber_restraint_values[key] = value
83 +
    # Restraint values - Amber NMR-style
84 +
    amber_restraint_values = get_restraint_values(restraint, phase, window)
113 85
114 86
    # Prepare AMBER NMR-style restraint
115 87
    atoms = "".join([iat1, iat2, iat3, iat4])

@@ -0,0 +1,42 @@
Loading
1 +
"""
2 +
Tests the dummy atoms utilities.
3 +
"""
4 +
5 +
import os
6 +
7 +
import parmed as pmd
8 +
import pytest
9 +
10 +
import paprika.dummy as dummy
11 +
from paprika.align import zalign
12 +
13 +
14 +
def test_extract_dummy():
15 +
    cb6 = pmd.load_file(
16 +
        os.path.join(os.path.dirname(__file__), "../data/cb6-but/vac.pdb")
17 +
    )
18 +
    aligned_cb6 = zalign(cb6, ":CB6", ":BUT")
19 +
    aligned_cb6 = dummy.add_dummy(aligned_cb6, residue_name="DM1", z=-6.0)
20 +
    aligned_cb6 = dummy.add_dummy(aligned_cb6, residue_name="DM2", z=-9.0)
21 +
    aligned_cb6 = dummy.add_dummy(aligned_cb6, residue_name="DM3", z=-11.2, y=2.2)
22 +
23 +
    dm_index = []
24 +
    for atom in aligned_cb6.topology.atoms():
25 +
        if atom.residue.name in ["DM1", "DM2", "DM3"]:
26 +
            dm_index.append(atom.index)
27 +
28 +
    assert len(dm_index) != 0
29 +
30 +
    dummy_dict = dummy.extract_dummy_atoms(aligned_cb6, serial=False)
31 +
32 +
    assert all(dummy_dict["DM1"]["pos"] == [0.0, 0.0, -6.0])
33 +
    assert all(dummy_dict["DM2"]["pos"] == [0.0, 0.0, -9.0])
34 +
    assert all(dummy_dict["DM3"]["pos"] == [0.0, 2.2, -11.2])
35 +
36 +
    assert dummy_dict["DM1"]["mass"] == 208.0
37 +
    assert dummy_dict["DM2"]["mass"] == 208.0
38 +
    assert dummy_dict["DM3"]["mass"] == 208.0
39 +
40 +
    assert dummy_dict["DM1"]["idx"] == dm_index[0]
41 +
    assert dummy_dict["DM2"]["idx"] == dm_index[1]
42 +
    assert dummy_dict["DM3"]["idx"] == dm_index[2]

@@ -186,3 +186,46 @@
Loading
186 186
                residue_name[0:3], atom_name, atom_type[0:2]
187 187
            )
188 188
        )
189 +
190 +
191 +
def extract_dummy_atoms(structure, resname=None, serial=True):
192 +
    """
193 +
    Extract information about dummy atoms from a parmed structure and
194 +
    returns the information as a dictionary.
195 +
196 +
    Parameters
197 +
    ----------
198 +
    structure : :class:`parmed.structure.Structure`
199 +
        The parmed structure object we want to extract from
200 +
    resname : list
201 +
        List of residue name for the dummy atoms (default: ["DM1", "DM2", "DM3"])
202 +
    serial : bool
203 +
        get indices in serial (starts from 1) or index (starts from 0)
204 +
205 +
    Returns
206 +
    -------
207 +
    dummy_atoms : dict
208 +
209 +
    Output example
210 +
    --------------
211 +
212 +
        main keys: {'DM1', 'DM2', 'DM3'}
213 +
        sub keys: 'pos'      - cartesian coordinates (x,y,z)
214 +
                  'idx'      - atom indices
215 +
                  'idx_type' - type of atom index (serial or index)
216 +
                  'mass'     - mass of dummy atom
217 +
    """
218 +
219 +
    if resname is None:
220 +
        resname = ["DM1", "DM2", "DM3"]
221 +
222 +
    dummy_atoms = {name: {} for name in resname}
223 +
224 +
    for dummy_atom in resname:
225 +
        residue = f":{dummy_atom}"
226 +
        dummy_atoms[dummy_atom]["pos"] = structure[residue].coordinates[0]
227 +
        dummy_atoms[dummy_atom]["mass"] = [atom.mass for atom in structure[residue].atoms][0]
228 +
        dummy_atoms[dummy_atom]["idx"] = utils.index_from_mask(structure, residue, amber_index=serial)[0]
229 +
        dummy_atoms[dummy_atom]["idx_type"] = "serial" if serial else "index"
230 +
231 +
    return dummy_atoms

@@ -11,6 +11,31 @@
Loading
11 11
logger = logging.getLogger(__name__)
12 12
13 13
14 +
def get_key(dct, value):
15 +
    """
16 +
    Get dictionary key given the value.
17 +
18 +
    NOTE: this function will return a list of keys if there are more than one key with the same value
19 +
          in the order they are found.
20 +
    """
21 +
    key = [key for key in dct if (dct[key] == value)]
22 +
23 +
    if len(key) > 1:
24 +
        logger.warning(
25 +
            "There more than one key with the same value. Please check if this is the desired output."
26 +
        )
27 +
28 +
    return key
29 +
30 +
31 +
def override_dict(dct, custom):
32 +
    """Overrides dictionary values from that of a custom dictionary."""
33 +
    for key, value in custom.items():
34 +
        if value is not None:
35 +
            logger.debug("Overriding {} = {}".format(key, value))
36 +
            dct[key] = value
37 +
38 +
14 39
def return_parmed_structure(filename):
15 40
    """
16 41
    Return a structure object from a filename.
@@ -34,6 +59,7 @@
Loading
34 59
        logger.error("Unable to load file: {}".format(filename))
35 60
    return structure
36 61
62 +
37 63
@lru_cache(maxsize=32)
38 64
def index_from_mask(structure, mask, amber_index=False):
39 65
    """
@@ -181,6 +207,7 @@
Loading
181 207
182 208
    return energies
183 209
210 +
184 211
def parse_mdout(file):
185 212
    """
186 213
    Return energies from an AMBER `mdout` file.

@@ -0,0 +1,407 @@
Loading
1 +
import logging
2 +
import os
3 +
4 +
import numpy as np
5 +
from parmed.structure import Structure as ParmedStructureClass
6 +
7 +
from paprika.dummy import extract_dummy_atoms
8 +
from paprika.restraints.utils import get_bias_potential_type, parse_window
9 +
from paprika.utils import get_key, return_parmed_structure
10 +
11 +
logger = logging.getLogger(__name__)
12 +
13 +
_PI_ = np.pi
14 +
15 +
16 +
class Plumed:
17 +
    """
18 +
    This class converts restraints generated in pAPRika DAT_restraints into Plumed restraints.
19 +
20 +
    # TODO: possibly change this module to use the python wrapper of Plumed.
21 +
22 +
    Example:
23 +
    -------
24 +
        >>> plumed = Plumed()
25 +
        >>> plumed.file_name = 'plumed.dat'
26 +
        >>> plumed.path = './windows'
27 +
        >>> plumed.window_list = window_list
28 +
        >>> plumed.restraint_list = restraint_list
29 +
        >>> plumed.dump_to_file()
30 +
31 +
    plumed.dat:
32 +
        UNITS LENGTH=A ENERGY=kcal/mol TIME=ns
33 +
        # Collective variables
34 +
        c1 : DISTANCE ATOMS=175,150 NOPBC
35 +
        c2 : ANGLE    ATOMS=176,175,150 NOPBC
36 +
        c3 : ANGLE    ATOMS=175,150,165 NOPBC
37 +
        # Bias potential
38 +
        RESTRAINT   ARG=c7  AT=  6.000 KAPPA=  10.00
39 +
        RESTRAINT   ARG=c8  AT=  3.142 KAPPA= 200.00
40 +
        RESTRAINT   ARG=c9  AT=  3.142 KAPPA= 200.00
41 +
    """
42 +
43 +
    @property
44 +
    def file_name(self):
45 +
        """
46 +
        str: The Plumed file name where the restraints will be written to.
47 +
        """
48 +
        return self._file_name
49 +
50 +
    @file_name.setter
51 +
    def file_name(self, value: str):
52 +
        self._file_name = value
53 +
54 +
    @property
55 +
    def window_list(self):
56 +
        """
57 +
        list: The list of APR windows
58 +
        """
59 +
        return self._window_list
60 +
61 +
    @window_list.setter
62 +
    def window_list(self, value: list):
63 +
        self._window_list = value
64 +
65 +
    @property
66 +
    def restraint_list(self):
67 +
        """
68 +
        list: The list of restraints to be converted to Plumed.
69 +
        """
70 +
        return self._restraint_list
71 +
72 +
    @restraint_list.setter
73 +
    def restraint_list(self, value: list):
74 +
        self._restraint_list = value
75 +
76 +
    @property
77 +
    def path(self):
78 +
        """
79 +
        os.PathLike: The parent directory that contains the simulation windows.
80 +
        """
81 +
        return self._path
82 +
83 +
    @path.setter
84 +
    def path(self, value: str):
85 +
        self._path = value
86 +
87 +
    @property
88 +
    def uses_legacy_k(self):
89 +
        """
90 +
        bool: Option to specify whether the force constant parsed into DAT_restraint()
91 +
        is Amber-style or Gromacs/NAMD-style. Amber-style force constants have their value
92 +
        multiplied by a factor of 1/2 whereas Gromacs/NAMD-style does not. Plumed follows
93 +
        the Gromacs/NAMD-style for the force constant and the equations below demonstrates
94 +
        this point.
95 +
96 +
            * Amber:  U = K x (x - x0)²    * Plumed: U = 1/2 x k x (x - x0)²
97 +
            --> K(Amber) = 1/2 k(Plumed)
98 +
99 +
        i.e. "uses_legacy_k" is set to True (default) the force constants will be multiplied by 2.
100 +
        """
101 +
        return self._uses_legacy_k
102 +
103 +
    @uses_legacy_k.setter
104 +
    def uses_legacy_k(self, value: bool):
105 +
        self._uses_legacy_k = value
106 +
107 +
    @property
108 +
    def units(self):
109 +
        """
110 +
        dict: Dictionary of units for Plumed as strings. The dictionary requires the key values
111 +
        of 'energy', 'length, 'time'.
112 +
        """
113 +
        return self._units
114 +
115 +
    @units.setter
116 +
    def units(self, value: dict):
117 +
        self._units = value
118 +
119 +
    def __init__(self):
120 +
        self._file_name = "plumed.dat"
121 +
        self._restraint_list = None
122 +
        self._window_list = None
123 +
        self._path = "./"
124 +
        self._uses_legacy_k = True
125 +
        self.k_factor = 1.0
126 +
        self._units = None
127 +
        self.header_line = None
128 +
        self.group_index = None
129 +
        self.group_atoms = None
130 +
131 +
    def _initialize(self):
132 +
        # Set factor for spring constant
133 +
        if self.uses_legacy_k:
134 +
            self.k_factor = 2.0
135 +
136 +
        # Check units
137 +
        # NOTE: We have to resort to strings until (if) we migrate all quantities to Pint/SimTK units
138 +
        if self.units is None:
139 +
            self.units = {"energy": "kcal/mol", "length": "A", "time": "ns"}
140 +
        _check_plumed_units(self.units)
141 +
142 +
        # header line
143 +
        self.header_line = (
144 +
            f"UNITS LENGTH={self.units['length']} "
145 +
            f"ENERGY={self.units['energy']} "
146 +
            f"TIME={self.units['time']}"
147 +
        )
148 +
149 +
    def dump_to_file(self):
150 +
151 +
        self._initialize()
152 +
153 +
        # Loop over APR windows
154 +
        for windows in self.window_list:
155 +
156 +
            window, phase = parse_window(windows)
157 +
158 +
            # Check if file exist and write header line
159 +
            with open(os.path.join(self.path, windows, self.file_name), "w") as file:
160 +
                file.write(self.header_line + "\n")
161 +
162 +
            cv_index = 1
163 +
            cv_dict = {}
164 +
            cv_lines = []
165 +
            bias_lines = []
166 +
167 +
            self.group_index = 1
168 +
            self.group_atoms = {}
169 +
170 +
            # Parse each restraint in the list
171 +
            for restraint in self.restraint_list:
172 +
                # Skip restraint if the target or force constant is not defined.
173 +
                # Example: wall restraints only used during the attach phase.
174 +
                try:
175 +
                    target = restraint.phase[phase]["targets"][window]
176 +
                    force_constant = (
177 +
                        restraint.phase[phase]["force_constants"][window]
178 +
                        * self.k_factor
179 +
                    )
180 +
                except TypeError:
181 +
                    continue
182 +
183 +
                # Convert list to comma-separated string
184 +
                atom_index = self._get_atom_indices(restraint)
185 +
                atom_string = ",".join(map(str, atom_index))
186 +
187 +
                # Determine collective variable type
188 +
                colvar_type = "distance"
189 +
                if len(atom_index) == 3:
190 +
                    colvar_type = "angle"
191 +
                    target *= _PI_ / 180.0
192 +
                elif len(atom_index) == 4:
193 +
                    colvar_type = "torsion"
194 +
                    target *= _PI_ / 180.0
195 +
196 +
                # Determine bias type for this restraint
197 +
                bias_type = get_bias_potential_type(restraint, phase, window)
198 +
199 +
                # Append cv strings to lists
200 +
                # The code below prevents duplicate cv definition.
201 +
                # While not necessary, it makes the plumed file cleaner.
202 +
                if not get_key(cv_dict, atom_string):
203 +
                    cv_key = f"c{cv_index}"
204 +
                    cv_dict[cv_key] = atom_string
205 +
206 +
                    cv_lines.append(
207 +
                        f"{cv_key}: {colvar_type.upper()} ATOMS={atom_string} NOPBC\n"
208 +
                    )
209 +
                    bias_lines.append(
210 +
                        f"{bias_type.upper()} ARG={cv_key} AT={target:.4f} KAPPA={force_constant:.2f}\n"
211 +
                    )
212 +
                else:
213 +
                    cv_key = get_key(cv_dict, atom_string)[0]
214 +
215 +
                    bias_lines.append(
216 +
                        f"{bias_type.upper()} ARG={cv_key} AT={target:.4f} KAPPA={force_constant:.2f}\n"
217 +
                    )
218 +
219 +
                # Increment cv index
220 +
                cv_index += 1
221 +
222 +
            # Write collective variables to file
223 +
            self._write_colvar_to_file(windows, cv_lines, bias_lines)
224 +
225 +
    def _write_colvar_to_file(self, window, cv_list, bias_list):
226 +
        with open(os.path.join(self.path, window, self.file_name), "a") as file:
227 +
            if len(self.group_atoms) != 0:
228 +
                file.write("# Centroid groups\n")
229 +
                for key, value in self.group_atoms.items():
230 +
                    file.write(f"{key}: COM ATOMS={value}\n")
231 +
232 +
            file.write("# Collective variables\n")
233 +
            for line in cv_list:
234 +
                file.write(line)
235 +
236 +
            file.write("# Bias potentials\n")
237 +
            for line in bias_list:
238 +
                file.write(line)
239 +
240 +
    def _get_atom_indices(self, restraint):
241 +
        # Check atom index setting
242 +
        index_shift = 0
243 +
        if not restraint.amber_index:
244 +
            index_shift = 1
245 +
            logger.debug("Atom indices starts from 0 --> shifting indices by 1.")
246 +
247 +
        # Collect DAT atom indices
248 +
        atom_index = []
249 +
250 +
        if not restraint.group1:
251 +
            atom_index.append(restraint.index1[0] + index_shift)
252 +
        else:
253 +
            igr1 = ""
254 +
            for index in restraint.index1:
255 +
                igr1 += "{},".format(index + index_shift)
256 +
257 +
            if not get_key(self.group_atoms, igr1):
258 +
                self.group_atoms[f"g{self.group_index}"] = igr1
259 +
                self.group_index += 1
260 +
261 +
            atom_index.append(get_key(self.group_atoms, igr1)[0])
262 +
263 +
        if not restraint.group2:
264 +
            atom_index.append(restraint.index2[0] + index_shift)
265 +
        else:
266 +
            igr2 = ""
267 +
            for index in restraint.index2:
268 +
                igr2 += "{},".format(index + index_shift)
269 +
270 +
            if not get_key(self.group_atoms, igr2):
271 +
                self.group_atoms[f"g{self.group_index}"] = igr2
272 +
                self.group_index += 1
273 +
274 +
            atom_index.append(get_key(self.group_atoms, igr2)[0])
275 +
276 +
        if restraint.index3 and not restraint.group3:
277 +
            atom_index.append(restraint.index3[0] + index_shift)
278 +
        elif restraint.group3:
279 +
            igr3 = ""
280 +
            for index in restraint.index3:
281 +
                igr3 += "{},".format(index + index_shift)
282 +
283 +
            if not get_key(self.group_atoms, igr3):
284 +
                self.group_atoms[f"g{self.group_index}"] = igr3
285 +
                self.group_index += 1
286 +
287 +
            atom_index.append(get_key(self.group_atoms, igr3)[0])
288 +
289 +
        if restraint.index4 and not restraint.group4:
290 +
            atom_index.append(restraint.index4[0] + index_shift)
291 +
        elif restraint.group4:
292 +
            igr4 = ""
293 +
            for index in restraint.index4:
294 +
                igr4 += "{},".format(index + index_shift)
295 +
296 +
            if not get_key(self.group_atoms, igr4):
297 +
                self.group_atoms[f"g{self.group_index}"] = igr4
298 +
                self.group_index += 1
299 +
300 +
            atom_index.append(get_key(self.group_atoms, igr4)[0])
301 +
302 +
        return atom_index
303 +
304 +
    def add_dummy_atoms_to_file(self, structure):
305 +
        # Extract dummy atoms
306 +
        for windows in self.window_list:
307 +
            if isinstance(structure, str):
308 +
                structure = return_parmed_structure(structure)
309 +
            elif isinstance(structure, ParmedStructureClass):
310 +
                pass
311 +
            else:
312 +
                raise Exception(
313 +
                    "add_dummy_atoms_to_file does not support the type associated with structure: "
314 +
                    + type(structure)
315 +
                )
316 +
317 +
            dummy_atoms = extract_dummy_atoms(structure, serial=True)
318 +
319 +
            # Write dummy atom info to plumed file
320 +
            plumed_file = os.path.join(self.path, windows, self.file_name)
321 +
            if os.path.isfile(plumed_file):
322 +
                with open(plumed_file, "a") as file:
323 +
                    _write_dummy_to_file(file, dummy_atoms)
324 +
            else:
325 +
                raise Exception(f"ERROR: '{plumed_file}' file does not exists!")
326 +
327 +
328 +
def _check_plumed_units(units):
329 +
    """
330 +
    Checks the specified units and makes sure that its supported.
331 +
    """
332 +
    if units["energy"] not in ["kj/mol", "kcal/mol"]:
333 +
        raise Exception(
334 +
            f"Specified unit for energy ({units['energy']}) is not supported."
335 +
        )
336 +
337 +
    if units["length"] not in ["nm", "A"]:
338 +
        raise Exception(
339 +
            f"Specified unit for length ({units['length']}) is not supported."
340 +
        )
341 +
342 +
    if units["time"] not in ["ps", "fs", "ns"]:
343 +
        raise Exception(f"Specified unit for time ({units['time']}) is not supported.")
344 +
345 +
346 +
def _write_dummy_to_file(file, dummy_atoms, kpos=100.0):
347 +
    """
348 +
    Append to the "plumed.dat" file the dummy atoms colvar definition and position restraints
349 +
350 +
    Parameters
351 +
    ----------
352 +
    file : class '_io.TextIOWrapper'
353 +
        The file object handle to save the plumed file.
354 +
    dummy_atoms : dict
355 +
        Dictionary containing information about the dummy atoms.
356 +
    kpos : float
357 +
        Spring constant used to restrain dummy atoms (kcal/mol/A^2).
358 +
359 +
    Output
360 +
    ------
361 +
    dm1: POSITION ATOM = 170 NOPBC
362 +
    dm2: POSITION ATOM = 171 NOPBC
363 +
    dm3: POSITION ATOM = 172 NOPBC
364 +
    RESTRAINT...
365 +
    ARG = dm1.x, dm1.y, dm1.z, dm2.x, dm2.y, dm2.z, dm3.x, dm3.y, dm3.z
366 +
    AT = 19.68, 20.3, 26.9, 19.68, 20.3, 23.9, 19.68, 22.5, 21.7
367 +
    KAPPA = 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0
368 +
    ...
369 +
    RESTRAINT
370 +
371 +
    """
372 +
373 +
    file.write("# Dummy Atoms\n")
374 +
    file.write(f"dm1: POSITION ATOM={dummy_atoms['DM1']['idx']} NOPBC\n")
375 +
    file.write(f"dm2: POSITION ATOM={dummy_atoms['DM2']['idx']} NOPBC\n")
376 +
    file.write(f"dm3: POSITION ATOM={dummy_atoms['DM3']['idx']} NOPBC\n")
377 +
378 +
    arg = "dm1.x,dm1.y,dm1.z," "dm2.x,dm2.y,dm2.z," "dm3.x,dm3.y,dm3.z,"
379 +
380 +
    at = (
381 +
        f"{dummy_atoms['DM1']['pos'][0]:0.3f},"
382 +
        f"{dummy_atoms['DM1']['pos'][1]:0.3f},"
383 +
        f"{dummy_atoms['DM1']['pos'][2]:0.3f},"
384 +
    )
385 +
    at += (
386 +
        f"{dummy_atoms['DM2']['pos'][0]:0.3f},"
387 +
        f"{dummy_atoms['DM2']['pos'][1]:0.3f},"
388 +
        f"{dummy_atoms['DM2']['pos'][2]:0.3f},"
389 +
    )
390 +
    at += (
391 +
        f"{dummy_atoms['DM3']['pos'][0]:0.3f},"
392 +
        f"{dummy_atoms['DM3']['pos'][1]:0.3f},"
393 +
        f"{dummy_atoms['DM3']['pos'][2]:0.3f},"
394 +
    )
395 +
396 +
    kappa = (
397 +
        f"{kpos:0.1f},{kpos:0.1f},{kpos:0.1f},"
398 +
        f"{kpos:0.1f},{kpos:0.1f},{kpos:0.1f},"
399 +
        f"{kpos:0.1f},{kpos:0.1f},{kpos:0.1f},"
400 +
    )
401 +
402 +
    file.write("RESTRAINT ...\n")
403 +
    file.write(f"ARG={arg}\n")
404 +
    file.write(f"AT={at}\n")
405 +
    file.write(f"KAPPA={kappa}\n")
406 +
    file.write("LABEL=dummy\n")
407 +
    file.write("... RESTRAINT\n")

@@ -1,3 +1,14 @@
Loading
1 +
import logging
2 +
3 +
import numpy as np
4 +
5 +
from paprika.utils import override_dict
6 +
7 +
logger = logging.getLogger(__name__)
8 +
9 +
_PI_ = np.pi
10 +
11 +
1 12
def parse_window(window):
2 13
    """
3 14
    Utility function to use a path to index a :class:`paprika.restraints.DAT_restraint` instance.
@@ -25,4 +36,125 @@
Loading
25 36
        raise Exception("Cannot determine the phase for this restraint.")
26 37
    window = int(window[1:])
27 38
28 -
    return window, phase

@@ -3,7 +3,6 @@
Loading
3 3
import numpy as np
4 4
import parmed as pmd
5 5
import pytraj as pt
6 -
7 6
from paprika import utils
8 7
9 8
logger = logging.getLogger(__name__)
@@ -307,8 +306,8 @@
Loading
307 306
                for phs in ["attach", "pull", "release"]:
308 307
                    for value in ["force_constants", "targets"]:
309 308
                        if (
310 -
                                self_dictionary["phase"][phs][value] is None
311 -
                                and other_dictionary["phase"][phs][value] is None
309 +
                            self_dictionary["phase"][phs][value] is None
310 +
                            and other_dictionary["phase"][phs][value] is None
312 311
                        ):
313 312
                            continue
314 313
                        else:
@@ -328,9 +327,10 @@
Loading
328 327
        restraint_dictionary: dict
329 328
            The restraint dictionary corresponding whose values will be set.
330 329
        method: str
331 -
            Which "method" will be used to set the remaining values. The "method" refers to which combination of restraint
332 -
            inputs are provided -- for example, initial value, final value, and number of windows or initial value, increment
333 -
            size, and number of increments -- which are defined in :meth:`paprika.restraints.DAT_restraint.initialize`.
330 +
            Which "method" will be used to set the remaining values. The "method" refers to which
331 +
            combination of restraint inputs are provided -- for example, initial value, final value,
332 +
            and number of windows or initial value, increment size, and number of increments -- which
333 +
            are defined in :meth:`paprika.restraints.DAT_restraint.initialize`.
334 334
335 335
            This appears to be a ``str`` but should be an ``int``.
336 336
@@ -345,39 +345,50 @@
Loading
345 345
        # Attach/Release, Force Constant Method 1
346 346
        if phase in ("a", "r") and method == "1":
347 347
            force_constants = np.linspace(
348 -
                restraint_dictionary["fc_initial"], restraint_dictionary["fc_final"],
349 -
                restraint_dictionary["num_windows"]
348 +
                restraint_dictionary["fc_initial"],
349 +
                restraint_dictionary["fc_final"],
350 +
                restraint_dictionary["num_windows"],
350 351
            )
351 352
352 353
        # Attach/Release, Force Constant Method 1a
353 354
        elif phase in ("a", "r") and method == "1a":
354 -
            force_constants = np.linspace(0.0, restraint_dictionary["fc_final"], restraint_dictionary["num_windows"])
355 +
            force_constants = np.linspace(
356 +
                0.0,
357 +
                restraint_dictionary["fc_final"],
358 +
                restraint_dictionary["num_windows"],
359 +
            )
355 360
356 361
        # Attach/Release, Force Constant Method 2
357 362
        elif phase in ("a", "r") and method == "2":
358 363
            force_constants = np.arange(
359 364
                restraint_dictionary["fc_initial"],
360 -
                restraint_dictionary["fc_final"]+restraint_dictionary["fc_increment"],
365 +
                restraint_dictionary["fc_final"] + restraint_dictionary["fc_increment"],
361 366
                restraint_dictionary["fc_increment"],
362 367
            )
363 368
364 369
        # Attach/Release, Force Constant Method 2a
365 370
        elif phase in ("a", "r") and method == "2a":
366 371
            force_constants = np.arange(
367 -
                0.0, restraint_dictionary["fc_final"]+restraint_dictionary["fc_increment"],
368 -
                restraint_dictionary["fc_increment"]
372 +
                0.0,
373 +
                restraint_dictionary["fc_final"] + restraint_dictionary["fc_increment"],
374 +
                restraint_dictionary["fc_increment"],
369 375
            )
370 376
371 377
        # Attach/Release, Force Constant Method 3
372 378
        elif phase in ("a", "r") and method == "3":
373 379
            force_constants = np.asarray(
374 -
                [fraction * restraint_dictionary["fc_final"] for fraction in restraint_dictionary["fraction_list"]]
380 +
                [
381 +
                    fraction * restraint_dictionary["fc_final"]
382 +
                    for fraction in restraint_dictionary["fraction_list"]
383 +
                ]
375 384
            )
376 385
377 386
        # Attach/Release, Force Constant Method 4
378 387
        elif phase in ("a", "r") and method == "4":
379 388
            fractions = np.arange(
380 -
                0, 1.0+restraint_dictionary["fraction_increment"], restraint_dictionary["fraction_increment"]
389 +
                0,
390 +
                1.0 + restraint_dictionary["fraction_increment"],
391 +
                restraint_dictionary["fraction_increment"],
381 392
            )
382 393
            force_constants = np.asarray(
383 394
                [fraction * restraint_dictionary["fc_final"] for fraction in fractions]
@@ -389,24 +400,32 @@
Loading
389 400
390 401
        # Attach/Release, Target Method
391 402
        if phase in ("a", "r"):
392 -
            targets = np.asarray([restraint_dictionary["target"]] * len(force_constants))
403 +
            targets = np.asarray(
404 +
                [restraint_dictionary["target"]] * len(force_constants)
405 +
            )
393 406
394 407
        # Pull, Target Method 1
395 408
        if phase == "p" and method == "1":
396 409
            targets = np.linspace(
397 -
                restraint_dictionary["target_initial"], restraint_dictionary["target_final"],
398 -
                restraint_dictionary["num_windows"]
410 +
                restraint_dictionary["target_initial"],
411 +
                restraint_dictionary["target_final"],
412 +
                restraint_dictionary["num_windows"],
399 413
            )
400 414
401 415
        # Pull, Target Method 1a
402 416
        elif phase == "p" and method == "1a":
403 -
            targets = np.linspace(0.0, restraint_dictionary["target_final"], restraint_dictionary["num_windows"])
417 +
            targets = np.linspace(
418 +
                0.0,
419 +
                restraint_dictionary["target_final"],
420 +
                restraint_dictionary["num_windows"],
421 +
            )
404 422
405 423
        # Pull, Target Method 2
406 424
        elif phase == "p" and method == "2":
407 425
            targets = np.arange(
408 426
                restraint_dictionary["target_initial"],
409 -
                restraint_dictionary["target_final"]+restraint_dictionary["target_increment"],
427 +
                restraint_dictionary["target_final"]
428 +
                + restraint_dictionary["target_increment"],
410 429
                restraint_dictionary["target_increment"],
411 430
            )
412 431
@@ -414,7 +433,8 @@
Loading
414 433
        elif phase == "p" and method == "2a":
415 434
            targets = np.arange(
416 435
                0.0,
417 -
                restraint_dictionary["target_final"]+restraint_dictionary["target_increment"],
436 +
                restraint_dictionary["target_final"]
437 +
                + restraint_dictionary["target_increment"],
418 438
                restraint_dictionary["target_increment"],
419 439
            )
420 440
@@ -430,10 +450,15 @@
Loading
430 450
        # Pull, Target Method 4
431 451
        elif phase == "p" and method == "4":
432 452
            fractions = np.arange(
433 -
                0, 1.0+restraint_dictionary["fraction_increment"], restraint_dictionary["fraction_increment"]
453 +
                0,
454 +
                1.0 + restraint_dictionary["fraction_increment"],
455 +
                restraint_dictionary["fraction_increment"],
434 456
            )
435 457
            targets = np.asarray(
436 -
                [fraction * restraint_dictionary["target_final"] for fraction in fractions]
458 +
                [
459 +
                    fraction * restraint_dictionary["target_final"]
460 +
                    for fraction in fractions
461 +
                ]
437 462
            )
438 463
439 464
        # Pull, Target Method 5
@@ -496,8 +521,8 @@
Loading
496 521
        targets = None
497 522
498 523
        if (
499 -
                self.attach["num_windows"] is not None
500 -
                and self.attach["fc_final"] is not None
524 +
            self.attach["num_windows"] is not None
525 +
            and self.attach["fc_final"] is not None
501 526
        ):
502 527
            if self.attach["fc_initial"] is not None:
503 528
                logger.debug("Attach, Method #1")
@@ -507,8 +532,8 @@
Loading
507 532
                force_constants, targets = self._calc_method("a", self.attach, "1a")
508 533
509 534
        elif (
510 -
                self.attach["fc_increment"] is not None
511 -
                and self.attach["fc_final"] is not None
535 +
            self.attach["fc_increment"] is not None
536 +
            and self.attach["fc_final"] is not None
512 537
        ):
513 538
            if self.attach["fc_initial"] is not None:
514 539
                logger.debug("Attach, Method #2")
@@ -518,15 +543,15 @@
Loading
518 543
                force_constants, targets = self._calc_method("a", self.attach, "2a")
519 544
520 545
        elif (
521 -
                self.attach["fraction_list"] is not None
522 -
                and self.attach["fc_final"] is not None
546 +
            self.attach["fraction_list"] is not None
547 +
            and self.attach["fc_final"] is not None
523 548
        ):
524 549
            logger.debug("Attach, Method #3")
525 550
            force_constants, targets = self._calc_method("a", self.attach, "3")
526 551
527 552
        elif (
528 -
                self.attach["fraction_increment"] is not None
529 -
                and self.attach["fc_final"] is not None
553 +
            self.attach["fraction_increment"] is not None
554 +
            and self.attach["fc_final"] is not None
530 555
        ):
531 556
            logger.debug("Attach, Method #4")
532 557
            force_constants, targets = self._calc_method("a", self.attach, "4")
@@ -563,8 +588,8 @@
Loading
563 588
            self.pull["target_initial"] = self.phase["attach"]["targets"][-1]
564 589
565 590
        if (
566 -
                self.pull["num_windows"] is not None
567 -
                and self.pull["target_final"] is not None
591 +
            self.pull["num_windows"] is not None
592 +
            and self.pull["target_final"] is not None
568 593
        ):
569 594
            if self.pull["target_initial"] is not None:
570 595
                logger.debug("Pull, Method #1")
@@ -574,8 +599,8 @@
Loading
574 599
                force_constants, targets = self._calc_method("p", self.pull, "1a")
575 600
576 601
        elif (
577 -
                self.pull["target_increment"] is not None
578 -
                and self.pull["target_final"] is not None
602 +
            self.pull["target_increment"] is not None
603 +
            and self.pull["target_final"] is not None
579 604
        ):
580 605
            if self.pull["target_initial"] is not None:
581 606
                logger.debug("Pull, Method #2")
@@ -585,15 +610,15 @@
Loading
585 610
                force_constants, targets = self._calc_method("p", self.pull, "2a")
586 611
587 612
        elif (
588 -
                self.pull["fraction_list"] is not None
589 -
                and self.pull["target_final"] is not None
613 +
            self.pull["fraction_list"] is not None
614 +
            and self.pull["target_final"] is not None
590 615
        ):
591 616
            logger.debug("Pull, Method #3")
592 617
            force_constants, targets = self._calc_method("p", self.pull, "3")
593 618
594 619
        elif (
595 -
                self.pull["fraction_increment"] is not None
596 -
                and self.pull["target_final"] is not None
620 +
            self.pull["fraction_increment"] is not None
621 +
            and self.pull["target_final"] is not None
597 622
        ):
598 623
            logger.debug("Pull, Method #4")
599 624
            force_constants, targets = self._calc_method("p", self.pull, "4")
@@ -644,8 +669,8 @@
Loading
644 669
                    self.release[key] = self.attach[key]
645 670
646 671
        if (
647 -
                self.release["num_windows"] is not None
648 -
                and self.release["fc_final"] is not None
672 +
            self.release["num_windows"] is not None
673 +
            and self.release["fc_final"] is not None
649 674
        ):
650 675
            if self.release["fc_initial"] is not None:
651 676
                logger.debug("Release, Method #1")
@@ -655,8 +680,8 @@
Loading
655 680
                force_constants, targets = self._calc_method("r", self.release, "1a")
656 681
657 682
        elif (
658 -
                self.release["fc_increment"] is not None
659 -
                and self.release["fc_final"] is not None
683 +
            self.release["fc_increment"] is not None
684 +
            and self.release["fc_final"] is not None
660 685
        ):
661 686
            if self.release["fc_initial"] is not None:
662 687
                logger.debug("Release, Method #2")
@@ -666,15 +691,15 @@
Loading
666 691
                force_constants, targets = self._calc_method("r", self.release, "2a")
667 692
668 693
        elif (
669 -
                self.release["fraction_list"] is not None
670 -
                and self.release["fc_final"] is not None
694 +
            self.release["fraction_list"] is not None
695 +
            and self.release["fc_final"] is not None
671 696
        ):
672 697
            logger.debug("Release, Method #3")
673 698
            force_constants, targets = self._calc_method("r", self.release, "3")
674 699
675 700
        elif (
676 -
                self.release["fraction_increment"] is not None
677 -
                and self.release["fc_final"] is not None
701 +
            self.release["fraction_increment"] is not None
702 +
            and self.release["fc_final"] is not None
678 703
        ):
679 704
            logger.debug("Release, Method #4")
680 705
            force_constants, targets = self._calc_method("r", self.release, "4")
@@ -742,12 +767,12 @@
Loading
742 767
743 768
744 769
def static_DAT_restraint(
745 -
        restraint_mask_list,
746 -
        num_window_list,
747 -
        ref_structure,
748 -
        force_constant,
749 -
        continuous_apr=True,
750 -
        amber_index=False,
770 +
    restraint_mask_list,
771 +
    num_window_list,
772 +
    ref_structure,
773 +
    force_constant,
774 +
    continuous_apr=True,
775 +
    amber_index=False,
751 776
):
752 777
    """
753 778
    Create a restraint whose value does not change during a calculation.
@@ -874,7 +899,9 @@
Loading
874 899
        logger.debug('All restraints are not "continuous_apr" style.')
875 900
        all_continuous_apr = False
876 901
    else:
877 -
        raise ValueError("All restraints must have the same setting for ``.continuous_apr``.")
902 +
        raise ValueError(
903 +
            "All restraints must have the same setting for ``.continuous_apr``."
904 +
        )
878 905
879 906
    window_list = []
880 907
    phases = ["attach", "pull", "release"]
@@ -902,27 +929,27 @@
Loading
902 929
903 930
                if phase == "attach" and all_continuous_apr:
904 931
                    window_list += [
905 -
                        phase[0]+str("{:03.0f}".format(val))
906 -
                        for val in np.arange(0, max_count-1, 1)
932 +
                        phase[0] + str("{:03.0f}".format(val))
933 +
                        for val in np.arange(0, max_count - 1, 1)
907 934
                    ]
908 935
                elif phase == "attach" and not all_continuous_apr:
909 936
                    window_list += [
910 -
                        phase[0]+str("{:03.0f}".format(val))
937 +
                        phase[0] + str("{:03.0f}".format(val))
911 938
                        for val in np.arange(0, max_count, 1)
912 939
                    ]
913 940
                elif phase == "pull":
914 941
                    window_list += [
915 -
                        phase[0]+str("{:03.0f}".format(val))
942 +
                        phase[0] + str("{:03.0f}".format(val))
916 943
                        for val in np.arange(0, max_count, 1)
917 944
                    ]
918 945
                elif phase == "release" and all_continuous_apr:
919 946
                    window_list += [
920 -
                        phase[0]+str("{:03.0f}".format(val))
947 +
                        phase[0] + str("{:03.0f}".format(val))
921 948
                        for val in np.arange(1, max_count, 1)
922 949
                    ]
923 950
                elif phase == "release" and not all_continuous_apr:
924 951
                    window_list += [
925 -
                        phase[0]+str("{:03.0f}".format(val))
952 +
                        phase[0] + str("{:03.0f}".format(val))
926 953
                        for val in np.arange(0, max_count, 1)
927 954
                    ]
928 955
        else:
@@ -938,6 +965,11 @@
Loading
938 965
                "phase.".format(phase)
939 966
            )
940 967
968 +
    # Check that each restraint have at least two atoms
969 +
    for restraint in restraint_list:
970 +
        if not restraint.index1 or not restraint.index2:
971 +
            raise Exception("There must be at least two atoms in a restraint.")
972 +
941 973
    logger.info("Restraints appear to be consistent")
942 974
943 975
    if create_window_list:
Files Coverage
paprika 77.22%
setup.py 16.08%
Project Totals (30 files) 74.04%
871.1
TRAVIS_OS_NAME=osx
871.2
TRAVIS_PYTHON_VERSION=3.6
TRAVIS_OS_NAME=linux
1
# Codecov configuration to make it a bit less noisy
2
coverage:
3
  status:
4
    patch: false
5
    project:
6
      default:
7
        threshold: 50%
8
comment:
9
  layout: "header"
10
  require_changes: false
11
  branches: null
12
  behavior: default
13
  flags: null
14
  paths: null
Sunburst
The inner-most circle is the entire project, moving away from the center are folders then, finally, a single file. The size and color of each slice is representing the number of statements and the coverage, respectively.
Icicle
The top section represents the entire project. Proceeding with folders and finally individual files. The size and color of each slice is representing the number of statements and the coverage, respectively.
Grid
Each block represents a single file in the project. The size and color of each block is represented by the number of statements and the coverage, respectively.
Loading