Skip to content

GitLab

Switch branch/tag
Find file Normal viewHistoryPermalink
TestUnitCell.py 2.77 KB
Newer Older
Jason Swails's avatar
Fix bug in unitcell.py (computeLengthsAndAngles), which resulted in a NameError.  
Jason Swails committed
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 
import unittest
import math
from simtk.unit import *
from simtk.openmm.app.internal.unitcell import computePeriodicBoxVectors
from simtk.openmm.app.internal.unitcell import computeLengthsAndAngles

def strip_units(x):
    if is_quantity(x): return x.value_in_unit_system(md_unit_system)
    return x

class TestAmberPrmtopFile(unittest.TestCase):

    """Test the AmberPrmtopFile.createSystem() method."""

    def testComputePBCVectors(self):
        """ Tests computing periodic box vectors """
        deg90 = 90 * degrees
        vecs = computePeriodicBoxVectors(1, 2, 3, deg90, deg90, deg90)
        a, b, c = vecs
        self.assertAlmostEqual(a[0]/nanometers, 1)
        self.assertAlmostEqual(a[1]/nanometers, 0)
        self.assertAlmostEqual(a[2]/nanometers, 0)
        self.assertAlmostEqual(b[0]/nanometers, 0)
        self.assertAlmostEqual(b[1]/nanometers, 2)
        self.assertAlmostEqual(b[2]/nanometers, 0)
        self.assertAlmostEqual(c[0]/nanometers, 0)
        self.assertAlmostEqual(c[1]/nanometers, 0)
        self.assertAlmostEqual(c[2]/nanometers, 3)

        # Make sure round-trip works
        la, lb, lc, al, be, ga = computeLengthsAndAngles(vecs)
        self.assertAlmostEqual(la, 1)
        self.assertAlmostEqual(lb, 2)
        self.assertAlmostEqual(lc, 3)
        self.assertAlmostEqual(al, math.pi / 2)
        self.assertAlmostEqual(be, math.pi / 2)
        self.assertAlmostEqual(ga, math.pi / 2)

        # Now test a truncated octahedron. Can't do a simple round-trip though,
        # due to the reduced form. So test the *second* round-trip, which should
        # yield the same measurements
        vecs = computePeriodicBoxVectors(4.24388485, 4.24388485, 4.24388485,
                109.4712190*degrees, 109.4712190*degrees, 109.4712190*degrees)
        la, lb, lc, al, be, ga = computeLengthsAndAngles(vecs)
        vecs2 = computePeriodicBoxVectors(la, lb, lc, al, be, ga)
        la2, lb2, lc2, al2, be2, ga2 = computeLengthsAndAngles(vecs2)

        # Now make sure that the round-trip worked
        self.assertAlmostEqual(strip_units(la), strip_units(la2))
        self.assertAlmostEqual(strip_units(lb), strip_units(lb2))
        self.assertAlmostEqual(strip_units(lc), strip_units(lc2))
        self.assertAlmostEqual(strip_units(al), strip_units(al2))
        self.assertAlmostEqual(strip_units(be), strip_units(be2))
        self.assertAlmostEqual(strip_units(ga), strip_units(ga2))

        # Check that the vectors are the same
        a1, a2, a3 = vecs
        b1, b2, b3 = vecs2
        for x, y in zip(a1, b1):
            self.assertAlmostEqual(strip_units(x), strip_units(y))
        for x, y in zip(a2, b2):
            self.assertAlmostEqual(strip_units(x), strip_units(y))
        for x, y in zip(a3, b3):
            self.assertAlmostEqual(strip_units(x), strip_units(y))