sire/tests/io/test_ambercmap.py at 3da5fd3052b2bd72c8a94e6a7ad3fe20c19fda40 · OpenBioSim/sire · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
import sire as sr

import pytest


def test_amber_cmap(tmpdir, amber_cmap):
    mols = amber_cmap

    cmap = mols[0].property("cmap")

    params = cmap.parameters()

    # there are 18 cmap parameters in this file
    assert len(params) == 18

    unique_params = []

    for param in cmap.parameters():
        if param.parameter() not in unique_params:
            unique_params.append(param.parameter())

    # there are 9 different terms from this file
    assert len(unique_params) == 9

    # write the file and re-read it
    dir = tmpdir.mkdir("test_amber_cmap")

    f = sr.save(mols, dir.join("output"), format="prm7")

    mols2 = sr.load(f)

    cmap2 = mols2[0].property("cmap")

    params2 = cmap2.parameters()

    assert len(params2) == 18

    unique_params2 = []

    for param in cmap2.parameters():
        if param.parameter() not in unique_params2:
            unique_params2.append(param.parameter())

    assert len(unique_params2) == 9

    # check that the parameters are the same
    for param in params:
        assert param in params2

    # write again - this should produce an identical file,
    # as we made sure to sort indexes and parameters in
    # a consistent way
    f2 = sr.save(mols2, dir.join("output2"), format="prm7")

    with open(f[0], "r") as f1:
        with open(f2[0], "r") as f2:
            for i, (line1, line2) in enumerate(zip(f1, f2)):
                # skip the first line since it includes a timestamp
                if i > 0:
                    assert line1 == line2


def test_amber_multichain_cmap(tmpdir, multichain_cmap):
    """Regression test for multi-chain CMAP write bug.

    When a topology has more than one molecule with CMAP terms, the write path
    previously applied a spurious '/= 3' to the per-molecule atom offset,
    producing wrong global atom indices for molecules 2, 3, … .  The re-parse
    of the generated text then raised:
      "there is a cmap between more than one different molecule"
    """
    mols = multichain_cmap

    # Collect per-molecule CMAP term counts indexed by position in the
    # molecule list. There must be at least two molecules with CMAP
    # to trigger the multi-chain bug.
    cmap_counts = {}
    for i, mol in enumerate(mols.molecules()):
        if mol.has_property("cmap"):
            cmap_counts[i] = len(mol.property("cmap").parameters())

    assert (
        len(cmap_counts) >= 2
    ), "Expected at least two molecules with CMAP terms in this topology"

    dir = tmpdir.mkdir("test_amber_multichain_cmap")

    # Write the topology back to file. (This is where the bug used to trigger.)
    f = sr.save(mols, dir.join("output"), format="prm7")

    # Reead the file back in.
    mols2 = sr.load(f)

    # Each molecule must still have the same number of CMAP terms after the
    # roundtrip.
    for i, count in cmap_counts.items():
        mol2 = mols2[i]
        assert mol2.has_property(
            "cmap"
        ), f"Molecule at index {i} lost its cmap property after roundtrip"
        count2 = len(mol2.property("cmap").parameters())
        assert (
            count2 == count
        ), f"Molecule at index {i}: CMAP count changed from {count} to {count2}"

    # Verify a second write also succeeds without error.
    sr.save(mols2, dir.join("output2"), format="prm7")


def test_amber_cmap_grotop(tmpdir, amber_cmap):
    """Testing reading and writing from amber to gromacs and back to amber."""
    mols = amber_cmap.clone()

    dir = tmpdir.mkdir("test_amber_cmap_grotop")

    # Save to a temporary file in GroTop format.
    f = sr.save(mols, dir.join("output"), format="GroTop")

    # Load the saved file.
    mols2 = sr.load(f, show_warnings=False)

    # Save back to prmtop format.
    f2 = sr.save(mols2, dir.join("output"), format="prmtop")

    # Load the saved file again.
    mols3 = sr.load(f2, show_warnings=False)

    # check that all of the cmap parameters are the same
    cmaps = mols[0].property("cmap").parameters()
    cmaps2 = mols2[0].property("cmap").parameters()
    cmaps3 = mols3[0].property("cmap").parameters()

    for cmap in cmaps:
        # find the same cmap in the other two molecules
        found = False

        for cmap2 in cmaps2:
            if (
                cmap.atom0() == cmap2.atom0()
                and cmap.atom1() == cmap2.atom1()
                and cmap.atom2() == cmap2.atom2()
                and cmap.atom3() == cmap2.atom3()
                and cmap.atom4() == cmap2.atom4()
            ):
                # make sure that the parameter is the same
                for val, val2 in zip(
                    cmap.parameter().values(), cmap2.parameter().values()
                ):
                    assert val == pytest.approx(val2, 1e-3)
                found = True

        assert found

        found = False

        for cmap3 in cmaps3:
            if (
                cmap.atom0() == cmap3.atom0()
                and cmap.atom1() == cmap3.atom1()
                and cmap.atom2() == cmap3.atom2()
                and cmap.atom3() == cmap3.atom3()
                and cmap.atom4() == cmap3.atom4()
            ):
                # make sure that the parameter is the same
                for val, val3 in zip(
                    cmap.parameter().values(), cmap3.parameter().values()
                ):
                    assert val == pytest.approx(val3, 1e-3)
                found = True

        assert found