1 38
import os
2 38
import shutil
3

4 38
import parmed as pmd
5 38
import pytest
6

7 38
from paprika import restraints
8 38
from paprika.restraints import amber
9 38
from paprika.restraints.restraints import create_window_list
10 38
from paprika.simulate import AMBER, GROMACS, NAMD
11 38
from paprika.tests import addons
12 38
from paprika.utils import parse_mden, parse_mdout
13

14

15 38
@pytest.fixture(scope="function", autouse=True)
16 38
def clean_files(directory="tmp"):
17
    # This happens before the test function call
18 38
    if os.path.isdir(directory):
19 0
        shutil.rmtree(directory)
20 38
    os.makedirs(directory)
21 38
    yield
22
    # This happens after the test function call
23 38
    shutil.rmtree(directory)
24

25

26 38
@addons.using_sander
27
# @addons.using_pmemd_cuda
28
def test_amber_single_window_gbmin(clean_files):
29
    # Distance restraint
30 38
    restraint = restraints.DAT_restraint()
31 38
    restraint.continuous_apr = True
32 38
    restraint.amber_index = True
33 38
    restraint.topology = os.path.join(
34
        os.path.dirname(__file__), "../data/k-cl/k-cl.pdb"
35
    )
36 38
    restraint.mask1 = ":K+"
37 38
    restraint.mask2 = ":Cl-"
38 38
    restraint.attach["target"] = 4.5
39 38
    restraint.attach["fraction_list"] = [0.00, 0.04, 0.181, 0.496, 1.000]
40 38
    restraint.attach["fc_final"] = 5.0
41 38
    restraint.pull["fc"] = restraint.attach["fc_final"]
42 38
    restraint.pull["target_initial"] = restraint.attach["target"]
43 38
    restraint.pull["target_final"] = 18.5
44 38
    restraint.pull["num_windows"] = 5
45 38
    restraint.initialize()
46

47 38
    windows_directory = os.path.join("tmp", "k-cl", "windows")
48 38
    window_list = create_window_list([restraint])
49

50 38
    for window in window_list:
51 38
        os.makedirs(os.path.join(windows_directory, window))
52 38
        with open(
53
            os.path.join(windows_directory, window, "restraints.in"), "a"
54
        ) as file:
55 38
            string = amber.amber_restraint_line(restraint, window)
56 38
            file.write(string)
57

58 38
    for window in window_list:
59 38
        if window[0] == "a":
60 38
            structure = pmd.load_file(
61
                os.path.join(os.path.dirname(__file__), "../data/k-cl/k-cl-sol.prmtop"),
62
                os.path.join(os.path.dirname(__file__), "../data/k-cl/k-cl-sol.rst7"),
63
                structure=True,
64
            )
65 38
            for atom in structure.atoms:
66 38
                if atom.name == "Cl-":
67 38
                    atom.xz = 2.65
68 38
            structure.save(
69
                os.path.join(windows_directory, window, "k-cl.prmtop"), overwrite=True
70
            )
71 38
            structure.save(
72
                os.path.join(windows_directory, window, "k-cl.rst7"), overwrite=True
73
            )
74 38
            structure.save(
75
                os.path.join(windows_directory, window, "k-cl.pdb"), overwrite=True
76
            )
77

78 38
        elif window[0] == "p":
79 38
            structure = pmd.load_file(
80
                os.path.join(os.path.dirname(__file__), "../data/k-cl/k-cl-sol.prmtop"),
81
                os.path.join(os.path.dirname(__file__), "../data/k-cl/k-cl-sol.rst7"),
82
                structure=True,
83
            )
84 38
            target_difference = (
85
                restraint.phase["pull"]["targets"][int(window[1:])]
86
                - restraint.phase["pull"]["targets"][0]
87
            )
88

89 38
            for atom in structure.atoms:
90 38
                if atom.name == "Cl-":
91 38
                    atom.xz += target_difference
92 38
            structure.save(
93
                os.path.join(windows_directory, window, "k-cl.prmtop"), overwrite=True
94
            )
95 38
            structure.save(
96
                os.path.join(windows_directory, window, "k-cl.rst7"), overwrite=True
97
            )
98 38
            structure.save(
99
                os.path.join(windows_directory, window, "k-cl.pdb"), overwrite=True
100
            )
101

102 38
    gbsim = AMBER()
103 38
    gbsim.path = os.path.join("tmp", "k-cl", "windows", "a003")
104 38
    gbsim.executable = "sander"
105 38
    gbsim.topology = "k-cl.prmtop"
106 38
    gbsim.prefix = "minimize"
107 38
    gbsim.coordinates = "k-cl.rst7"
108 38
    gbsim.config_gb_min()
109 38
    gbsim.cntrl["maxcyc"] = 1
110 38
    gbsim.cntrl["ncyc"] = 1
111 38
    gbsim.run()
112

113 38
    gbsim.config_gb_md()
114 38
    gbsim.prefix = "md"
115 38
    gbsim.coordinates = "minimize.rst7"
116 38
    gbsim.cntrl["nstlim"] = 1
117 38
    gbsim.cntrl["ntwe"] = 1
118 38
    gbsim.cntrl["ntpr"] = 1
119 38
    gbsim.cntrl["ig"] = 777
120 38
    gbsim.run()
121

122 38
    mden = parse_mden(os.path.join("tmp", "k-cl", "windows", "a003", "md.mden"))
123

124 38
    assert pytest.approx(mden["Bond"][0]) == 0
125 38
    assert pytest.approx(mden["Angle"][0]) == 0
126 38
    assert pytest.approx(mden["Dihedral"][0]) == 0
127 38
    assert pytest.approx(mden["V14"][0]) == 0
128 38
    assert pytest.approx(mden["E14"][0]) == 0
129

130 38
    assert pytest.approx(mden["VDW"][0], 0.1) == 25956.13225
131 38
    assert pytest.approx(mden["Ele"][0], 0.1) == -18828.99631
132 38
    assert pytest.approx(mden["Total"][0], 0.1) == 7127.13594
133

134

135 38
def test_amber_minimization(clean_files):
136 38
    simulation = AMBER()
137 38
    simulation.path = os.path.join("tmp")
138

139 38
    shutil.copy(
140
        os.path.join(os.path.dirname(__file__), "../data/k-cl/k-cl.prmtop"), "tmp"
141
    )
142 38
    shutil.copy(
143
        os.path.join(os.path.dirname(__file__), "../data/k-cl/k-cl.rst7"), "tmp"
144
    )
145

146 38
    simulation.executable = "sander"
147 38
    simulation.restraint_file = None
148

149 38
    simulation.prefix = "minimize"
150 38
    simulation.topology = "k-cl.prmtop"
151 38
    simulation.coordinates = "k-cl.rst7"
152

153 38
    simulation.config_gb_min()
154
    # Turn off GB for now.
155 38
    simulation.cntrl["igb"] = 0
156 38
    simulation.cntrl["ntb"] = 0
157

158 38
    simulation.run()
159

160 38
    mdout = parse_mdout(os.path.join("tmp", "minimize.out"))
161

162 38
    assert pytest.approx(mdout["Bond"][-1]) == 0
163 38
    assert pytest.approx(mdout["Angle"][-1]) == 0
164 38
    assert pytest.approx(mdout["Dihedral"][-1]) == 0
165 38
    assert pytest.approx(mdout["V14"][-1]) == 0
166 38
    assert pytest.approx(mdout["E14"][-1]) == 0
167

168 38
    assert pytest.approx(mdout["VDW"][0], 0.1) == 6.5734
169 38
    assert pytest.approx(mdout["Ele"][0], 0.1) == -211.7616
170

171

172 38
def test_gromacs_config(clean_files):
173
    # Test PBC
174 38
    simulation = GROMACS()
175 38
    simulation.path = os.path.join("tmp")
176 38
    simulation.prefix = "test"
177 38
    simulation.config_pbc_md(
178
        ensemble=GROMACS.Ensemble.NPT,
179
        integrator=GROMACS.Integrator.VelocityVerlet,
180
        thermostat=GROMACS.Thermostat.NoseHoover,
181
        barostat=GROMACS.Barostat.ParrinelloRahman,
182
    )
183 38
    simulation.control["nsteps"] = 1250000
184 38
    simulation._write_input_file()
185

186 38
    f = open(os.path.join(simulation.path, simulation.prefix + ".mdp"))
187 38
    lines = f.readlines()
188

189
    # Check if critical keys are written
190 38
    assert "nsteps" in "".join(lines)
191 38
    assert "tcoupl" in "".join(lines)
192 38
    assert "pcoupl" in "".join(lines)
193

194
    # Check if the values are correct
195 38
    for line in lines:
196 38
        if line.startswith("nsteps"):
197 38
            assert int(line.split()[-1]) == 1250000
198

199 38
        elif line.startswith("tcoupl"):
200 38
            assert line.split()[-1] == "nose-hoover"
201

202 38
        elif line.startswith("pcoupl") and not line.startswith("pcoupltype"):
203 38
            assert line.split()[-1] == "Parrinello-Rahman"
204

205 38
        elif line.startswith("pcoupltype"):
206 38
            assert line.split()[-1] == "isotropic"
207

208 38
        elif line.startswith("integrator"):
209 38
            assert line.split()[-1] == "md-vv"
210

211
    # Test for tc-groups
212 38
    simulation.tc_groups = ["CB8", "AMT", "HOH"]
213 38
    simulation._write_input_file()
214

215 38
    f = open(os.path.join(simulation.path, simulation.prefix + ".mdp"))
216 38
    lines = f.readlines()
217

218 38
    for line in lines:
219 38
        if line.startswith("tc-grps"):
220 38
            assert line.split()[2] == "CB8"
221 38
            assert line.split()[3] == "AMT"
222 38
            assert line.split()[4] == "HOH"
223

224 38
        elif line.startswith("tau-t"):
225 0
            assert len(line.split()) == 5
226

227 38
        elif line.startswith("ref-t"):
228 0
            assert len(line.split()) == 5
229

230
    # Test vacuum
231 38
    simulation = GROMACS()
232 38
    simulation.path = os.path.join("tmp")
233 38
    simulation.prefix = "test"
234 38
    simulation.config_vac_md(
235
        integrator=GROMACS.Integrator.LangevinDynamics,
236
        thermostat=GROMACS.Thermostat.VelocityRescaling,
237
    )
238 38
    simulation._write_input_file()
239

240 38
    f = open(os.path.join(simulation.path, simulation.prefix + ".mdp"))
241 38
    lines = f.readlines()
242

243 38
    assert "tcoupl" not in "".join(lines)
244

245 38
    for line in lines:
246 38
        if line.startswith("integrator"):
247 38
            assert line.split()[-1] == "sd"
248

249 38
        elif line.startswith("tau-t"):
250 0
            assert float(line.split()[-1]) == 0.1
251

252 38
        elif line.startswith("rcoulomb"):
253 38
            assert float(line.split()[-1]) == 333.3
254

255 38
        elif line.startswith("rvdw"):
256 38
            assert float(line.split()[-1]) == 333.3
257

258 38
        elif line.startswith("DispCorr"):
259 38
            assert line.split()[-1] == "no"
260

261
    # Test min
262 38
    simulation = GROMACS()
263 38
    simulation.path = os.path.join("tmp")
264 38
    simulation.prefix = "test"
265 38
    simulation.config_pbc_min(
266
        optimizer=GROMACS.Optimizer.ConjugateGradient,
267
    )
268 38
    simulation._write_input_file()
269

270 38
    f = open(os.path.join(simulation.path, simulation.prefix + ".mdp"))
271 38
    lines = f.readlines()
272

273 38
    for line in lines:
274 38
        if line.startswith("integrator"):
275 38
            assert line.split()[-1] == "cg"
276

277

278 38
def test_namd_config(clean_files):
279
    # Test PBC
280 38
    simulation = NAMD()
281 38
    simulation.topology = "k-cl.prmtop"
282 38
    simulation.coordinates = "k-cl.rst7"
283 38
    simulation.checkpoint = "equilibration"
284 38
    simulation.path = os.path.join("tmp")
285 38
    simulation.prefix = "test"
286 38
    simulation.config_pbc_md(
287
        ensemble=NAMD.Ensemble.NPT,
288
        thermostat=NAMD.Thermostat.LoweAnderson,
289
        barostat=NAMD.Barostat.Berendsen,
290
    )
291 38
    simulation._write_input_file()
292

293 38
    f = open(os.path.join(simulation.path, simulation.prefix + ".conf"))
294 38
    lines = f.readlines()
295

296 38
    assert "ambercoor" in "".join(lines)
297 38
    assert "parmfile" in "".join(lines)
298 38
    assert "readexclusions" in "".join(lines)
299 38
    assert "scnb" in "".join(lines)
300 38
    assert "bincoordinates" in "".join(lines)
301 38
    assert "binvelocities" in "".join(lines)
302 38
    assert "extendedSystem" in "".join(lines)
303 38
    assert "loweAndersen" in "".join(lines)
304

305 38
    for line in lines:
306 38
        if line.startswith("amber") and not line.startswith("ambercoor"):
307 38
            assert line.split()[-1] == "yes"
308

309 38
        elif line.startswith("readexclusions"):
310 38
            assert line.split()[-1] == "yes"
311

312 38
        elif line.startswith("scnb"):
313 38
            assert float(line.split()[-1]) == 2.0
314

315 38
        elif line.startswith("loweAndersenCutoff"):
316 38
            assert float(line.split()[-1]) == 2.7
317

318 38
        elif line.startswith("BerendsenPressureCompressibility"):
319 38
            assert float(line.split()[-1]) == 4.57e-5
320

321 38
        elif line.startswith("bincoordinates"):
322 38
            assert line.split()[-1] == simulation.checkpoint + ".coor"
323

324 38
        elif line.startswith("binvelocities"):
325 38
            assert line.split()[-1] == simulation.checkpoint + ".vel"
326

327 38
        elif line.startswith("extendedSystem"):
328 38
            assert line.split()[-1] == simulation.checkpoint + ".xsc"
329

330
    # Test vac
331 38
    simulation = NAMD()
332 38
    simulation.topology = "k-cl.prmtop"
333 38
    simulation.coordinates = "k-cl.rst7"
334 38
    simulation.checkpoint = "equilibration"
335 38
    simulation.path = os.path.join("tmp")
336 38
    simulation.prefix = "test"
337 38
    simulation.config_vac_md(
338
        thermostat=NAMD.Thermostat.Langevin,
339
    )
340 38
    simulation.control["run"] = 500000
341 38
    simulation._write_input_file()
342

343 38
    f = open(os.path.join(simulation.path, simulation.prefix + ".conf"))
344 38
    lines = f.readlines()
345

346 38
    assert "GBIS" in "".join(lines)
347

348 38
    for line in lines:
349 38
        if line.startswith("PME"):
350 38
            assert line.split()[-1] == "off"
351

352 38
        elif line.startswith("GBIS"):
353 38
            assert line.split()[-1] == "off"
354

355 38
        elif line.startswith("cutoff"):
356 38
            assert float(line.split()[-1]) == 999.0
357

358 38
        elif line.startswith("run"):
359 38
            assert int(line.split()[-1]) == 500000
360

361
    # Test GB
362 38
    simulation = NAMD()
363 38
    simulation.topology = "k-cl.prmtop"
364 38
    simulation.coordinates = "k-cl.rst7"
365 38
    simulation.checkpoint = "equilibration"
366 38
    simulation.path = os.path.join("tmp")
367 38
    simulation.prefix = "test"
368 38
    simulation.config_gb_md(
369
        thermostat=NAMD.Thermostat.Langevin,
370
    )
371 38
    simulation.implicit["ionConcentration"] = 0.15
372 38
    simulation._write_input_file()
373

374 38
    f = open(os.path.join(simulation.path, simulation.prefix + ".conf"))
375 38
    lines = f.readlines()
376

377 38
    assert "GBIS" in "".join(lines)
378 38
    assert "solventDielectric" in "".join(lines)
379 38
    assert "ionConcentration" in "".join(lines)
380

381 38
    for line in lines:
382 38
        if line.startswith("PME"):
383 38
            assert line.split()[-1] == "off"
384

385 38
        elif line.startswith("GBIS"):
386 38
            assert line.split()[-1] == "on"
387

388 38
        elif line.startswith("ionConcentration"):
389 38
            assert float(line.split()[-1]) == 0.15
390

391 38
        elif line.startswith("cutoff"):
392 38
            assert float(line.split()[-1]) == 999.0
393

394 38
        elif line.startswith("run"):
395 38
            assert int(line.split()[-1]) == 5000

Read our documentation on viewing source code .

Loading