def test_structure(self):
        quartz = self.quartz.structure
        np.testing.assert_array_equal(quartz.lattice.matrix,
                                      [[4.913400, 0, 0],
                                       [-2.456700, 4.255129, 0],
                                       [0, 0, 5.405200]])
        self.assertEqual(quartz.formula, "Si3 O6")
        self.assertNotIn("molecule-ID", self.quartz.atoms.columns)

        ethane = self.ethane.structure
        np.testing.assert_array_equal(ethane.lattice.matrix,
                                      np.diag([10.0] * 3))
        lbounds = np.array(self.ethane.box.bounds)[:, 0]
        coords = self.ethane.atoms[["x", "y", "z"]].values - lbounds
        np.testing.assert_array_equal(ethane.cart_coords, coords)
        np.testing.assert_array_equal(ethane.site_properties["charge"],
                                      self.ethane.atoms["q"])
        tatb = self.tatb.structure
        frac_coords = tatb.frac_coords[381]
        real_frac_coords = frac_coords - np.floor(frac_coords)
        np.testing.assert_array_almost_equal(real_frac_coords,
                                             [0.01553397,
                                              0.71487872,
                                              0.14134139])

        co = Structure.from_spacegroup(194,
                                       Lattice.hexagonal(2.50078, 4.03333),
                                       ["Co"], [[1/3, 2/3, 1/4]])
        ld_co = LammpsData.from_structure(co)
        self.assertEqual(ld_co.structure.composition.reduced_formula, "Co")
        ni = Structure.from_spacegroup(225, Lattice.cubic(3.50804),
                                       ["Ni"], [[0, 0, 0]])
        ld_ni = LammpsData.from_structure(ni)
        self.assertEqual(ld_ni.structure.composition.reduced_formula, "Ni")

    def from_file(name, filename, data_filename=None, is_forcefield=False):
        """
        Read in the input settings from json file as ordereddict. Also
        reads in the datafile if provided.
        Note: with monty.serialization.loadfn the order of paramters in the
        json file is not preserved

        Args:
            filename (string): name of the file with the lamps control
                paramters
            data_filename (string): path to the data file
            is_forcefield (bool): whether the data file has forcefield and
                topology info in it.

        Returns:
            DictLammpsInput
        """
        with open(filename) as f:
            config_dict = json.load(f, object_pairs_hook=OrderedDict)
        lammps_data = None
        if data_filename:
            if is_forcefield:
                lammps_data = LammpsForceFieldData.from_file(data_filename)
            else:
                lammps_data = LammpsData.from_file(data_filename)
        return DictLammpsInput(name, config_dict, lammps_data)

    def from_file(name, filename, lammps_data=None, data_filename="in.data",
                  user_lammps_settings={}, is_forcefield=False):
        """
        Read in the input settings from json file as ordereddict.
        Note: with monty.serialization.loadfn the order of paramters in the
        json file is not preserved

        Args:
            filename (string): name of the file with the lamps control
                paramters
            lammps_data (string/LammpsData/LammpsForceFieldData): path to the
                data file or an appropriate object
            data_filename (string): name of the the lammps data file
            user_lammps_settings (dict): User lammps settings
            is_forcefield (bool): whether the data file has forcefield and
                topology info in it. This is required only if lammps_data is
                a path to the data file instead of a data object

        Returns:
            DictLammpsInput
        """
        with open(filename) as f:
            config_dict = json.load(f, object_pairs_hook=OrderedDict)
        lammps_data = lammps_data
        if isinstance(lammps_data, six.string_types):
            if is_forcefield:
                lammps_data = LammpsForceFieldData.from_file(lammps_data)
            else:
                lammps_data = LammpsData.from_file(lammps_data)
        return DictLammpsInput(name, config_dict, lammps_data=lammps_data,
                               data_filename=data_filename,
                               user_lammps_settings=user_lammps_settings)

    def from_file(cls, name, input_template, user_settings,
                  lammps_data=None, data_filename="in.data",
                  is_forcefield=False):
        """
        Returns LammpsInputSet from  input file template and input data.

        Args:
            name (str)
            input_template (string): path to the input template file.
            user_settings (dict): User lammps settings, the keys must
                correspond to the keys in the template.
            lammps_data (string/LammpsData/LammpsForceFieldData): path to the
                data file or an appropriate object
            data_filename (string): name of the the lammps data file.
            is_forcefield (bool): whether the data file has forcefield and
                topology info in it. This is required only if lammps_data is
                a path to the data file instead of a data object.

        Returns:
            LammpsInputSet
        """
        user_settings["data_file"] = data_filename
        lammps_input = LammpsInput.from_file(input_template, user_settings)
        if isinstance(lammps_data, six.string_types):
            if is_forcefield:
                lammps_data = LammpsForceFieldData.from_file(lammps_data)
            else:
                lammps_data = LammpsData.from_file(lammps_data)
        return cls(name, lammps_input, lammps_data=lammps_data,
                   data_filename=data_filename)

 def setUpClass(cls):
     h2o_coords = [[9.626, 6.787, 12.673],
                   [9.626, 8.420, 12.673],
                   [10.203, 7.604, 12.673]]
     molecule = Molecule(["H", "H", "O"], h2o_coords)
     box_size = [[0.0, 10.0], [0.0, 10.0], [0.0, 10.0]]
     cls.lammps_data = LammpsData.from_structure(molecule, box_size)

    def test_from_ff_and_topologies(self):
        mass = OrderedDict()
        mass["H"] = 1.0079401
        mass["O"] = 15.999400
        nonbond_coeffs = [[0.00774378, 0.98], [0.1502629, 3.1169]]
        topo_coeffs = {"Bond Coeffs": [{"coeffs": [176.864, 0.9611],
                                        "types": [("H", "O")]}],
                       "Angle Coeffs": [{"coeffs": [42.1845, 109.4712],
                                         "types": [("H", "O", "H")]}]}
        ff = ForceField(mass.items(), nonbond_coeffs, topo_coeffs)
        with gzip.open(os.path.join(test_dir, "topologies_ice.json.gz")) as f:
            topo_dicts = json.load(f)
        topologies = [Topology.from_dict(d) for d in topo_dicts]
        box_bounds = [[-0.75694412, 44.165558],
                      [0.38127473, 47.066074],
                      [0.17900842, 44.193867]]
        ice = LammpsData.from_ff_and_topologies(ff=ff, topologies=topologies,
                                                box_bounds=box_bounds)
        atoms = ice.atoms
        bonds = ice.topology["Bonds"]
        angles = ice.topology["Angles"]
        np.testing.assert_array_equal(atoms.index.values,
                                      np.arange(1, len(atoms) + 1))
        np.testing.assert_array_equal(bonds.index.values,
                                      np.arange(1, len(bonds) + 1))
        np.testing.assert_array_equal(angles.index.values,
                                      np.arange(1, len(angles) + 1))

        i = random.randint(0, len(topologies) - 1)
        sample = topologies[i]
        in_atoms = ice.atoms[ice.atoms["molecule-ID"] == i + 1]
        np.testing.assert_array_equal(in_atoms.index.values,
                                      np.arange(3 * i + 1, 3 * i + 4))
        np.testing.assert_array_equal(in_atoms["type"].values, [2, 1, 1])
        np.testing.assert_array_equal(in_atoms["q"].values, sample.charges)
        np.testing.assert_array_equal(in_atoms[["x", "y", "z"]].values,
                                      sample.sites.cart_coords)
        broken_topo_coeffs = {"Bond Coeffs": [{"coeffs": [176.864, 0.9611],
                                               "types": [("H", "O")]}],
                              "Angle Coeffs": [{"coeffs": [42.1845, 109.4712],
                                                "types": [("H", "H", "H")]}]}
        broken_ff = ForceField(mass.items(), nonbond_coeffs,
                               broken_topo_coeffs)
        ld_woangles = LammpsData.from_ff_and_topologies(ff=broken_ff,
                                                        topologies=[sample],
                                                        box_bounds=box_bounds)
        self.assertNotIn("Angles", ld_woangles.topology)

 def test_write_file(self):
     filename1 = "test1.data"
     self.ethane.write_file(filename=filename1)
     c2h6 = LammpsData.from_file(filename1)
     pd.testing.assert_frame_equal(c2h6.masses, self.ethane.masses)
     pd.testing.assert_frame_equal(c2h6.atoms, self.ethane.atoms)
     ff_kw = random.sample(self.ethane.force_field.keys(), 1)[0]
     pd.testing.assert_frame_equal(c2h6.force_field[ff_kw],
                                   self.ethane.force_field[ff_kw], ff_kw)
     topo_kw = random.sample(self.ethane.topology.keys(), 1)[0]
     pd.testing.assert_frame_equal(c2h6.topology[topo_kw],
                                   self.ethane.topology[topo_kw], topo_kw)
     filename2 = "test2.data"
     self.virus.write_file(filename=filename2)
     v = LammpsData.from_file(filename2, atom_style="angle")
     pd.testing.assert_frame_equal(v.force_field["PairIJ Coeffs"],
                                   self.virus.force_field["PairIJ Coeffs"])

 def __init__(self, data_file, trajectory_file, log_file="log.lammps"):
     self.data_file = os.path.abspath(data_file)
     self.trajectory_file = os.path.abspath(trajectory_file)
     self.log_file = os.path.abspath(log_file)
     self.log = LammpsLog(log_file)
     self.lammps_data = LammpsData.from_file(self.data_file)
     self._set_mol_masses_and_charges()
     self._parse_trajectory()

 def __init__(self, data_file, trajectory_file, log_file="log.lammps",
              is_forcefield=False):
     self.data_file = data_file
     self.trajectory_file = trajectory_file
     self.log_file = log_file
     self.lammps_log = LammpsLog(log_file)
     if is_forcefield:
         self.lammps_data = LammpsForceFieldData.from_file(data_file)
     else:
         self.lammps_data = LammpsData.from_file(data_file)
     self._set_mol_masses_and_charges()
     self._parse_trajectory()

    def test_md(self):
        s = Structure.from_spacegroup(225, Lattice.cubic(3.62126),
                                      ["Cu"], [[0, 0, 0]])
        ld = LammpsData.from_structure(s, atom_style="atomic")
        ff = "\n".join(["pair_style eam", "pair_coeff * * Cu_u3.eam"])
        md = LammpsRun.md(data=ld, force_field=ff, temperature=1600.0,
                          nsteps=10000)
        md.write_inputs(output_dir="md")
        with open(os.path.join("md", "in.md")) as f:
            md_script = f.read()
        script_string = """# Sample input script template for MD

# Initialization

units           metal
atom_style      atomic

# Atom definition

read_data       md.data
#read_restart    md.restart

# Force field settings (consult official document for detailed formats)

pair_style eam
pair_coeff * * Cu_u3.eam

# Create velocities
velocity        all create 1600.0 142857 mom yes rot yes dist gaussian

# Ensemble constraints
#fix             1 all nve
fix             1 all nvt temp 1600.0 1600.0 0.1
#fix             1 all npt temp 1600.0 1600.0 0.1 iso $pressure $pressure 1.0

# Various operations within timestepping
#fix             ...
#compute         ...

# Output settings
#thermo_style    custom ...  # control the thermo data type to output
thermo          100  # output thermo data every N steps
#dump            1 all atom 100 traj.*.gz  # dump a snapshot every 100 steps

# Actions
run             10000
"""
        self.assertEqual(md_script, script_string)
        self.assertTrue(os.path.exists(os.path.join("md", "md.data")))

    def test_write_lammps_inputs(self):
        # script template
        with open(os.path.join(test_dir, "kappa.txt")) as f:
            kappa_template = f.read()
        kappa_settings = {"method": "heat"}
        write_lammps_inputs(output_dir="heat", script_template=kappa_template,
                            settings=kappa_settings)
        with open(os.path.join("heat", "in.lammps")) as f:
            kappa_script = f.read()
        fix_hot = re.search(r"fix\s+hot\s+all\s+([^\s]+)\s+", kappa_script)
        # placeholders supposed to be filled
        self.assertEqual(fix_hot.group(1), "heat")
        fix_cold = re.search(r"fix\s+cold\s+all\s+([^\s]+)\s+", kappa_script)
        self.assertEqual(fix_cold.group(1), "heat")
        lattice = re.search(r"lattice\s+fcc\s+(.*)\n", kappa_script)
        # parentheses not supposed to be filled
        self.assertEqual(lattice.group(1), "${rho}")
        pair_style = re.search(r"pair_style\slj/cut\s+(.*)\n", kappa_script)
        self.assertEqual(pair_style.group(1), "${rc}")

        with open(os.path.join(test_dir, "in.peptide")) as f:
            peptide_script = f.read()
        # copy data file
        src = os.path.join(test_dir, "data.quartz")
        write_lammps_inputs(output_dir="path", script_template=peptide_script,
                            data=src)
        dst = os.path.join("path", "data.peptide")
        self.assertTrue(filecmp.cmp(src, dst, shallow=False))
        # write data file from obj
        obj = LammpsData.from_file(src, atom_style="atomic")
        write_lammps_inputs(output_dir="obj", script_template=peptide_script,
                            data=obj)
        obj_read = LammpsData.from_file(os.path.join("obj", "data.peptide"),
                                        atom_style="atomic")
        pd.testing.assert_frame_equal(obj_read.masses, obj.masses)
        pd.testing.assert_frame_equal(obj_read.atoms, obj.atoms)

    def from_structure(cls, structure, input_template, user_settings,
                       data_filename="in.data", name="basic"):
        """
        Returns inputset from structure

        Args:
            structure (Structure/Molecule):
            input_template (string): path to the input template file.
            user_settings (dict): User lammps settings, the keys must
                correspond to the keys in the template.
            data_filename (string): name of the the lammps data file.

        Returns:
            LammpsInputSet
        """

        lammps_data = LammpsData.from_structure(structure)
        return cls.from_file(name, input_template, user_settings,
                             lammps_data=lammps_data, data_filename=data_filename)

 def test_json_dict(self):
     encoded = json.dumps(self.ethane.as_dict())
     decoded = json.loads(encoded)
     c2h6 = LammpsData.from_dict(decoded)
     pd.testing.assert_frame_equal(c2h6.masses, self.ethane.masses)
     pd.testing.assert_frame_equal(c2h6.atoms, self.ethane.atoms)
     ff = self.ethane.force_field
     key, target_df = random.sample(ff.items(), 1)[0]
     self.assertIsNone(
         pd.testing.assert_frame_equal(c2h6.force_field[key], target_df,
                                       check_dtype=False),
         key
     )
     topo = self.ethane.topology
     key, target_df = random.sample(topo.items(), 1)[0]
     self.assertIsNone(
         pd.testing.assert_frame_equal(c2h6.topology[key], target_df),
         key
     )

 def test_from_structure(self):
     latt = Lattice.monoclinic(9.78746, 4.75058, 8.95892, 115.9693)
     structure = Structure.from_spacegroup(15, latt, ["Os", "O", "O"],
                                           [[0, 0.25583, 0.75],
                                            [0.11146, 0.46611, 0.91631],
                                            [0.11445, 0.04564, 0.69518]])
     velocities = np.random.randn(20, 3) * 0.1
     structure.add_site_property("velocities", velocities)
     ld = LammpsData.from_structure(structure=structure,
                                    ff_elements=["O", "Os", "Na"])
     i = random.randint(0, 19)
     a = latt.matrix[0]
     va = velocities[i].dot(a) / np.linalg.norm(a)
     self.assertAlmostEqual(va, ld.velocities.loc[i + 1, "vx"])
     self.assertAlmostEqual(velocities[i, 1],
                            ld.velocities.loc[i + 1, "vy"])
     np.testing.assert_array_almost_equal(ld.masses["mass"],
                                          [22.989769, 190.23, 15.9994])
     np.testing.assert_array_equal(ld.atoms["type"], [2] * 4 + [3] * 16)

    def from_file(cls, name, input_template, user_settings,
                  lammps_data=None, data_filename="in.data"):
        """
        Returns LammpsInputSet from  input file template and input data.

        Args:
            name (str)
            input_template (string): path to the input template file.
            user_settings (dict): User lammps settings, the keys must
                correspond to the keys in the template.
            lammps_data (string/LammpsData): path to the
                data file or an appropriate object
            data_filename (string): name of the the lammps data file.

        Returns:
            LammpsInputSet
        """
        user_settings["data_file"] = data_filename
        lammps_input = LammpsInput.from_file(input_template, user_settings)
        if isinstance(lammps_data, six.string_types):
            lammps_data = LammpsData.from_file(lammps_data)
        return cls(name, lammps_input, lammps_data=lammps_data,
                   data_filename=data_filename)

 def test_from_file(self):
     # general tests
     pep = self.peptide
     # header stats and Nos. of columns
     self.assertEqual(pep.masses.shape, (14, 1))
     self.assertEqual(pep.atoms.shape, (2004, 9))
     self.assertListEqual(list(pep.atoms.columns),
                          ["molecule-ID", "type", "q", "x", "y", "z",
                           "nx", "ny", "nz"])
     topo = pep.topology
     self.assertEqual(topo["Bonds"].shape, (1365, 3))
     self.assertEqual(topo["Angles"].shape, (786, 4))
     self.assertEqual(topo["Dihedrals"].shape, (207, 5))
     self.assertEqual(topo["Impropers"].shape, (12, 5))
     ff = pep.force_field
     self.assertEqual(ff["Pair Coeffs"].shape, (14, 4))
     self.assertEqual(ff["Bond Coeffs"].shape, (18, 2))
     self.assertEqual(ff["Angle Coeffs"].shape, (31, 4))
     self.assertEqual(ff["Dihedral Coeffs"].shape, (21, 4))
     self.assertEqual(ff["Improper Coeffs"].shape, (2, 2))
     # header box
     np.testing.assert_array_equal(pep.box_bounds,
                                   [[36.840194, 64.211560],
                                    [41.013691, 68.385058],
                                    [29.768095, 57.139462]])
     # body
     self.assertEqual(pep.masses.at[7, "mass"], 12.0110)
     self.assertEqual(ff["Pair Coeffs"].at[9, "coeff3"], 0.152100)
     self.assertEqual(ff["Bond Coeffs"].at[5, "coeff2"], 1.430000)
     self.assertEqual(ff["Angle Coeffs"].at[21, "coeff2"], 120.000000)
     self.assertEqual(ff["Dihedral Coeffs"].at[10, "coeff1"], 0.040000)
     self.assertEqual(ff["Improper Coeffs"].at[2, "coeff1"], 20.000000)
     self.assertEqual(pep.atoms.at[29, "molecule-ID"], 1)
     self.assertEqual(pep.atoms.at[29, "type"], 7)
     self.assertEqual(pep.atoms.at[29, "q"], -0.020)
     self.assertAlmostEqual(pep.atoms.at[29, "x"], 42.96709)
     self.assertEqual(pep.atoms.at[1808, "molecule-ID"], 576)
     self.assertEqual(pep.atoms.at[1808, "type"], 14)
     self.assertAlmostEqual(pep.atoms.at[1808, "y"], 58.64352)
     self.assertEqual(pep.atoms.at[1808, "nx"], -1)
     self.assertAlmostEqual(pep.velocities.at[527, "vz"], -0.010889)
     self.assertEqual(topo["Bonds"].at[47, "type"], 8)
     self.assertEqual(topo["Bonds"].at[47, "atom2"], 54)
     self.assertEqual(topo["Bonds"].at[953, "atom1"], 1384)
     self.assertEqual(topo["Angles"].at[105, "type"], 19)
     self.assertEqual(topo["Angles"].at[105, "atom3"], 51)
     self.assertEqual(topo["Angles"].at[376, "atom2"], 772)
     self.assertEqual(topo["Dihedrals"].at[151, "type"], 14)
     self.assertEqual(topo["Dihedrals"].at[151, "atom4"], 51)
     self.assertEqual(topo["Impropers"].at[4, "atom4"], 32)
     # class 2 and comments
     ethane = self.ethane
     self.assertEqual(ethane.masses.shape, (2, 1))
     self.assertEqual(ethane.atoms.shape, (8, 9))
     class2 = ethane.force_field
     self.assertEqual(class2["Pair Coeffs"].shape, (2, 2))
     self.assertEqual(class2["Bond Coeffs"].shape, (2, 4))
     self.assertEqual(class2["Angle Coeffs"].shape, (2, 4))
     self.assertEqual(class2["Dihedral Coeffs"].shape, (1, 6))
     self.assertEqual(class2["Improper Coeffs"].shape, (2, 2))
     self.assertEqual(class2["BondBond Coeffs"].at[2, "coeff3"], 1.1010)
     self.assertEqual(class2["BondAngle Coeffs"].at[2, "coeff4"], 1.1010)
     self.assertEqual(class2["AngleAngle Coeffs"].at[2, "coeff6"],
                      107.6600)
     self.assertEqual(class2["AngleAngle Coeffs"].at[2, "coeff6"],
                      107.6600)
     self.assertEqual(class2["AngleAngleTorsion Coeffs"].at[1, "coeff3"],
                      110.7700)
     self.assertEqual(class2["EndBondTorsion Coeffs"].at[1, "coeff8"],
                      1.1010)
     self.assertEqual(class2["MiddleBondTorsion Coeffs"].at[1, "coeff4"],
                      1.5300)
     self.assertEqual(class2["BondBond13 Coeffs"].at[1, "coeff3"], 1.1010)
     self.assertEqual(class2["AngleTorsion Coeffs"].at[1, "coeff8"],
                      110.7700)
     # box_tilt and another atom_style
     quartz = self.quartz
     np.testing.assert_array_equal(quartz.box_tilt, [-2.456700, 0.0, 0.0])
     self.assertListEqual(list(quartz.atoms.columns),
                          ["type", "x", "y", "z"])
     self.assertAlmostEqual(quartz.atoms.at[7, "x"], 0.299963)
     # PairIJ Coeffs section
     virus = self.virus
     pairij = virus.force_field["PairIJ Coeffs"]
     self.assertEqual(pairij.at[7, "id1"], 3)
     self.assertEqual(pairij.at[7, "id2"], 3)
     self.assertEqual(pairij.at[7, "coeff2"], 2.1)
     # sort_id
     nvt = LammpsData.from_file(filename=os.path.join(test_dir,
                                                      "nvt.data"),
                                sort_id=True)
     atom_id = random.randint(1, 648)
     self.assertEqual(nvt.atoms.loc[atom_id].name, atom_id)

    def test_disassemble(self):
        # general tests
        c = LammpsData.from_file(os.path.join(test_dir, "crambin.data"))
        _, c_ff, topos = c.disassemble()
        mass_info = [('N1', 14.0067), ('H1', 1.00797), ('C1', 12.01115),
                     ('H2', 1.00797), ('C2', 12.01115), ('O1', 15.9994),
                     ('C3', 12.01115), ('O2', 15.9994), ('H3', 1.00797),
                     ('C4', 12.01115), ('N2', 14.0067), ('C5', 12.01115),
                     ('S1', 32.064), ('C6', 12.01115), ('N3', 14.0067),
                     ('C7', 12.01115), ('C8', 12.01115), ('C9', 12.01115),
                     ('O3', 15.9994)]
        self.assertListEqual(c_ff.mass_info, mass_info)
        np.testing.assert_array_equal(c_ff.nonbond_coeffs,
                                      c.force_field["Pair Coeffs"].values)
        base_kws = ["Bond", "Angle", "Dihedral", "Improper"]
        for kw in base_kws:
            ff_kw = kw + " Coeffs"
            i = random.randint(0, len(c_ff.topo_coeffs[ff_kw]) - 1)
            sample_coeff = c_ff.topo_coeffs[ff_kw][i]
            np.testing.\
                assert_array_equal(sample_coeff["coeffs"],
                                   c.force_field[ff_kw].iloc[i].values,
                                   ff_kw)
        topo = topos[-1]
        atoms = c.atoms[c.atoms["molecule-ID"] == 46]
        np.testing.assert_array_equal(topo.sites.cart_coords,
                                      atoms[["x", "y", "z"]])
        np.testing.assert_array_equal(topo.charges, atoms["q"])
        atom_labels = [m[0] for m in mass_info]
        self.assertListEqual(topo.sites.site_properties["ff_map"],
                             [atom_labels[i - 1] for i in atoms["type"]])
        shift = min(atoms.index)
        for kw in base_kws:
            ff_kw = kw + " Coeffs"
            ff_coeffs = c_ff.topo_coeffs[ff_kw]
            topo_kw = kw + "s"
            topos_df = c.topology[topo_kw]
            topo_df = topos_df[topos_df["atom1"] >= shift]
            topo_arr = topo_df.drop("type", axis=1).values
            np.testing.assert_array_equal(topo.topologies[topo_kw],
                                          topo_arr - shift, topo_kw)
            sample_topo = random.sample(list(topo_df.itertuples(False, None)),
                                        1)[0]
            topo_type_idx = sample_topo[0] - 1
            topo_type = tuple([atom_labels[i - 1] for i in
                               atoms.loc[sample_topo[1:], "type"]])

            self.assertIn(topo_type, ff_coeffs[topo_type_idx]["types"], ff_kw)
        # test no guessing element and pairij as nonbond coeffs
        v = self.virus
        _, v_ff, _ = v.disassemble(guess_element=False)
        self.assertDictEqual(v_ff.maps["Atoms"],
                             dict(Qa1=1, Qb1=2, Qc1=3, Qa2=4))
        pairij_coeffs = v.force_field["PairIJ Coeffs"].drop(["id1", "id2"],
                                                            axis=1)
        np.testing.assert_array_equal(v_ff.nonbond_coeffs,
                                      pairij_coeffs.values)
        # test class2 ff
        _, e_ff, _ = self.ethane.disassemble()
        e_topo_coeffs = e_ff.topo_coeffs
        for k in ["BondBond Coeffs", "BondAngle Coeffs"]:
            self.assertIn(k, e_topo_coeffs["Angle Coeffs"][0], k)
        for k in ["MiddleBondTorsion Coeffs", "EndBondTorsion Coeffs",
                  "AngleTorsion Coeffs", "AngleAngleTorsion Coeffs",
                  "BondBond13 Coeffs"]:
            self.assertIn(k, e_topo_coeffs["Dihedral Coeffs"][0], k)
        self.assertIn("AngleAngle Coeffs",
                      e_topo_coeffs["Improper Coeffs"][0])

 def test_from_file(self):
     peptide = self.peptide
     # header stats
     self.assertEqual(len(peptide.atoms), 2004)
     topology = peptide.topology
     self.assertEqual(len(topology["Bonds"]), 1365)
     self.assertEqual(len(topology["Angles"]), 786)
     self.assertEqual(len(topology["Dihedrals"]), 207)
     self.assertEqual(len(topology["Impropers"]), 12)
     self.assertEqual(len(peptide.masses), 14)
     force_field = peptide.force_field
     self.assertEqual(len(force_field["Pair Coeffs"]), 14)
     self.assertEqual(len(force_field["Bond Coeffs"]), 18)
     self.assertEqual(len(force_field["Angle Coeffs"]), 31)
     self.assertEqual(len(force_field["Dihedral Coeffs"]), 21)
     self.assertEqual(len(force_field["Improper Coeffs"]), 2)
     # header box
     np.testing.assert_array_equal(peptide.box_bounds,
                                   [[36.840194, 64.211560],
                                    [41.013691, 68.385058],
                                    [29.768095, 57.139462]])
     # body  last line of each section
     self.assertDictEqual(peptide.masses[-1], {"id": 14, "mass": 1.0100})
     self.assertDictEqual(force_field["Pair Coeffs"][-1],
                          {"id": 14, "coeffs": [0.046000, 0.400014,
                                                0.046000, 0.400014]})
     self.assertDictEqual(force_field["Bond Coeffs"][-1],
                          {"id": 18, "coeffs": [450.000000, 0.957200]})
     self.assertDictEqual(force_field["Angle Coeffs"][-1],
                          {"id": 31, "coeffs": [55.000000, 104.520000,
                                                0.000000, 0.000000]})
     self.assertDictEqual(force_field["Dihedral Coeffs"][-1],
                          {"id": 21, "coeffs": [0.010000, 3, 0, 1.000000]})
     for c in force_field["Dihedral Coeffs"][-1]["coeffs"][1:2]:
         self.assertIsInstance(c, int)
     self.assertDictEqual(force_field["Improper Coeffs"][-1],
                          {"id": 2, "coeffs": [20.000000, 0.000000]})
     self.assertDictEqual(peptide.atoms[-1],
                          {"id": 2004, "molecule-ID": 641, "type": 14,
                           "q": 0.417, "x": 56.55074, "y": 49.75049,
                           "z": 48.61854, "nx": 1, "ny": 1, "nz": 1})
     self.assertDictEqual(peptide.velocities[-1],
                          {"id": 2004, "velocity": [-0.010076,
                                                    -0.005729,
                                                    -0.026032]})
     self.assertDictEqual(topology["Bonds"][-1],
                          {"id": 1365, "type": 18, "bond": [2002, 2003]})
     self.assertDictEqual(topology["Angles"][-1],
                          {"id": 786, "type": 31,
                           "angle": [2003, 2002, 2004]})
     self.assertDictEqual(topology["Dihedrals"][-1],
                          {"id": 207, "type": 3,
                           "dihedral": [64, 79, 80, 82]})
     self.assertDictEqual(topology["Impropers"][-1],
                          {"id": 12, "type": 2,
                           "improper": [79, 64, 80, 81]})
     # box_tilt and another atom_style
     quartz = self.quartz
     np.testing.assert_array_equal(quartz.box_tilt, [-2.456700, 0.0, 0.0])
     self.assertDictEqual(quartz.atoms[-1],
                          {"id": 9, "type": 2, "x": 1.375998,
                           "y": -1.140800, "z": -2.443511})
     # sort_id
     nvt = LammpsData.from_file(filename=os.path.join(test_dir,
                                                      "nvt.data"),
                                sort_id=True)
     atom_id = random.randint(1, 648)
     self.assertEqual(nvt.atoms[atom_id - 1]["id"], atom_id)
     # PairIJ Coeffs section
     virus = LammpsData.from_file(filename=os.path.join(test_dir,
                                                        "virus.data"),
                                  atom_style="angle")
     n = len(virus.masses)
     pairij = virus.force_field["PairIJ Coeffs"]
     self.assertEqual(len(pairij), n * (n + 1) / 2)
     self.assertDictEqual(pairij[-1],
                          {"id1": 4, "id2": 4, "coeffs": [1, 1, 1.1225]})

 def setUpClass(cls):
     polymer_chain = Molecule.from_file(os.path.join(test_dir, "polymer_chain.xyz"))
     box_size = [[0.0, 20.0], [0.0, 20.0], [0.0, 20.0]]
     cls.lammps_data = LammpsData.from_structure(polymer_chain, box_size)

 def test_from_file(self):
     self.lammps_data.write_data_file(os.path.join(test_dir, "lammps_data.dat"))
     lammps_data = LammpsData.from_file(os.path.join(test_dir, "lammps_data.dat"))
     self.assertEqual(str(lammps_data), str(self.lammps_data))