lammpsdata.py

"""Read and write LAMMPS data files."""

from __future__ import division

import copy
import datetime
import os
import warnings
from itertools import count
from pathlib import Path

import numpy as np
import unyt as u
from unyt.array import allclose_units

import gmso
from gmso.abc.abstract_site import Molecule
from gmso.core.angle import Angle
from gmso.core.atom import Atom
from gmso.core.atom_type import AtomType
from gmso.core.bond import Bond
from gmso.core.box import Box
from gmso.core.dihedral import Dihedral
from gmso.core.element import element_by_mass
from gmso.core.improper import Improper
from gmso.core.topology import Topology
from gmso.core.views import PotentialFilters
from gmso.formats.formats_registry import loads_as, saves_as
from gmso.lib.potential_templates import PotentialTemplateLibrary
from gmso.utils.compatibility import check_compatibility
from gmso.utils.conversions import convert_kelvin_to_energy_units
from gmso.utils.sorting import (
    reindex_molecules,
    sort_by_types,
    sort_connection_members,
)
from gmso.utils.units import LAMMPS_UnitSystems, write_out_parameter_and_units

pfilter = PotentialFilters.UNIQUE_SORTED_NAMES


# TODO: Write in header of each potential type any conversions that happened
# TODO: write in file header the source of the xml
@saves_as(".lammps", ".lammpsdata", ".data")
def write_lammpsdata(
    top,
    filename,
    atom_style="full",
    unit_style="real",
    strict_potentials=False,
    strict_units=False,
    lj_cfactorsDict=None,
):
    """Output a LAMMPS data file.

    Outputs a LAMMPS data file in the 'full' atom style format.
    Assumes use of 'real' units.
    See http://lammps.sandia.gov/doc/atom_style.html for more information on atom styles.

    Parameters
    ----------
    Topology : `Topology`
        A Topology Object
    filename : str
        Path of the output file
    atom_style : str, optional, default='full'
        Defines the style of atoms to be saved in a LAMMPS data file.
        The following atom styles are currently supported: 'full', 'atomic', 'charge', 'molecular'
        see http://lammps.sandia.gov/doc/atom_style.html for more information on atom styles.
    unit_style : str, optional, default='real'
        Can be any of "real", "lj", "metal", "si", "cgs", "electron", "micro", "nano". Otherwise
        an error will be thrown. These are defined in _unit_style_factory. See
        https://docs.lammps.org/units.html for LAMMPS documentation.
    strict_potentials : bool, optional, default False
        Tells the writer how to treat conversions. If False, then check for conversions
        to usable potential styles found in default_parameterMaps. If True, then error if
        potentials are not compatible.
    strict_units : bool, optional, default False
        Tells the writer how to treat unit conversions. If False, then check for conversions
        to unit styles defined in _unit_style_factory. If True, then error if parameter units
        do not match.
    lj_cfactorsDict : (None, dict), optional, default None
        If using unit_style="lj" only, can pass a dictionary with keys of ("mass", "energy",
        "length", "charge"), or any combination of these, and they will be used to non-
        dimensionalize all values in the topology. If any key is not passed, default values
        will be pulled from the topology (see _default_lj_val). These are the largest: sigma,
        epsilon, atomtype.mass, and atomtype.charge from the topology.

    Notes
    -----
    See http://lammps.sandia.gov/doc/2001/data_format.html for a full description of the LAMMPS data format.
    This is a work in progress, as only a subset of everything LAMMPS supports is currently available.
    However, please raise issues as the current writer has been set up to eventually grow to support
    all LAMMPS styles.

    Some of this function has been adopted from `mdtraj`'s support of the LAMMPSTRJ trajectory format.
    See https://github.com/mdtraj/mdtraj/blob/master/mdtraj/formats/lammpstrj.py for details.

    """
    if atom_style not in ["full", "atomic", "molecular", "charge"]:
        raise ValueError(
            'Atom style "{}" is invalid or is not currently supported'.format(
                atom_style
            )
        )

    if unit_style not in [
        "real",
        "lj",
        "metal",
        "si",
        "cgs",
        "electron",
        "micro",
        "nano",
    ]:
        raise ValueError(
            'Unit style "{}" is invalid or is not currently supported'.format(
                unit_style
            )
        )
    if unit_style != "lj" and lj_cfactorsDict:
        raise ValueError("lj_cfactorsDict argument is only used if unit_style is lj.")
    base_unyts = LAMMPS_UnitSystems(unit_style)
    default_parameterMaps = {  # TODO: sites are not checked currently because gmso
        # doesn't store pair potential eqn the same way as the connections.
        "impropers": [
            "HarmonicImproperPotential",
            "HarmonicTorsionPotential",
            "PeriodicTorsionPotential",
        ],
        "dihedrals": ["OPLSTorsionPotential", "PeriodicTorsionPotential"],
        "angles": ["LAMMPSHarmonicAnglePotential"],
        "bonds": ["LAMMPSHarmonicBondPotential"],
        # "sites":"LennardJonesPotential",
        # "sites":"CoulombicPotential"
    }

    # TODO: Use strict_x, (e.g. x=bonds) to validate what topology attrs to convert
    if not strict_potentials:
        _try_default_potential_conversions(top, default_parameterMaps)
    potentialsMap = _validate_potential_compatibility(top)
    potential_typesDict = {}
    for potential in potentialsMap:
        pot_container = potential.__class__.__name__
        potStr = pot_container.lower() + "s"
        potStr = potStr[:-5] + "_" + potStr[-5:]
        if not potential_typesDict.get(potStr):
            potential_typesDict[potStr] = {potentialsMap[potential]}
        else:
            potential_typesDict[potStr].add(potentialsMap[potential])

    dihedral_parser = _identify_dihedral_parser(top, potential_typesDict)
    improper_parser = _identify_improper_parser(top, potential_typesDict)

    if strict_units:
        _validate_unit_compatibility(top, base_unyts)
    else:
        if base_unyts and unit_style != "lj":
            lj_cfactorsDict = None
        else:  # LJ unit styles
            if lj_cfactorsDict is None:
                lj_cfactorsDict = {}
            source_factorsList = list(lj_cfactorsDict.keys())
            defaultsList = ["length", "energy", "mass", "charge"]
            for source_factor in defaultsList + source_factorsList:
                if source_factor not in defaultsList:
                    raise ValueError(
                        f"Conversion factor {source_factor} is not used. Pleas only provide some of {defaultsList}"
                    )
                if lj_cfactorsDict.get(source_factor):
                    continue
                default_val_from_topology = _default_lj_val(top, source_factor)
                lj_cfactorsDict[source_factor] = lj_cfactorsDict.get(
                    source_factor, default_val_from_topology
                )

    reindex_molecules(top)  # reset the topology molecule index to match with lammps
    path = Path(filename)
    if not path.parent.exists():
        msg = "Provided path to file that does not exist"
        raise FileNotFoundError(msg)

    with open(path, "w") as out_file:
        _write_header(out_file, top, atom_style, dihedral_parser)
        _write_box(out_file, top, base_unyts, lj_cfactorsDict)
        all_ordered_typesDict = {}
        if top.is_fully_typed():
            _write_atomtypes(out_file, top, base_unyts, lj_cfactorsDict)
            _write_pairtypes(out_file, top, base_unyts, lj_cfactorsDict)
            if top.bond_types:
                sorted_bondsList = _write_bondtypes(
                    out_file, top, base_unyts, lj_cfactorsDict
                )
                all_ordered_typesDict["bonds"] = sorted_bondsList
            if top.angle_types:
                sorted_anglesList = _write_angletypes(
                    out_file, top, base_unyts, lj_cfactorsDict
                )
                all_ordered_typesDict["angles"] = sorted_anglesList
            if top.dihedral_types:
                sorted_dihedralsList = (
                    _write_dihedraltypes(  # return a list of dihedraltypes
                        out_file,
                        top,
                        base_unyts,
                        dihedral_parser,
                        lj_cfactorsDict,
                    )
                )
                all_ordered_typesDict["dihedrals"] = sorted_dihedralsList
            if top.improper_types:
                sorted_impropersList = _write_impropertypes(
                    out_file, top, base_unyts, improper_parser, lj_cfactorsDict
                )
                all_ordered_typesDict["impropers"] = sorted_impropersList

        _write_site_data(out_file, top, atom_style, base_unyts, lj_cfactorsDict)
        for conn in ["bonds", "angles", "dihedrals", "impropers"]:
            connIter = getattr(top, conn)
            conn_typesList = all_ordered_typesDict.get(conn)
            if connIter and conn_typesList:
                _write_conn_data(out_file, top, conn, conn_typesList)


@loads_as(".lammps", ".lammpsdata", ".data")
def read_lammpsdata(
    filename,
    atom_style="full",
    unit_style="real",
):
    """Read in a lammps data file as a GMSO topology.

    Parameters
    ----------
    filename : str
        LAMMPS data file
    atom_style : str, optional, default='full'
        Inferred atom style defined by LAMMPS, be certain that this is provided
        accurately.
    unit_style : str, optional, default='real
        LAMMPS unit style used for writing the datafile. Can be "real", "lj",
        "metal", "si", "cgs", "electron", "micro", "nano".

    Returns
    -------
    top : GMSO Topology
        A GMSO Topology object

    Notes
    -----
    See http://lammps.sandia.gov/doc/2001/data_format.html for a full description of the LAMMPS data format.

    This is a work in progress, as only several atom styles, potential styles, and unit styes are currently supported.

    Currently supporting the following atom styles: 'full'

    Currently supporting the following unit styles: 'real', "real", "lj", "metal", "si", "cgs",
    "electron", "micro", "nano".

    Currently supporting the following potential styles: 'lj'
    Currently supporting the following bond styles: 'harmonic'
    Currently supporting the following angle styles: 'harmonic'
    Currently supporting the following dihedral styles: 'opls'
    Currently supporting the following improper styles: 'harmonic'

    """
    top = Topology()

    # Validate 'atom_style'
    if atom_style not in ["full"]:
        raise ValueError(
            'Atom Style "{}" is invalid or is not currently supported'.format(
                atom_style
            )
        )

    # Validate 'unit_style'
    if unit_style not in [
        "real",
        "lj",
        "metal",
        "si",
        "cgs",
        "electron",
        "micro",
        "nano",
    ]:
        raise ValueError(
            'Unit Style "{}" is invalid or is not currently supported'.format(
                unit_style
            )
        )
    base_unyts = LAMMPS_UnitSystems(unit_style)

    # Parse box information
    _get_box_coordinates(filename, base_unyts, top)
    # Parse atom type information
    top, type_list = _get_ff_information(filename, base_unyts, top)
    # Parse atom information
    _get_atoms(filename, top, base_unyts, type_list)
    # Parse connection (bonds, angles, dihedrals, impropers) information
    # TODO: Add more atom styles
    if atom_style in ["full"]:
        _get_connection(filename, top, base_unyts, connection_type="bond")
        _get_connection(filename, top, base_unyts, connection_type="angle")
        _get_connection(filename, top, base_unyts, connection_type="dihedral")
        _get_connection(filename, top, base_unyts, connection_type="improper")

    top.update_topology()

    return top


def get_units(base_unyts, dimension):
    """Get u.Unit for specific LAMMPS unit style with given dimension."""
    # Need separate angle units for harmonic force constant and angle
    if base_unyts.usystem.name == "lj":
        if dimension == "angle":
            return u.radian
        return u.dimensionless

    if dimension == "angle_eq":
        return u.degree  # LAMMPS specifies different units for some angles, such as equilibrium angles

    return u.Unit(base_unyts.usystem[dimension], registry=base_unyts.reg)


def _get_connection(filename, topology, base_unyts, connection_type):
    """Parse connection types."""
    # TODO: check for other connection types besides the defaults
    with open(filename, "r") as lammps_file:
        types = False
        for i, line in enumerate(lammps_file):
            if connection_type in line.split():
                n_connection_types = int(line.split()[0])
                types = True
            if connection_type.capitalize() in line.split():
                break
    if types is False:
        return topology
    templates = PotentialTemplateLibrary()
    connection_type_lines = open(filename, "r").readlines()[
        i + 2 : i + n_connection_types + 2
    ]
    connection_type_list = list()
    for line in connection_type_lines:
        if connection_type == "bond":
            template_potential = templates["LAMMPSHarmonicBondPotential"]
            # Multiply 'k' by 2 since LAMMPS includes 1/2 in the term
            conn_params = {
                "k": float(line.split()[1])
                * get_units(base_unyts, "energy")
                / get_units(base_unyts, "length") ** 2
                * 2,
                "r_eq": float(line.split()[2]) * get_units(base_unyts, "length"),
            }
            name = template_potential.name
            expression = template_potential.expression
            variables = template_potential.independent_variables
            c_type = getattr(gmso, "BondType")(
                name=name,
                parameters=conn_params,
                expression=expression,
                independent_variables=variables,
            )
        elif connection_type == "angle":
            template_potential = templates["LAMMPSHarmonicAnglePotential"]
            # Multiply 'k' by 2 since LAMMPS includes 1/2 in the term
            conn_params = {
                "k": float(line.split()[1])
                * get_units(base_unyts, "energy")
                / get_units(base_unyts, "angle") ** 2
                * 2,
                "theta_eq": float(line.split()[2]) * get_units(base_unyts, "angle_eq"),
            }
            name = template_potential.name
            expression = template_potential.expression
            variables = template_potential.independent_variables
            c_type = getattr(gmso, "AngleType")(
                name=name,
                parameters=conn_params,
                expression=expression,
                independent_variables=variables,
            )
        elif connection_type == "dihedral":
            template_potential = templates["OPLSTorsionPotential"]
            conn_params = {
                "k1": float(line.split()[1]) * get_units(base_unyts, "energy"),
                "k2": float(line.split()[2]) * get_units(base_unyts, "energy"),
                "k3": float(line.split()[3]) * get_units(base_unyts, "energy"),
                "k4": float(line.split()[4]) * get_units(base_unyts, "energy"),
            }
            name = template_potential.name
            expression = template_potential.expression
            variables = template_potential.independent_variables
            c_type = getattr(gmso, "DihedralType")(
                name=name,
                parameters=conn_params,
                expression=expression,
                independent_variables=variables,
            )
        elif connection_type == "improper":
            template_potential = templates["HarmonicImproperPotential"]
            conn_params = {
                "k": float(line.split()[2])
                * get_units(base_unyts, "energy")
                / get_units(base_unyts, "energy") ** 2
                * 2,
                "phi_eq": float(line.split()[3]) * get_units(base_unyts, "angle_eq"),
            }
            name = template_potential.name
            expression = template_potential.expression
            variables = template_potential.independent_variables
            c_type = getattr(gmso, "ImproperType")(
                name=name,
                parameters=conn_params,
                expression=expression,
                independent_variables=variables,
            )

        connection_type_list.append(c_type)

    with open(filename, "r") as lammps_file:
        for i, line in enumerate(lammps_file):
            if connection_type + "s" in line.split():
                n_connections = int(line.split()[0])
            if connection_type.capitalize() + "s" in line.split():
                break
    connection_lines = open(filename, "r").readlines()[i + 2 : i + n_connections + 2]
    # Determine number of sites to generate
    if connection_type == "bond":
        n_sites = 2
    elif connection_type == "angle":
        n_sites = 3
    else:
        n_sites = 4
    for i, line in enumerate(connection_lines):
        site_list = list()
        for j in range(n_sites):
            site = topology.sites[int(line.split()[j + 2]) - 1]
            site_list.append(site)
        ctype = copy.copy(connection_type_list[int(line.split()[1]) - 1])
        ctype.member_types = tuple(map(lambda x: x.atom_type.name, site_list))
        ctype.member_classes = ctype.member_types
        if connection_type == "bond":
            connection = Bond(
                connection_members=site_list,
                bond_type=ctype,
            )
        elif connection_type == "angle":
            connection = Angle(
                connection_members=site_list,
                angle_type=ctype,
            )
        elif connection_type == "dihedral":
            connection = Dihedral(
                connection_members=site_list,
                dihedral_type=ctype,
            )
        elif connection_type == "improper":
            connection = Improper(
                connection_members=site_list,
                improper_type=ctype,
            )
        topology.add_connection(connection)

    return topology


def _get_atoms(filename, topology, base_unyts, type_list):
    """Parse the atom information in the LAMMPS data file."""
    with open(filename, "r") as lammps_file:
        for i, line in enumerate(lammps_file):
            if "atoms" in line.split():
                n_atoms = int(line.split()[0])
            if "Atoms" in line.split():
                break
    atom_lines = open(filename, "r").readlines()[i + 2 : i + n_atoms + 2]
    for line in atom_lines:
        atom_line = line.split()
        atom_type = atom_line[2]
        charge = u.unyt_quantity(float(atom_line[3]), get_units(base_unyts, "charge"))
        coord = u.unyt_array(
            [float(atom_line[4]), float(atom_line[5]), float(atom_line[6])]
        ) * get_units(base_unyts, "length")
        site = Atom(
            charge=charge,
            position=coord,
            atom_type=copy.deepcopy(type_list[int(atom_type) - 1]),  # 0-index
            molecule=Molecule(
                name=atom_line[1], number=int(atom_line[1]) - 1
            ),  # 0-index
        )
        element = element_by_mass(site.atom_type.mass.value)
        site.name = element.name if element else site.atom_type.name
        site.element = element
        topology.add_site(site)

    return topology


def _get_box_coordinates(filename, base_unyts, topology):
    """Parse box information."""
    with open(filename, "r") as lammps_file:
        for line in lammps_file:
            if "xlo" in line.split():
                break
        x_line = line.split()
        y_line = lammps_file.readline().split()
        z_line = lammps_file.readline().split()

        x = float(x_line[1]) - float(x_line[0])
        y = float(y_line[1]) - float(y_line[0])
        z = float(z_line[1]) - float(z_line[0])

        # Check if box is triclinic
        tilts = lammps_file.readline().split()
        if "xy" in tilts:
            xy = float(tilts[0])
            xz = float(tilts[1])
            yz = float(tilts[2])

            xhi = float(x_line[1]) - np.max([0.0, xy, xz, xy + xz])
            xlo = float(x_line[0]) - np.min([0.0, xy, xz, xy + xz])
            yhi = float(y_line[1]) - np.max([0.0, yz])
            ylo = float(y_line[0]) - np.min([0.0, yz])
            zhi = float(z_line[1])
            zlo = float(z_line[0])

            lx = xhi - xlo
            ly = yhi - ylo
            lz = zhi - zlo

            c = np.sqrt(lz**2 + xz**2 + yz**2)
            b = np.sqrt(ly**2 + xy**2)
            a = lx

            alpha = np.arccos((yz * ly + xy * xz) / (b * c))
            beta = np.arccos(xz / c)
            gamma = np.arccos(xy / b)

            # Box Information
            lengths = u.unyt_array([a, b, c], get_units(base_unyts, "length"))
            angles = u.unyt_array([alpha, beta, gamma], get_units(base_unyts, "angle"))
            topology.box = Box(lengths, angles)
        else:
            # Box Information
            lengths = u.unyt_array([x, y, z], get_units(base_unyts, "length"))
            topology.box = Box(lengths)

        return topology


def _get_ff_information(filename, base_unyts, topology):
    """Parse atom-type information."""
    with open(filename, "r") as lammps_file:
        types = False
        for i, line in enumerate(lammps_file):
            if "atom" in line:
                n_atomtypes = int(line.split()[0])
                types = True
            elif "Masses" in line:
                break
    if types is False:
        return topology
    mass_lines = open(filename, "r").readlines()[i + 2 : i + n_atomtypes + 2]
    type_list = list()
    for line in mass_lines:
        atom_type = AtomType(
            name=line.split()[0],
            mass=float(line.split()[1]) * get_units(base_unyts, "mass"),
        )
        type_list.append(atom_type)

    with open(filename, "r") as lammps_file:
        for i, line in enumerate(lammps_file):
            if "Pair" in line:
                break
    # Need to figure out if we're going have mixing rules printed out
    # Currently only reading in LJ params
    warn_ljcutBool = False
    pair_lines = open(filename, "r").readlines()[i + 2 : i + n_atomtypes + 2]
    for i, pair in enumerate(pair_lines):
        if len(pair.split()) == 3:
            type_list[i].parameters["sigma"] = float(pair.split()[2]) * get_units(
                base_unyts, "length"
            )
            type_list[i].parameters["epsilon"] = float(pair.split()[1]) * get_units(
                base_unyts, "energy"
            )
        elif len(pair.split()) == 4:
            warn_ljcutBool = True

    if warn_ljcutBool:
        warnings.warn(
            "Currently not reading in LJ cutoff values."
            "These should be specified in the engine run files."
        )

    return topology, type_list


def _accepted_potentials():
    """List of accepted potentials that LAMMPS can support."""
    templates = PotentialTemplateLibrary()
    lennard_jones_potential = templates["LennardJonesPotential"]
    harmonic_bond_potential = templates["LAMMPSHarmonicBondPotential"]
    harmonic_angle_potential = templates["LAMMPSHarmonicAnglePotential"]
    periodic_torsion_potential = templates["PeriodicTorsionPotential"]
    harmonic_improper_potential = templates["HarmonicImproperPotential"]
    opls_torsion_potential = templates["OPLSTorsionPotential"]
    accepted_potentialsList = [
        lennard_jones_potential,
        harmonic_bond_potential,
        harmonic_angle_potential,
        periodic_torsion_potential,
        harmonic_improper_potential,
        opls_torsion_potential,
    ]
    return accepted_potentialsList


def _validate_potential_compatibility(top):
    """Check compatability of topology object potentials with LAMMPSDATA format."""
    pfilter = PotentialFilters.UNIQUE_EXPRESSION
    pot_types = check_compatibility(
        top, _accepted_potentials(), site_pfilter=pfilter, conn_pfilter=pfilter
    )
    return pot_types


def _validate_unit_compatibility(top, base_unyts):
    """Check compatability of topology object units with LAMMPSDATA format."""
    for attribute in ["sites", "bonds", "angles", "dihedrals", "impropers"]:
        if attribute == "sites":
            atype = "atom_types"
        else:
            atype = attribute[:-1] + "_types"
        parametersList = [
            (parameter, name)
            for attr_type in getattr(top, atype)
            for name, parameter in attr_type.parameters.items()
        ]
        for parameter, name in parametersList:
            assert np.isclose(
                float(base_unyts.convert_parameter(parameter, n_decimals=6, name=name)),
                parameter.value,
                atol=1e-3,
            ), (
                f"Units System {base_unyts.usystem} is not compatible with {atype} with value {parameter}"
            )


def _write_header(out_file, top, atom_style, dihedral_parser):
    """Write Lammps file header."""
    out_file.write(
        "{} written by {} at {} using the GMSO LAMMPS Writer\n\n\n".format(
            top.name if top.name is not None else "Topology",
            os.environ.get("USER"),
            str(datetime.datetime.now()),
        )
    )
    out_file.write("{:d} atoms\n".format(top.n_sites))
    if atom_style in ["full", "molecular"]:
        out_file.write("{:d} bonds\n".format(top.n_bonds))
        out_file.write("{:d} angles\n".format(top.n_angles))
        if dihedral_parser in [parse_opls_style_dihedral]:  # no layered dihedrals
            n_dihedrals = top.n_dihedrals
        elif dihedral_parser in [parse_charmm_style_dihedral]:  # layered dihedrals
            n_dihedrals = 0
            for dihedral in top.dihedrals:
                param = next(iter(dihedral.dihedral_type.parameters.values()))
                if isinstance(param, u.unyt_quantity):
                    n_dihedrals += 1
                else:
                    n_dihedrals += len(param)
        elif dihedral_parser is None:
            n_dihedrals = 0
        out_file.write("{:d} dihedrals\n".format(n_dihedrals))
        out_file.write("{:d} impropers\n\n".format(top.n_impropers))

    # TODO: allow users to specify filter_by syntax
    out_file.write("{:d} atom types\n".format(len(top.atom_types(filter_by=pfilter))))
    if top.n_bonds > 0 and atom_style in ["full", "molecular"]:
        out_file.write(
            "{:d} bond types\n".format(len(top.bond_types(filter_by=pfilter)))
        )
    if top.n_angles > 0 and atom_style in ["full", "molecular"]:
        out_file.write(
            "{:d} angle types\n".format(len(top.angle_types(filter_by=pfilter)))
        )
    if top.n_dihedrals > 0 and atom_style in ["full", "molecular"]:
        unique_dtypes = top.dihedral_types(filter_by=pfilter)
        nkeys = len(next(iter(unique_dtypes)).parameters.keys())
        nparams = 0  # write out the total number of found for dihedrals
        for potential in unique_dtypes:
            for param in potential.parameters.values():
                paramList = param.tolist()
                if isinstance(paramList, float):
                    nparams += 1
                else:
                    for _ in param.tolist():
                        nparams += 1
        ntypes = int(
            nparams / nkeys
        )  # allows us to count multiples for ones stored in a single object
        out_file.write("{:d} dihedral types\n".format(ntypes))
    if top.n_impropers > 0 and atom_style in ["full", "molecular"]:
        out_file.write(
            "{:d} improper types\n".format(len(top.improper_types(filter_by=pfilter)))
        )

    out_file.write("\n")


def _write_box(out_file, top, base_unyts, cfactorsDict):
    """Write GMSO Topology box to LAMMPS file."""
    if allclose_units(
        top.box.angles,
        u.unyt_array([90, 90, 90], "degree"),
        rtol=1e-5,
        atol=1e-8,
    ):
        box_lengths = [
            float(base_unyts.convert_parameter(top.box.lengths[i], cfactorsDict))
            for i in range(3)
        ]
        for i, dim in enumerate(["x", "y", "z"]):
            out_file.write(
                "{0:.6f} {1:.6f} {2}lo {2}hi\n".format(0, box_lengths[i], dim)
            )
        out_file.write("0.000000 0.000000 0.000000 xy xz yz\n")
    else:
        box_lengths = [
            float(base_unyts.convert_parameter(top.box.lengths[i], cfactorsDict))
            for i in range(3)
        ]
        vectors = (box_lengths * top.box.get_unit_vectors().T).T

        xhi = vectors[0][0]
        yhi = vectors[1][1]
        zhi = vectors[2][2]
        xy = vectors[1][0]
        xz = vectors[2][0]
        yz = vectors[2][1]
        xlo = u.unyt_array(0, xy.units)
        ylo = u.unyt_array(0, xy.units)
        zlo = u.unyt_array(0, xy.units)

        xlo_bound = xlo + u.unyt_array(np.min([0.0, xy, xz, xy + xz]), xy.units)
        xhi_bound = xhi + u.unyt_array(np.max([0.0, xy, xz, xy + xz]), xy.units)
        ylo_bound = ylo + u.unyt_array(np.min([0.0, yz]), xy.units)
        yhi_bound = yhi + u.unyt_array(np.max([0.0, yz]), xy.units)
        zlo_bound = zlo
        zhi_bound = zhi

        out_file.write(
            "{0:.6f} {1:.6f} xlo xhi\n".format(xlo_bound.value, xhi_bound.value)
        )
        out_file.write(
            "{0:.6f} {1:.6f} ylo yhi\n".format(ylo_bound.value, yhi_bound.value)
        )
        out_file.write(
            "{0:.6f} {1:.6f} zlo zhi\n".format(zlo_bound.value, zhi_bound.value)
        )
        out_file.write(
            "{0:.6f} {1:.6f} {2:.6f} xy xz yz\n".format(xy.value, xz.value, yz.value)
        )


def _write_atomtypes(out_file, top, base_unyts, cfactorsDict):
    """Write out atomtypes in GMSO topology to LAMMPS file."""
    out_file.write("\nMasses\n")
    out_file.write(f"#\tmass ({base_unyts.usystem['mass']})\n")
    atypesView = sorted(top.atom_types(filter_by=pfilter), key=lambda x: x.name)
    for atom_type in atypesView:
        out_file.write(
            "{:d}\t{}\t# {}\n".format(
                atypesView.index(atom_type) + 1,
                base_unyts.convert_parameter(atom_type.mass, cfactorsDict),
                atom_type.name,
            )
        )


def _write_pairtypes(out_file, top, base_unyts, cfactorsDict):
    """Write out pair interaction to LAMMPS file."""
    # TODO: Handling of modified cross-interactions is not considered from top.pairpotential_types
    # Pair coefficients
    test_atomtype = top.sites[0].atom_type
    out_file.write(f"\nPair Coeffs # {test_atomtype.expression}\n")
    nb_style_orderTuple = (
        "epsilon",
        "sigma",
    )  # this will vary with new pair styles
    param_labels = [
        write_out_parameter_and_units(
            key,
            convert_kelvin_to_energy_units(test_atomtype.parameters[key], "kJ"),
            base_unyts,
        )
        for key in nb_style_orderTuple
    ]
    out_file.write("#\t" + "\t".join(param_labels) + "\n")
    sorted_atomtypes = sorted(top.atom_types(filter_by=pfilter), key=lambda x: x.name)
    for idx, param in enumerate(sorted_atomtypes):
        out_file.write(
            "{}\t{:7}\t\t{:7}\t\t# {}\n".format(
                idx + 1,
                *[
                    base_unyts.convert_parameter(
                        convert_kelvin_to_energy_units(param.parameters[key], "kJ"),
                        cfactorsDict,
                        n_decimals=5,
                    )
                    for key in nb_style_orderTuple
                ],
                param.name,
            )
        )


def _write_bondtypes(out_file, top, base_unyts, cfactorsDict):
    """Write out bonds to LAMMPS file."""
    # TODO: Use any accepted lammps styles (only takes harmonic now)
    test_bondtype = top.bonds[0].bond_type
    out_file.write(f"\nBond Coeffs #{test_bondtype.name}\n")
    bond_style_orderTuple = ("k", "r_eq")
    param_labels = [
        write_out_parameter_and_units(
            key,
            convert_kelvin_to_energy_units(test_bondtype.parameters[key], "kJ"),
            base_unyts,
        )
        for key in bond_style_orderTuple
    ]

    out_file.write("#\t" + "\t".join(param_labels) + "\n")
    bond_types = list(top.bond_types(filter_by=pfilter))
    bond_types.sort(key=lambda x: sorted(x.member_types))
    for idx, bond_type in enumerate(bond_types):
        member_types = sorted([bond_type.member_types[0], bond_type.member_types[1]])
        out_file.write(
            "{}\t{:7}\t{:7}\t\t# {}\t{}\n".format(
                idx + 1,
                *[
                    base_unyts.convert_parameter(
                        convert_kelvin_to_energy_units(bond_type.parameters[key], "kJ"),
                        cfactorsDict,
                        n_decimals=6,
                    )
                    for key in bond_style_orderTuple
                ],
                *member_types,
            )
        )
    return bond_types


def _write_angletypes(out_file, top, base_unyts, cfactorsDict):
    """Write out angles to LAMMPS file."""
    # TODO: Use any accepted lammps parameters, only harmonic now
    test_angletype = top.angles[0].angle_type
    out_file.write(f"\nAngle Coeffs #{test_angletype.name}\n")
    angle_style_orderTuple = (
        "k",
        "theta_eq",
    )  # this will vary with new angle styles
    param_labels = [
        write_out_parameter_and_units(
            key,
            convert_kelvin_to_energy_units(test_angletype.parameters[key], "kJ"),
            base_unyts,
        )
        for key in angle_style_orderTuple
    ]
    out_file.write("#\t" + "\t".join(param_labels) + "\n")
    angle_types = list(top.angle_types(filter_by=pfilter))
    angle_types.sort(
        key=lambda x: (
            x.member_types[1],
            min(x.member_types[::2]),
            max(x.member_types[::2]),
        )
    )
    for idx, angle_type in enumerate(angle_types):
        out_file.write(
            "{}\t{:7}\t{:7}\t#{:11s}\t{:11s}\t{:11s}\n".format(
                idx + 1,
                *[
                    base_unyts.convert_parameter(
                        convert_kelvin_to_energy_units(
                            angle_type.parameters[key], "kJ"
                        ),
                        cfactorsDict,
                        n_decimals=6,
                        name=key,
                    )
                    for key in angle_style_orderTuple
                ],
                *angle_type.member_types,
            )
        )
    return angle_types


def _write_dihedraltypes(out_file, top, base_unyts, parser, cfactorsDict):
    """Write out dihedrals to LAMMPS file."""
    test_dihedraltype = top.dihedrals[0].dihedral_type
    out_file.write(f"\nDihedral Coeffs #{test_dihedraltype.name}\n")
    param_labels0 = parser(test_dihedraltype)  # tuple (paramsList, params_namesList)

    if isinstance(
        param_labels0[0][0], list
    ):  # check for parsing out multiple instances from the dihedral
        param_labels = [
            write_out_parameter_and_units(
                name, convert_kelvin_to_energy_units(param, "kJ"), base_unyts
            )
            for param, name in zip(param_labels0[0][0], param_labels0[1])
        ]
    else:
        param_labels = [
            write_out_parameter_and_units(
                name, convert_kelvin_to_energy_units(param, "kJ"), base_unyts
            )
            for param, name in zip(param_labels0[0], param_labels0[1])
        ]
    out_file.write("#\t" + "\t".join(param_labels) + "\n")
    indexList = list(top.dihedral_types(filter_by=pfilter))
    index_membersList = [
        (dihedral_type, sort_by_types(dihedral_type)) for dihedral_type in indexList
    ]
    index_membersList.sort(key=lambda x: ([x[1][i] for i in [1, 2, 0, 3]]))
    # handle variable lengths for parameters
    base_msg = "{}\t"  # handles index
    end_msg = "# {}\t{}\t{}\t{}\n"

    if (
        parser.__name__ == "parse_opls_style_dihedral"
    ):  # one opls parameter per dihedral type
        dihedral_typesList = []
        for idx, (dihedral_type, members) in enumerate(index_membersList):
            param_labels = parser(dihedral_type)
            variable_msg = "{:8}\t" * len(param_labels[1])
            full_msg = base_msg + variable_msg + end_msg
            out_file.write(
                full_msg.format(
                    idx + 1,
                    *[
                        base_unyts.convert_parameter(
                            convert_kelvin_to_energy_units(parameter, "kJ"),
                            cfactorsDict,
                            n_decimals=6,
                            name=parameterStr,
                        )
                        for parameter, parameterStr in zip(*parser(dihedral_type))
                    ],
                    *members,
                )
            )
            dihedral_typesList.append(dihedral_type)

    elif parser.__name__ == "parse_charmm_style_dihedral":
        ndecimalsDict = {"k": 6, "n": 0, "phi_eq": 0, "weights": 1}
        idx = 0
        dihedral_typesList = []
        for dihedral_type, members in index_membersList:
            parameter_termList, parameterStrList = parser(dihedral_type)
            variable_msg = "{:8}\t" * len(parameterStrList)
            full_msg = base_msg + variable_msg + end_msg
            for parameter_terms in parameter_termList:  # list of params on each line
                out_file.write(
                    full_msg.format(
                        idx + 1,
                        *[
                            base_unyts.convert_parameter(
                                convert_kelvin_to_energy_units(parameter, "kJ"),
                                cfactorsDict,
                                n_decimals=ndecimalsDict[parameterStr],
                                name=parameterStr,
                            )
                            for parameter, parameterStr in zip(
                                parameter_terms, parameterStrList
                            )
                        ],
                        *members,
                    )
                )
                dihedral_typesList.append(
                    dihedral_type
                )  # add dihedral type multiple times if it is layered
                idx += 1
    return dihedral_typesList


def parse_opls_style_dihedral(dihedral_type):
    """Take a dihedral type and list parameters as expected in lammps outputs."""
    parametersList = []
    namesList = ["k1", "k2", "k3", "k4"]
    for k in namesList:
        parametersList.append(dihedral_type.parameters[k])

    return parametersList, namesList


def parse_charmm_style_dihedral(dihedral_type, weightsArray=None):
    """Take a dihedral type and list parameters as expected in lammps outputs."""
    kArray = dihedral_type.parameters["k"].flatten()
    nArray = dihedral_type.parameters["n"].flatten()
    phi_eqArray = dihedral_type.parameters["phi_eq"].flatten()
    if not weightsArray:  # used for amber forcefield weights are 0
        weightsArray = np.zeros(kArray.size) * u.dimensionless
    allParamsList = []
    for a, b, c, d in zip(kArray, nArray, phi_eqArray, weightsArray):
        allParamsList.append([a, b, c, d])
    return allParamsList, ["k", "n", "phi_eq", "weights"]


def _write_impropertypes(out_file, top, base_unyts, parser, cfactorsDict):
    """Write out impropers to LAMMPS file."""
    test_impropertype = top.impropers[0].improper_type
    out_file.write(f"\nImproper Coeffs #{test_impropertype.name}\n")
    param_labels0 = parser(test_impropertype)  # tuple (paramsList, params_namesList)

    if isinstance(
        param_labels0[0][0], list
    ):  # check for parsing out multiple instances from the dihedral
        param_labels = [
            write_out_parameter_and_units(
                name, convert_kelvin_to_energy_units(param, "kJ"), base_unyts
            )
            for param, name in zip(param_labels0[0][0], param_labels0[1])
        ]
    else:
        param_labels = [
            write_out_parameter_and_units(
                name, convert_kelvin_to_energy_units(param, "kJ"), base_unyts
            )
            for param, name in zip(param_labels0[0], param_labels0[1])
        ]
    out_file.write("#\t" + "\t".join(param_labels) + "\n")
    indexList = list(top.improper_types(filter_by=pfilter))
    index_membersList = [
        (improper_type, sort_by_types(improper_type)) for improper_type in indexList
    ]
    index_membersList.sort(key=lambda x: ([x[1][i] for i in [0, 1, 2, 3]]))
    # handle variable lengths for parameters
    base_msg = "{}\t"  # handles index
    end_msg = "# {}\t{}\t{}\t{}\n"

    if (
        parser.__name__ == "parse_cvff_style_improper"
        or "parse_harmonic_style_improper"
    ):  # one cvff parameter per dihedral type
        ndecimalsDict = {"k": 6, "n": 0, "phi_eq": 0}
        idx = 0
        for improper_type, members in index_membersList:
            param_labels = parser(improper_type)
            variable_msg = "{:8}\t" * len(param_labels[1])
            full_msg = base_msg + variable_msg + end_msg
            out_file.write(
                full_msg.format(
                    idx + 1,
                    *[
                        base_unyts.convert_parameter(
                            convert_kelvin_to_energy_units(parameter, "kJ"),
                            cfactorsDict,
                            n_decimals=ndecimalsDict[parameterStr],
                            name=parameterStr,
                        )
                        for parameter, parameterStr in zip(*parser(improper_type))
                    ],
                    *members,
                )
            )
            idx += 1
    return index_membersList  # cvff is not layered, so no added to list


def parse_cvff_style_improper(improper_type):
    """Take a dihedral type and list parameters as expected in lammps outputs."""
    parametersList = []
    namesList = ["k", "n", "phi_eq"]
    for k in namesList:
        parametersList.append(improper_type.parameters[k])
    return parametersList, namesList


def parse_harmonic_style_improper(improper_type):
    """Take a dihedral type and list parameters as expected in lammps outputs."""
    parametersList = []
    namesList = ["k", "phi_eq"]
    for k in namesList:
        parametersList.append(improper_type.parameters[k])
    return parametersList, namesList


def _write_site_data(out_file, top, atom_style, base_unyts, cfactorsDict):
    """Write atomic positions and charges to LAMMPS file.."""
    out_file.write(f"\nAtoms #{atom_style}\n\n")
    if atom_style == "atomic":
        atom_line = "{index:d}\t{type_index:d}\t{x:.8}\t{y:.8}\t{z:.8}\n"
    elif atom_style == "charge":
        atom_line = "{index:d}\t{type_index:d}\t{charge:.8}\t{x:.8}\t{y:.8}\t{z:.8}\n"
    elif atom_style == "molecular":
        atom_line = (
            "{index:d}\t{moleculeid:d}\t{type_index:d}\t{x:.8}\t{y:.8}\t{z:.8}\n"
        )
    elif atom_style == "full":
        atom_line = "{index:d}\t{moleculeid:d}\t{type_index:d}\t{charge:.8}\t{x:.8}\t{y:.8}\t{z:.8}\n"

    unique_sorted_typesList = sorted(
        top.atom_types(filter_by=pfilter), key=lambda x: x.name
    )
    for i, site in enumerate(top.sites):
        out_file.write(
            atom_line.format(
                index=i + 1,
                moleculeid=site.molecule.number + 1,  # index is 0-based in GMSO
                type_index=unique_sorted_typesList.index(site.atom_type) + 1,
                charge=base_unyts.convert_parameter(
                    site.charge,
                    cfactorsDict,
                    n_decimals=6,
                ),
                x=base_unyts.convert_parameter(
                    site.position[0],
                    cfactorsDict,
                    n_decimals=6,
                ),
                y=base_unyts.convert_parameter(
                    site.position[1], cfactorsDict, n_decimals=6
                ),
                z=base_unyts.convert_parameter(
                    site.position[2], cfactorsDict, n_decimals=6
                ),
            )
        )


def _angle_order_sorter(angle_typesList):
    return [angle_typesList[i] for i in [1, 0, 2]]


def _dihedral_order_sorter(dihedral_typesList):
    return [dihedral_typesList[i] for i in [1, 2, 0, 3]]


def _improper_order_sorter(improper_typesList):
    return [improper_typesList[i] for i in [0, 1, 2, 3]]


sorting_funcDict = {
    "bonds": None,
    "angles": _angle_order_sorter,
    "dihedrals": _dihedral_order_sorter,
    "impropers": _improper_order_sorter,
}


def _write_conn_data(out_file, top, connStr, sorted_typesList):
    """Write all connections to LAMMPS datafile."""
    out_file.write(f"\n{connStr.capitalize()}\n\n")

    i = 0
    for conn in getattr(top, connStr):
        ctype_members = sort_by_types(getattr(conn, connStr[:-1] + "_type"))
        indexList = [
            ind
            for ind, ele in zip(count(), sorted_typesList)
            if sort_by_types(ele) == ctype_members
        ]
        for index in indexList:
            typeStr = f"{i + 1:<6d}\t{index + 1:<6d}\t"
            sorted_membersList = sort_connection_members(conn, sort_by="index", top=top)
            indexStr = "\t".join(
                [
                    str(top.get_index(member) + 1).ljust(6)
                    for member in sorted_membersList
                ]
            )
            out_file.write(typeStr + indexStr + "\n")
            i += 1


def _try_default_potential_conversions(top, potentialsDict):
    """Take a topology a convert all potentials to the style in potentialDict."""
    for pot_container in potentialsDict:
        containerStr = pot_container[:-1] + "_types"
        if getattr(top, containerStr):
            for potential in potentialsDict[pot_container]:
                top.convert_potential_styles({pot_container: potential})
        elif getattr(top, pot_container):
            raise AttributeError(
                f"Missing parameters in {pot_container} for {top.get_untyped(pot_container)}"
            )


def _default_lj_val(top, source):
    """Generate default lj non-dimensional values from topology."""
    if source == "length":
        return copy.deepcopy(
            max(list(map(lambda x: x.parameters["sigma"], top.atom_types)))
        )
    elif source == "energy":
        return copy.deepcopy(
            max(list(map(lambda x: x.parameters["epsilon"], top.atom_types)))
        )
    elif source == "mass":
        return copy.deepcopy(max(list(map(lambda x: x.mass, top.atom_types))))
    elif source == "charge":
        return copy.deepcopy(max(list(map(lambda x: x.charge, top.atom_types))))
    else:
        raise ValueError(
            f"Provided {source} for default LJ cannot be found in the topology."
        )


def _identify_dihedral_parser(top, potential_typesDict):
    if not getattr(top, "dihedral_types"):
        return None
    # This is where dihedral_parser should get found
    parserDict = {
        "PeriodicTorsionPotential": parse_charmm_style_dihedral,
        "OPLSTorsionPotential": parse_opls_style_dihedral,
    }
    assert (
        len(potential_typesDict["dihedral_types"]) == 1
    )  # only allowing one potential type atm
    dihedralparser = parserDict[potential_typesDict["dihedral_types"].pop()]
    return dihedralparser


def _identify_improper_parser(top, potential_typesDict):
    if not getattr(top, "improper_types"):
        return None
    # This is where improper_parser should be stored
    parserDict = {
        "PeriodicTorsionPotential": parse_cvff_style_improper,
        "HarmonicTorsionPotential": parse_harmonic_style_improper,
        "HarmonicImproperPotential": parse_harmonic_style_improper,
    }
    assert (
        len(potential_typesDict["improper_types"]) == 1
    )  # only allowing one potential type atm
    improper_parser = parserDict[potential_typesDict["improper_types"].pop()]
    return improper_parser