Adding of the MMTK and Scientific tests

BuildServer/Unix/tests.sh

@@ -17,7 +17,21 @@ if [ $status -ne 0 ]; then
 	${MDANSE_SOURCE_DIR}/BuildServer/Unix/clean.sh
 	exit $status
 fi
-cd ../..
+cd ${MDANSE_SOURCE_DIR}
+
+#############################
+# DEPENDENCIES TESTS
+#############################
+echo -e "${BLUE}""Performing dependencies tests""${NORMAL}"
+cd Tests/DependeciesTests
+${PYTHONEXE} AllTests.py
+status=$?
+if [ $status -ne 0 ]; then
+	echo -e "${RED}" "One or several unit tests failed"
+	${MDANSE_SOURCE_DIR}/BuildServer/Unix/clean.sh
+	exit $status
+fi
+cd ${MDANSE_SOURCE_DIR}
 
 #############################
 # FUCNTIONAL TESTS
@@ -32,4 +46,5 @@ if [ $status -ne 0 ]; then
 	echo -e "${RED}" "One or several functional tests failed"
 	${MDANSE_SOURCE_DIR}/BuildServer/Unix/clean.sh
 	exit $status
-fi
\ No newline at end of file
+fi
+cd ${MDANSE_SOURCE_DIR}
\ No newline at end of file

Tests/DependenciesTests/AllTests.py

+# Run the complete test suite
+#
+# Written by Konrad Hinsen.
+#
+
+import unittest
+
+import mmtk_basic_tests
+import mmtk_energy_tests
+import mmtk_enm_tests
+import mmtk_internal_coordinate_tests
+import mmtk_normal_mode_tests
+import mmtk_pickle_tests
+import mmtk_restraint_tests
+import mmtk_subspace_tests
+import mmtk_trajectory_tests
+import mmtk_universe_tests
+
+
+import scientific_clustering_tests
+import scientific_intepolation_tests
+import scientific_signal_tests
+import scientific_vector_tests
+
+
+def suite():
+    test_suite = unittest.TestSuite()
+    test_suite.addTests(scientific_clustering_tests.suite())
+    test_suite.addTests(scientific_intepolation_tests.suite())
+    test_suite.addTests(scientific_signal_tests.suite())
+    test_suite.addTests(scientific_vector_tests.suite())
+    test_suite.addTests(mmtk_basic_tests.suite())
+    test_suite.addTests(mmtk_energy_tests.suite())
+    test_suite.addTests(mmtk_enm_tests.suite())
+    test_suite.addTests(mmtk_internal_coordinate_tests.suite())
+    test_suite.addTests(mmtk_normal_mode_tests.suite())
+    test_suite.addTests(mmtk_pickle_tests.suite())
+    test_suite.addTests(mmtk_restraint_tests.suite())
+    test_suite.addTests(mmtk_subspace_tests.suite())
+    test_suite.addTests(mmtk_trajectory_tests.suite())
+    test_suite.addTests(mmtk_universe_tests.suite())
+    return test_suite
+
+if __name__ == '__main__':
+    unittest.TextTestRunner(verbosity=2).run(suite())
+

Tests/DependenciesTests/mmtk_basic_tests.py

+# Basic tests
+#
+# Written by Konrad Hinsen
+#
+
+import numpy
+
+import unittest
+import MMTK
+from MMTK.Proteins import Protein
+
+class GroupOfAtomTest(unittest.TestCase):
+
+    """
+    Test the methods of Collections.GroupOfAtoms
+    """
+
+    def setUp(self):
+        self.molecule = MMTK.Molecule('water')
+        self.results = {}
+        self.results['numberOfAtoms'] = 3
+
+    def test_numberOfAtoms(self):
+        self.assertEqual(self.molecule.numberOfAtoms(),
+                         self.results['numberOfAtoms'])
+
+
+class SuperpositionTest(unittest.TestCase):
+
+    """
+    Test the findTransformation method
+    """
+
+    def test_rotation_translation(self):
+        for m in [MMTK.Molecule('water'), Protein('bala1')]:
+            universe = MMTK.InfiniteUniverse()
+            universe.addObject(m)
+            ref_conf = universe.copyConfiguration()
+            universe.translateBy(MMTK.Vector(0.1, -1.3, 1.2))
+            universe.rotateAroundOrigin(MMTK.Vector(0., 1., 1.), 0.7)
+            tr, rms = universe.findTransformation(ref_conf)
+            self.assert_(abs(rms) < 1.e-5)
+            universe.applyTransformation(tr)
+            self.assert_(universe.rmsDifference(ref_conf) < 1.e-5)
+            axis, angle = tr.rotation().axisAndAngle()
+            self.assert_(abs(angle - 0.7) < 1.e-5)
+            self.assert_(abs(abs(numpy.cos(axis.angle(MMTK.Vector(0., 1., 1.))))-1)
+                         < 1.e-5)
+
+def suite():
+    loader = unittest.TestLoader()
+    s = unittest.TestSuite()
+    s.addTest(loader.loadTestsFromTestCase(GroupOfAtomTest))
+    s.addTest(loader.loadTestsFromTestCase(SuperpositionTest))
+    return s
+
+if __name__ == '__main__':
+    unittest.main()

Tests/DependenciesTests/mmtk_energy_tests.py

+# Energy tests
+#
+# Written by Konrad Hinsen
+#
+
+import unittest
+from mmtk_subsets import SubsetTest
+import numpy
+
+from MMTK import *
+from MMTK.MoleculeFactory import MoleculeFactory
+from MMTK.ForceFields import Amber99ForceField, LennardJonesForceField
+from MMTK_forcefield import NonbondedList
+from MMTK.Random import randomPointInBox
+from MMTK.Geometry import SCLattice
+from MMTK.Utility import pairs
+from Scientific.Geometry import ex, ey, ez
+from cStringIO import StringIO
+import itertools
+
+factory = MoleculeFactory()
+factory.createGroup('dihedral_test')
+
+factory.addAtom('dihedral_test', 'C1', 'C')
+factory.addAtom('dihedral_test', 'C2', 'C')
+factory.addAtom('dihedral_test', 'C3', 'C')
+factory.addAtom('dihedral_test', 'C4', 'C')
+
+factory.addBond('dihedral_test', 'C1', 'C2')
+factory.addBond('dihedral_test', 'C2', 'C3')
+factory.addBond('dihedral_test', 'C3', 'C4')
+
+factory.setAttribute('dihedral_test', 'C1.amber_atom_type', 'D1')
+factory.setAttribute('dihedral_test', 'C2.amber_atom_type', 'D2')
+factory.setAttribute('dihedral_test', 'C3.amber_atom_type', 'D3')
+factory.setAttribute('dihedral_test', 'C4.amber_atom_type', 'D4')
+factory.setAttribute('dihedral_test', 'C1.amber_charge', 0.)
+factory.setAttribute('dihedral_test', 'C2.amber_charge', 0.)
+factory.setAttribute('dihedral_test', 'C3.amber_charge', 0.)
+factory.setAttribute('dihedral_test', 'C4.amber_charge', 0.)
+
+factory.setPosition('dihedral_test', 'C1', Vector(1., 0., 0.))
+factory.setPosition('dihedral_test', 'C2', Vector(0., 0., 0.))
+factory.setPosition('dihedral_test', 'C3', Vector(0., 0., 1.))
+factory.setPosition('dihedral_test', 'C4', Vector(1., 0., 1.))
+
+
+def sorted_tuple(pair):
+    if pair[0] < pair[1]:
+        return (pair[0], pair[1])
+    else:
+        return (pair[1], pair[0])
+
+
+class DihedralTest(unittest.TestCase):
+
+    def setUp(self):
+        self.universe = InfiniteUniverse()
+        self.universe.addObject(factory.retrieveMolecule('dihedral_test'))
+        self.mod_template = """Amber parameters
+
+MASS
+D1 12.0
+D2 12.0
+D3 12.0
+D4 12.0
+
+BOND
+D1-D2    0.0    1.000
+D2-D3    0.0    1.000
+D3-D4    0.0    1.000
+
+ANGL
+D1-D2-D3      0.0      90.0
+D2-D3-D4      0.0      90.0
+
+DIHEDRAL
+D1-D2-D3-D4   1%15.6f%15.6f%15.6f
+
+NONB
+  D1          1.0000    0.0000   0.00000
+  D2          1.0000    0.0000   0.00000
+  D3          1.0000    0.0000   0.00000
+  D4          1.0000    0.0000   0.00000
+"""
+
+    def _gradientTest(self, delta = 0.0001):
+        e0, grad = self.universe.energyAndGradients()
+        atoms = self.universe.atomList()
+        for a in atoms:
+            num_grad = []
+            for v in [ex, ey, ez]:
+                x = a.position()
+                a.setPosition(x+delta*v)
+                eplus = self.universe.energy()
+                a.setPosition(x-delta*v)
+                eminus = self.universe.energy()
+                a.setPosition(x)
+                num_grad.append(0.5*(eplus-eminus)/delta)
+            self.assert_((Vector(num_grad)-grad[a]).length() < 1.e-2)
+
+    def _forceConstantTest(self, delta = 0.0001):
+        e0, grad0, fc = self.universe.energyGradientsAndForceConstants()
+        atoms = self.universe.atomList()
+        for a1, a2 in itertools.chain(itertools.izip(atoms, atoms),
+                                      pairs(atoms)):
+            num_fc = []
+            for v in [ex, ey, ez]:
+                x = a1.position()
+                a1.setPosition(x+delta*v)
+                e_plus, grad_plus = self.universe.energyAndGradients()
+                a1.setPosition(x-delta*v)
+                e_minus, grad_minus = self.universe.energyAndGradients()
+                a1.setPosition(x)
+                num_fc.append(0.5*(grad_plus[a2]-grad_minus[a2])/delta)
+            num_fc = numpy.array([a.array for a in num_fc])
+            diff = numpy.fabs(numpy.ravel(num_fc-fc[a1, a2].array))
+            error = numpy.maximum.reduce(diff)
+            self.assert_(error < 5.e-2)
+            
+    def _dihedralTerm(self, n, phase, V):
+
+        mod_file = self.mod_template % \
+                   (V/(Units.kcal/Units.mol), phase/Units.deg, n)
+        ff = Amber99ForceField(mod_files=[StringIO(mod_file)])
+        self.universe.setForceField(ff)
+
+        param = self.universe.energyEvaluatorParameters()
+        i1, i2, i3, i4, n_test, phase_test, V_test = \
+            param['cosine_dihedral_term'][0]
+        self.assertEqual(n_test, n)
+        # The accuracy is no better than five digits because the
+        # parameters pass through a text representation.
+        self.assertAlmostEqual(phase_test, phase, 5)
+        self.assertAlmostEqual(V_test, V, 5)
+
+        two_pi = 2.*numpy.pi
+        m = self.universe[0]
+        for angle in numpy.arange(0., two_pi, 0.1):
+            m.C4.setPosition(Vector(numpy.cos(angle), numpy.sin(angle), 1.))
+            e = self.universe.energyTerms()['cosine dihedral angle']
+            da = self.universe.dihedral(m.C1, m.C2, m.C3, m.C4)
+            e_ref = V*(1.+numpy.cos(n*angle-phase))
+            self.assertAlmostEqual(angle % two_pi, da % two_pi, 14)
+            self.assertAlmostEqual(e, e_ref, 5)
+            self._gradientTest()
+            self._forceConstantTest()
+
+    def test_dihedral(self):
+        for n in [1, 2, 3, 4]:
+            for phase in numpy.arange(0., 6., 0.5):
+                for V in [-10., 2.]:
+                    self._dihedralTerm(n, phase, V)
+
+
+
+class NonbondedListTest:
+
+    def test_nonbondedList(self):
+
+        self.universe.configuration()
+        atoms = self.universe.atomList()
+        atom_indices = numpy.array([a.index for a in self.universe.atomList()])
+        empty = numpy.zeros((0, 2), numpy.int)
+
+        for cutoff in [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.]:
+
+            nblist = NonbondedList(empty, empty, atom_indices,
+                                   self.universe._spec, cutoff)
+            nblist.update(self.universe.configuration().array)
+            distances = nblist.pairDistances()
+            pairs1 = nblist.pairIndices()
+            pairs1 = [sorted_tuple(pairs1[i]) for i in range(len(pairs1))
+                      if distances[i] < cutoff]
+            pairs1.sort(lambda a, b: cmp(a[0], b[0]) or cmp(a[1], b[1]))
+
+            pairs2 = []
+            for i in range(len(atoms)):
+                for j in range(i+1, len(atoms)):
+                    d = self.universe.distance(atoms[i], atoms[j])
+                    if d < cutoff:
+                        pairs2.append(sorted_tuple((atoms[i].index,
+                                                    atoms[j].index)))
+            pairs2.sort(lambda a, b: cmp(a[0], b[0]) or cmp(a[1], b[1]))
+
+            self.assertEqual(pairs1, pairs2)
+
+
+class InfiniteUniverseNonbondedListTest(unittest.TestCase,
+                                        NonbondedListTest):
+    def setUp(self):
+        self.universe = InfiniteUniverse()
+        for i in range(100):
+            p = randomPointInBox(1.)
+            self.universe.addObject(Atom('C', position = p))
+
+class OrthorhombicUniverseNonbondedListTest(unittest.TestCase,
+                                            NonbondedListTest):
+    def setUp(self):
+        self.universe = OrthorhombicPeriodicUniverse((1.5, 0.8, 1.1))
+        for i in range(100):
+            p = self.universe.boxToRealCoordinates(randomPointInBox(1.))
+            self.universe.addObject(Atom('C', position = p))
+
+class ParallelepipedicUniverseNonbondedListTest(unittest.TestCase,
+                                                NonbondedListTest):
+    def setUp(self):
+        a = Vector(1.5, 0.3, 0.)
+        b = Vector(-0.1, 0.8, 0.2)
+        c = Vector(0., -0.2, 1.1)
+        self.universe = ParallelepipedicPeriodicUniverse((a, b, c))
+        for i in range(100):
+            p = self.universe.boxToRealCoordinates(randomPointInBox(1.))
+            self.universe.addObject(Atom('C', position = p))
+
+class LennardJonesSubsetTest(unittest.TestCase,
+                             SubsetTest):
+
+    def setUp(self):
+        self.universe = OrthorhombicPeriodicUniverse((2., 2., 2.),
+                                                     LennardJonesForceField())
+        for point in SCLattice(0.5, 4):
+            self.universe.addObject(Atom('Ar', position=point))
+        self.subset1 = Collection(self.universe.atomList()[:10])
+        self.subset2 = Collection(self.universe.atomList()[20:30])
+
+            
+def suite():
+    loader = unittest.TestLoader()
+    s = unittest.TestSuite()
+    s.addTest(loader.loadTestsFromTestCase(DihedralTest))
+    s.addTest(loader.loadTestsFromTestCase(InfiniteUniverseNonbondedListTest))
+    s.addTest(loader.loadTestsFromTestCase(OrthorhombicUniverseNonbondedListTest))
+    s.addTest(loader.loadTestsFromTestCase(ParallelepipedicUniverseNonbondedListTest))
+    s.addTest(loader.loadTestsFromTestCase(LennardJonesSubsetTest))
+    return s
+
+
+if __name__ == '__main__':
+    unittest.main()

Tests/DependenciesTests/mmtk_enm_tests.py

+# Elastic network model tests
+#
+# Written by Konrad Hinsen
+#
+
+import unittest
+from MMTK import *
+from MMTK.ForceFields import CalphaForceField, DeformationForceField, \
+                             AnisotropicNetworkForceField
+from MMTK.Proteins import Protein
+from mmtk_subsets import SubsetTest
+
+class CalphaFFSubsetTest(unittest.TestCase,
+                         SubsetTest):
+
+    def setUp(self):
+        self.universe = InfiniteUniverse(CalphaForceField())
+        protein = Protein('insulin_calpha')
+        self.universe.addObject(protein)        
+        self.subset1 = protein[0]
+        self.subset2 = protein[1]
+
+
+class DeformationFFSubsetTest(unittest.TestCase,
+                              SubsetTest):
+
+    def setUp(self):
+        self.universe = InfiniteUniverse(DeformationForceField(cutoff=5.))
+        protein = Protein('insulin_calpha')
+        self.universe.addObject(protein)        
+        self.subset1 = protein[0]
+        self.subset2 = protein[1]
+
+
+class ANMFFSubsetTest(unittest.TestCase,
+                      SubsetTest):
+
+    def setUp(self):
+        self.universe = InfiniteUniverse(AnisotropicNetworkForceField(cutoff=5.))
+        protein = Protein('insulin_calpha')
+        self.universe.addObject(protein)        
+        self.subset1 = protein[0]
+        self.subset2 = protein[1]
+
+
+def suite():
+    loader = unittest.TestLoader()
+    s = unittest.TestSuite()
+    s.addTest(loader.loadTestsFromTestCase(CalphaFFSubsetTest))
+    s.addTest(loader.loadTestsFromTestCase(DeformationFFSubsetTest))
+    s.addTest(loader.loadTestsFromTestCase(ANMFFSubsetTest))
+    return s
+
+
+if __name__ == '__main__':
+    unittest.main()

Tests/DependenciesTests/mmtk_internal_coordinate_tests.py

+# Internal coordinate tests
+#
+# Written by Konrad Hinsen
+#
+
+import unittest
+
+import numpy
+from MMTK.Proteins import Protein
+from MMTK.InternalCoordinates import BondLength, BondAngle, DihedralAngle
+
+class ICTest(unittest.TestCase):
+
+    def setUp(self):
+        # Pick a proline residue because of its cycle
+        self.protein = Protein("insulin.pdb")
+        self.chain = self.protein[1]
+        self.res = self.chain[27]
+
+    def test_phipsi(self):
+        phi, psi = self.res.phiPsi()
+        phi_o = self.res.phiAngle()
+        psi_o = self.res.psiAngle()
+        self.assertAlmostEqual(phi, phi_o.getValue(), 10)
+        self.assertAlmostEqual(psi, psi_o.getValue(), 10)
+
+    def test_changephi(self):
+        phi = self.res.phiAngle()
+        current_phi = phi.getValue()
+        for dphi in [-0.2, -0.1, 0.1, 0.2]:
+            phi.setValue(current_phi+dphi)
+            self.assertAlmostEqual(phi.getValue()  % (2.*numpy.pi),
+                                   (current_phi+dphi) % (2.*numpy.pi), 10)
+    def test_changebond(self):
+        b = BondLength(self.res.peptide.C_alpha, self.res.peptide.C)
+        current_b = b.getValue()
+        for db in [-0.02, -0.01, 0.01, 0.02]:
+            b.setValue(current_b + db)
+            self.assertAlmostEqual(b.getValue(), current_b + db)
+
+    def test_changeangle(self):
+        a = BondAngle(self.res.peptide.N, self.res.peptide.C_alpha,
+                      self.res.peptide.C)
+        current_a = a.getValue()
+        for da in [-0.2, -0.1, 0.1, 0.2]:
+            a.setValue(current_a+da)
+            self.assertAlmostEqual(a.getValue()  % (2.*numpy.pi),
+                                   (current_a+da) % (2.*numpy.pi), 10)
+
+    def test_cyclic(self):
+        self.assertRaises(ValueError, BondLength,
+                          self.res.peptide.N, self.res.peptide.C_alpha)
+        self.assertRaises(ValueError, BondAngle,
+                          self.res.peptide.N, self.res.peptide.C_alpha,
+                          self.res.sidechain.C_beta)
+        self.assertRaises(ValueError, DihedralAngle,
+                          self.res.peptide.N, self.res.peptide.C_alpha,
+                          self.res.sidechain.C_beta, self.res.sidechain.C_gamma)
+            
+def suite():
+    return unittest.TestLoader().loadTestsFromTestCase(ICTest)
+
+if __name__ == '__main__':
+    unittest.main()

Tests/DependenciesTests/mmtk_normal_mode_tests.py

+# Normal mode tests
+#
+# Written by Konrad Hinsen
+#
+
+import unittest
+import MMTK
+from MMTK.Proteins import Protein
+from MMTK.ForceFields import Amber99ForceField
+from MMTK.Subspace import RigidMotionSubspace, PairDistanceSubspace
+from MMTK import NormalModes
+
+class WaterTest(unittest.TestCase):
+
+    """
+    Test NormalModes.EnergeticModes
+    and NormalModes.VibrationalModes
+    """
+
+    def setUp(self):
+        self.universe = MMTK.InfiniteUniverse(Amber99ForceField())
+        self.universe.water = MMTK.Molecule('water')
+
+    def test_energeticModes(self):
+        emodes = NormalModes.EnergeticModes(self.universe)
+        fc = emodes.force_constants[6:]
+        self.assertAlmostEqual(fc[0], 757849.3957485439, 5)
+        self.assertAlmostEqual(fc[1], 1041551.2240706938, 5)
+        self.assertAlmostEqual(fc[2], 1388251.2, 5)
+        for i in range(len(emodes)):
+            mi = emodes.rawMode(i)
+            norm_sq = mi.dotProduct(mi)
+            self.assertAlmostEqual(norm_sq, 1.)
+            for j in range(i+1, len(emodes)):
+                overlap = mi.dotProduct(emodes.rawMode(j))
+                self.assert_(overlap < 1.e-15)
+        self.assertAlmostEqual(emodes[6].norm(), 0.0025656740328985168)
+        self.assertAlmostEqual(emodes[7].norm(), 0.0021885627126988836)
+        self.assertAlmostEqual(emodes[8].norm(), 0.0018956532696964624)
+        f = emodes.fluctuations()
+        self.assertAlmostEqual(f[self.universe.water.O],  3.359215370401356e-06)
+        self.assertAlmostEqual(f[self.universe.water.H1], 2.061890142153454e-06)
+        self.assertAlmostEqual(f[self.universe.water.H2], 2.061890142153454e-06)
+        af = emodes.anisotropicFluctuations()
+        self.assertAlmostEqual(f[self.universe.water.O],
+                               af[self.universe.water.O].trace())
+        self.assertAlmostEqual(f[self.universe.water.H1],
+                               af[self.universe.water.H1].trace())
+        self.assertAlmostEqual(f[self.universe.water.H2],
+                               af[self.universe.water.H2].trace())
+
+    def test_vibrationalModes(self):
+        vmodes = NormalModes.VibrationalModes(self.universe)
+        freq = vmodes.frequencies[6:]
+        self.assertAlmostEqual(freq[0], 87.220841117866954)
+        self.assertAlmostEqual(freq[1], 112.01238171677888)
+        self.assertAlmostEqual(freq[2], 181.05242207954848)
+        for i in range(len(vmodes)):
+            mi = vmodes.rawMode(i)
+            norm_sq = mi.dotProduct(mi)
+            self.assertAlmostEqual(norm_sq, 1.)
+            for j in range(i+1, len(vmodes)):
+                overlap = mi.dotProduct(vmodes.rawMode(j))
+                self.assert_(overlap < 1.e-15)
+        self.assertAlmostEqual(vmodes[6].norm(), 0.0038454773367577063)
+        self.assertAlmostEqual(vmodes[7].norm(), 0.0030509249071616175)
+        self.assertAlmostEqual(vmodes[8].norm(), 0.001953454891033823)
+        f = vmodes.fluctuations()
+        self.assertAlmostEqual(f[self.universe.water.O], 8.0141737611250569e-08)
+        self.assertAlmostEqual(f[self.universe.water.H1], 6.9381391949153065e-06)
+        self.assertAlmostEqual(f[self.universe.water.H2], 6.9381391949153023e-06)
+        af = vmodes.anisotropicFluctuations()
+        self.assertAlmostEqual(f[self.universe.water.O],
+                               af[self.universe.water.O].trace())
+        self.assertAlmostEqual(f[self.universe.water.H1],
+                               af[self.universe.water.H1].trace())