slochower / pAPRika
Showing 2 of 5 files from the diff.

@@ -395,9 +395,14 @@
Loading
395 395
        self._waters_added_history = [0]
396 396
        self._write_save_lines = True
397 397
398 -
    def build(self):
398 +
    def build(self, clean_files: bool = True):
399 399
        """
400 400
        Build the ``TLeap`` system.
401 +
402 +
        Parameters
403 +
        ----------
404 +
        clean_files : bool, optional
405 +
            Whether to delete log files after completion.
401 406
        """
402 407
403 408
        log.debug("Running tleap.build() in {}".format(self.output_path))
@@ -426,6 +431,14 @@
Loading
426 431
        else:
427 432
            self.solvate()
428 433
434 +
        # Cleanup
435 +
        if clean_files:
436 +
            os.remove(os.path.join(self.output_path, self.output_prefix + ".tleap.in"))
437 +
            os.remove(os.path.join(self.output_path, "leap.log"))
438 +
439 +
            if self.water_model["frcmod"] == "bind3p":
440 +
                os.remove(os.path.join(self.output_path, "frcmod.bind3p"))
441 +
429 442
    def filter_template(self):
430 443
        """
431 444
        Filter out any ``template_lines`` that may interfere with solvation.
@@ -468,6 +481,23 @@
Loading
468 481
469 482
        self.template_lines = filtered_lines
470 483
484 +
        # Add the water leaprc library to template_lines
485 +
        if self.pbc_type and self.water_model:
486 +
            # We need to source the water model after all other modules are loaded
487 +
            # because some modules will overwrite the water model definition.
488 +
            # Example case: source leaprc.protein.ff14SB contains tip3p parameters so
489 +
            # we need to load this first before the water frcmod file.
490 +
            insert_index = 0
491 +
            for idx, line in enumerate(self.template_lines):
492 +
                if "source" in line:
493 +
                    insert_index = idx + 1
494 +
            self.template_lines.insert(
495 +
                insert_index, f"loadamberparams frcmod.{self.water_model['frcmod']}"
496 +
            )
497 +
            self.template_lines.insert(
498 +
                insert_index, f"source leaprc.water.{self.water_model['library']}"
499 +
            )
500 +
471 501
    def write_input(self):
472 502
        """
473 503
        Write a ``TLeap`` input file based on ``template_lines`` and other things we have set.
@@ -480,12 +510,32 @@
Loading
480 510
            except OSError:
481 511
                raise
482 512
483 -
        with open(file_path, "w") as f:
484 -
            # load the water leaprc library
485 -
            if self.pbc_type and self.water_model:
486 -
                f.write(f"source leaprc.water.{self.water_model['library']}\n")
487 -
                f.write(f"loadamberparams frcmod.{self.water_model['frcmod']}\n")
513 +
        # If Bind3P, write frcmod.bind3p file
514 +
        if self.water_model["frcmod"] == "bind3p":
515 +
            frcmod_path = os.path.join(self.output_path, "frcmod.bind3p")
516 +
            with open(frcmod_path, "w") as f:
517 +
                f.write(
518 +
                    """\
519 +
This is the additional/replacement parameter set for Bind3P water
520 +
MASS
521 +
OW    16.0
522 +
HW     1.008   0.000
523 +
524 +
BOND
525 +
OW-HW  553.0    0.9572      Bind3P water
526 +
HW-HW  553.0    1.5136      Bind3P water
527 +
528 +
ANGLE
488 529
530 +
DIHE
531 +
532 +
NONBON
533 +
  OW          1.7577  0.1818             Bind3P water model
534 +
  HW          0.0000  0.0000             Bind3P water model
535 +
"""
536 +
                )
537 +
538 +
        with open(file_path, "w") as f:
489 539
            for line in self.template_lines:
490 540
                f.write(line + "\n")
491 541
@@ -1392,11 +1442,11 @@
Loading
1392 1442
            The water model to use, models supported are ["spc", "opc", "tip3p", "tip4p"].
1393 1443
            Strings are case-insensitive.
1394 1444
        model_type: str
1395 -
            The particular type of the water model, default is None and the key in parenthesis is the water box
1396 -
            used in TLeap. Strings are case-insensitive.
1445 +
            The particular type of the water model, default is ``None`` and the key in parenthesis is the water box
1446 +
            used in ``TLeap``. Strings are case-insensitive.
1397 1447
            * spc: None (SPCBOX), "flexible" (SPCFWBOX), "quantum" (QSPCFWBOX)
1398 1448
            * opc: None (OPCBOX), "three-point" (OPC3BOX)
1399 -
            * tip3p: None (TIP3PBOX), "flexible" (TIP3PFBOX), "force-balance" (FB3BOX)
1449 +
            * tip3p: None (TIP3PBOX), "flexible" (TIP3PFBOX), "force-balance" (FB3BOX), "bind3p" (TIP3PBOX)
1400 1450
            * tip4p: None (TIP4PBOX), "ewald" (TIP4PEWBOX), "force-balance" (FB4BOX)
1401 1451
        """
1402 1452
        if model.lower() not in ["spc", "opc", "tip3p", "tip4p"]:
@@ -1435,9 +1485,11 @@
Loading
1435 1485
1436 1486
        if model.lower() == "tip3p":
1437 1487
            library = "tip3p"
1488 +
            water_box = "TIP3PBOX"
1438 1489
            if model_type is None:
1439 1490
                frcmod = "tip3p"
1440 -
                water_box = "TIP3PBOX"
1491 +
            elif model_type.lower() == "bind3p":
1492 +
                frcmod = "bind3p"
1441 1493
            elif model_type.lower() == "force-balance":
1442 1494
                frcmod = "tip3pfb"
1443 1495
                water_box = "FB3BOX"

@@ -10,6 +10,8 @@
Loading
10 10
import numpy as np
11 11
import parmed as pmd
12 12
import pytest
13 +
import simtk.unit as unit
14 +
from simtk.openmm import NonbondedForce
13 15
14 16
from paprika.build.align import zalign
15 17
from paprika.build.system import TLeap
@@ -102,6 +104,50 @@
Loading
102 104
    assert int(grepped_waters) == waters
103 105
104 106
107 +
def test_solvation_bind3p(clean_files):
108 +
    logger.debug("Solvating with 500 Bind3P waters in a truncated octahedron...")
109 +
    sys = TLeap()
110 +
    sys.output_path = "tmp"
111 +
    sys.loadpdb_file = os.path.join(
112 +
        os.path.dirname(__file__), "../data/cb6-but/cb6-but.pdb"
113 +
    )
114 +
    sys.target_waters = 2000
115 +
    sys.output_prefix = "solvate"
116 +
    sys.pbc_type = PBCBox.cubic
117 +
    sys.set_water_model("tip3p", model_type="bind3p")
118 +
    sys.template_lines = [
119 +
        "source leaprc.gaff",
120 +
        "loadamberparams ../paprika/data/cb6-but/cb6.frcmod",
121 +
        "CB6 = loadmol2 ../paprika/data/cb6-but/cb6.mol2",
122 +
        "BUT = loadmol2 ../paprika/data/cb6-but/but.mol2",
123 +
        "model = loadpdb ../paprika/data/cb6-but/cb6-but.pdb",
124 +
    ]
125 +
    sys.build()
126 +
127 +
    prmtop = os.path.join("./tmp/solvate.prmtop")
128 +
    inpcrd = os.path.join("./tmp/solvate.rst7")
129 +
    structure = pmd.load_file(prmtop, inpcrd, structure=True)
130 +
    water_index = [
131 +
        atom.index
132 +
        for atom in structure.topology.atoms()
133 +
        if atom.residue.name == "WAT" and atom.residue.index == 10
134 +
    ]
135 +
    system = structure.createSystem()
136 +
    nonbonded = [
137 +
        force for force in system.getForces() if isinstance(force, NonbondedForce)
138 +
    ][0]
139 +
    oxygen = nonbonded.getParticleParameters(water_index[0])
140 +
    assert (
141 +
        pytest.approx(oxygen[1].value_in_unit(unit.angstrom) - 3.1319, abs=1e-3) == 0.0
142 +
    )
143 +
    assert (
144 +
        pytest.approx(
145 +
            oxygen[2].value_in_unit(unit.kilocalorie_per_mole) - 0.1818, abs=1e-3
146 +
        )
147 +
        == 0.0
148 +
    )
149 +
150 +
105 151
@pytest.mark.slow
106 152
def test_solvation_spatial_size(clean_files):
107 153
    """ Test that we can solvate CB6-BUT with an buffer size in Angstroms. """
@@ -214,7 +260,7 @@
Loading
214 260
    sys.neutralize = False
215 261
    sys.pbc_type = PBCBox.rectangular
216 262
    sys.add_ions = ["NA", "0.150M", "CL", "0.150M", "K", "0.100m", "BR", "0.100m"]
217 -
    sys.build()
263 +
    sys.build(clean_files=False)
218 264
219 265
    # Molarity Check
220 266
    obs_num_na = sp.check_output(
@@ -344,7 +390,7 @@
Loading
344 390
    sys.output_prefix = "multi"
345 391
    sys.pbc_type = None
346 392
    sys.neutralize = False
347 -
    sys.build()
393 +
    sys.build(clean_files=False)
348 394
349 395
    with open(f"{temporary_directory}/multi.tleap.in", "r") as f:
350 396
        lines = f.read()
Files Coverage
paprika 75.14%
Project Totals (38 files) 75.14%
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