Compare 8fec922 ... +0 ... 62beac9


@@ -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])

@@ -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,76 +40,48 @@
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,"
53 53
        igr1 = ""
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,"
64 61
        igr2 = ""
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])

@@ -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

@@ -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

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Learn more Showing 5 files with coverage changes found.

Changes in paprika/tests/test_restraints.py
+195
Loading file...
New file paprika/restraints/plumed.py
New
Loading file...
New file paprika/tests/test_dummy.py
New
Loading file...
Changes in paprika/utils.py
-1
+1
Loading file...
Changes in paprika/restraints/utils.py
+28
+4
Loading file...
Files Coverage
paprika -2.68% 77.22%
setup.py 16.08%
Project Totals (30 files) 74.04%
Loading