#133 Add support for Plumed-based restraints

Tests the restraints utilities.
"""

import pytest
import os

import pytest

from paprika.restraints.restraints import *
from paprika.restraints.utils import *

logger = logging.getLogger(__name__)


@@ -17,7 +19,9 @@

    rest1.amber_index = True
    rest1.continuous_apr = False
    rest1.auto_apr = False
    rest1.topology = os.path.join(os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb")
    rest1.topology = os.path.join(
        os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb"
    )
    rest1.mask1 = ":CB6@O,O2,O4,O6,O8,O10"
    rest1.mask2 = ":BUT@C3"
    rest1.attach["target"] = 3.0

@@ -73,7 +77,9 @@

    rest2.amber_index = True
    rest2.continuous_apr = False
    rest2.auto_apr = False
    rest2.topology = os.path.join(os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb")
    rest2.topology = os.path.join(
        os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb"
    )
    rest2.mask1 = ":CB6@O,O2,O4,O6,O8,O10"
    rest2.mask2 = ":BUT@C3"
    rest2.mask3 = ":BUT@C"

@@ -131,7 +137,9 @@

    rest3.amber_index = True
    rest3.continuous_apr = False
    rest3.auto_apr = True
    rest3.topology = os.path.join(os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb")
    rest3.topology = os.path.join(
        os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb"
    )
    rest3.mask1 = ":CB6@O2"
    rest3.mask2 = ":CB6@O"
    rest3.mask3 = ":BUT@C3"

@@ -188,7 +196,9 @@

    rest4.amber_index = True
    rest4.continuous_apr = False
    rest4.auto_apr = False
    rest4.topology = os.path.join(os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb")
    rest4.topology = os.path.join(
        os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb"
    )
    rest4.mask1 = ":CB6@O2"
    rest4.mask2 = ":CB6@O"
    rest4.mask3 = ":BUT@C3"

@@ -243,7 +253,9 @@

    rest5.amber_index = True
    rest5.continuous_apr = False
    rest5.auto_apr = False
    rest5.topology = os.path.join(os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb")
    rest5.topology = os.path.join(
        os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb"
    )
    rest5.mask1 = ":CB6@O,O2,O4,O6,O8,O10"
    rest5.mask2 = ":BUT@C*"
    rest5.attach["target"] = 0.0

@@ -295,7 +307,9 @@

    rest6.amber_index = True
    rest6.continuous_apr = False
    rest6.auto_apr = False
    rest6.topology = os.path.join(os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb")
    rest6.topology = os.path.join(
        os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb"
    )
    rest6.mask1 = ":CB6@O,O2,O4,O6,O8,O10"
    rest6.mask2 = ":BUT@C*"
    rest6.attach["target"] = 0.0

@@ -352,7 +366,9 @@

    rest7.amber_index = True
    rest7.continuous_apr = True
    rest7.auto_apr = False
    rest7.topology = os.path.join(os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb")
    rest7.topology = os.path.join(
        os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb"
    )
    rest7.mask1 = ":1@O,O1,:BUT@H1"
    rest7.mask2 = ":CB6@N"
    rest7.attach["target"] = 0.0

@@ -400,7 +416,9 @@

    rest8.amber_index = True
    rest8.continuous_apr = False
    rest8.auto_apr = False
    rest8.topology = os.path.join(os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb")
    rest8.topology = os.path.join(
        os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb"
    )
    rest8.mask1 = ":CB6@O"
    rest8.mask2 = ":BUT@C3"
    rest8.attach["target"] = 0.0

@@ -429,7 +447,9 @@

    rest9.amber_index = True
    rest9.continuous_apr = False
    rest9.auto_apr = False
    rest9.topology = os.path.join(os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb")
    rest9.topology = os.path.join(
        os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb"
    )
    rest9.mask1 = ":CB6@O"
    rest9.mask2 = ":BUT@C3"
    rest9.pull["fc"] = 3.0

@@ -458,7 +478,9 @@

    rest10.amber_index = True
    rest10.continuous_apr = False
    rest10.auto_apr = False
    rest10.topology = os.path.join(os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb")
    rest10.topology = os.path.join(
        os.path.dirname(__file__), "../data/cb6-but/cb6-but-notcentered.pdb"
    )
    rest10.mask1 = ":CB6@O"
    rest10.mask2 = ":BUT@C3"
    rest10.release["target"] = 0.0

@@ -487,4 +509,274 @@

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

import logging

from paprika.restraints.utils import parse_window
from paprika.restraints.utils import get_restraint_values, parse_window

logger = logging.getLogger(__name__)


@@ -40,76 +40,48 @@

    """

    window, phase = parse_window(window)
    if restraint.phase[phase]["force_constants"] is None and restraint.phase[phase]["targets"] is None:
    if (
        restraint.phase[phase]["force_constants"] is None
        and restraint.phase[phase]["targets"] is None
    ):
        return ""

    if not restraint.index1:
        iat1 = " "
        raise Exception("There must be at least two atoms in a restraint.")
    elif not restraint.group1:
    if not restraint.group1:
        iat1 = "{},".format(restraint.index1[0])
    else:
        iat1 = "-1,"
        igr1 = ""
        for index in restraint.index1:
            igr1 += "{},".format(index)

    if not restraint.index2:
        iat2 = " "
        raise Exception("There must be at least two atoms in a restraint.")
    elif not restraint.group2:
    if not restraint.group2:
        iat2 = "{},".format(restraint.index2[0])
    else:
        iat2 = "-1,"
        igr2 = ""
        for index in restraint.index2:
            igr2 += "{},".format(index)

    if not restraint.index3:
        iat3 = ""
    elif not restraint.group3:
    iat3 = ""
    if restraint.index3 and not restraint.group3:
        iat3 = "{},".format(restraint.index3[0])
    else:
    elif restraint.group3:
        iat3 = "-1,"
        igr3 = ""
        for index in restraint.index3:
            igr3 += "{},".format(index)

    if not restraint.index4:
        iat4 = ""
    elif not restraint.group4:
    iat4 = ""
    if restraint.index4 and not restraint.group4:
        iat4 = "{},".format(restraint.index4[0])
    else:
    elif restraint.group4:
        iat4 = "-1,"
        igr4 = ""
        for index in restraint.index4:
            igr4 += "{},".format(index)


    # Set upper/lower bounds depending on whether distance, angle, or torsion
    lower_bound = 0.0
    upper_bound = 999.0

    if restraint.mask3 and not restraint.mask4:
        upper_bound = 180.0

    if restraint.mask3 and restraint.mask4:
        lower_bound = restraint.phase[phase]["targets"][window] - 180.0
        upper_bound = restraint.phase[phase]["targets"][window] + 180.0

    amber_restraint_values = {
        "r1": lower_bound,
        "r2": restraint.phase[phase]["targets"][window],
        "r3": restraint.phase[phase]["targets"][window],
        "r4": upper_bound,
        "rk2": restraint.phase[phase]["force_constants"][window],
        "rk3": restraint.phase[phase]["force_constants"][window],
    }

    for key, value in restraint.custom_restraint_values.items():
        if value is not None:
            logger.debug("Overriding {} = {}".format(key, value))
            amber_restraint_values[key] = value
    # Restraint values - Amber NMR-style
    amber_restraint_values = get_restraint_values(restraint, phase, window)

    # Prepare AMBER NMR-style restraint
    atoms = "".join([iat1, iat2, iat3, iat4])

@@ -2,13 +2,14 @@

import os
import subprocess as sp
from collections import OrderedDict
from sys import platform

logger = logging.getLogger(__name__)


class Simulation(object):
    """
    A wrapper that can be used to set AMBER simulation parameters .
    A wrapper that can be used to set AMBER simulation parameters.
    """

    @property

@@ -34,13 +35,13 @@

        self._executable = value

    @property
    def CUDA_VISIBLE_DEVICES(self):
    def gpu_devices(self):
        """A wrapper around the environmental variable ``CUDA_VISIBLE_DEVICES``."""
        return self._CUDA_VISIBLE_DEVICES
        return self._gpu_devices

    @CUDA_VISIBLE_DEVICES.setter
    def CUDA_VISIBLE_DEVICES(self, value):
        self._CUDA_VISIBLE_DEVICES = value
    @gpu_devices.setter
    def gpu_devices(self, value):
        self._gpu_devices = value

    @property
    def phase(self):

@@ -75,13 +76,44 @@


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

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

    @restraint_file.setter
    def restraint_file(self, value):
        self._restraint_file = value

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

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

    @plumed_file.setter
    def plumed_file(self, value: str):
        self._plumed_file = value

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

    @plumed_kernel_library.setter
    def plumed_kernel_library(self, value: str):
        self._plumed_kernel_library = value

    @property
    def converged(self) -> bool:
        """bool: Whether the simulation is converged.

@@ -105,14 +137,14 @@

    @prefix.setter
    def prefix(self, new_prefix):
        self._prefix = new_prefix
        self.input = new_prefix+".in"
        self.inpcrd = new_prefix+".inpcrd"
        self.ref = new_prefix+".inpcrd"
        self.output = new_prefix+".out"
        self.restart = new_prefix+".rst7"
        self.mdinfo = new_prefix+".mdinfo"
        self.mdcrd = new_prefix+".nc"
        self.mden = new_prefix+".mden"
        self.input = new_prefix + ".in"
        self.inpcrd = new_prefix + ".inpcrd"
        self.ref = new_prefix + ".inpcrd"
        self.output = new_prefix + ".out"
        self.restart = new_prefix + ".rst7"
        self.mdinfo = new_prefix + ".mdinfo"
        self.mdcrd = new_prefix + ".nc"
        self.mden = new_prefix + ".mden"

    @property
    def cntrl(self):

@@ -208,24 +240,26 @@

        # I/O
        self._path = "."
        self._executable = "sander"
        self._CUDA_VISIBLE_DEVICES = None
        self._gpu_devices = None
        self._phase = None
        self._window = None
        self._topology = "prmtop"
        self._restraint_file = "restraints.in"
        self._restraint_file = None
        self._plumed_file = None
        self._plumed_kernel_library = None
        self.title = "PBC MD Simulation"
        self.converged = False

        # File names
        self._prefix = "md"
        self.input = self._prefix+".in"
        self.inpcrd = self._prefix+".inpcrd"
        self.ref = self._prefix+".inpcrd"
        self.output = self._prefix+".out"
        self.restart = self._prefix+".rst7"
        self.mdinfo = self._prefix+".mdinfo"
        self.mdcrd = self._prefix+".nc"
        self.mden = self._prefix+".mden"
        self.input = self._prefix + ".in"
        self.inpcrd = self._prefix + ".inpcrd"
        self.ref = self._prefix + ".inpcrd"
        self.output = self._prefix + ".out"
        self.restart = self._prefix + ".rst7"
        self.mdinfo = self._prefix + ".mdinfo"
        self.mdcrd = self._prefix + ".nc"
        self.mden = self._prefix + ".mden"

        # Input file cntrl settings (Default = NTP)
        self._cntrl = OrderedDict()

@@ -373,7 +407,13 @@


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

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

        f.write(" /\n")

    def _amber_write_input_file(self):

@@ -394,11 +434,14 @@

            if self.cntrl["nmropt"] == 1:
                if self.wt is not None:
                    for line in self.wt:
                        f.write(" "+line+"\n")
                        f.write(" " + line + "\n")
                f.write(" &wt type = 'END', /\n")

                # Specify Amber NMR file
                if self.restraint_file is not None:
                    f.write("DISANG = {}\n".format(self.restraint_file))
                    f.write("LISTOUT = POUT\n\n")

            if self.group is not None:
                f.write("{:s}".format(self.group))


@@ -425,13 +468,13 @@

                # maxcyc
                ncyc = self.cntrl["ncyc"]
                maxcyc = self.cntrl["maxcyc"]
                burn_in = int(float(ncyc)+0.20 * (float(maxcyc)-float(ncyc)))
                burn_in = int(float(ncyc) + 0.20 * (float(maxcyc) - float(ncyc)))
                # If the burn_in value is nuts, then just set it to zero
                if burn_in < 0 or burn_in >= maxcyc:
                    burn_in = 0
                # Set an end_soft value that is 75% of way between ncyc and
                # maxcyc
                end_soft = int(float(ncyc)+0.60 * (float(maxcyc)-float(ncyc)))
                end_soft = int(float(ncyc) + 0.60 * (float(maxcyc) - float(ncyc)))
                self.wt = [
                    "&wt type = 'NB', istep1=0, istep2={:.0f}, value1 = 0.0, value2=0.0, IINC=50, /".format(
                        burn_in

@@ -441,6 +484,12 @@

                    ),
                ]

            # Check restraints file
            if self.restraint_file and self.plumed_file:
                raise Exception(
                    "Cannot use both NMR-style and Plumed-style restraints at the same time."
                )

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


@@ -450,7 +499,7 @@

                logger.info("Running MD at {}".format(self.path))

            # Create executable list for subprocess
            exec_list = self.executable.split()+["-O", "-p", self.topology]
            exec_list = self.executable.split() + ["-O", "-p", self.topology]
            if self.ref is not None:
                exec_list += ["-ref", self.ref]
            exec_list += [

@@ -470,30 +519,60 @@

            if self.mden is not None:
                exec_list += ["-e", self.mden]

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

            # Execute
            if self.CUDA_VISIBLE_DEVICES:
                amber_output = sp.Popen(
                    exec_list,
                    cwd=self.path,
                    stdout=sp.PIPE,
                    stderr=sp.PIPE,
                    env=dict(
                        os.environ, CUDA_VISIBLE_DEVICES=str(self.CUDA_VISIBLE_DEVICES)
                    ),
            # Add CUDA_VISIBLE_DEVICES variable to env dict
            if self.gpu_devices is not None:
                os.environ = dict(
                    os.environ, CUDA_VISIBLE_DEVICES=str(self.gpu_devices)
                )

            # Set the Plumed kernel library path
            if self.plumed_kernel_library is None:
                # using "startswith" as recommended from https://docs.python.org/3/library/sys.html#sys.platform
                if platform.startswith("linux"):
                    self.plumed_kernel_library = os.path.join(
                        os.environ["CONDA_PREFIX"], "lib", "libplumedKernel.so"
                    )
                elif platform.startswith("darwin"):
                    self.plumed_kernel_library = os.path.join(
                        os.environ["CONDA_PREFIX"], "lib", "libplumedKernel.dylib"
                    )
                else:
                    logger.warning(
                        f"Plumed is not supported with the '{platform}' platform."
                    )

            # Add Plumed kernel library to env dict
            if os.path.isfile(self.plumed_kernel_library):
                os.environ = dict(os.environ, PLUMED_KERNEL=self.plumed_kernel_library)
            else:
                amber_output = sp.Popen(
                    exec_list, cwd=self.path, stdout=sp.PIPE, stderr=sp.PIPE
                logger.warning(
                    f"Plumed kernel library {self.plumed_kernel_library} does not exist, "
                    "PLUMED_KERNEL environment variable is not set."
                )

            amber_output = amber_output.stdout.read().splitlines()
            # Execute
            amber_output = sp.Popen(
                exec_list,
                cwd=self.path,
                stdout=sp.PIPE,
                stderr=sp.PIPE,
                env=os.environ,
            )

            amber_stdout = amber_output.stdout.read().splitlines()
            amber_stderr = amber_output.stderr.read().splitlines()

            # Report any stdout/stderr which are output from execution
            if amber_output:
                logger.info("STDOUT/STDERR received from AMBER execution")
                for line in amber_output:
            if amber_stdout:
                logger.info("STDOUT received from AMBER execution")
                for line in amber_stdout:
                    logger.info(line)

            if amber_stderr:
                logger.info("STDERR received from AMBER execution")
                for line in amber_stderr:
                    logger.info(line)

            # Check completion status

@@ -0,0 +1,42 @@

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

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

→

Loading file...

Changes in paprika/restraints/utils.py

→

+28

Loading file...

74.04% 50.00% -0.01%

add exception test for restraint files

jeff231li

#133 eec6940

CI Passed

Hiding 2 contexual commits

-2

-3

73.97% 100.00% +0.26%

add more cases for get_bias_potential_type

jeff231li

#133 cabf1f6

CI Passed

+38

73.71% ø -0.06%

fix info about plumed kernel library

jeff231li

#133 c9a70b5

CI Passed

Hiding 3 contexual commits

+25

+16

73.78% 59.26% +0.24%

moved restraint checking outside of restraints/plumed.py & restraints/amber.py

jeff231li

#133 689decf

CI Passed

-7

-11

#133 2a16292

Hiding 2 contexual commits

+35

-1

+36

#133 e5718ca

Hiding 7 contexual commits

+1 Files

+193

+219

-26

#133 cb43afd

-40

#133 82a211b

Hiding 1 contexual commits

-1

71.26% ø +0.02%

add tutorial for plumed

jeff231li

#133 a26735c

CI Passed

Hiding 1 contexual commits

-1

50475e2

Hiding 7 contexual commits

+1 Files

+242

+10

+232

Pull Request Base Commit

76.10% ø ø

Merge pull request #132 from slochower/fix-docs

jeff231li

8fec922

CI Passed

Files		•	•	•	Coverage
paprika	+486	+291	ø	+195	^-2.68% 77.22%
setup.py	ø	ø	ø	ø	16.08%
Project Totals (30 files)	3,825	2,832	0	993	74.04%

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

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):

#133 Add support for Plumed-based restraints

No flags found

Show failed CI and contexual commits `c`

33 Commits

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

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

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

Navigation	Overlay
`t` Navigate files	`h` Toggle hits
`y` Change url to tip of branch	`m` Toggle misses
`b / v` Jump to prev/next hit line	`p` Toggle partial
`z / x` Jump to prev/next missed or partial line	`1..9` Toggle flags
`shift + o` Open current page in GitHub	`a` Toggle all on
`/ or ?` Show keyboard shortcuts dialog	`c` Toggle context lines or commits

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

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

#133 Add support for Plumed-based restraints

No flags found

Show failed CI and contexual commits c 33 Commits

Show failed CI and contexual commits `c`

33 Commits