convert_parmed.py

"""Module support for converting to/from ParmEd objects."""

import copy
import warnings
from operator import attrgetter, itemgetter

import numpy as np
import unyt as u
from symengine import expand

import gmso
from gmso.core.element import element_by_atomic_number, element_by_symbol
from gmso.core.views import PotentialFilters, get_parameters

pfilter = PotentialFilters.UNIQUE_PARAMETERS
from gmso.exceptions import GMSOError
from gmso.lib.potential_templates import PotentialTemplateLibrary
from gmso.utils.io import has_parmed, import_

if has_parmed:
    pmd = import_("parmed")

lib = PotentialTemplateLibrary()


def from_parmed(structure, refer_type=True):
    """Convert a parmed.Structure to a gmso.Topology.

    Convert a parametrized or un-parametrized parmed.Structure object to a topology.Topology.
    Specifically, this method maps Structure to Topology and Atom to Site.
    This method can only convert AtomType, BondType AngleType, DihedralType, and
    ImproperType.

    Parameters
    ----------
    structure : parmed.Structure
        parmed.Structure instance that need to be converted.
    refer_type : bool, optional, default=True
        Whether or not to transfer AtomType, BondType, AngleType,
        DihedralType, and ImproperType information

    Returns
    -------
    top : gmso.Topology
    """
    msg = "Provided argument is not a Parmed Structure"
    assert isinstance(structure, pmd.Structure), msg

    top = gmso.Topology(name=structure.title)
    site_map = dict()

    if np.all(structure.box):
        # add gmso box from structure
        top.box = gmso.Box(
            (structure.box[0:3] * u.angstrom).in_units(u.nm),
            angles=u.degree * structure.box[3:6],
        )
    top.combining_rule = structure.combining_rule

    # Consolidate parmed atomtypes and relate topology atomtypes
    if refer_type:
        pmd_top_atomtypes = _atom_types_from_pmd(structure)

    ind_res = _check_independent_residues(structure)
    for residue in structure.residues:
        for atom in residue.atoms:
            # add atom to sites in gmso
            element = (
                element_by_atomic_number(atom.element) if atom.element else None
            )
            if residue.number == -1:  # use default value of 0 in GMSO
                residue_number = 0
            else:
                residue_number = residue.number - 1
            site = gmso.Atom(
                name=atom.name,
                charge=atom.charge * u.elementary_charge,
                position=[atom.xx, atom.xy, atom.xz] * u.angstrom,
                atom_type=None,
                residue=(residue.name, residue_number),
                element=element,
            )
            site.molecule = (residue.name, residue_number) if ind_res else None
            site.atom_type = (
                copy.deepcopy(pmd_top_atomtypes[atom.atom_type])
                if refer_type and isinstance(atom.atom_type, pmd.AtomType)
                else None
            )
            if site.atom_type:
                site.atom_type.charge = atom.charge * u.elementary_charge
            site_map[atom] = site
            top.add_site(site)

    harmonicbond_potential = lib["HarmonicBondPotential"]
    name = harmonicbond_potential.name
    expression = harmonicbond_potential.expression
    variables = harmonicbond_potential.independent_variables
    for bond in structure.bonds:
        # Generate bonds and harmonic parameters
        # If typed, assumed to be harmonic bonds
        top_connection = gmso.Bond(
            connection_members=_sort_bond_members(
                top, site_map, *attrgetter("atom1", "atom2")(bond)
            )
        )
        if refer_type and isinstance(bond.type, pmd.BondType):
            conn_params = {
                "k": (2 * bond.type.k * u.Unit("kcal / (angstrom**2 * mol)")),
                "r_eq": bond.type.req * u.angstrom,
            }
            _add_conn_type_from_pmd(
                connStr="BondType",
                pmd_conn=bond,
                gmso_conn=top_connection,
                conn_params=conn_params,
                name=name,
                expression=expression,
                variables=variables,
            )
        top.add_connection(top_connection, update_types=False)

    harmonicangle_potential = lib["HarmonicAnglePotential"]
    name = harmonicangle_potential.name
    expression = harmonicangle_potential.expression
    variables = harmonicangle_potential.independent_variables
    for angle in structure.angles:
        # Generate angles and harmonic parameters
        # If typed, assumed to be harmonic angles
        top_connection = gmso.Angle(
            connection_members=_sort_angle_members(
                top, site_map, *attrgetter("atom1", "atom2", "atom3")(angle)
            )
        )
        if refer_type and isinstance(angle.type, pmd.AngleType):
            conn_params = {
                "k": (2 * angle.type.k * u.Unit("kcal / (radian**2 * mol)")),
                "theta_eq": (angle.type.theteq * u.degree),
            }
            _add_conn_type_from_pmd(
                connStr="AngleType",
                pmd_conn=angle,
                gmso_conn=top_connection,
                conn_params=conn_params,
                name=name,
                expression=expression,
                variables=variables,
            )
        top.add_connection(top_connection, update_types=False)

    periodic_torsion_potential = lib["PeriodicTorsionPotential"]
    name_proper = periodic_torsion_potential.name
    expression_proper = periodic_torsion_potential.expression
    variables_proper = periodic_torsion_potential.independent_variables
    periodic_imp_potential = lib["PeriodicImproperPotential"]
    name_improper = periodic_imp_potential.name
    expression_improper = periodic_imp_potential.expression
    variables_improper = periodic_imp_potential.independent_variables
    for dihedral in structure.dihedrals:
        # Generate dihedrals and impropers from structure.dihedrals
        # If typed, assumed to be periodic
        if dihedral.improper:
            top_connection = gmso.Improper(
                connection_members=_sort_improper_members(
                    top,
                    site_map,
                    *attrgetter("atom1", "atom2", "atom3", "atom4")(dihedral),
                )
            )
            if refer_type and isinstance(dihedral.type, pmd.DihedralType):
                conn_params = {
                    "k": (dihedral.type.phi_k * u.Unit("kcal / mol")),
                    "phi_eq": (dihedral.type.phase * u.degree),
                    "n": dihedral.type.per * u.dimensionless,
                }
                _add_conn_type_from_pmd(
                    connStr="ImproperType",
                    pmd_conn=dihedral,
                    gmso_conn=top_connection,
                    conn_params=conn_params,
                    name=name_improper,
                    expression=expression_improper,
                    variables=variables_improper,
                )
        else:
            top_connection = gmso.Dihedral(
                connection_members=_sort_dihedral_members(
                    top,
                    site_map,
                    *attrgetter("atom1", "atom2", "atom3", "atom4")(dihedral),
                )
            )
            if refer_type and isinstance(dihedral.type, pmd.DihedralType):
                conn_params = {
                    "k": (dihedral.type.phi_k * u.Unit("kcal / mol")),
                    "phi_eq": (dihedral.type.phase * u.degree),
                    "n": dihedral.type.per * u.dimensionless,
                }
                _add_conn_type_from_pmd(
                    connStr="DihedralType",
                    pmd_conn=dihedral,
                    gmso_conn=top_connection,
                    conn_params=conn_params,
                    name=name_proper,
                    expression=expression_proper,
                    variables=variables_proper,
                )
        top.add_connection(top_connection, update_types=False)

    ryckaert_bellemans_torsion_potential = lib[
        "RyckaertBellemansTorsionPotential"
    ]
    name = ryckaert_bellemans_torsion_potential.name
    expression = ryckaert_bellemans_torsion_potential.expression
    variables = ryckaert_bellemans_torsion_potential.independent_variables
    for rb_torsion in structure.rb_torsions:
        # Generate dihedrals from structure rb_torsions
        # If typed, assumed to be ryckaert bellemans torsions
        top_connection = gmso.Dihedral(
            connection_members=_sort_dihedral_members(
                top,
                site_map,
                *attrgetter("atom1", "atom2", "atom3", "atom4")(rb_torsion),
            )
        )
        if refer_type and isinstance(rb_torsion.type, pmd.RBTorsionType):
            conn_params = {
                "c0": (rb_torsion.type.c0 * u.Unit("kcal/mol")),
                "c1": (rb_torsion.type.c1 * u.Unit("kcal/mol")),
                "c2": (rb_torsion.type.c2 * u.Unit("kcal/mol")),
                "c3": (rb_torsion.type.c3 * u.Unit("kcal/mol")),
                "c4": (rb_torsion.type.c4 * u.Unit("kcal/mol")),
                "c5": (rb_torsion.type.c5 * u.Unit("kcal/mol")),
            }
            _add_conn_type_from_pmd(
                connStr="DihedralType",
                pmd_conn=rb_torsion,
                gmso_conn=top_connection,
                conn_params=conn_params,
                name=name,
                expression=expression,
                variables=variables,
            )
        top.add_connection(top_connection, update_types=False)

    periodic_torsion_potential = lib["HarmonicTorsionPotential"]
    name = periodic_torsion_potential.name
    expression = periodic_torsion_potential.expression
    variables = periodic_torsion_potential.independent_variables
    for improper in structure.impropers:
        # Generate impropers from structure impropers
        # If typed, assumed to be harmonic torsions
        top_connection = gmso.Improper(
            connection_members=_sort_improper_members(
                top,
                site_map,
                *attrgetter("atom3", "atom2", "atom1", "atom4")(improper),
            )
        )
        if refer_type and isinstance(improper.type, pmd.ImproperType):
            conn_params = {
                "k": (improper.type.psi_k * u.kcal / (u.mol * u.radian**2)),
                "phi_eq": (improper.type.psi_eq * u.degree),
            }
            _add_conn_type_from_pmd(
                connStr="ImproperType",
                pmd_conn=improper,
                gmso_conn=top_connection,
                conn_params=conn_params,
                name=name,
                expression=expression,
                variables=variables,
            )
        top.add_connection(top_connection, update_types=False)

    top.update_topology()
    return top


def _atom_types_from_pmd(structure):
    """Convert ParmEd atomtypes to GMSO AtomType.

    This function take in a Parmed Structure, iterate through its
    atom's atom_type, create a corresponding GMSO.AtomType, and
    finally return a dictionary containing all pairs of pmd.AtomType
    and GMSO.AtomType

    Parameter
    ----------
        structure: pmd.Structure
            Parmed Structure that needed to be converted.

    Return
    ------
        pmd_top_atomtypes : dict
            A dictionary linking a pmd.AtomType object to its
            corresponding GMSO.AtomType object.
    """
    unique_atom_types = [
        atom.atom_type
        for atom in structure.atoms
        if isinstance(atom.atom_type, pmd.AtomType)
    ]
    pmd_top_atomtypes = {}
    for atom_type in unique_atom_types:
        if atom_type.atomic_number:
            element = element_by_atomic_number(atom_type.atomic_number).symbol
        else:
            element = atom_type.name
        top_atomtype = gmso.AtomType(
            name=atom_type.name,
            charge=atom_type.charge * u.elementary_charge,
            tags={"element": element},
            expression="4*epsilon*((sigma/r)**12 - (sigma/r)**6)",
            parameters={
                "sigma": atom_type.sigma * u.angstrom,
                "epsilon": atom_type.epsilon * u.Unit("kcal / mol"),
            },
            independent_variables={"r"},
            mass=atom_type.mass,
        )
        pmd_top_atomtypes[atom_type] = top_atomtype
    return pmd_top_atomtypes


def _sort_bond_members(top, site_map, atom1, atom2):
    return sorted(
        [site_map[atom1], site_map[atom2]], key=lambda x: top.get_index(x)
    )


def _sort_angle_members(top, site_map, atom1, atom2, atom3):
    sorted_angles = sorted(
        [site_map[atom1], site_map[atom3]], key=lambda x: top.get_index(x)
    )
    return (sorted_angles[0], site_map[atom2], sorted_angles[1])


# function to check reversibility of dihedral type
rev_dih_order = lambda top, site_map, x, y: top.get_index(
    site_map[x]
) > top.get_index(site_map[y])


def _sort_dihedral_members(top, site_map, atom1, atom2, atom3, atom4):
    if rev_dih_order(top, site_map, atom2, atom3):
        return itemgetter(atom4, atom3, atom2, atom1)(site_map)
    return itemgetter(atom1, atom2, atom3, atom4)(site_map)


def _sort_improper_members(top, site_map, atom1, atom2, atom3, atom4):
    sorted_impropers = sorted(
        [site_map[atom2], site_map[atom3], site_map[atom4]],
        key=lambda x: top.get_index(x),
    )
    return (site_map[atom1], *sorted_impropers)


def _add_conn_type_from_pmd(
    connStr, pmd_conn, gmso_conn, conn_params, name, expression, variables
):
    """Create a GMSO connection type and add to the conneciton object.

    This function creates the connection type object and add it to the
    connection object provided.

    Parameters
    ----------
    connStr : str
        The name of the connection type. Accepted values include
        "BondType", "AngleType", "DihedralType", and "ImproperType".
    pmd_conn : pmd.Bond/Angle/Dihedral/Improper
        The parmed connection object.
    gmso_conn : gmso.Bond/Angle/Dihedral/Improper
        The GMSO connection object.
    conn_params : dict
        The potential expression parameters in dictionary form.
    name : str
        Name of the potential form.
    expression :  expression
        The potential expression form.
    variables : dict
        The independent variables.
    """
    try:
        member_types = list(
            map(lambda x: x.atom_type.name, gmso_conn.connection_members)
        )
    except AttributeError:
        member_types = list(
            map(lambda x: f"{x}: {x.atom_type})", gmso_conn.connection_members)
        )
        raise AttributeError(
            f"Parmed structure is missing atomtypes. One of the atomtypes in \
            {member_types} is missing a type from the ParmEd structure.\
            Try using refer_type=False to not look for a parameterized structure."
        )
    try:
        get_classes = lambda x: (
            x.atom_type.atomclass if x.atom_type.atomclass else x.atom_type.name
        )
        member_classes = list(map(get_classes, gmso_conn.connection_members))
    except AttributeError:
        member_classes = list(
            map(
                lambda x: f"{x}: {x.atom_type.name})",
                gmso_conn.connection_members,
            )
        )
    top_conntype = getattr(gmso, connStr)(
        name=name,
        parameters=conn_params,
        expression=expression,
        independent_variables=variables,
        member_types=member_types,
        member_classes=member_classes,
    )
    conntypeStr = connStr.lower()[:-4] + "_type"
    setattr(gmso_conn, conntypeStr, top_conntype)


def to_parmed(top, refer_type=True):
    """Convert a gmso.topology.Topology to a parmed.Structure.

    At this point we only assume a three level structure for topology
    Topology - Molecule - Residue - Sites, which transform to three level of
    Parmed Structure - Residue - Atoms (gmso Molecule level will be skipped).

    Parameters
    ----------
    top : topology.Topology
        topology.Topology instance that need to be converted
    refer_type : bool, optional, default=True
        Whether or not to transfer AtomType, BondType, AngleTye,
        and DihedralType information

    Returns
    -------
    structure : parmed.Structure
    """
    # Sanity check
    msg = "Provided argument is not a topology.Topology."
    assert isinstance(top, gmso.Topology)

    # Set up Parmed structure and define general properties
    structure = pmd.Structure()
    structure.title = top.name
    structure.box = (
        np.concatenate(
            (
                top.box.lengths.to("angstrom").value,
                top.box.angles.to("degree").value,
            )
        )
        if top.box
        else None
    )

    # Maps
    atom_map = dict()  # Map site to atom
    bond_map = dict()  # Map top's bond to structure's bond
    angle_map = dict()  # Map top's angle to strucutre's angle
    dihedral_map = dict()  # Map top's dihedral to structure's dihedral

    # Set up unparametrized system
    # Build up atom
    for site in top.sites:
        if site.element:
            atomic_number = site.element.atomic_number
        else:
            atomic_number = 0
        pmd_atom = pmd.Atom(
            atomic_number=atomic_number,
            name=site.name,
            mass=site.mass.to(u.amu).value if site.mass else None,
            charge=(
                site.charge.to(u.elementary_charge).value
                if site.charge
                else None
            ),
        )
        pmd_atom.xx, pmd_atom.xy, pmd_atom.xz = site.position.to(
            "angstrom"
        ).value

        # Add atom to structure
        if site.residue:
            structure.add_atom(
                pmd_atom,
                resname=site.residue.name,
                resnum=site.residue.number + 1,
            )
        else:
            structure.add_atom(pmd_atom, resname="RES", resnum=-1)
        atom_map[site] = pmd_atom

    # "Claim" all of the item it contains and subsequently index all of its item
    structure.residues.claim()

    # Create and add bonds to Parmed structure
    for bond in top.bonds:
        site1, site2 = bond.connection_members
        pmd_bond = pmd.Bond(atom_map[site1], atom_map[site2])
        structure.bonds.append(pmd_bond)
        bond_map[bond] = pmd_bond

    # Create and add angles to Parmed structure
    for angle in top.angles:
        site1, site2, site3 = angle.connection_members
        pmd_angle = pmd.Angle(atom_map[site1], atom_map[site2], atom_map[site3])
        structure.angles.append(pmd_angle)
        angle_map[angle] = pmd_angle

    # Create and add dihedrals to Parmed structure

    for dihedral in top.dihedrals:
        site1, site2, site3, site4 = dihedral.connection_members
        pmd_dihedral = pmd.Dihedral(
            atom_map[site1], atom_map[site2], atom_map[site3], atom_map[site4]
        )
        if dihedral.connection_type and expand(
            dihedral.connection_type.expression
        ) == expand(
            "c0 * cos(phi)**0 + "
            + "c1 * cos(phi)**1 + "
            + "c2 * cos(phi)**2 + "
            + "c3 * cos(phi)**3 + "
            + "c4 * cos(phi)**4 + "
            + "c5 * cos(phi)**5"
        ):
            structure.rb_torsions.append(pmd_dihedral)
        else:
            structure.dihedrals.append(pmd_dihedral)
        dihedral_map[dihedral] = pmd_dihedral

    # Set up structure for Connection Type conversion
    if refer_type:
        # Need to add a warning if Topology does not have types information
        if top.atom_types:
            _atom_types_from_gmso(top, structure, atom_map)
        if top.bond_types:
            _bond_types_from_gmso(top, structure, bond_map)
        if top.angle_types:
            _angle_types_from_gmso(top, structure, angle_map)
        if top.dihedral_types:
            _dihedral_types_from_gmso(top, structure, dihedral_map)

    return structure


def _check_independent_residues(structure):
    """Check to see if residues will constitute independent graphs."""
    # Copy from foyer forcefield.py
    for res in structure.residues:
        atoms_in_residue = set([*res.atoms])
        bond_partners_in_residue = [
            item
            for sublist in [atom.bond_partners for atom in res.atoms]
            for item in sublist
        ]
        # Handle the case of a 'residue' with no neighbors
        if not bond_partners_in_residue:
            continue
        if set(atoms_in_residue) != set(bond_partners_in_residue):
            return False
    return True


def _atom_types_from_gmso(top, structure, atom_map):
    """Convert gmso.Topology AtomType to parmed.Structure AtomType.

    This function will first check the AtomType expression of Topology and make sure it match with the one default in Parmed.
    After that, it would start atomtyping and parametrizing this part of the structure.

    Parameters
    ----------
    top : topology.Topology
        The topology that need to be converted
    structure: parmed.Structure
        The destination parmed Structure
    """
    # Maps
    atype_map = dict()
    for atom_type in top.atom_types(
        filter_by=PotentialFilters.UNIQUE_NAME_CLASS
    ):
        msg = "Atom type {} expression does not match Parmed AtomType default expression".format(
            atom_type.name
        )
        assert expand(atom_type.expression) == expand(
            "4*epsilon*(-sigma**6/r**6 + sigma**12/r**12)"
        ), msg
        # Extract Topology atom type information
        atype_name = atom_type.name
        # Convert charge to elementary_charge
        atype_charge = float(atom_type.charge.to("Coulomb").value) / (
            1.6 * 10 ** (-19)
        )
        atype_sigma = float(atom_type.parameters["sigma"].to("angstrom").value)
        atype_epsilon = float(
            atom_type.parameters["epsilon"].to("kcal/mol").value
        )
        if atom_type.mass:
            atype_mass = float(atom_type.mass.to("amu").value)
        else:
            atype_mass = float(
                element_by_symbol(atom_type.name).mass.to("amu").value
            )
        atype_atomic_number = getattr(
            element_by_symbol(atom_type.name), "atomic_number", None
        )
        atype_rmin = atype_sigma * 2 ** (1 / 6) / 2  # to rmin/2
        # Create unique Parmed AtomType object
        atype = pmd.AtomType(
            atype_name,
            None,
            atype_mass,
            atype_atomic_number,
            atype_charge,
        )
        atype.set_lj_params(atype_epsilon, atype_rmin)
        # Type map to match AtomType to its name
        atype_map[atype_name] = atype

    for site in top.sites:
        # Assign atom_type to atom
        pmd_atom = atom_map[site]
        pmd_atom.type = site.atom_type.name
        pmd_atom.atom_type = atype_map[site.atom_type.name]


def _bond_types_from_gmso(top, structure, bond_map):
    """Convert gmso.Topology BondType to parmed.Structure BondType.

    This function will first check the BondType expression of Topology and make sure it match with the one default in Parmed.
    After that, it would start atomtyping and parametrizing this part of the structure.

    Parameters
    ----------
    top : topology.Topology
        The topology that need to be converted
    structure: parmed.Structure
        The destination parmed Structure
    """
    btype_map = dict()
    for bond_type in top.bond_types(filter_by=pfilter):
        msg = "Bond type {} expression does not match Parmed BondType default expression".format(
            bond_type.name
        )
        assert expand(bond_type.expression) == expand(
            "0.5 * k * (r-r_eq)**2"
        ), msg
        # Extract Topology bond_type information
        btype_k = 0.5 * float(
            bond_type.parameters["k"].to("kcal / (angstrom**2 * mol)").value
        )
        btype_r_eq = float(bond_type.parameters["r_eq"].to("angstrom").value)
        # Create unique Parmed BondType object
        btype = pmd.BondType(btype_k, btype_r_eq)
        # Type map to match Topology BondType parameters with Parmed BondType
        btype_map[get_parameters(bond_type)] = btype
        # Add BondType to structure.bond_types
        structure.bond_types.append(btype)

    for bond in top.bonds:
        # Assign bond_type to bond
        pmd_bond = bond_map[bond]
        pmd_bond.type = btype_map[get_parameters(bond.bond_type)]
    structure.bond_types.claim()


def _angle_types_from_gmso(top, structure, angle_map):
    """Convert gmso.Topology AngleType to parmed.Structure AngleType.

    This function will first check the AngleType expression of Topology and make sure it match with the one default in Parmed.
    After that, it would start atomtyping and parametrizing the structure.

    Parameters
    ----------
    top : topology.Topology
        The topology that need to be converted
    structure: parmed.Structure
        The destination parmed Structure
    """
    agltype_map = dict()
    for angle_type in top.angle_types(filter_by=pfilter):
        msg = "Angle type {} expression does not match Parmed AngleType default expression".format(
            angle_type.name
        )
        assert expand(angle_type.expression) == expand(
            "0.5 * k * (theta-theta_eq)**2"
        ), msg
        # Extract Topology angle_type information
        agltype_k = 0.5 * float(
            angle_type.parameters["k"].to("kcal / (radian**2 * mol)").value
        )
        agltype_theta_eq = float(
            angle_type.parameters["theta_eq"].to("degree").value
        )
        # Create unique Parmed AngleType object
        agltype = pmd.AngleType(agltype_k, agltype_theta_eq)
        # Type map to match Topology AngleType with Parmed AngleType
        #
        for key, value in agltype_map.items():
            if value == agltype:
                agltype = value
                break
        agltype_map[get_parameters(angle_type)] = agltype
        # Add AngleType to structure.angle_types
        if agltype not in structure.angle_types:
            structure.angle_types.append(agltype)

    for angle in top.angles:
        # Assign angle_type to angle
        pmd_angle = angle_map[angle]
        pmd_angle.type = agltype_map[get_parameters(angle.connection_type)]
    structure.angle_types.claim()


def _dihedral_types_from_gmso(top, structure, dihedral_map):
    """Convert gmso.Topology DihedralType to parmed.Structure DihedralType.

    This function will first check the DihedralType expression of Topology and
    make sure it match with the one default in Parmed.
    After that, it would start atomtyping and parametrizing the structure.

    Parameters
    ----------
    top : topology.Topology
        The topology that need to be converted
    structure: parmed.Structure
        The destination parmed Structure
    """
    dtype_map = dict()
    for dihedral_type in top.dihedral_types(filter_by=pfilter):
        msg = "Dihedral type {} expression does not match Parmed DihedralType default expressions (Periodics, RBTorsions)".format(
            dihedral_type.name
        )
        if expand(dihedral_type.expression) == expand(
            "k * (1 + cos(n * phi - phi_eq))**2"
        ):
            dtype_k = float(dihedral_type.parameters["k"].to("kcal/mol").value)
            dtype_phi_eq = float(
                dihedral_type.parameters["phi_eq"].to("degrees").value
            )
            dtype_n = float(dihedral_type.parameters["n"].value)
            # Create unique Parmed DihedralType object
            dtype = pmd.DihedralType(dtype_k, dtype_n, dtype_phi_eq)
            # Add DihedralType to structure.dihedral_types
            structure.dihedral_types.append(dtype)
        elif expand(dihedral_type.expression) == expand(
            "c0 * cos(phi)**0 + "
            + "c1 * cos(phi)**1 + "
            + "c2 * cos(phi)**2 + "
            + "c3 * cos(phi)**3 + "
            + "c4 * cos(phi)**4 + "
            + "c5 * cos(phi)**5"
        ):
            dtype_c0 = float(
                dihedral_type.parameters["c0"].to("kcal/mol").value
            )
            dtype_c1 = float(
                dihedral_type.parameters["c1"].to("kcal/mol").value
            )
            dtype_c2 = float(
                dihedral_type.parameters["c2"].to("kcal/mol").value
            )
            dtype_c3 = float(
                dihedral_type.parameters["c3"].to("kcal/mol").value
            )
            dtype_c4 = float(
                dihedral_type.parameters["c4"].to("kcal/mol").value
            )
            dtype_c5 = float(
                dihedral_type.parameters["c5"].to("kcal/mol").value
            )
            # Create unique DihedralType object
            dtype = pmd.RBTorsionType(
                dtype_c0,
                dtype_c1,
                dtype_c2,
                dtype_c3,
                dtype_c4,
                dtype_c5,
                list=structure.rb_torsion_types,
            )
            # Add RBTorsionType to structure.rb_torsion_types
            structure.rb_torsion_types.append(dtype)
            # dtype._idx = len(structure.rb_torsion_types) - 1
        else:
            raise GMSOError("msg")
        dtype_map[get_parameters(dihedral_type)] = dtype

    for dihedral in top.dihedrals:
        pmd_dihedral = dihedral_map[dihedral]
        pmd_dihedral.type = dtype_map[get_parameters(dihedral.connection_type)]
    structure.dihedral_types.claim()
    structure.rb_torsions.claim()