def testMol3SingleWrite(self):
     """ Tests writing mol3 file of single ResidueTemplate """
     mol2 = formats.Mol2File.parse(get_fn('tripos9.mol2'))
     formats.Mol2File.write(mol2, get_fn('tripos9.mol3', written=True),
                            mol3=True)
     self.assertTrue(diff_files(get_fn('tripos9.mol3', written=True),
                                get_saved_fn('tripos9.mol3')))

 def test_box_from_system(self):
     """ Tests loading box from System """
     parm = load_file(get_fn('solv2.parm7'), get_fn('solv2.rst7'))
     system = parm.createSystem(nonbondedMethod=app.PME,
                                nonbondedCutoff=8*u.angstroms)
     top = openmm.load_topology(parm.topology, system)
     np.testing.assert_allclose(parm.box, top.box)

 def setUp(self):
     self.lfh=sys.stderr
     self.verbose=False
     self.pathPdbxDataFile = get_fn("1kip.cif")
     self.pathOutputFile = get_fn("testOutputDataFile.cif", written=True)
     if not os.path.exists(get_fn('writes')):
         os.makedirs(get_fn('writes'))

 def testAdd(self):
     """ Tests combining AmberParm instances """
     parm1 = readparm.AmberParm(get_fn('phenol.prmtop'))
     parm2 = readparm.AmberParm(get_fn('biphenyl.prmtop'))
     comb = parm1 + parm2
     self.assertEqual(len(comb.atoms), len(parm1.atoms) + len(parm2.atoms))
     for a1, a2 in zip(comb.atoms, parm1.atoms + parm2.atoms):
         self.assertEqual(a1.name, a2.name)
         self.assertEqual(a1.mass, a2.mass)
         self.assertEqual(a1.charge, a2.charge)
         self.assertEqual(a1.radii, a2.radii)
     self.assertEqual(len(comb.residues), len(parm1.residues) + len(parm2.residues))
     for r1, r2 in zip(comb.residues, parm1.residues + parm2.residues):
         self.assertEqual(len(r1), len(r2))
         self.assertEqual(r1.name, r2.name)
         self.assertEqual(r1.chain, r2.chain)
     # In-place now
     parm1 += parm2
     self.assertEqual(len(parm1.atoms), len(comb.atoms))
     for a1, a2 in zip(comb.atoms, parm1.atoms):
         self.assertEqual(a1.name, a2.name)
         self.assertEqual(a1.mass, a2.mass)
         self.assertEqual(a1.charge, a2.charge)
         self.assertEqual(a1.radii, a2.radii)
     self.assertEqual(len(parm1.residues), len(comb.residues))
     for r1, r2 in zip(comb.residues, parm1.residues):
         self.assertEqual(len(r1), len(r2))
         self.assertEqual(r1.name, r2.name)
         self.assertEqual(r1.chain, r2.chain)

    def test_geometric_combining_rule_energy(self):
        """ Tests converting geom. comb. rule energy from Gromacs to Amber """
        top = load_file(os.path.join(get_fn('05.OPLS'), 'topol.top'),
                        xyz=os.path.join(get_fn('05.OPLS'), 'conf.gro'))
        self.assertEqual(top.combining_rule, 'geometric')
        del top.rb_torsions[:]
        parm = load_file(get_saved_fn('opls.parm7'),
                         xyz=os.path.join(get_fn('05.OPLS'), 'conf.gro'))
        self.assertEqual(parm.combining_rule, 'geometric')
        self.assertFalse(parm.has_NBFIX())

        sysg = top.createSystem()
        sysa = parm.createSystem()

        cong = mm.Context(sysg, mm.VerletIntegrator(0.001), CPU)
        cona = mm.Context(sysa, mm.VerletIntegrator(0.001), CPU)

        cong.setPositions(top.positions)
        cona.setPositions(top.positions)

        self._check_energies(top, cong, parm, cona)

        # Make an NBFIX
        self.assertFalse(parm.has_NBFIX())
        parm.parm_data['LENNARD_JONES_ACOEF'][-4] = 10.0
        self.assertTrue(parm.has_NBFIX())
        parm.createSystem()

 def testAmberMdcrdNumpy(self):
     """ Test writing ASCII trajectory file passing numpy arrays """
     import numpy as np
     box = np.asarray([15, 15, 15])
     Mdcrd = asciicrd.AmberMdcrd
     crd = Mdcrd(get_fn('testc.mdcrd', written=True), natom=15, hasbox=False,
                 mode='w', title='Test file')
     coorddata = np.arange(45).reshape((15,3))
     crd.add_coordinates(coorddata)
     crd.add_coordinates(coorddata+1)
     crd.add_coordinates(coorddata+2)
     crd.add_coordinates(coorddata+3)
     crd.add_coordinates(coorddata+4)
     crd.close()
     crd = Mdcrd(get_fn('testcb.mdcrd', written=True), natom=18, hasbox=True,
                 mode='w', title='Test file')
     coorddata = np.arange(54).reshape((18,3))
     crd.add_coordinates(coorddata)
     crd.add_box(box)
     crd.add_coordinates(coorddata+1)
     crd.add_box(box)
     crd.add_coordinates(coorddata+2)
     crd.add_box(box)
     crd.add_coordinates(coorddata+3)
     crd.add_box(box)
     crd.add_coordinates(coorddata+4)
     crd.add_box(box)
     crd.close()
     self._check_written_mdcrds(box)

 def _check_written_restarts(self, box):
     # Now try to read them and verify the information (keep in mind that the
     # restart velocities are scaled down then back up, so you'll need to use
     # assertAlmostEqual in this case).
     rst = readparm.Rst7.open(get_fn('testc.rst7', written=True))
     self.assertFalse(rst.hasbox)
     self.assertFalse(rst.hasvels)
     for x1, x2 in zip(rst.coordinates, range(27)):
         self.assertEqual(x1, x2)
     rst = readparm.Rst7.open(get_fn('testcb.rst7', written=True))
     self.assertTrue(rst.hasbox)
     self.assertFalse(rst.hasvels)
     for x1, x2 in zip(rst.coordinates, range(21)):
         self.assertEqual(x1, x2)
     for x1, x2 in zip(rst.box, box):
         self.assertEqual(x1, x2)
     rst = readparm.Rst7.open(get_fn('testcv.rst7', written=True))
     self.assertTrue(rst.hasvels)
     self.assertFalse(rst.hasbox)
     for x1, x2 in zip(rst.coordinates, range(60)):
         self.assertEqual(x1, x2)
     for x1, x2 in zip(rst.vels, reversed(range(60))):
         self.assertAlmostEqual(x1, x2, places=5)
     rst = readparm.Rst7.open(get_fn('testcvb.rst7', written=True))
     self.assertTrue(rst.hasvels)
     self.assertTrue(rst.hasbox)
     for x1, x2 in zip(rst.coordinates, range(45)):
         self.assertEqual(x1, x2)
     for x1, x2 in zip(rst.vels, reversed(range(45))):
         self.assertAlmostEqual(x1, x2, places=5)
     for x1, x2 in zip(rst.box, box):
         self.assertEqual(x1, x2)

    def test_dppc(self):
        """ Tests non-standard Gromacs force fields and nonbonded exceptions """
        # We know what we're doing
        warnings.filterwarnings('ignore', category=GromacsWarning)
        top = load_file(os.path.join(get_fn('12.DPPC'), 'topol.top'),
                        xyz=os.path.join(get_fn('12.DPPC'), 'conf.gro'))
        self.assertEqual(top.combining_rule, 'lorentz')

        # Create the system and context, then calculate the energy decomposition
        system = top.createSystem()
        context = mm.Context(system, mm.VerletIntegrator(0.001), CPU)
        context.setPositions(top.positions)
        energies = energy_decomposition(top, context, nrg=u.kilojoules_per_mole)

        # Compare with Lee-Ping's answers.
        self.assertAlmostEqual(energies['bond'], 0)
        self.assertAlmostEqual(energies['angle'], 1405.7354199, places=4)
        self.assertAlmostEqual(energies['dihedral'], 236.932663255, places=4)
        self.assertAlmostEqual(energies['improper'], 33.201541811, places=4)
        self.assertAlmostEqual(energies['rb_torsion'], 428.0550599, places=4)
        self.assertRelativeEqual(energies['nonbonded'], -16432.8092955, places=4)
        gmxfrc = get_forces_from_xvg(os.path.join(get_fn('12.DPPC'), 'force.xvg'))
        ommfrc = context.getState(getForces=True).getForces().value_in_unit(
                    u.kilojoules_per_mole/u.nanometer)
        zero_ep_frc(ommfrc, top)
        max_diff = get_max_diff(gmxfrc, ommfrc)
        self.assertLess(max_diff, 5)

    def test_round_trip(self):
        """ Test ParmEd -> OpenMM round trip with Gromacs system """
        # Use DPPC to get RB-torsions tested. Also check that it initially fails
        # with the CustomNonbondedForce
        warnings.filterwarnings('ignore', category=GromacsWarning)
        top = load_file(os.path.join(get_fn('12.DPPC'), 'topol.top'),
                        xyz=os.path.join(get_fn('12.DPPC'), 'conf.gro'))
        self.assertEqual(top.combining_rule, 'lorentz')

        system = top.createSystem()
        def bad_system():
            return openmm.load_topology(top.topology, system).createSystem()
        warnings.filterwarnings('ignore', category=OpenMMWarning)
        self.assertTrue(
                openmm.load_topology(top.topology, system).unknown_functional
        )
        self.assertRaises(ParmedError, bad_system)
        for i in range(len(system.getForces())):
            if isinstance(system.getForce(i), mm.CustomNonbondedForce):
                system.removeForce(i)
                break
        system2 = openmm.load_topology(top.topology, system).createSystem()
        con1 = mm.Context(system, mm.VerletIntegrator(0.001), CPU)
        con2 = mm.Context(system2, mm.VerletIntegrator(0.001), CPU)
        con1.setPositions(top.positions)
        con2.setPositions(top.positions)
        ene1 = energy_decomposition(top, con1)
        ene2 = energy_decomposition(top, con2)
        self.assertAlmostEqual(ene1['bond'], ene2['bond'])
        self.assertAlmostEqual(ene1['angle'], ene2['angle'])
        self.assertAlmostEqual(ene1['dihedral'], ene2['dihedral'])
        self.assertAlmostEqual(ene1['improper'], ene2['improper'])
        self.assertAlmostEqual(ene1['rb_torsion'], ene2['rb_torsion'])
        self.assertAlmostEqual(ene1['nonbonded'], ene2['nonbonded'])

    def test_jac(self):
        """ Tests the JAC benchmark Gromacs system nrg and force (no PBC) """
        # Load the top and gro files
        top = load_file(os.path.join(get_fn('07.DHFR-Liquid-NoPBC'), 'topol.top'),
                        xyz=os.path.join(get_fn('07.DHFR-Liquid-NoPBC'), 'conf.gro'))
        self.assertEqual(top.combining_rule, 'lorentz')

        # Create the system and context, then calculate the energy decomposition
        system = top.createSystem()
        context = mm.Context(system, mm.VerletIntegrator(0.001), CPU)
        context.setPositions(top.positions)
        energies = energy_decomposition(top, context, nrg=u.kilojoules_per_mole)

        # Compare with Lee-Ping's answers. Make sure we zero-out forces for
        # virtual sites, since OMM doesn't do this and Gromacs does.
        self.assertAlmostEqual(energies['bond'], 35.142565, places=3)
        self.assertAlmostEqual(energies['angle'], 3735.514669, places=3)
        self.assertAlmostEqual(energies['dihedral'], 7277.741635, delta=0.002)
        self.assertRelativeEqual(energies['nonbonded'], -288718.981405, places=4)
        gmxfrc = get_forces_from_xvg(
                os.path.join(get_fn('07.DHFR-Liquid-NoPBC'), 'force.xvg'))
        ommfrc = context.getState(getForces=True).getForces().value_in_unit(
                    u.kilojoules_per_mole/u.nanometer)
        zero_ep_frc(ommfrc, top)
        max_diff = get_max_diff(gmxfrc, ommfrc)
        self.assertLess(max_diff, 0.5)

    def test_jac_pme(self):
        """ Tests the JAC benchmark Gromacs system nrg and force (PME) """
        # Load the top and gro files
        top = load_file(os.path.join(get_fn('09.DHFR-PME'), 'topol.top'),
                        xyz=os.path.join(get_fn('09.DHFR-PME'), 'conf.gro'))
        self.assertEqual(top.combining_rule, 'lorentz')

        # Create the system and context, then calculate the energy decomposition
        system = top.createSystem(nonbondedMethod=app.PME,
                                  constraints=app.HBonds,
                                  nonbondedCutoff=0.9*u.nanometers,
                                  ewaldErrorTolerance=1.0e-5)
        context = mm.Context(system, mm.VerletIntegrator(0.001), Reference)
        context.setPositions(top.positions)
        energies = energy_decomposition(top, context, nrg=u.kilojoules_per_mole)

        # Compare with Lee-Ping's answers. Make sure we zero-out forces for
        # virtual sites, since OMM doesn't do this and Gromacs does.
        self.assertAlmostEqual(energies['bond'], 1628.54739, places=3)
        self.assertAlmostEqual(energies['angle'], 3604.58751, places=3)
        self.assertAlmostEqual(energies['dihedral'], 6490.00844, delta=0.002)
        self.assertRelativeEqual(energies['nonbonded'], 23616.457584, delta=0.002)
        gmxfrc = get_forces_from_xvg(
                os.path.join(get_fn('09.DHFR-PME'), 'force.xvg'))
        ommfrc = context.getState(getForces=True).getForces().value_in_unit(
                    u.kilojoules_per_mole/u.nanometer)
        zero_ep_frc(ommfrc, top)
        max_diff = get_max_diff(gmxfrc, ommfrc)
        self.assertLess(max_diff, 5)

    def test_small_double_peptide(self):
        """ Test interacting peptides Gromacs system nrg and frc (no PBC) """
        # Load the top and gro files
        top = load_file(os.path.join(get_fn('03.AlaGlu'), 'topol.top'),
                        xyz=os.path.join(get_fn('03.AlaGlu'), 'conf.gro'))
        self.assertEqual(top.combining_rule, 'lorentz')

        # create the system and context, then calculate the energy decomposition
        system = top.createSystem()
        context = mm.Context(system, mm.VerletIntegrator(0.001), CPU)
        context.setPositions(top.positions)
        energies = energy_decomposition(top, context, nrg=u.kilojoules_per_mole)

        # Compare with Lee-Ping's answers. Make sure we zero-out forces for
        # virtual sites, since OMM doesn't do this and Gromacs does.
        self.assertAlmostEqual(energies['bond'], 1.5307999, places=4)
        self.assertAlmostEqual(energies['angle'], 6.5804488, places=4)
        self.assertAlmostEqual(energies['dihedral'], 80.379714, places=4)
        self.assertAlmostEqual(energies['nonbonded'], -275.142487, places=3)
        gmxfrc = get_forces_from_xvg(os.path.join(get_fn('03.AlaGlu'), 'force.xvg'))
        ommfrc = context.getState(getForces=True).getForces().value_in_unit(
                    u.kilojoules_per_mole/u.nanometer)
        zero_ep_frc(ommfrc, top)
        max_diff = get_max_diff(gmxfrc, ommfrc)
        self.assertLess(max_diff, 0.05)

def test_component_for_duck_typing():
    view = NGLWidget()
    traj = pt.load(nv.datafiles.PDB)
    view.add_component(get_fn('tz2.pdb'))
    view.add_component(get_fn('tz2_2.pdb.gz'))
    view.add_trajectory(nv.PyTrajTrajectory(traj))
    view.component_0.add_representation('cartoon')
    
    c0 = view[0]
    c1 = view[1]
    assert hasattr(view, 'component_0')
    assert hasattr(view, 'component_1')
    assert hasattr(view, 'trajectory_0')
    assert hasattr(view.trajectory_0, 'n_frames')
    assert hasattr(view.trajectory_0, 'get_coordinates')
    assert hasattr(view.trajectory_0, 'get_structure_string')

    c0.show()
    c0.hide()

    view.remove_component(c0.id)
    assert not hasattr(view, 'component_2')

    # negative indexing
    assert view[-1]._index == c1._index

 def testMol3SingleWriteStruct(self):
     """ Tests writing mol3 file of single-residue Structure """
     mol2 = formats.Mol2File.parse(get_fn('tripos9.mol2'), structure=True)
     formats.Mol2File.write(mol2, get_fn('tripos9struct.mol3', written=True),
                            mol3=True)
     self.assertTrue(diff_files(get_fn('tripos9struct.mol3', written=True),
                                get_saved_fn('tripos9struct.mol3')))

 def tearDown(self):
     try:
         for f in os.listdir(get_fn('writes')):
             os.unlink(get_fn(f, written=True))
         os.rmdir(get_fn('writes'))
     except OSError:
         pass

 def testWriteCharmm27Top(self):
     """ Tests writing a Gromacs topology file with CHARMM 27 FF """
     top = load_file(get_fn('1aki.charmm27.top'))
     GromacsTopologyFile.write(top,
             get_fn('1aki.charmm27.top', written=True))
     top2 = load_file(get_fn('1aki.charmm27.top', written=True))
     self._charmm27_checks(top)

 def testAmberMdcrd(self):
     """ Test writing ASCII trajectory file """
     box = [15, 15, 15]
     Mdcrd = asciicrd.AmberMdcrd
     crd = Mdcrd(get_fn('testc.mdcrd', written=True), natom=15, hasbox=False,
                 mode='w', title='Test file')
     crd.add_coordinates(list(range(45)))
     crd.add_coordinates([x+1 for x in range(45)])
     crd.add_coordinates([x+2 for x in range(45)])
     crd.add_coordinates([x+3 for x in range(45)])
     crd.add_coordinates([x+4 for x in range(45)])
     crd.close()
     crd = Mdcrd(get_fn('testcb.mdcrd', written=True), natom=18, hasbox=True,
                 mode='w', title='Test file')
     crd.add_coordinates(list(range(54)))
     crd.add_box(box)
     crd.add_coordinates([x+1 for x in range(54)])
     crd.add_box(box)                          
     crd.add_coordinates([x+2 for x in range(54)])
     crd.add_box(box)                          
     crd.add_coordinates([x+3 for x in range(54)])
     crd.add_box(box)                          
     crd.add_coordinates([x+4 for x in range(54)])
     crd.add_box(box)
     crd.close()
     self._check_written_mdcrds(box)

 def testParmSetParsing(self):
     """ Tests parsing a set of Amber parameter files """
     params = parameters.AmberParameterSet(
             os.path.join(get_fn('parm'), 'parm99.dat'),
             os.path.join(get_fn('parm'), 'frcmod.ff99SB'),
             os.path.join(get_fn('parm'), 'frcmod.parmbsc0'),
     )
     self.assertGreater(_num_unique_types(params.atom_types), 0)
     self.assertGreater(_num_unique_types(params.bond_types), 0)
     self.assertGreater(_num_unique_types(params.angle_types), 0)
     self.assertGreater(_num_unique_types(params.dihedral_types), 0)
     self.assertGreater(_num_unique_types(params.improper_periodic_types), 0)
     # Check that parameters were properly overridden. parm99.dat defines
     # C-N-CT-C torsion as follows:
     #
     # C -N -CT-C    1    0.850       180.000          -2.
     # C -N -CT-C    1    0.800         0.000           1.
     #
     # whereas ff99SB defines that torsion as:
     #
     # C -N -CT-C    1    0.00          0.0            -4.
     # C -N -CT-C    1    0.42          0.0            -3.
     # C -N -CT-C    1    0.27          0.0            -2.
     # C -N -CT-C    1    0.00          0.0             1.
     #
     # Since ff99SB is loaded last, this should be the one that is stored
     self.assertEqual(len(params.dihedral_types[('C','N','CT','C')]), 4)
     self.assertEqual(params.dihedral_types[('C','N','CT','C')][0],
                      topologyobjects.DihedralType(0, 4, 0, 1.2, 2.0))
     self.assertEqual(params.dihedral_types[('C','N','CT','C')][1],
                      topologyobjects.DihedralType(0.42, 3, 0, 1.2, 2.0))
     self.assertEqual(params.dihedral_types[('C','N','CT','C')][2],
                      topologyobjects.DihedralType(0.27, 2, 0, 1.2, 2.0))
     self.assertEqual(params.dihedral_types[('C','N','CT','C')][3],
                      topologyobjects.DihedralType(0, 1, 0, 1.2, 2.0))

 def testBadFileUsage(self):
     """ Check that illegal file usage results in desired exceptions """
     Restart = asciicrd.AmberAsciiRestart
     Mdcrd = asciicrd.AmberMdcrd
     box = [10, 10, 10, 90, 90, 90]
     rst = Restart(get_fn('testc.rst7', written=True), 'w', natom=9,
                   hasbox=True)
     def assign(obj, stmnt):
         rst = crd = obj
         exec(stmnt)
     try:
         self.assertRaises(ValueError,
                           lambda: assign(rst, 'rst.coordinates=range(20)'))
         self.assertRaises(RuntimeError,
                           lambda: assign(rst, 'rst.box=[10]*3+[90]*3'))
         rst.coordinates = list(range(27))
         rst.box = box
         self.assertRaises(RuntimeError, lambda:
                           assign(rst, 'rst.velocities=list(range(27))'))
     finally:
         rst.close()
     crd = Mdcrd(get_fn('testc.mdcrd', written=True), natom=15, hasbox=True,
                 mode='w', title='Test file')
     s = 'list(range(45))'
     s2 = 'list(range(42))'
     try:
         crd.add_coordinates(eval(s))
         self.assertRaises(RuntimeError,
                           lambda: assign(crd, 'crd.add_coordinates(%s)'%s))
         crd.add_box([10, 10, 10])
         self.assertRaises(ValueError,
                           lambda: assign(crd, 'crd.add_coordinates(%s)'%s2))
     finally:
         crd.close()

 def _check_written_restarts(self, box):
     # Now try to read them and verify the information (keep in mind that the
     # restart velocities are scaled down then back up, so you'll need to use
     # assertAlmostEqual in this case).
     rst = readparm.Rst7.open(get_fn('testc.rst7', written=True))
     self.assertFalse(rst.hasbox)
     self.assertFalse(rst.hasvels)
     np.testing.assert_equal(rst.coordinates.flatten(), list(range(27)))
     rst = readparm.Rst7.open(get_fn('testcb.rst7', written=True))
     self.assertTrue(rst.hasbox)
     self.assertFalse(rst.hasvels)
     np.testing.assert_equal(rst.coordinates.flatten(), list(range(21)))
     np.testing.assert_equal(rst.box.flatten(), box)
     rst = readparm.Rst7.open(get_fn('testcv.rst7', written=True))
     self.assertTrue(rst.hasvels)
     self.assertFalse(rst.hasbox)
     np.testing.assert_equal(rst.coordinates,
             np.arange(60).reshape(rst.coordinates.shape))
     np.testing.assert_allclose(rst.velocities,
             np.array(list(reversed(range(60)))).reshape(rst.velocities.shape))
     rst = readparm.Rst7.open(get_fn('testcvb.rst7', written=True))
     self.assertTrue(rst.hasvels)
     self.assertTrue(rst.hasbox)
     np.testing.assert_equal(rst.coordinates,
             np.arange(45).reshape(rst.coordinates.shape))
     np.testing.assert_allclose(rst.velocities,
             np.array(list(reversed(range(45)))).reshape(rst.velocities.shape))
     np.testing.assert_equal(rst.box, box)