Added a bunch of tests...

packages/Qpyl/common.py

@@ -83,7 +83,7 @@ def backup_file(filename):
         while os.path.lexists(os.path.join(di, backup_filename)):
             backup_filename = "#%s.%d#" % (fn, i)
             i += 1
-            if i > 100:
+            if i > 10:
                 logger.warning("You have more than 100 backed up files... "
                                "Cleaning time...")
         shutil.copy2(filename, os.path.join(di, backup_filename))

packages/Qpyl/core/qdyninp.py

@@ -110,7 +110,8 @@ def logical(x):
                         ("atom_restraints", list),
                         ("sequence_restraints", list),
                         ("distance_restraints", list),
-                        ("angle_restraints", list)
+                        ("angle_restraints", list),
+                        ("wall_restraints", list)
                       ]) 
 
 class QDynInputError(Exception):

packages/Qpyl/core/qpotential.py

@@ -193,7 +193,7 @@ def torsion_energy(phi, fc, periodicity, npaths, phi0):
     """
     phi = math.pi/180.0 * phi   # degrees to radians
     phi0 = math.pi/180.0 * phi0   # degrees to radians
-    return fc/npaths * (1 + math.cos(periodicity*phi - phi0))
+    return 1.0 * fc/npaths * (1 + math.cos(periodicity*phi - phi0))
 
 def improper_energy_periodic(phi, fc, periodicity, phi0):
     """

packages/Qpyl/qanalysis.py

@@ -483,6 +483,10 @@ def __init__(self, logfiles, timeunit="ps", stepsize=None):
             print colt, np.mean( [ x for j,x in enumerate(temps[i]) if temps[0][j] >= midpoint ] )
 
         """
+
+        logger.warning("QAnalyseDyns is under development. "
+                       "Use at your own risk!")
+
         if not logfiles:
             raise QAnalyseDynError("No logfiles given")
         self.analysed = []

qscripts-cli/q_analysedyns.py

@@ -34,11 +34,17 @@
 import argparse
 import os
 from collections import OrderedDict as ODict
+import logging
 
 from Qpyl.qanalysis import QAnalyseDyns, QAnalyseDynError
-from Qpyl.common import backup_file
+from Qpyl.common import backup_file, SpecialFormatter
 from Qpyl import plotdata
 
+logger = logging.getLogger('Qpyl')
+logger.setLevel(logging.INFO)
+handler = logging.StreamHandler(sys.stdout)
+handler.setFormatter(SpecialFormatter())
+logger.addHandler(handler)
 
 parser = argparse.ArgumentParser(description="""
 Script for extracting temperatures and energies from QDyn outputs.

tests/Qpyl/common_test.py

@@ -0,0 +1,38 @@
+#########################
+# py.test test functions
+#########################
+
+
+import pytest
+import shutil
+import os
+from Qpyl.common import np, backup_file
+
+def is_close(a, b, rel_tol=1e-09, abs_tol=0.0):
+    return abs(a-b) <= max(rel_tol * max(abs(a), abs(b)), abs_tol)
+
+
+def test_backupfile(tmpdir):
+    shutil.copy("data/qdyn.inp.1", tmpdir.dirname)
+    fn = backup_file(os.path.join(tmpdir.dirname, "qdyn.inp.1"))
+    inpstr = open("data/qdyn.inp.1").read()
+    inpstr2 = open(os.path.join(tmpdir.dirname, fn)).read()
+    assert inpstr == inpstr2
+
+
+class TestStats:
+    def test_mean(self):
+        vals = [i**2 for i in range(1, 21)]
+        assert is_close(np.mean(vals), 143.5)
+
+    def test_median(self):
+        vals = [i**2 for i in range(1, 21)]
+        assert is_close(np.median(vals), 110.5)
+
+    def test_std(self):
+        vals = [i**2 for i in range(1, 21)]
+        assert is_close(np.std(vals), 127.9023064686)
+
+    def test_sem(self):
+        vals = [i**2 for i in range(1, 21)]
+        assert is_close(np.std_error(vals), 28.5998251743)

tests/Qpyl/core/qcalc_test.py

@@ -6,6 +6,9 @@
 
 from Qpyl.core.qcalc import QCalcOutput
 
+def is_close(a, b, rel_tol=1e-09, abs_tol=0.0):
+    return abs(a-b) <= max(rel_tol * max(abs(a), abs(b)), abs_tol)
+
 
 class TestQFepOutput:
 
@@ -24,19 +27,19 @@ def test_header(self, qco1):
 
     def test_gc1(self, qco1):
         lj, el = qco1.results["gc"].get_columns(["E_LJ", "E_EL"])
-        assert abs(sum(lj) - -10.40) < 1e-7
-        assert abs(sum(el) - -212.96) < 1e-7
+        assert is_close(sum(lj), -10.40)
+        assert is_close(sum(el), -212.96)
 
     def test_gc2(self, qco1):
         els = dict(qco1.results["gc"].get_rows(columns=["Residue", "E_EL"]))
-        assert abs(els[35] - 24.27) < 1e-7
-        assert abs(els[227] - -13.65) < 1e-7
-        assert abs(els[286] - 22.21) < 1e-7
+        assert is_close(els[35], 24.27)
+        assert is_close(els[227], -13.65)
+        assert is_close(els[286], 22.21)
 
     def test_calcs(self, qco2):
-        assert abs(sum(qco2.results["1"].get_columns()[1]) - 14.189316) < 1e-7
-        assert abs(sum(qco2.results["2"].get_columns()[1]) - 16.44) < 1e-7
-        assert abs(sum(qco2.results["3"].get_columns()[1]) - 32.81) < 1e-7
-        assert abs(sum(qco2.results["4"].get_columns()[1]) - 179.94) < 1e-7
-        assert abs(sum(qco2.results["5"].get_columns()[1]) - 1183.01) < 1e-7
+        assert is_close(sum(qco2.results["1"].get_columns()[1]), 14.189316)
+        assert is_close(sum(qco2.results["2"].get_columns()[1]), 16.44)
+        assert is_close(sum(qco2.results["3"].get_columns()[1]), 32.81)
+        assert is_close(sum(qco2.results["4"].get_columns()[1]), 179.94)
+        assert is_close(sum(qco2.results["5"].get_columns()[1]), 1183.01)
 

tests/Qpyl/core/qdyninp_test.py

@@ -0,0 +1,61 @@
+#########################
+# py.test test functions
+#########################
+
+
+import pytest
+from Qpyl.core.qdyninp import QDynInput, QDynInputError
+
+class TestQDynInput: 
+
+    def test_parse_input(self):
+        # parse and write back, see if it's the same
+        qdis = open("data/qdyn.inp.1").read()
+        assert QDynInput(qdis).get_string() == qdis
+
+    def test_update_parms(self):
+        # update the params, see what happens
+        qdis = open("data/qdyn.inp.1").read()
+        qdi = QDynInput(qdis, parameters = {"md": {"temperature":100}})
+        assert qdi.parameters["md"]["temperature"] == 100
+        qdi.update(input_string=qdis)
+        assert int(qdi.parameters["md"]["temperature"]) == 300
+        qdi.update(parameters = {"md": {"temperature":100}})
+        assert qdi.parameters["md"]["temperature"] == 100
+
+    def test_parse_input_fail(self):
+        # fail on typos
+        qdis = open("data/qdyn.inp.1").read()
+        with pytest.raises(QDynInputError):
+            QDynInput(qdis.replace("temperature", "tmperature"))
+
+        with pytest.raises(QDynInputError):
+            QDynInput(qdis.replace("[lambdas]\n", ""))
+
+        with pytest.raises(QDynInputError):
+            QDynInput(qdis.replace("off", "whatisthis?"))
+
+        with pytest.raises(QDynInputError):
+            QDynInput(qdis.replace("300", "300_"))
+
+        # FEP and lambda must exist at the same time
+        with pytest.raises(QDynInputError):
+            qdi = QDynInput(qdis)
+            del qdi.parameters["lambdas"]
+            qdi.get_string()
+
+        # random_seed and initial_vel must exist at the same time
+        with pytest.raises(QDynInputError):
+            qdi = QDynInput(qdis)
+            del qdi.parameters["md"]["random_seed"]
+            qdi.get_string()
+
+        # energy file and interval must be both specified
+        with pytest.raises(QDynInputError):
+            qdi = QDynInput(qdis)
+            del qdi.parameters["intervals"]["energy"]
+            qdi.get_string()
+
+
+
+

tests/Qpyl/core/qfep_test.py

@@ -6,6 +6,9 @@
 
 from Qpyl.core.qfep import QFepOutput, QFepOutputError
 
+def is_close(a, b, rel_tol=1e-09, abs_tol=0.0):
+    return abs(a-b) <= max(rel_tol * max(abs(a), abs(b)), abs_tol)
+
 
 class TestQFepOutput:
 
@@ -46,10 +49,10 @@ def test_part0(self, qfo1):
 
     def test_part1(self, qfo1):
         sum_dg = sum(qfo1.part1.data.get_columns(["dG"])[0])
-        assert abs(sum_dg - 1651.77) < 1e-7
+        assert is_close(sum_dg, 1651.77)
 
     def test_dgfep1(self, qfo1):
-        assert abs(qfo1.part1.dg - -10.204) < 1e-7
+        assert is_close(qfo1.part1.dg, -10.204)
 
     def test_lra1(self, qfo1):
         # LRA of EQbond
@@ -60,7 +63,7 @@ def test_lra1(self, qfo1):
         lras = qfo1.part0.calc_lra(1.0, 0.0)
         print lras
         lra_eqbond = lras.get_columns(["LRA"])[0][1]
-        assert abs(lra_eqbond - -6.7350) < 1e-7
+        assert is_close(lra_eqbond, -6.7350)
 
     def test_lra2(self, qfo2):
         lras = qfo2.part0.calc_lra(1.0, 0.0)
@@ -71,48 +74,48 @@ def test_lra2(self, qfo2):
         assert lras == expected
 
     def test_dgs1(self, qfo2):
-        assert abs(qfo2.part3.dga - 23.98) < 1e-7
-        assert abs(qfo2.part3.dg0 - -21.98) < 1e-7
+        assert is_close(qfo2.part3.dga, 23.98)
+        assert is_close(qfo2.part3.dg0, -21.98)
 
     def test_dgs2(self, qfo2):
         sub_calc = qfo2.sub_calcs["ex_full_49"]
-        assert abs(sub_calc.part3.dga - 23.28) < 1e-7
-        assert abs(sub_calc.part3.dg0 - -23.45) < 1e-7
+        assert is_close(sub_calc.part3.dga, 23.28)
+        assert is_close(sub_calc.part3.dg0, -23.45)
 
     def test_dgs3(self, qfo2):
         sub_calc = qfo2.sub_calcs["ex_full_200"]
-        assert abs(sub_calc.part3.dga - 25.82) < 1e-7
-        assert abs(sub_calc.part3.dg0 - -15.95) < 1e-7
+        assert is_close(sub_calc.part3.dga, 25.82)
+        assert is_close(sub_calc.part3.dg0, -15.95)
 
     def test_dgs4(self, qfo2):
         sub_calc = qfo2.sub_calcs["ex_full_221"]
-        assert abs(sub_calc.part3.dga - 9.77) < 1e-7
-        assert abs(sub_calc.part3.dg0 - -52.08) < 1e-7
+        assert is_close(sub_calc.part3.dga, 9.77)
+        assert is_close(sub_calc.part3.dg0, -52.08)
 
     def test_dgs5(self, qfo2):
         sub_calc = qfo2.sub_calcs["ex_full_49_200_221"]
-        assert abs(sub_calc.part3.dga - 10.69) < 1e-7
-        assert abs(sub_calc.part3.dg0 - -47.67) < 1e-7
+        assert is_close(sub_calc.part3.dga, 10.69)
+        assert is_close(sub_calc.part3.dg0, -47.67)
 
     def test_dgs6(self, qfo2):
         sub_calc = qfo2.sub_calcs["ex_el_49_200_221"]
-        assert abs(sub_calc.part3.dga - 10.64) < 1e-7
-        assert abs(sub_calc.part3.dg0 - -47.71) < 1e-7
+        assert is_close(sub_calc.part3.dga, 10.64)
+        assert is_close(sub_calc.part3.dg0, -47.71)
 
     def test_dgs7(self, qfo2):
         sub_calc = qfo2.sub_calcs["ex_vdw_49_200_221"]
-        assert abs(sub_calc.part3.dga - 24.04) < 1e-7
-        assert abs(sub_calc.part3.dg0 - -21.94) < 1e-7
+        assert is_close(sub_calc.part3.dga, 24.04)
+        assert is_close(sub_calc.part3.dg0, -21.94)
 
     def test_dgs8(self, qfo2):
         sub_calc = qfo2.sub_calcs["QCP"]
-        assert abs(sub_calc.part3.dga - 23.31) < 1e-7
-        assert abs(sub_calc.part3.dg0 - -22.12) < 1e-7
+        assert is_close(sub_calc.part3.dga, 23.31)
+        assert is_close(sub_calc.part3.dg0, -22.12)
 
     def test_dgs9(self, qfo2):
         sub_calc = qfo2.sub_calcs["QCP_mass"]
-        assert abs(sub_calc.part3.dga - 23.60) < 1e-7
-        assert abs(sub_calc.part3.dg0 - -21.95) < 1e-7
+        assert is_close(sub_calc.part3.dga, 23.60)
+        assert is_close(sub_calc.part3.dg0, -21.95)
 
 
 

tests/Qpyl/core/qlibrary_test.py

@@ -3,12 +3,14 @@
 #########################
 
 
-import re
 import pytest
 
 from Qpyl.core.qlibrary import QLib, QLibError
 from Qpyl.core.qstructure import QStruct
 
+def is_close(a, b, rel_tol=1e-09, abs_tol=0.0):
+    return abs(a-b) <= max(rel_tol * max(abs(a), abs(b)), abs_tol)
+
 class TestQ:
     def test_read_write_lib(self):
         qlib = QLib("amber")
@@ -36,7 +38,7 @@ def test_atom(self):
         trp_CG = qlib.residue_dict["TRP"].atoms[7]
         assert trp_CG.name == "CG"
         assert trp_CG.atom_type == "Cstar"
-        assert abs(trp_CG.charge - -0.1415) < 1e-7
+        assert is_close(trp_CG.charge, -0.1415)
         assert trp_CG.residue.name == "TRP"
 
     def test_residue(self):
@@ -48,7 +50,7 @@ def test_residue(self):
         assert asp.connections == ["head N", "tail C"]
         hoh = qlib.residue_dict["HOH"]
         assert int(hoh.info["solvent"]) == 1
-        assert abs(float(hoh.info["density"]) - 0.0335) < 1e-7
+        assert is_close(float(hoh.info["density"]), 0.0335)
 
 
     def test_rounding(self):
@@ -58,15 +60,15 @@ def test_rounding(self):
         old_charge = prc.atoms[0].charge
         prc.atoms[0].charge -= 1e-7
         prc.rescale(prc.charge_groups[0], 1)
-        assert abs(prc.atoms[0].charge - old_charge) < 1e-7
+        assert is_close(prc.atoms[0].charge, old_charge)
 
     def test_rescale(self):
         qlib = QLib("oplsaa")
         qlib.read_lib("data/prc.lib")
         prc = qlib.residue_dict["PRC"]
         prc.atoms[0].charge -= 0.5
         prc.rescale(prc.charge_groups[0], 1)
-        assert abs(prc.atoms[0].charge - -0.7256169) < 1e-7
+        assert is_close(prc.atoms[0].charge, -0.72561699999)
 
     def test_rescale_badatom_fail(self):
         qlib = QLib("oplsaa")
@@ -208,7 +210,7 @@ def test_read_ffld(self):
         assert len(qlib.residue_dict["NMA"].atoms) == 6
         ash = qlib.residue_dict["ASH"]
         assert ash.atoms[1].atom_type == "CT1_C1_224"
-        assert abs(ash.atoms[1].charge - 0.14) < 1e-7
+        assert is_close(ash.atoms[1].charge, 0.14)
         assert "tail C" in ash.connections
         assert "head N" in ash.connections