text: import datetime
  import os
  from collections import OrderedDict
  from warnings import warn
  
  import numpy as np
  from parmed.periodic_table import Element
  from parmed.utils.io import genopen
  
  from mbuild.box import Box
  from mbuild.compound import Compound
  from mbuild.utils.conversion import (
      RB_to_CHARMM,
      base10_to_base16_alph_num,
      base10_to_base26_alph,
      base10_to_base52_alph,
      base10_to_base62_alph_num,
  )
  from mbuild.utils.sorting import natural_sort
  from mbuild.utils.specific_ff_to_residue import specific_ff_to_residue
  
  
  def _get_bond_type_key(
      bond, sigma_conversion_factor, epsilon_conversion_factor
  ):
      """Get the bond_type key for a bond
  
      Parameters
      ----------
      bond : mbuild.compound.Compound
          The bond information from the mbuild.compound.Compound
      sigma_conversion_factor : float or int
          The sigma conversion factor
      epsilon_conversion_factor : float or int
          The epsilon conversion factor
  
      Returns
      ----------
      tuple, (bond_k_constant, bond_bo_length, bond_atom_1_and_2_types_tuple,
          bond_atom_1_residue_name, bond_atom_2_residue_name)
          bond_k_constant : float
              Harmonic bond k-constant or bond energy scaling factor
          bond_bo_length : float
              Harmonic bonds equilibrium length
          bond_atom_1_and_2_types_tuple : tuple
              A sorted tuple ofstrings for the bonded atom types for atoms 1 and 2.
          bond_atom_1_residue_name : str
              The residue name for atom 1 in the bond.
          bond_atom_2_residue_name : str
              The residue name for atom 2 in the bond.
      """
      bond_k_constant = round(
          bond.type.k
          * (sigma_conversion_factor ** 2 / epsilon_conversion_factor),
          8,
      )
      bond_bo_length = round(bond.type.req / sigma_conversion_factor, 8)
      bond_atom_1_and_2_types_tuple = tuple(
          sorted((bond.atom1.type, bond.atom2.type))
      )
      bond_atom_1_residue_name = bond.atom1.residue.name
      bond_atom_2_residue_name = bond.atom2.residue.name
  
      return (
          bond_k_constant,
          bond_bo_length,
          bond_atom_1_and_2_types_tuple,
          bond_atom_1_residue_name,
          bond_atom_2_residue_name,
      )
  
  
  def _get_angle_type_key(
      angle, sigma_conversion_factor, epsilon_conversion_factor
  ):
      """Get the angle_type key for an harmonic angle
  
      Parameters
      ----------
      angle : parmed.topologyobjects.Angle
          The angle information from the parmed.topologyobjects.Angle
      sigma_conversion_factor : float or int
          The sigma conversion factor
      epsilon_conversion_factor : float or int
          The epsilon conversion factor
  
      Returns
      ----------
      tuple, (angle_k_constant, angle_theta_o, angle_center_atom_type_2, angle_end_atom_types_1_and_3_tuple,
          angle_residue_atom_1, angle_residue_atom_2, angle_residue_atom_3)
          angle_k_constant : float
              Harmonic angle k-constant or bond energy scaling factor
          angle_theta_o : float
              Harmonic equilbrium angle between the atoms
          angle_center_atom_type_2 : str
              The center atom type for the angle (atom type 2)
          angle_end_atom_types_1_and_3_tuple : tuple
              A sorted tuple of atom types (strings) for the end angle atoms
              (atoms 1 and 3).
          angle_atom_1_residue_name : str
              The residue name for atom 1 in the angle (end atom).
          angle_atom_2_residue_name : str
              The residue name for atom 2 in the angle (center atom).
          angle_atom_3_residue_name : str
              The residue name for atom 3 in the angle (end atom).
      """
  
      angle_k_constant = round(
          angle.type.k
          * (sigma_conversion_factor ** 2 / epsilon_conversion_factor),
          8,
      )
      angle_theta_o = round(angle.type.theteq, 8)
      angle_center_atom_type_2 = angle.atom2.type
      angle_end_atom_types_1_and_3_tuple = tuple(
          sorted((angle.atom1.type, angle.atom3.type))
      )
      angle_residue_atom_1 = angle.atom1.residue.name
      angle_residue_atom_2 = angle.atom2.residue.name
      angle_residue_atom_3 = angle.atom3.residue.name
  
      return (
          angle_k_constant,
          angle_theta_o,
          angle_center_atom_type_2,
          angle_end_atom_types_1_and_3_tuple,
          angle_residue_atom_1,
          angle_residue_atom_2,
          angle_residue_atom_3,
      )
  
  
  def _get_dihedral_rb_torsion_key(dihedral, epsilon_conversion_factor):
      """Get the dihedral_type key for a Ryckaert-Bellemans (RB) dihedrals/torsions
  
      Parameters
      ----------
      dihedral : parmed.topologyobjects.Dihedral
          The dihedral information from the parmed.topologyobjects.Angle
      epsilon_conversion_factor : float or int
          The epsilon conversion factor
  
      Returns
      ----------
      tuple, (dihed_type_RB_c0,  dihed_type_RB_c1, dihed_type_RB_c2, dihed_type_RB_c3, dihed_type_RB_c4,
          dihed_type_RB_c5, dihed_type_scee, dihed_type_scnb, dihed_atom_1_type, dihed_atom_2_type,
          dihed_atom_3_type, dihed_atom_4_type, dihed_atom_1_res_type, dihed_atom_2_res_type,
          dihed_atom_3_res_type,  dihed_atom_4_res_type)
          dihed_type_RB_c0 : float
              Ryckaert-Bellemans (RB) dihedrals/torsions C0 constant.
          dihed_type_RB_c1 : float
              Ryckaert-Bellemans (RB) dihedrals/torsions C1 constant.
          dihed_type_RB_c2 : float
              Ryckaert-Bellemans (RB) dihedrals/torsions C2 constant.
          dihed_type_RB_c3 : float
              Ryckaert-Bellemans (RB) dihedrals/torsions C3 constant.
          dihed_type_RB_c4 : float
              Ryckaert-Bellemans (RB) dihedrals/torsions C4 constant.
          dihed_type_RB_c5 : float
              Ryckaert-Bellemans (RB) dihedrals/torsions C5 constant.
          dihed_type_scee : float, default = 1.0
              The 1-4 electrostatic scaling factor
          dihed_type_scnb : float, default = 1.0
              The 1-4 Lennard-Jones scaling factor.
          dihed_atom_1_type : str
              The atom type for atom number 1 in the dihedral
          dihed_atom_2_type : str
              The atom type for atom number 2 in the dihedral
          dihed_atom_3_type : str
              The atom type for atom number 3 in the dihedral
          dihed_atom_4_type : str
              The atom type for atom number 4 in the dihedral
          dihed_atom_1_res_type : str
              The residue name for atom number 1 in the dihedral
          dihed_atom_2_res_type : str
              The residue name for atom number 2 in the dihedral
          dihed_atom_3_res_type : str
              The residue name for atom number 3 in the dihedral
          dihed_atom_4_res_type : str
              The residue name for atom number 4 in the dihedral
      """
  
      lj_unit = 1 / epsilon_conversion_factor
  
      dihed_type_RB_c0 = round(dihedral.type.c0 * lj_unit, 8)
      dihed_type_RB_c1 = round(dihedral.type.c1 * lj_unit, 8)
      dihed_type_RB_c2 = round(dihedral.type.c2 * lj_unit, 8)
      dihed_type_RB_c3 = round(dihedral.type.c3 * lj_unit, 8)
      dihed_type_RB_c4 = round(dihedral.type.c4 * lj_unit, 8)
      dihed_type_RB_c5 = round(dihedral.type.c5 * lj_unit, 8)
  
      dihed_type_scee = round(dihedral.type.scee, 4)
      dihed_type_scnb = round(dihedral.type.scnb, 4)
  
      dihed_atom_1_type = dihedral.atom1.type
      dihed_atom_2_type = dihedral.atom2.type
      dihed_atom_3_type = dihedral.atom3.type
      dihed_atom_4_type = dihedral.atom4.type
  
      dihed_atom_1_res_type = dihedral.atom1.residue.name
      dihed_atom_2_res_type = dihedral.atom2.residue.name
      dihed_atom_3_res_type = dihedral.atom3.residue.name
      dihed_atom_4_res_type = dihedral.atom4.residue.name
  
      return (
          dihed_type_RB_c0,
          dihed_type_RB_c1,
          dihed_type_RB_c2,
          dihed_type_RB_c3,
          dihed_type_RB_c4,
          dihed_type_RB_c5,
          dihed_type_scee,
          dihed_type_scnb,
          dihed_atom_1_type,
          dihed_atom_2_type,
          dihed_atom_3_type,
          dihed_atom_4_type,
          dihed_atom_1_res_type,
          dihed_atom_2_res_type,
          dihed_atom_3_res_type,
          dihed_atom_4_res_type,
      )
  
  
  def _get_improper_type_key(improper, epsilon_conversion_factor):
      """Get the improper_type key for the harmonic improper
  
      Parameters
      ----------
      improper : parmed.topologyobjects.Dihedral
          The improper information from the parmed.topologyobjects.Angle
      epsilon_conversion_factor : float or int
          The epsilon conversion factor
  
      Returns
      ----------
      tuple, (improper_k_constant, improper_psi_o, improper_atom_1_type,
          (improper_atom_2_type, improper_atom_3_type,  improper_atom_4_type), improper_atom_1_res_type,
          (improper_atom_2_res_type, improper_atom_3_res_type, improper_atom_4_res_type)
          improper_k_constant : float
              Harmonic k-constant or bond energy scaling factor
          improper_psi_o : float
              Harmonic equilbrium improper angle
          improper_atom_1_type : str
              The atom type for atom number 1 in the dihedral
          improper_atom_2_type : str
              The atom type for atom number 2 in the dihedral
          improper_atom_3_type : str
              The atom type for atom number 3 in the dihedral
          improper_atom_4_type : str
              The atom type for atom number 4 in the dihedral
          improper_atom_1_res_type : str
              The residue name for atom number 1 in the dihedral
          improper_atom_2_res_type : str
              The residue name for atom number 2 in the dihedral
          improper_atom_3_res_type : str
              The residue name for atom number 3 in the dihedral
          improper_atom_4_res_type : str
              The residue name for atom number 4 in the dihedral
      """
      lj_unit = 1 / epsilon_conversion_factor
  
      improper_k_constant = round(improper.type.psi_k * lj_unit, 8)
      improper_psi_o = round(improper.type.psi_eq, 8)
      improper_atom_1_type = improper.atom1.type
      improper_atom_2_type = improper.atom2.type
      improper_atom_3_type = improper.atom3.type
      improper_atom_4_type = improper.atom4.type
      improper_atom_1_res_type = improper.atom1.residue.name
      improper_atom_2_res_type = improper.atom2.residue.name
      improper_atom_3_res_type = improper.atom3.residue.name
      improper_atom_4_res_type = improper.atom4.residue.name
  
      return (
          improper_k_constant,
          improper_psi_o,
          improper_atom_1_type,
          (improper_atom_2_type, improper_atom_3_type, improper_atom_4_type),
          improper_atom_1_res_type,
          (
              improper_atom_2_res_type,
              improper_atom_3_res_type,
              improper_atom_4_res_type,
          ),
      )
  
  
  def _get_unique_bond_types(
      structure, sigma_conversion_factor, epsilon_conversion_factor
  ):
      """Get the unique bond types for a structure in a dictionary
  
      Parameters
      ----------
      structure : parmed.structure.Structure
         This is a parmed stucture input (parmed.structure.Structure)
      sigma_conversion_factor : float or int
          The sigma conversion factor
      epsilon_conversion_factor : float or int
          The epsilon conversion factor
  
      Returns
      ----------
      bond_key_dict : dict, {(float, float, (str, str), str, str) : unique_number}
          This provides a way to uniquely number the harmonic bond types
          by providing all the bond parameters as a key and the
          unique number as the value. An example of the dict is below:
          {(bond_k_constant, bond_bo_length, (bond_atom_type_1, bond_atom_type_2),
            bond_atom_1_residue_name, bond_atom_2_residue_name ) : 1,
           ...,
           (bond_k_constant, bond_bo_length, (bond_atom_type_1, bond_atom_type_2),
            bond_atom_1_residue_name, bond_atom_2_residue_name ) : n
          }
      """
  
      unique_bond_set = set()
      for bond in structure.bonds:
          unique_bond_set.add(
              _get_bond_type_key(
                  bond, sigma_conversion_factor, epsilon_conversion_factor
              )
          )
  
      bond_key_dict = {
          bond_key: i + 1 for i, bond_key in enumerate(unique_bond_set)
      }
  
      return bond_key_dict
  
  
  def _get_unique_angle_types(
      structure, sigma_conversion_factor, epsilon_conversion_factor
  ):
      """Get the unique angle types for a structure and return a dictionary
  
      Parameters
      ----------
      structure : parmed.structure.Structure
         This is a parmed stucture input (parmed.structure.Structure)
      sigma_conversion_factor : float or int
          The sigma conversion factor
      epsilon_conversion_factor : float or int
          The epsilon conversion factor
  
      Returns
      ----------
      angle_key_dict : dict, {(float, float, str, (str, str), str, str, str) : unique_number}
          This provides a way to uniquely number the harmonic angle types
          by providing all the angle parameters as a key and the
          unique number as the value. An example of the dict is below:
          {(angle_k_constant, angle_theta_o, angle_center_atom_type_2,
            (angle_end_atom_type_1, angle_end_atom_type_3),
             angle_residue_atom_1, angle_residue_atom_2, angle_residue_atom_3
           ) : 1,
           ...,
           (angle_k_constant, angle_theta_o, angle_center_atom_type_2,
           (angle_end_atom_type_1, angle_end_atom_type_3),
            angle_residue_atom_1, angle_residue_atom_2, angle_residue_atom_3) : n
          }
      """
  
      unique_angle_set = set()
      for angle in structure.angles:
          unique_angle_set.add(
              _get_angle_type_key(
                  angle, sigma_conversion_factor, epsilon_conversion_factor
              )
          )
  
      angle_key_dict = {
          angle_key: i + 1 for i, angle_key in enumerate(unique_angle_set)
      }
  
      return angle_key_dict
  
  
  def _get_unique_rb_torsion_types(structure, epsilon_conversion_factor):
      """Get the unique rb torsion types for a structure and return a dictionary
  
      Parameters
      ----------
      structure : parmed.structure.Structure
         This is a parmed stucture input (parmed.structure.Structure)
      epsilon_conversion_factor : float or int
          The epsilon conversion factor
  
      Returns
      ----------
      dihed_key_dict : dict, {(float, float, float, float, float, float, float, float,
                               str, str, str, str, str, str, str, str) : unique_number}
          This provides a way to uniquely number the Ryckaert-Bellemans (RB)
          dihedral types by providing all the dihedral parameters as a key and the
          unique number as the value. An example of the dict is below:
          {(dihed_type_RB_c0, dihed_type_RB_c1, dihed_type_RB_c2,
            dihed_type_RB_c3, dihed_type_RB_c4, dihed_type_RB_c5,
            dihed_type_scee, dihed_type_scnb,
            dihed_atom_1_type, dihed_atom_2_type,
            dihed_atom_3_type, dihed_atom_4_type,
            dihed_atom_1_res_type, dihed_atom_2_res_type,
            dihed_atom_3_res_type, dihed_atom_4_res_type
           ) : 1,
           ...,
          (dihed_type_RB_c0, dihed_type_RB_c1, dihed_type_RB_c2,
            dihed_type_RB_c3, dihed_type_RB_c4, dihed_type_RB_c5,
            dihed_type_scee, dihed_type_scnb,
            dihed_atom_1_type, dihed_atom_2_type,
            dihed_atom_3_type, dihed_atom_4_type,
            dihed_atom_1_res_type, dihed_atom_2_res_type,
            dihed_atom_3_res_type, dihed_atom_4_res_type
           ) : n
          }
      """
      unique_dihedral_set = set()
      for dihedral in structure.rb_torsions:
          unique_dihedral_set.add(
              _get_dihedral_rb_torsion_key(dihedral, epsilon_conversion_factor)
          )
  
      dihed_key_dict = {
          dihed_key: i + 1 for i, dihed_key in enumerate(unique_dihedral_set)
      }
  
      return dihed_key_dict
  
  
  def _get_unique_improper_types(structure, epsilon_conversion_factor):
      """Get the unique improper types for a structure  and return a dictionary
  
      Parameters
      ----------
      structure : parmed.structure.Structure
         This is a parmed stucture input (parmed.structure.Structure)
      epsilon_conversion_factor : float or int
          The epsilon conversion factor
  
      Returns
      ----------
      improper_key_dict : dict, {(float, float, str, (str, str, str), str, (str, str, str)) : unique_number}
          This provides a way to uniquely number the harmonic improper
          types by providing all the improper parameters as a key and the
          unique number as the value. An example of the dict is below:
          {(improper_k_constant, improper_psi_o,
            improper_atom_1_type, (improper_atom_2_type,
            improper_atom_3_type, improper_atom_4_type),
            improper_atom_1_res_type, (improper_atom_2_res_type,
            improper_atom_3_res_type, improper_atom_4_res_type)
           ) : 1,
           ...,
           (improper_k_constant, improper_psi_o,
            improper_atom_1_type, (improper_atom_2_type,
            improper_atom_3_type, improper_atom_4_type),
            improper_atom_1_res_type, (improper_atom_2_res_type,
            improper_atom_3_res_type, improper_atom_4_res_type)
           ) : n
          }
      """
      unique_improper_set = set()
      for improper in structure.impropers:
          unique_improper_set.add(
              _get_improper_type_key(improper, epsilon_conversion_factor)
          )
  
      improper_key_dict = {
          improper_key: i + 1
          for i, improper_key in enumerate(unique_improper_set)
      }
  
      return improper_key_dict
  
  
  def _get_bond_types(
      structure, sigma_conversion_factor, epsilon_conversion_factor
  ):
      """Get a list of unique bond types that are used to create the Charmm style
      parameter file (i.e. force field file).  This also checks that the reverse bond order
      is considered the same unique bond type.
  
      Parameters
      ----------
      structure : parmed.Structure
          The parmed structure
      sigma_conversion_factor : float or int
          The sigma conversion factor
      epsilon_conversion_factor : float or int
          The epsilon conversion factor
  
      Returns
      ----------
      bond_types : list
          The bond types by number in the structure
      unique_bond_types : OrderedDict, ((float, float, (str, str), str, str), unique_number)
          This provides the unique harmonic bond types, numbering, and the data values
          so it can easily be extracted.  An example of the OrderedDict is below:
          OrderedDict([(bond_k_constant, bond_bo_length, (bond_atom_type_1, bond_atom_type_2),
                        bond_atom_1_residue_name, bond_atom_2_residue_name), 1),
                       ...,
                      (bond_k_constant, bond_bo_length, (bond_atom_type_1, bond_atom_type_2),
                       bond_atom_1_residue_name, bond_atom_2_residue_name), n)])
      """
  
      unique_bond_types = _get_unique_bond_types(
          structure, sigma_conversion_factor, epsilon_conversion_factor
      )
  
      bond_types = [
          unique_bond_types[
              _get_bond_type_key(
                  bond, sigma_conversion_factor, epsilon_conversion_factor
              )
          ]
          for bond in structure.bonds
      ]
  
      unique_bond_check_dict = {}
      for i_value_bond, i_key_bond in unique_bond_types.items():
          i_value_duplicated = False
          for j_value_bond, j_key_bond in unique_bond_types.items():
              j_value_bond_reorder = (
                  j_value_bond[0],
                  j_value_bond[1],
                  j_value_bond[2][0],
                  j_value_bond[2][0],
                  j_value_bond[3],
                  j_value_bond[4],
              )
  
              if i_value_bond == j_value_bond_reorder:
                  i_value_duplicated = True
                  if i_value_bond[2][0] > j_value_bond[2][0]:
                      unique_bond_check_dict.update(
                          {j_value_bond: len(unique_bond_check_dict)}
                      )
                  else:
                      unique_bond_check_dict.update(
                          {i_value_bond: len(unique_bond_check_dict)}
                      )
  
              if i_value_duplicated is False:
                  unique_bond_check_dict.update(
                      {i_value_bond: len(unique_bond_check_dict)}
                  )
  
      unique_bond_types = OrderedDict(
          [(y, x) for y, x in unique_bond_check_dict.items()]
      )
  
      return bond_types, unique_bond_types
  
  
  def _get_angle_types(
      structure,
      sigma_conversion_factor,
      epsilon_conversion_factor,
      use_urey_bradleys=False,
  ):
      """
      Get a list of unique angle types that are used to create the Charmm style
      parameter file (i.e. force field file).  This also checks that the alternately
      ordered angle types are considered the same unique angle type.
  
      Parameters
      ----------
      structure : parmed.Structure
          The parmed structure
      sigma_conversion_factor : float or int
          The sigma conversion factor
      epsilon_conversion_factor : float or int
          The epsilon conversion factor
      use_urey_bradleys : bool
          The option that Urey-Bradleys are included in the angles
  
      Returns
      ----------
      angle_types : list
          The angle types by number in the structure
      unique_angle_types : OrderedDict, ((float, float, (str, str), str, str), unique_number)
          This provides the unique harmonic angles types, numbering, and the data values
          so it can easily be extracted.  An example of the OrderedDict is below:
          OrderedDict([(angle_k_constant, angle_theta_o, angle_center_atom_type_2,
                       (angle_end_atom_type_1, angle_end_atom_type_3),
                        angle_residue_atom_1, angle_residue_atom_2, angle_residue_atom_3), 1,
                       ...,
                       (angle_k_constant, angle_theta_o, angle_center_atom_type_2,
                       (angle_end_atom_type_1, angle_end_atom_type_3),
                       angle_residue_atom_1, angle_residue_atom_2, angle_residue_atom_3), n])
      """
  
      if use_urey_bradleys:
          print_warn_text = (
              "WARNING: Urey-Bradleys are not available in the current "
              "version of this psf, pdb, and GOMC writer."
          )
          warn(print_warn_text)
          return None, None
      else:
          unique_angle_types = _get_unique_angle_types(
              structure, sigma_conversion_factor, epsilon_conversion_factor
          )
  
          angle_types = [
              unique_angle_types[
                  _get_angle_type_key(
                      angle, sigma_conversion_factor, epsilon_conversion_factor
                  )
              ]
              for angle in structure.angles
          ]
  
      unique_angle_check_dict = {}
      for i_value_ang, i_key_ang in unique_angle_types.items():
          i_value_duplicated = False
          for j_value_ang, j_key_ang in unique_angle_types.items():
              j_value_ang_reorder = (
                  j_value_ang[0],
                  j_value_ang[1],
                  j_value_ang[2],
                  j_value_ang[3][0],
                  j_value_ang[3][1],
                  j_value_ang[4],
                  j_value_ang[5],
                  j_value_ang[6],
              )
  
              if i_value_ang == j_value_ang_reorder:
                  i_value_duplicated = True
                  if i_value_ang[2] > j_value_ang[2]:
                      unique_angle_check_dict.update(
                          {j_value_ang: len(unique_angle_check_dict)}
                      )
                  else:
                      unique_angle_check_dict.update(
                          {i_value_ang: len(unique_angle_check_dict)}
                      )
  
          if not i_value_duplicated:
              unique_angle_check_dict.update(
                  {i_value_ang: len(unique_angle_check_dict)}
              )
  
      unique_angle_types = OrderedDict(
          [(y, x) for y, x in unique_angle_check_dict.items()]
      )
  
      return angle_types, unique_angle_types
  
  
  def _get_dihedral_types(
      structure, use_rb_torsions, use_dihedrals, epsilon_conversion_factor
  ):
      """
      Get a list of unique dihedral types that are used to create the Charmm style
      parameter file (i.e. force field file).  This also checks that the alternately
      ordered dihedral types are considered the same unique dihedral type.
  
      Parameters
      ----------
      structure : parmed.Structure
          The parmed structure
      use_rb_torsions : bool
          The Ryckaert-Bellemans (RB) dihedrals/torsions
      use_dihedrals : bool
          The CHARMM style dihedrals style equations.
      epsilon_conversion_factor : float or int
          The epsilon conversion factor
  
      Returns
      ----------
      dihedral_types : list
          The dihedral types by number in the structure
      unique_dihedral_types : OrderedDict, ([(float, float, float, float, float, float, float, float,
                                              str, str, str, str, str, str, str, str), unique_number])
          This provides the unique dihedrals types, numbering, and the data values
          so it can easily be extracted.  An example of the OrderedDict is below:
          OrderedDict([(dihed_type_RB_c0, dihed_type_RB_c1, dihed_type_RB_c2,
                        dihed_type_RB_c3, dihed_type_RB_c4, dihed_type_RB_c5,
                        dihed_type_scee, dihed_type_scnb,
                        dihed_atom_1_type, dihed_atom_2_type,
                        dihed_atom_3_type, dihed_atom_4_type,
                        dihed_atom_1_res_type, dihed_atom_2_res_type,
                        dihed_atom_3_res_type, dihed_atom_4_res_type), 1,
                       ...,
                       (dihed_type_RB_c0, dihed_type_RB_c1, dihed_type_RB_c2,
                        dihed_type_RB_c3, dihed_type_RB_c4, dihed_type_RB_c5,
                        dihed_type_scee, dihed_type_scnb,
                        dihed_atom_1_type, dihed_atom_2_type,
                        dihed_atom_3_type, dihed_atom_4_type,
                        dihed_atom_1_res_type, dihed_atom_2_res_type,
                        dihed_atom_3_res_type, dihed_atom_4_res_type), n])
      """
      if use_rb_torsions:
          unique_dihedral_types = _get_unique_rb_torsion_types(
              structure, epsilon_conversion_factor
          )
  
          dihedral_types = [
              unique_dihedral_types[
                  _get_dihedral_rb_torsion_key(
                      dihedral, epsilon_conversion_factor
                  )
              ]
              for dihedral in structure.rb_torsions
          ]
  
      elif use_dihedrals:
          print_warn_text = (
              "WARNING: Using the charmm style and impropers is not "
              "available in the current version of this psf, pdb, and GOMC writer."
          )
          warn(print_warn_text)
          return None, None
  
      unique_dihedral_check_dict = OrderedDict()
      for i_value_dihed, i_key_dihed in unique_dihedral_types.items():
          i_value_duplicated = False
          for j_value_dihed, j_key_dihed in unique_dihedral_types.items():
              j_value_dihed_reorder = (
                  j_value_dihed[0],
                  j_value_dihed[1],
                  j_value_dihed[2],
                  j_value_dihed[3],
                  j_value_dihed[4],
                  j_value_dihed[5],
                  j_value_dihed[6],
                  j_value_dihed[7],
                  j_value_dihed[11],
                  j_value_dihed[10],
                  j_value_dihed[9],
                  j_value_dihed[8],
                  j_value_dihed[15],
                  j_value_dihed[14],
                  j_value_dihed[13],
                  j_value_dihed[12],
              )
  
              if i_value_dihed == j_value_dihed_reorder:
                  i_value_duplicated = True
                  if i_value_dihed[8] > j_value_dihed[8]:
                      unique_dihedral_check_dict.update(
                          {j_value_dihed: len(unique_dihedral_check_dict) + 1}
                      )
                  else:
                      unique_dihedral_check_dict.update(
                          {i_value_dihed: len(unique_dihedral_check_dict) + 1}
                      )
          if i_value_duplicated is False:
              unique_dihedral_check_dict.update(
                  {i_value_dihed: len(unique_dihedral_check_dict) + 1}
              )
  
      unique_dihedral_types = OrderedDict(
          [(y, x) for y, x in unique_dihedral_check_dict.items()]
      )
  
      return dihedral_types, unique_dihedral_types
  
  
  def _get_impropers(structure, epsilon_conversion_factor):
      """
      Get a list of unique improper types that are used to create the Charmm style
      parameter file (i.e. force field file).  This also checks that the alternately
      ordered improper types are considered the same unique improper type.
  
      Parameters
      ----------
      structure : parmed.Structure
          The parmed structure
      epsilon_conversion_factor : float or int
          The epsilon conversion factor
  
      Returns
      ----------
      improper_types : list
          The improper types by number in the structure
      unique_improper_types : OrderedDict, ([(float, float, str, (str, str, str), str, (str, str, str)), unique_number])
          This provides the unique improper types, numbering, and the data values
          so it can easily be extracted.  An example of the OrderedDict is below:
          OrderedDict([(improper_k_constant, improper_psi_o,
            improper_atom_1_type, (improper_atom_2_type,
            improper_atom_3_type, improper_atom_4_type),
            improper_atom_1_res_type, (improper_atom_2_res_type,
            improper_atom_3_res_type, improper_atom_4_res_type)
           ), 1,
           ...,
           (improper_k_constant, improper_psi_o,
            improper_atom_1_type, (improper_atom_2_type,
            improper_atom_3_type, improper_atom_4_type),
            improper_atom_1_res_type, (improper_atom_2_res_type,
            improper_atom_3_res_type, improper_atom_4_res_type)
           ), n ])
      """
      unique_improper_types = _get_unique_improper_types(
          structure, epsilon_conversion_factor
      )
      improper_types = [
          unique_improper_types[
              _get_improper_type_key(improper, epsilon_conversion_factor)
          ]
          for improper in structure.impropers
      ]
  
      unique_improper_check_dict = OrderedDict()
      for i_value_improper, i_key_improper in unique_improper_types.items():
          i_value_duplicated = False
  
          i_value_improper_k_constant = i_value_improper[0]
          i_value_improper_psi_o = i_value_improper[1]
  
          i_value_impr_1_atoms_reorder = i_value_improper[2]
          i_value_impr_atoms_2_3_4_reorder = sorted(
              set(
                  i_value_improper[3][0],
                  i_value_improper[3][1],
                  i_value_improper[3][2],
              )
          )
  
          i_value_impr_1_res_reorder = i_value_improper[4]
          i_value_impr_res_2_3_4_reorder = sorted(
              set(
                  i_value_improper[5][0],
                  i_value_improper[5][1],
                  i_value_improper[5][2],
              )
          )
          i_improper_reformed = (
              i_value_improper_k_constant,
              i_value_improper_psi_o,
              i_value_impr_1_atoms_reorder,
              i_value_impr_atoms_2_3_4_reorder,
              i_value_impr_1_res_reorder,
              i_value_impr_res_2_3_4_reorder,
          )
  
          for j_value_improper, j_key_improper in unique_improper_types.items():
              j_value_improper_k_constant = j_value_improper[0]
              j_value_improper_psi_o = j_value_improper[1]
  
              j_value_impr_1_atoms_reorder = j_value_improper[2]
              j_value_impr_atoms_2_3_4_reorder = sorted(
                  set(
                      j_value_improper[3][0],
                      j_value_improper[3][1],
                      j_value_improper[3][2],
                  )
              )
              j_value_impr_1_res_reorder = j_value_improper[4]
              j_value_impr_res_2_3_4_reorder = sorted(
                  set(
                      j_value_improper[5][0],
                      j_value_improper[5][1],
                      j_value_improper[5][2],
                  )
              )
  
              j_improper_reformed = (
                  j_value_improper_k_constant,
                  j_value_improper_psi_o,
                  j_value_impr_1_atoms_reorder,
                  j_value_impr_atoms_2_3_4_reorder,
                  j_value_impr_1_res_reorder,
                  j_value_impr_res_2_3_4_reorder,
              )
  
              if i_improper_reformed == j_improper_reformed:
                  i_value_duplicated = True
                  if i_value_improper[2] > j_value_improper[2]:
                      unique_improper_check_dict.update(
                          {j_value_improper: len(unique_improper_check_dict) + 1}
                      )
                  else:
                      unique_improper_check_dict.update(
                          {i_value_improper: len(unique_improper_check_dict) + 1}
                      )
          if i_value_duplicated is False:
              unique_improper_check_dict.update(
                  {i_value_improper: len(unique_improper_check_dict) + 1}
              )
  
      unique_improper_types = OrderedDict(
          [(y, x) for y, x in unique_improper_check_dict.items()]
      )
  
      return improper_types, unique_improper_types
  
  
  def unique_atom_naming(
      structure, residue_id_list, residue_names_list, bead_to_atom_name_dict=None
  ):
      """
      Generates unique atom/bead names for each molecule, which is required for some
      simulation types (Example: special Monte Carlo moves)
  
      Parameters
      ----------
      structure : compound object
      residue_id_list : list, in sequential order
              The residue ID for every atom in the system
      residue_names_list : list, in sequential order
          The atom names for every atom in the system
      bead_to_atom_name_dict: dictionary ; optional, default =None
          For all atom names/elements/beads with 2 or less digits, this converts
          the atom name in the GOMC psf and pdb files to a unique atom name,
          provided they do not exceed 3844 atoms (62^2) of the same name/element/bead
          per residue. For all atom names/elements/beads with 3 digits, this converts
          the atom name in the GOMC psf and pdb files to a unique atom name,
          provided they do not exceed 62 of the same name/element pre residue.
          Example dictionary: {'_CH3':'C', '_CH2':'C', '_CH':'C', '_HC':'C'}
  
      Returns
       ----------
      unique_individual_atom_names_dict : dictionary
          All the unique atom names comno_piled into a dictionary.
      individual_atom_names_list : list, in sequential  order
          The atom names for every atom in the system
      missing_bead_to_atom_name : list, in sequential  order
          The bead names of any atoms beads that did not have a name specificed to them
          via the bead_to_atom_name_dict
      """
      unique_individual_atom_names_dict = {}
      individual_atom_names_list = []
      missing_bead_to_atom_name = []
      for i, atom in enumerate(structure.atoms):
          interate_thru_names = True
          j = 0
          while interate_thru_names is True:
              j = j + 1
              if str(atom.name)[:1] == "_":
                  if (
                      bead_to_atom_name_dict is not None
                      and (str(atom.name) in bead_to_atom_name_dict) is True
                  ):
                      if len(bead_to_atom_name_dict[str(atom.name)]) > 2:
                          text_to_write = (
                              "ERROR: only enter atom names that have 2 or less digits"
                              + " in the Bead to atom naming dictionary (bead_to_atom_name_dict)."
                          )
                          warn(text_to_write)
                          return None, None, None
                      else:
                          atom_name_value = bead_to_atom_name_dict[str(atom.name)]
                          no_digits_atom_name = 2
                  else:
                      missing_bead_to_atom_name.append(1)
                      atom_name_value = "BD"
                      no_digits_atom_name = 2
              elif len(str(atom.name)) > 2:
                  if len(str(atom.name)) == 3:
                      no_digits_atom_name = 1
                      atom_name_value = atom.name
                  else:
                      text_to_write = (
                          "ERROR: atom numbering will not work propery at"
                          + " the element has more than 4 charaters"
                      )
                      warn(text_to_write)
                      return None, None, None
              else:
                  no_digits_atom_name = 2
                  atom_name_value = atom.name
              atom_name_iteration = str(atom_name_value) + str(
                  base10_to_base62_alph_num(j)
              )
              atom_res_no_resname_atomname_iteration = (
                  str(residue_id_list[i])
                  + "_"
                  + str(residue_names_list[i])
                  + "_"
                  + atom_name_iteration
              )
  
              if (
                  unique_individual_atom_names_dict.get(
                      str(atom_res_no_resname_atomname_iteration)
                  )
                  is None
              ):
                  unique_individual_atom_names_dict.update(
                      {atom_res_no_resname_atomname_iteration: i + 1}
                  )
                  interate_thru_names = False
                  individual_atom_names_list.append(
                      str(atom_name_value)
                      + str(
                          str(base10_to_base62_alph_num(j))[-no_digits_atom_name:]
                      )
                  )
  
      if sum(missing_bead_to_atom_name) > 0:
          warn(
              "NOTE: All bead names were not found in the Bead to atom naming dictionary (bead_to_atom_name_dict) "
          )
  
      return [
          unique_individual_atom_names_dict,
          individual_atom_names_list,
          missing_bead_to_atom_name,
      ]
  
  
  def _lengths_angles_to_vectors(lengths, angles, precision=6):
      """Converts the length and angles into CellBasisVectors
  
      Parameters
      ----------
      lengths : list-like, shape=(3,), dtype=float
          Lengths of the edges of the box (user chosen units).
      angles : list-like, shape=(3,), dtype=float, default=None
          Angles (in degrees) that define the tilt of the edges of the box. If
          None is given, angles are assumed to be [90.0, 90.0, 90.0]. These are
          also known as alpha, beta, gamma in the crystallography community.
      precision : int, optional, default=6
          Control the precision of the floating point representation of box
          attributes. If none provided, the default is 6 decimals.
  
      Returns
      -------
      box_vectors: numpy.ndarray, [[float, float, float], [float, float, float], [float, float, float]]
          Three (3) sets vectors for box 0 each with 3 float values, which represent
          the vectors for the Charmm-style systems (units are the same as entered for lengths)
  
      """
  
      (a, b, c) = lengths
  
      (alpha, beta, gamma) = np.deg2rad(angles)
      cos_a = np.clip(np.cos(alpha), -1.0, 1.0)
      cos_b = np.clip(np.cos(beta), -1.0, 1.0)
      cos_g = np.clip(np.cos(gamma), -1.0, 1.0)
  
      sin_a = np.clip(np.sin(alpha), -1.0, 1.0)
      sin_b = np.clip(np.sin(beta), -1.0, 1.0)
      sin_g = np.clip(np.sin(gamma), -1.0, 1.0)
      a_vec = np.asarray([a, 0.0, 0.0])
  
      b_x = b * cos_g
      b_y = b * sin_g
      b_vec = np.asarray([b_x, b_y, 0.0])
  
      c_x = c * cos_b
      c_cos_y_term = (cos_a - (cos_b * cos_g)) / sin_g
      c_y = c * c_cos_y_term
      c_z = c * np.sqrt(1 - np.square(cos_b) - np.square(c_cos_y_term))
      c_vec = np.asarray([c_x, c_y, c_z])
      box_vectors = np.asarray((a_vec, b_vec, c_vec))
      box_vectors.reshape(3, 3)
      # still leaves some floating values in some cases
      box_vectors = np.around(box_vectors, decimals=precision)
  
      return box_vectors
  
  
  def _check_fixed_bonds_angles_lists(
      gomc_fix_bonds_and_or_angles,
      gomc_fix_bonds_and_or_angles_selection,
      residues,
  ):
      """Check the GOMC fixed bonds and angles lists for input errors.
  
      Parameters
      ----------
      gomc_fix_bonds_and_or_angles : list of strings, [str, ..., str]
          A list of the residues (i.e., molecules since GOMC currently considers a
          a whole molecule as a residue) to have their bonds and/or angles held
          rigid/fixed for the GOMC simulation engine.
          The `gomc_fix_bonds_angles`, `gomc_fix_bonds`, `gomc_fix_angles` are the only possible
          variables from the `Charmm` object to be entered.
          In GOMC, the residues currently are the same for every bead or atom in
          the molecules. Therefore, when the residue is selected, the whole molecule
          is selected.
      gomc_fix_bonds_and_or_angles_selection : str
          The name of the variable that is used but formatted as a string, which is fed
          to the error and information outputs. The
          `gomc_fix_bonds_angles`, `gomc_fix_bonds`, `gomc_fix_angles` are the only possible
          variables from the `Charmm` object to be entered.
          Whichever variable you choose, the variable name is just input as a
          string here. For example, if `gomc_fix_bonds_and_or_angles` is equal to
          gomc_fix_bonds_angles, then this should be 'gomc_fix_bonds_angles'
          (i.e., `gomc_fix_bonds_and_or_angles_selection` = 'gomc_fix_bonds_angles').
      residues : list, [str, ..., str]
          Labels of unique residues in the Compound. Residues are assigned by
          checking against Compound.name.  Only supply residue names as 4 character
          strings, as the residue names are truncated to 4 characters to fit in the
          psf and pdb file.
  
      Returns
      -------
      Provides a ValueError or TypeError if the input is not correct.
      """
  
      if gomc_fix_bonds_and_or_angles is not None and not isinstance(
          gomc_fix_bonds_and_or_angles, list
      ):
          print_error_message = (
              "ERROR: Please ensure the residue names in the ({}) variable "
              "are in a list.".format(gomc_fix_bonds_and_or_angles_selection)
          )
          raise TypeError(print_error_message)
  
      if isinstance(gomc_fix_bonds_and_or_angles, list):
          for gomc_fix_i in gomc_fix_bonds_and_or_angles:
              if gomc_fix_i not in residues:
                  print_error_message = (
                      "ERROR: Please ensure that all the residue names in the "
                      "{} list are also in the residues list.".format(
                          gomc_fix_bonds_and_or_angles_selection
                      )
                  )
                  raise ValueError(print_error_message)
              elif not isinstance(gomc_fix_i, str):
                  print_error_message = "ERROR: Please enter a fix_res_bonds list with only string values."
                  raise TypeError(print_error_message)
              else:
                  print(
                      "INFORMATION: The following residues will have these fixed parameters: "
                      + "gomc_fix_bonds = {}".format(gomc_fix_bonds_and_or_angles)
                  )
  
  
  # Currently the NBFIX is disabled as since only the OPLS and TRAPPE force fields are currently supported
  class Charmm:
      """Generates a Charmm object that is required to produce the Charmm style parameter
      (force field), PDB, PSF files, which are usable in the GOMC and NAMD engines.
      Additionally, this Charmm object is also used in generating the GOMC control file.
  
      The units for the GOMC data files.
          * Mw = g/mol
          * charge = e
          * Harmonic bonds : Kb = kcal/mol, b0 = Angstroms
          * Harmonic angles : Ktheta = kcal/mole/rad**2 , Theta0 = degrees
          * Dihedral angles: Ktheta = kcal/mole, n = interger (unitless), delta = degrees
          * Improper angles (currently unavailable) : TBD
          * LJ-NONBONDED : epsilon = kcal/mol, Rmin/2 = Angstroms
          * Mie-NONBONDED (currently unavailable): epsilon = K, sigma = Angstroms, n = interger (unitless)
          * Buckingham-NONBONDED (currently unavailable): epsilon = K, sigma = Angstroms, n = interger (unitless)
          * LJ-NBFIX (currently unavailable) : epsilon = kcal/mol, Rmin = Angstroms
          * Mie-NBFIX (currently unavailable) : same as Mie-NONBONDED
          * Buckingham-NBFIX (currently unavailable) : same as Buckingham-NONBONDED
  
      Note: units are the same as the NAMD units and the LAMMPS real units.  The atom style
      is the same as the lammps 'full' atom style format.
  
      Parameters
      ----------
      structure_box_0 : mbuild Compound object (mbuild.Compound) or mbuild Box object (mbuild.Box);
          If the structure has atoms/beads it must be an mbuild Compound.
          If the structure is empty it must be and mbuild Box object.
          Note: If 1 structures are provided (i.e., only structure_box_0),
          it must be an mbuild Compound.
          Note: If 2 structures are provided,
          only 1 structure can be an empty box (i.e., either structure_box_0 or structure_box_1)
      filename_box_0 : str
          The file name of the output file for structure_box_0.  Note: the extension should
          not be provided, as multiple extension (.pdb and .psf) are added to this name.
      structure_box_1 : mbuild Compound object (mbuild.Compound) or mbuild Box object (mbuild.Box), default = None;
          If the structure has atoms/beads it must be an mbuild Compound.
          Note: When running a GEMC or GCMC simulation the box 1 stucture should be input
          here.  Otherwise, there is no guarantee that any of the atom type and force field
          information will all work together correctly with box 0, if it is built separately.
          Note: If 2 structures are provided, only 1 structure can be an empty box
          (i.e., either structure_box_0 or structure_box_1).
      filename_box_1 : str , default = None
          The file name of the output file for structure_box_1 (Ex: for GCMC or GEMC simulations
          which have mulitiple simulation boxes).  Note: the extension should
          not be provided, as multiple extension (.pdb and .psf) are added to this name.
          Note: When running a GEMC or GCMC simulation the box 1 stucture should be input
          here.  Otherwise, there is no guarantee that any of the atom type and force field
          information will all work together correctly with box 0, if it is built separately.
      non_bonded_type : str, default = 'LJ' (i.e., Lennard-Jones )
          Specify the type of non-bonded potential for the GOMC force field files.
          Note: Currently, on the 'LJ' potential is supported.
      residues : list, [str, ..., str]
          Labels of unique residues in the Compound. Residues are assigned by
          checking against Compound.name.  Only supply residue names as 4 character
          strings, as the residue names are truncated to 4 characters to fit in the
          psf and pdb file.
      forcefield_selection : str or dictionary, default = None
          Apply a forcefield to the output file by selecting a force field XML file with
          its path or by using the standard force field name provided the `foyer` package.
          Note: to write the NAMD/GOMC force field, pdb, and psf files, the
          residues and forcefields must be provided in a str or
          dictionary.  If a dictionary is provided all residues must
          be specified to a force field.
          * Example dict for FF file: {'ETH' : 'oplsaa.xml', 'OCT': 'path_to_file/trappe-ua.xml'}
  
          * Example str for FF file: 'path_to file/trappe-ua.xml'
  
          * Example dict for standard FF names : {'ETH' : 'oplsaa', 'OCT': 'trappe-ua'}
  
          * Example str for standard FF names: 'trappe-ua'
  
          * Example of a mixed dict with both : {'ETH' : 'oplsaa', 'OCT': 'path_to_file/'trappe-ua.xml'}
      detect_forcefield_style : boolean, default = True
          If True, format NAMD/GOMC/LAMMPS parameters based on the contents of
          the parmed structure_box_0 and structure_box_1.
      gomc_fix_bonds_angles : list, default = None
          When list of residues is provided, the selected residues will have
          their bonds and angles fixed in the GOMC engine.  This is specifically
          for the GOMC engine and it changes the residue's bond constants (Kbs)
          and angle constants (Kthetas) values to 999999999999 in the
          FF file (i.e., the .inp file).
      bead_to_atom_name_dict : dict, optional, default =None
          For all atom names/elements/beads with 2 or less digits, this converts
          the atom name in the GOMC psf and pdb files to a unique atom name,
          provided they do not exceed 3844 atoms (62^2) of the same name/element/bead
          per residue. For all atom names/elements/beads with 3 digits, this converts
          the atom name in the GOMC psf and pdb files to a unique atom name,
          provided they do not exceed 62 of the same name/element pre residue.
  
          * Example dictionary: {'_CH3':'C', '_CH2':'C', '_CH':'C', '_HC':'C'}
  
          * Example name structure: {atom_type: first_part_pf atom name_without_numbering}
  
      fix_residue : list  or None, default = None
          Changes occcur in the pdb file only.
          When residues are listed here, all the atoms in the residue are
          fixed and can not move via setting the Beta values in the PDB
          file to 1.00.
          If neither fix_residue or fix_residue_in_box lists a
          residue or both equal None, then the Beta values for all the atoms
          in the residue are free to move in the simulation and Beta values
          in the PDB file is set to 0.00
      fix_residue_in_box : list  or None, default = None
          Changes occcur in the pdb file only.
          When residues are listed here, all the atoms in the residue
          can move within the box but cannot be transferred between boxes
          via setting the Beta values in the PDB file to 2.00.
          If neither fix_residue or fix_residue_in_box lists a
          residue or both equal None, then the Beta values for all the atoms
          in the residue are free to move in the simulation and Beta values
          in the PDB file is set to 0.00
      ff_filename : str, default =None
          If a string, it will write the  force field files that work in
          GOMC and NAMD structures.
      reorder_res_in_pdb_psf : bool, default =False
          If False, the order of of the atoms in the pdb file is kept in
          its original order, as in the Compound sent to the writer.
          If True, the order of the atoms is reordered based on their
          residue names in the 'residues' list that was entered.
  
      Attributes
      ----------
      input_error : bool
          This error is typically incurred from an error in the user's input values.
          However, it could also be due to a bug, provided the user is inputting
          the data as this Class intends.
      structure_box_0 : mbuild.compound.Compound
          The mbuild Compound for the input box 0
      structure_box_1 : mbuild.compound.Compound or None, default = None
          The mbuild Compound for the input box 1
      filename_box_0 : str
          The file name of the output file for structure_box_0.  Note: the extension should
          not be provided, as multiple extension (.pdb and .psf) are added to this name.
      filename_box_1 : str or None , default = None
          The file name of the output file for structure_box_1.  Note: the extension should
          not be provided, as multiple extension (.pdb and .psf) are added to this name.
          (i.e., either structure_box_0 or structure_box_1).
      non_bonded_type : str, default = 'LJ' (i.e., Lennard-Jones )
          Specify the type of non-bonded potential for the GOMC force field files.
          Note: Currently, on the 'LJ' potential is supported.
      residues : list, [str, ..., str]
          Labels of unique residues in the Compound. Residues are assigned by
          checking against Compound.name.  Only supply residue names as 4 character
          strings, as the residue names are truncated to 4 characters to fit in the
          psf and pdb file.
      forcefield_selection : str or dictionary, default = None
          Apply a forcefield to the output file by selecting a force field XML file with
          its path or by using the standard force field name provided the `foyer` package.
          Note: to write the NAMD/GOMC force field, pdb, and psf files, the
          residues and forcefields must be provided in a str or
          dictionary.  If a dictionary is provided all residues must
          be specified to a force field.
  
          * Example dict for FF file: {'ETH' : 'oplsaa.xml', 'OCT': 'path_to_file/trappe-ua.xml'}
  
          * Example str for FF file: 'path_to file/trappe-ua.xml'
  
          * Example dict for standard FF names : {'ETH' : 'oplsaa', 'OCT': 'trappe-ua'}
  
          * Example str for standard FF names: 'trappe-ua'
  
          * Example of a mixed dict with both : {'ETH' : 'oplsaa', 'OCT': 'path_to_file/'trappe-ua.xml'}
  
      detect_forcefield_style : bool, default = True
          If True, format NAMD/GOMC/LAMMPS parameters based on the contents of
          the parmed structure_box_0 and structure_box_1
      gomc_fix_bonds_angles : list, default = None
          When list of residues is provided, the selected residues will have
          their bonds and angles fixed and will ignore the relative bond energies and
          related angle energies in the GOMC engine. Note that GOMC
          does not sample bond stretching. This is specifically
          for the GOMC engine and it changes the residue's bond constants (Kbs)
          and angle constants (Kthetas) values to 999999999999 in the
          FF file (i.e., the .inp file).
          If the residues are listed in either the gomc_fix_angles or the gomc_fix_bonds_angles
          lists, the angles will be fixed for that residue.
          If the residues are listed in either the gomc_fix_bonds or the gomc_fix_bonds_angles
          lists, the bonds will be fixed for that residue.
          NOTE if this option is utilized it may cause issues if using the FF file in NAMD.
      gomc_fix_bonds : list, default = None
          When list of residues is provided, the selected residues will have their
          relative bond energies ignored in the GOMC engine. Note that GOMC
          does not sample bond stretching. This is specifically
          for the GOMC engine and it changes the residue's bond constants (Kbs)
          values to 999999999999 in the FF file (i.e., the .inp file).
          If the residues are listed in either the gomc_fix_bonds or the gomc_fix_bonds_angles
          lists, the relative bond energy will be ignored.
          NOTE if this option is utilized it may cause issues if using the FF file in NAMD.
      gomc_fix_angles : list, default = None
          When list of residues is provided, the selected residues will have
          their angles fixed and will ignore the related angle energies in the GOMC engine.
          This is specifically for the GOMC engine and it changes the residue's angle
          constants (Kthetas) values to 999999999999 in the FF file (i.e., the .inp file),
          which fixes the angles and ignores related angle energy.
          If the residues are listed in either the gomc_fix_angles or the gomc_fix_bonds_angles
          lists, the angles will be fixed and the related angle energy will be ignored
          for that residue.
          NOTE if this option is utilized it may cause issues if using the FF file in NAMD.
      bead_to_atom_name_dict : dict, optional, default =None
          For all atom names/elements/beads with 2 or less digits, this converts
          the atom name in the GOMC psf and pdb files to a unique atom name,
          provided they do not exceed 3844 atoms (62^2) of the same name/element/bead
          per residue. For all atom names/elements/beads with 3 digits, this converts
          the atom name in the GOMC psf and pdb files to a unique atom name,
          provided they do not exceed 62 of the same name/element pre residue.
  
          * Example dictionary: {'_CH3':'C', '_CH2':'C', '_CH':'C', '_HC':'C'}
  
          * Example name structure: {atom_type: first_part_pf atom name_without_numbering}
  
      fix_residue : list  or None, default = None
          Changes occcur in the pdb file only.
          When residues are listed here, all the atoms in the residue are
          fixed and can not move via setting the Beta values in the PDB
          file to 1.00.
          If neither fix_residue or fix_residue_in_box lists a
          residue or both equal None, then the Beta values for all the atoms
          in the residue are free to move in the simulation and Beta values
          in the PDB file is set to 0.00
      fix_residue_in_box : list  or None, default = None
          Changes occcur in the pdb file only.
          When residues are listed here, all the atoms in the residue
          can move within the box but cannot be transferred between boxes
          via setting the Beta values in the PDB file to 2.00.
          If neither fix_residue or fix_residue_in_box lists a
          residue or both equal None, then the Beta values for all the atoms
          in the residue are free to move in the simulation and Beta values
          in the PDB file is set to 0.00
      ff_filename : str, default =None
          If a string, it will write the  force field files that work in
          GOMC and NAMD structures.
      reorder_res_in_pdb_psf : bool, default =False
          If False, the order of of the atoms in the pdb file is kept in
          its original order, as in the Compound sent to the writer.
          If True, the order of the atoms is reordered based on their
          residue names in the 'residues' list that was entered.
      box_0 : Box
          The Box class that contains the attributes Lx, Ly, Lz for the length
          of the box 0 (units in nanometers (nm)). It also contains the xy, xz, and yz Tilt factors
          needed to displace an orthogonal box's xy face to its
          parallelepiped structure for box 0.
      box_1 : Box
          The Box class that contains the attributes Lx, Ly, Lz for the length
          of the box 1 (units in nanometers (nm)). It also contains the xy, xz, and yz Tilt factors
          needed to displace an orthogonal box's xy face to its
          parallelepiped structure for box 0.
      box_0_vectors : numpy.ndarray, [[float, float, float], [float, float, float], [float, float, float]]
          Three (3) sets vectors for box 0 each with 3 float values, which represent
          the vectors for the Charmm-style systems (units in Angstroms (Ang))
      box_1_vectors : numpy.ndarray, [[float, float, float], [float, float, float], [float, float, float]]
          Three (3) sets vectors for box 1 each with 3 float values, which represent
          the vectors for the Charmm-style systems (units in Angstroms (Ang))
      structure_box_0_ff : parmed.structure.Structure
          The box 0 structure (structure_box_0) after all the provided
          force fields are applied.
      structure_box_1_ff : parmed.structure.Structure
          The box 0 structure (structure_box_0) after all the provided
          force fields are applied. This only exists if the box 1 structure
          (structure_box_1) is provided.
      coulomb14scalar_dict_box_0 : dict
          The residue/moleclues (key) of box 0 and their corresponding
          coulombic 1-4 scalers (value).  Note: NAMD and GOMC can only have one (1)
          value for the coulombic 1-4 scalers, as they both only accept a
          single value in the NAMD and GOMC control files.
      coulomb14scalar_dict_box_1 : dict
          The residue/moleclues  (key) of box 1 and their corresponding
          coulombic 1-4 scalers (value).  Note: NAMD and GOMC can only have one (1)
          value for the coulombic 1-4 scalers, as they both only accept a
          single value in the NAMD and GOMC control files.
          This only exists if the box 1 structure (structure_box_1) is provided.
      LJ14scalar_dict_box_0 : dict
          The residue/moleclues (key) of box 0 and their corresponding
          Lennard-Jones (LJ) 1-4 scalers (value).  Note: NAMD and GOMC can have
          multiple values for the LJ 1-4 scalers, since they are provided as an
          individual input for each atom type in the force field (.inp) file.
      LJ14scalar_dict_box_1 : dict
          The residue/moleclues (key) of box 1 and their corresponding
          Lennard-Jones (LJ) 1-4 scalers (value).  Note: NAMD and GOMC can have
          multiple values for the LJ 1-4 scalers, since they are provided as an
          individual input for each atom type in the force field (.inp) file.
          This only exists if the box 1 structure (structure_box_1) is provided.
      residues_applied_list_box_0 : list
          The residues in box 0 that were found and had the force fields applied to them.
      residues_applied_list_box_1 : list
          The residues in box 1 that were found and had the force fields applied to them.
          This only exists if the box 1 structure (structure_box_1) is provided.
      boxes_for_simulation : int, [0, 1]
          The number of boxes used when writing the Charmm object and force fielding
          the system. If only box 0 is provided, the value is 0. If box 0 and box 1
          are provided, the value is 1.
      epsilon_dict : dict {str: float or int}
          The uniquely numbered atom type (key) and it's non-bonded epsilon
          coefficient (value).
      sigma_dict : dict {str: float or int}
          The uniquely numbered atom type (key) and it's non-bonded sigma
          coefficient (value).
      LJ_1_4_dict : dict {str: float or int}
          The uniquely numbered atom type (key) and it's non-bonded 1-4
          Lennard-Jones, LJ, scaling factor (value).
      coul_1_4 : float or int
          The non-bonded 1-4 coulombic scaling factor, which is the
          same for all the residues/molecules, regardless if
          differenct force fields are utilized.  Note: if
          1-4 coulombic scaling factor is not the same for all
          molecules the Charmm object will fail with an error.
      combined_1_4_coul_dict_per_residue : dict, {str: float or int}
          The residue name/molecule (key) and it's non-bonded 1-4 coulombic
          scaling factor (value).
      forcefield_dict : dict
          The uniquely numbered atom type (key) with it's corresponding
          foyer atom typing and residue name.  The residue name is added
          to provide a distinction between other residues with the same
          atom types.  This allows the CHARMM force field to fix the
          bonds and angles specific residues without effecting other
          residues with the same atom types.
      all_individual_atom_names_list : list
          A list of all the atom names for the combined structures
          (box 0 and box 1 (if supplied)), in order.
      all_residue_names_List : list
          A list of all the residue names for the combined structures
          (box 0 and box 1 (if supplied)), in order.
      max_residue_no : int
          The maximum number that the residue number will count to
          before restarting the counting back to 1, which is predetermined
          by the PDB format. This is a constant, which equals 9999
      max_resname_char : int
          The maximum number of characters allowed in the residue name,
          which is predetermined by the PDB format. This is a constant,
          which equals 4.
      all_res_unique_atom_name_dict : dict, {str : [str, ..., str]}
          A dictionary that provides the residue names (keys) and a list
          of the unique atom names in the residue (value), for the
          combined structures (box 0 and box 1 (if supplied)).
  
      Notes
      -----
      Impropers, Urey-Bradleys, and NBFIX are not currenly supported.
      Currently the NBFIX is disabled as since only the OPLS and TRAPPE force fields are supported.
      OPLS and CHARMM forcefield styles are supported (without impropers),
      AMBER forcefield styles are NOT supported.
  
      The atom typing is currently provided via a base 52 numbering (capital and lowercase lettering).
      This base 52 numbering allows for (52)^4 unique atom types.
  
      Unique atom names are provided if the system do not exceed 3844 atoms (62^2) of the same
      name/bead per residue (base 62 numbering). For all atom names/elements with 3 or less digits,
      this converts the atom name in the GOMC psf and pdb files to a unique atom name,
      provided they do not exceed 62 of the same name/element pre residue.
  
      Generating an empty box (i.e., pdb and psf files):
      Single Box system: Enter residues = [], but the accompanying structure (structure_box_0)
      must be an empty mb.Box. However, when doing this, the forcefield_selection
      must be supplied, or it will provide an error
      (i.e., forcefield_selection can not be equal to None).
      Dual Box System: Enter an empty mb.Box structure for either structure_box_0 or
      structure_box_1.
  
      In this current FF/psf/pdb writer, a residue type is essentially a molecule type.
      Therefore, it can only correctly write systems where every bead/atom in the molecule
      has the same residue name, and the residue name is specific to that molecule type.
      For example: a protein molecule with many residue names is not currently supported,
      but is planned to be supported in the future.
      """
  
      def __init__(
          self,
          structure_box_0,
          filename_box_0,
          structure_box_1=None,
          filename_box_1=None,
          non_bonded_type="LJ",
          forcefield_selection=None,
          residues=None,
          detect_forcefield_style=True,
          gomc_fix_bonds_angles=None,
          gomc_fix_bonds=None,
          gomc_fix_angles=None,
          bead_to_atom_name_dict=None,
          fix_residue=None,
          fix_residue_in_box=None,
          ff_filename=None,
          reorder_res_in_pdb_psf=False,
      ):
  
          # set all input variables to the class
          self.structure_box_0 = structure_box_0
          self.filename_box_0 = filename_box_0
          self.structure_box_1 = structure_box_1
          self.filename_box_1 = filename_box_1
          self.non_bonded_type = non_bonded_type
          self.forcefield_selection = forcefield_selection
          self.residues = residues
          self.detect_forcefield_style = detect_forcefield_style
          self.gomc_fix_bonds_angles = gomc_fix_bonds_angles
          self.gomc_fix_bonds = gomc_fix_bonds
          self.gomc_fix_angles = gomc_fix_angles
          self.bead_to_atom_name_dict = bead_to_atom_name_dict
          self.fix_residue = fix_residue
          self.fix_residue_in_box = fix_residue_in_box
          self.ff_filename = ff_filename
          self.reorder_res_in_pdb_psf = reorder_res_in_pdb_psf
  
          # value to check for errors, with  self.input_error = True or False. Set to False initally
          self.input_error = False
  
          if not isinstance(self.structure_box_0, (Compound, Box)):
              self.input_error = True
              print_error_message = (
                  "ERROR: The structure_box_0 expected to be of type: "
                  "{} or {}, received: {}".format(
                      type(Compound()),
                      type(Box(lengths=[1, 1, 1])),
                      type(structure_box_0),
                  )
              )
              raise TypeError(print_error_message)
  
          if self.structure_box_1 is not None and not isinstance(
              self.structure_box_1, (Compound, Box)
          ):
              self.input_error = True
              print_error_message = (
                  "ERROR: The structure_box_1 expected to be of type: "
                  "{} or {}, received: {}".format(
                      type(Compound()),
                      type(Box(lengths=[1, 1, 1])),
                      type(structure_box_1),
                  )
              )
              raise TypeError(print_error_message)
  
          if isinstance(self.structure_box_0, Box) and isinstance(
              self.structure_box_1, Box
          ):
              self.input_error = True
              print_error_message = (
                  "ERROR: Both structure_box_0 and structure_box_0 are empty Boxes {}. "
                  "At least 1 structure must be an mbuild compound {} with 1 "
                  "or more atoms in it".format(
                      type(Box(lengths=[1, 1, 1])), type(Compound())
                  )
              )
              raise TypeError(print_error_message)
  
          if self.structure_box_1 is None and not isinstance(
              self.structure_box_0, Compound
          ):
              self.input_error = True
              print_error_message = (
                  "ERROR: Only 1 structure is provided and it can not be an empty mbuild Box {}. "
                  "it must be an mbuild compound {} with at least 1 "
                  "or more atoms in it.".format(
                      type(Box(lengths=[1, 1, 1])), type(Compound())
                  )
              )
              raise TypeError(print_error_message)
  
          if not isinstance(self.residues, list):
              self.input_error = True
              print_error_message = "ERROR: Please enter the residues list (residues) in a list format."
              raise TypeError(print_error_message)
  
          if isinstance(self.residues, list):
              for each_residue in self.residues:
                  if not isinstance(each_residue, str):
                      self.input_error = True
                      print_error_message = "ERROR: Please enter a residues list (residues) with only string values."
                      raise TypeError(print_error_message)
  
          if self.residues is None:
              self.input_error = True
              print_error_message = (
                  "ERROR: Please enter the residues list (residues)"
              )
              raise TypeError(print_error_message)
          if not isinstance(self.filename_box_0, str):
              self.input_error = True
              print_error_message = (
                  "ERROR: Please enter the filename_box_0 as a string."
              )
              raise TypeError(print_error_message)
  
          unique_residue_test_name_list = []
          for res_m in range(0, len(self.residues)):
              if self.residues[res_m] not in unique_residue_test_name_list:
                  unique_residue_test_name_list.append(self.residues[res_m])
          if len(unique_residue_test_name_list) != len(self.residues):
              self.input_error = True
              print_error_message = "ERROR: Please enter the residues list (residues) that has only unique residue names."
              raise ValueError(print_error_message)
  
          if self.filename_box_1 is not None and not isinstance(
              self.filename_box_1, str
          ):
              self.input_error = True
              print_error_message = (
                  "ERROR: Please enter the filename_box_1 as a string."
              )
              raise TypeError(print_error_message)
  
          if self.ff_filename is not None:
              if not isinstance(self.ff_filename, str):
                  self.input_error = True
                  print_error_message = "ERROR: Please enter GOMC force field name (ff_filename) as a string."
                  raise TypeError(print_error_message)
              if isinstance(self.ff_filename, str):
                  extension_ff_name = os.path.splitext(self.ff_filename)[-1]
                  if extension_ff_name == "":
                      self.ff_filename = self.ff_filename + ".inp"
                  elif extension_ff_name == ".inp":
                      self.ff_filename = self.ff_filename + ""
                  elif extension_ff_name != ".inp":
                      self.input_error = True
                      print_error_message = (
                          "ERROR: Please enter GOMC force field name without an "
                          "extention or the .inp extension."
                      )
                      raise ValueError(print_error_message)
  
          if self.forcefield_selection is not None:
              print(
                  "write_gomcdata: forcefield_selection = "
                  + str(self.forcefield_selection)
                  + ", "
                  + "residues = "
                  + str(self.residues)
              )
              if not isinstance(self.forcefield_selection, (dict, str)):
                  self.input_error = True
                  print_error_message = (
                      "ERROR: The force field selection (forcefield_selection) "
                      "is not a string or a dictionary with all the residues specified "
                      'to a force field. -> String Ex: "path/trappe-ua.xml" or Ex: "trappe-ua" '
                      "Otherise provided a dictionary with all the residues specified "
                      "to a force field "
                      '->Dictionary Ex: {"Water" : "oplsaa", "OCT": "path/trappe-ua.xml"}, '
                      "Note: the file path must be specified the force field file if "
                      "a standard foyer force field is not used."
                  )
                  raise TypeError(print_error_message)
  
              if isinstance(self.forcefield_selection, str):
                  ff_name = self.forcefield_selection
                  self.forcefield_selection = {}
                  for i in range(0, len(self.residues)):
                      self.forcefield_selection.update(
                          {self.residues[i]: ff_name}
                      )
                  print(
                      "FF forcefield_selection = "
                      + str(self.forcefield_selection)
                  )
  
          elif self.forcefield_selection is None:
              self.input_error = True
              print_error_message = "ERROR: Please enter the forcefield_selection as it was not provided."
              raise TypeError(print_error_message)
  
          if self.residues is not None and not isinstance(self.residues, list):
              self.input_error = True
              print_error_message = (
                  "ERROR:  Please enter the residues (residues) in a list format"
              )
              raise TypeError(print_error_message)
  
          _check_fixed_bonds_angles_lists(
              self.gomc_fix_bonds_angles, "gomc_fix_bonds_angles", self.residues
          )
          _check_fixed_bonds_angles_lists(
              self.gomc_fix_bonds, "gomc_fix_bonds", self.residues
          )
          _check_fixed_bonds_angles_lists(
              self.gomc_fix_angles, "gomc_fix_angles", self.residues
          )
  
          if self.fix_residue is not None and not isinstance(
              self.fix_residue, list
          ):
              self.input_error = True
              print_error_message = (
                  "ERROR: Please enter the fix_residue in a list format"
              )
              raise TypeError(print_error_message)
  
          if isinstance(self.fix_residue, list):
              for fix_residue_q in self.fix_residue:
                  if fix_residue_q not in self.residues:
                      self.input_error = True
                      print_error_message = (
                          "Error: Please ensure that all the residue names in the fix_residue "
                          "list are also in the residues list."
                      )
                      raise ValueError(print_error_message)
  
          if self.fix_residue_in_box is not None and not isinstance(
              self.fix_residue_in_box, list
          ):
              self.input_error = True
              print_error_message = (
                  "ERROR: Please enter the fix_residue_in_box in a list format."
              )
              raise TypeError(print_error_message)
  
          if isinstance(self.fix_residue_in_box, list):
              for fix_residue_in_box_q in self.fix_residue_in_box:
                  if fix_residue_in_box_q not in self.residues:
                      self.input_error = True
                      print_error_message = (
                          "Error: Please ensure that all the residue names in the "
                          "fix_residue_in_box list are also in the residues list."
                      )
                      raise ValueError(print_error_message)
  
          if self.bead_to_atom_name_dict is not None and not isinstance(
              self.bead_to_atom_name_dict, dict
          ):
              self.input_error = True
              print_error_message = (
                  "ERROR: Please enter the a bead type to atom in the dictionary "
                  "(bead_to_atom_name_dict) so GOMC can properly evaluate the unique atom names"
              )
              raise TypeError(print_error_message)
  
          if isinstance(self.bead_to_atom_name_dict, dict):
              dict_list = []
              for key in self.bead_to_atom_name_dict.keys():
                  dict_list.append(key)
  
              for dict_lis_i in dict_list:
                  if not isinstance(dict_lis_i, str) or not isinstance(
                      self.bead_to_atom_name_dict[dict_lis_i], str
                  ):
                      print_error_message = "ERROR: Please enter the bead_to_atom_name_dict with only string inputs."
                      raise TypeError(print_error_message)
  
          print("******************************")
          print("")
  
          self.sub_1_for_base_52 = 1
  
          if self.structure_box_1:
              self.boxes_for_simulation = 2
          else:
              self.boxes_for_simulation = 1
  
          # write the Force fields
          self.combined_1_4_lj_dict_per_residue = {}
          self.combined_1_4_coul_dict_per_residue = {}
  
          if self.structure_box_1:
  
              print(
                  "GOMC FF writing each residues FF as a group for structure_box_0"
              )
              [
                  self.structure_box_0_ff,
                  self.coulomb14scalar_dict_box_0,
                  self.LJ14scalar_dict_box_0,
                  self.residues_applied_list_box_0,
              ] = specific_ff_to_residue(
                  self.structure_box_0,
                  forcefield_selection=self.forcefield_selection,
                  residues=self.residues,
                  reorder_res_in_pdb_psf=self.reorder_res_in_pdb_psf,
                  boxes_for_simulation=self.boxes_for_simulation,
              )
  
              print(
                  "GOMC FF writing each residues FF as a group for  structure_box_1"
              )
              [
                  self.structure_box_1_ff,
                  self.coulomb14scalar_dict_box_1,
                  self.LJ14scalar_dict_box_1,
                  self.residues_applied_list_box_1,
              ] = specific_ff_to_residue(
                  self.structure_box_1,
                  forcefield_selection=self.forcefield_selection,
                  residues=self.residues,
                  reorder_res_in_pdb_psf=self.reorder_res_in_pdb_psf,
                  boxes_for_simulation=self.boxes_for_simulation,
              )
  
              self.structure_box_0_and_1_ff = (
                  self.structure_box_0_ff + self.structure_box_1_ff
              )
              self.combined_1_4_lj_dict_per_residue.update(
                  self.LJ14scalar_dict_box_0
              )
              self.combined_1_4_lj_dict_per_residue.update(
                  self.LJ14scalar_dict_box_1
              )
              self.combined_1_4_coul_dict_per_residue.update(
                  self.coulomb14scalar_dict_box_0
              )
              self.combined_1_4_coul_dict_per_residue.update(
                  self.coulomb14scalar_dict_box_1
              )
  
              self.residues_applied_list_box_0_and_1 = (
                  self.residues_applied_list_box_0
              )
              for res_iter in self.residues_applied_list_box_1:
                  if res_iter not in self.residues_applied_list_box_0:
                      self.residues_applied_list_box_0_and_1.append(res_iter)
  
              for res_iter_0_1 in self.residues_applied_list_box_0_and_1:
                  if res_iter_0_1 not in self.residues:
                      self.input_error = True
                      print_error_message = "ERROR: All the residues were not used from the forcefield_selection "
                      "string or dictionary.  There may be residues below other specified "
                      "residues in the mbuild.Compound hierarchy.  If so, the residues "
                      "acquire the residue's force fields, which is at the top of the "
                      "hierarchy.  Alternatively, residues that are not in the structure "
                      r"may have been specified."
                      raise ValueError(print_error_message)
  
              for res_iter_0_1 in self.residues:
                  if res_iter_0_1 not in self.residues_applied_list_box_0_and_1:
                      self.input_error = True
                      print_error_message = (
                          "ERROR: All the residues were not used from the forcefield_selection "
                          "string or dictionary.  There may be residues below other specified "
                          "residues in the mbuild.Compound hierarchy.  If so, the residues "
                          "acquire the residue's force fields, which is at the top of the "
                          "hierarchy.  Alternatively, residues that are not in the structure "
                          "may have been specified."
                      )
                      raise ValueError(print_error_message)
  
              total_charge = sum(
                  [atom.charge for atom in self.structure_box_0_ff]
              )
              if round(total_charge, 4) != 0.0:
                  warn(
                      "System is not charge neutral for structure_box_0. "
                      "Total charge is {}.".format(total_charge)
                  )
  
              total_charge = sum(
                  [atom.charge for atom in self.structure_box_1_ff]
              )
              if round(total_charge, 4) != 0.0:
                  warn(
                      "System is not charge neutral for structure_box_1. "
                      "Total charge is {}.".format(total_charge)
                  )
  
              total_charge = sum(
                  [atom.charge for atom in self.structure_box_0_and_1_ff]
              )
              if round(total_charge, 4) != 0.0:
                  warn(
                      "System is not charge neutral for structure_0_and_1. "
                      "Total charge is {}.".format(total_charge)
                  )
  
          else:
              print(
                  "GOMC FF writing each residues FF as a group for structure_box_0"
              )
              [
                  self.structure_box_0_ff,
                  self.coulomb14scalar_dict_box_0,
                  self.LJ14scalar_dict_box_0,
                  self.residues_applied_list_box_0,
              ] = specific_ff_to_residue(
                  self.structure_box_0,
                  forcefield_selection=self.forcefield_selection,
                  residues=self.residues,
                  reorder_res_in_pdb_psf=self.reorder_res_in_pdb_psf,
                  boxes_for_simulation=self.boxes_for_simulation,
              )
  
              self.combined_1_4_lj_dict_per_residue.update(
                  self.LJ14scalar_dict_box_0
              )
              self.combined_1_4_coul_dict_per_residue.update(
                  self.coulomb14scalar_dict_box_0
              )
  
              for res_iter_0 in self.residues_applied_list_box_0:
                  if res_iter_0 not in self.residues:
                      self.input_error = True
                      print_error_message = (
                          "ERROR: All the residues were not used from the forcefield_selection "
                          "string or dictionary.  There may be residues below other specified "
                          "residues in the mbuild.Compound hierarchy.  If so, the residues "
                          "acquire the residue's force fields, which is at the top of the "
                          "hierarchy.  Alternatively, residues that are not in the structure "
                          "may have been specified."
                      )
                      raise ValueError(print_error_message)
  
              for res_iter_0 in self.residues:
                  if res_iter_0 not in self.residues_applied_list_box_0:
                      self.input_error = True
                      print_error_message = (
                          "ERROR: All the residues were not used from the forcefield_selection "
                          "string or dictionary.  There may be residues below other specified "
                          "residues in the mbuild.Compound hierarchy.  If so, the residues "
                          "acquire the residue's force fields, which is at the top of the "
                          "hierarchy.  Alternatively, residues that are not in the structure "
                          "may have been specified."
                      )
                      raise ValueError(print_error_message)
  
              total_charge = sum(
                  [atom.charge for atom in self.structure_box_0_ff]
              )
              if round(total_charge, 4) != 0.0:
                  warn(
                      "System is not charge neutral for structure_box_0. "
                      "Total charge is {}.".format(total_charge)
                  )
  
          print(
              "forcefield type from compound = " + str(self.forcefield_selection)
          )
          print(
              "coulomb14scale from compound = "
              + str(self.combined_1_4_coul_dict_per_residue)
          )
          print(
              "lj14scale from compound = "
              + str(self.combined_1_4_lj_dict_per_residue)
          )
  
          # lock the atom_style and unit_style for GOMC. Can be inserted into variables
          # once more functionality is built in
          self.atom_style = "full"
          self.unit_style = "real"
          # functional form type default.  Can be inserted into variables once more functionality is built in
          use_rb_torsions = True
          use_dihedrals = False
          use_urey_bradleys = False
  
          # Convert coordinates to real or other units (real only current option)
          if self.unit_style == "real":
              self.sigma_conversion_factor = 1
              self.epsilon_conversion_factor = 1
              self.mass_conversion_factor = 1
  
          if self.structure_box_1:
              self.types = np.array(
                  [
                      atom.type + "_" + str(atom.residue.name)
                      for atom in self.structure_box_0_and_1_ff.atoms
                  ]
              )
  
          else:
              self.types = np.array(
                  [
                      atom.type + "_" + str(atom.residue.name)
                      for atom in self.structure_box_0_ff.atoms
                  ]
              )
  
          self.unique_types = list(set(self.types))
          self.unique_types.sort(key=natural_sort)
  
          if self.structure_box_1:
              self.masses = (
                  np.array(
                      [atom.mass for atom in self.structure_box_0_and_1_ff.atoms]
                  )
                  / self.mass_conversion_factor
              )
              self.mass_dict = dict(
                  [
                      (self.unique_types.index(atom_type) + 1, mass)
                      for atom_type, mass in zip(self.types, self.masses)
                  ]
              )
  
          else:
              self.masses = (
                  np.array([atom.mass for atom in self.structure_box_0_ff.atoms])
                  / self.mass_conversion_factor
              )
              self.mass_dict = dict(
                  [
                      (self.unique_types.index(atom_type) + 1, mass)
                      for atom_type, mass in zip(self.types, self.masses)
                  ]
              )
  
          # added an index so the atom types can be converted to numbers as the type name is to long for insertion into
          # the pdb and psf files
          self.atom_types_to_index_value_dict = dict(
              [
                  (
                      self.unique_types[self.unique_types.index(atom_type)],
                      self.unique_types.index(atom_type),
                  )
                  for atom_type, mass in zip(self.types, self.masses)
              ]
          )
  
          # normalize by sigma
          self.box_0 = Box(
              lengths=np.array(
                  [
                      (0.1 * val) / self.sigma_conversion_factor
                      for val in self.structure_box_0_ff.box[0:3]
                  ]
              ),
              angles=self.structure_box_0_ff.box[3:6],
          )
  
          # create box 0 vector list and convert from nm to Ang and round to 6 decimals.
          # note mbuild standard lengths are in nm, so round to 6+1 = 7 then mutlipy by 10
          box_0_lengths_ang = (
              self.box_0.lengths[0] * 10,
              self.box_0.lengths[1] * 10,
              self.box_0.lengths[2] * 10,
          )
          self.box_0_vectors = _lengths_angles_to_vectors(
              box_0_lengths_ang, self.box_0.angles, precision=6
          )
  
          # Internally use nm
          if self.structure_box_1:
              self.box_1 = Box(
                  lengths=np.array(
                      [
                          (0.1 * val) / self.sigma_conversion_factor
                          for val in self.structure_box_1_ff.box[0:3]
                      ]
                  ),
                  angles=self.structure_box_1_ff.box[3:6],
              )
  
              # create box 1 vector list and convert from nm to Ang and round to 6 decimals.
              # note mbuild standard lengths are in nm, so round to 6+1 = 7 then mutlipy by 10
              box_1_lengths_ang = (
                  self.box_1.lengths[0] * 10,
                  self.box_1.lengths[1] * 10,
                  self.box_1.lengths[2] * 10,
              )
              self.box_1_vectors = _lengths_angles_to_vectors(
                  box_1_lengths_ang, self.box_1.angles, precision=6
              )
  
          # if self.structure_box_1 != None:
          if self.structure_box_1:
              self.structure_selection = self.structure_box_0_and_1_ff
          else:
              self.structure_selection = self.structure_box_0_ff
  
          # Non-Bonded forces
          epsilons = (
              np.array([atom.epsilon for atom in self.structure_selection.atoms])
              / self.epsilon_conversion_factor
          )
          sigmas = (
              np.array([atom.sigma for atom in self.structure_selection.atoms])
              / self.sigma_conversion_factor
          )
          forcefields = [
              atom.type + "_" + atom.residue.name
              for atom in self.structure_selection.atoms
          ]
          residues_all_list = [
              atom.residue.name for atom in self.structure_selection.atoms
          ]
  
          self.epsilon_dict = dict(
              [
                  (self.unique_types.index(atom_type), epsilon)
                  for atom_type, epsilon in zip(self.types, epsilons)
              ]
          )
          self.sigma_dict = dict(
              [
                  (self.unique_types.index(atom_type), sigma)
                  for atom_type, sigma in zip(self.types, sigmas)
              ]
          )
          self.LJ_1_4_dict = dict(
              [
                  (
                      self.unique_types.index(atom_type),
                      self.combined_1_4_lj_dict_per_residue[residues_all_list],
                  )
                  for atom_type, residues_all_list in zip(
                      self.types, residues_all_list
                  )
              ]
          )
          self.forcefield_dict = dict(
              [
                  (self.unique_types.index(atom_type), forcefield)
                  for atom_type, forcefield in zip(self.types, forcefields)
              ]
          )
  
          # ensure all 1,4-coulombic scaling factors are the same
          coul_1_4_List = []
          for p in self.combined_1_4_coul_dict_per_residue.values():
              coul_1_4_List.append(p)
          coul_1_4_set = set(coul_1_4_List)
          if len(coul_1_4_set) > 1:
              self.input_error = True
              print_error_message = (
                  "ERROR: There are multiple 1,4-coulombic scaling factors "
                  "GOMC will only accept a singular input for the 1,4-coulombic "
                  "scaling factors."
              )
              raise ValueError(print_error_message)
          else:
              self.coul_1_4 = list(coul_1_4_set)[0]
  
          # get all the unique atom name to check for the MEMC move in the gomc_conf_writer
          self.all_individual_atom_names_list = []
          self.all_residue_names_List = []
          if self.structure_box_1:
              list_of_structures = [
                  self.structure_box_0_ff,
                  self.structure_box_1_ff,
              ]
              stuct_only = [self.structure_box_0_ff, self.structure_box_1_ff]
          else:
              list_of_structures = [self.structure_box_0_ff]
              stuct_only = [self.structure_box_0_ff]
  
          for q_i in range(0, len(list_of_structures)):
              stuct_only_iteration = stuct_only[q_i]
  
              # caluculate the atom name and unique atom names
              residue_data_list = []
              residue_names_list = []
              for k, atom in enumerate(stuct_only_iteration.atoms):
                  residue_data_list.append(str(atom.residue))
                  residue_names_list.append(atom.residue.name)
  
              unique_residue_data_dict = {}
              unique_residue_data_list = []
              residue_data_name_list = []
  
              for m, residue in enumerate(stuct_only_iteration.residues):
                  unique_residue_data_list.append(
                      str(stuct_only_iteration.residues[m])
                  )
                  unique_residue_data_dict.update(
                      {unique_residue_data_list[m]: m + 1}
                  )
                  residue_data_name_list.append(
                      stuct_only_iteration.residues[m].name
                  )
  
              self.max_residue_no = 9999
              self.max_resname_char = 4
  
              res_no_chain_iter_corrected = []
              residue_id_list = []
              residue_id_adder_fixed_struct_wo_bonds = (
                  0  # for example zeolite used as fixed atoms wo bonds
              )
              for f, PSF_atom_iteration_0 in enumerate(
                  stuct_only_iteration.atoms
              ):
  
                  if f > 0:
                      if (
                          PSF_atom_iteration_0.residue.chain
                          == previous_residue_chain
                          and len(PSF_atom_iteration_0.bonds) == 0
                      ):
                          residue_id_adder_fixed_struct_wo_bonds += 1
  
                  previous_residue_chain = PSF_atom_iteration_0.residue.chain
  
                  residue_id_int = int(
                      unique_residue_data_dict[residue_data_list[f]]
                      + residue_id_adder_fixed_struct_wo_bonds
                  )
  
                  res_id_adder = int(
                      (residue_id_int % self.max_residue_no) % self.max_residue_no
                  )
                  if int(res_id_adder) == 0:
                      res_no_iteration_corrected = int(self.max_residue_no)
                  else:
                      res_no_iteration_corrected = res_id_adder
  
                  res_no_chain_iter_corrected.append(res_no_iteration_corrected)
                  residue_id_list.append(residue_id_int)
  
              # This converts the atom name in the GOMC psf and pdb files to unique atom names
              [
                  unique_individual_atom_names_dict_iter,
                  individual_atom_names_list_iter,
                  missing_bead_to_atom_name_iter,
              ] = unique_atom_naming(
                  stuct_only_iteration,
                  residue_id_list,
                  residue_names_list,
                  bead_to_atom_name_dict=self.bead_to_atom_name_dict,
              )
  
              if q_i == 0:
                  self.all_individual_atom_names_list = (
                      individual_atom_names_list_iter
                  )
  
                  self.all_residue_names_List = residue_names_list
              else:
  
                  self.all_individual_atom_names_list = (
                      self.all_individual_atom_names_list
                      + individual_atom_names_list_iter
                  )
  
                  self.all_residue_names_List = (
                      self.all_residue_names_List + residue_names_list
                  )
  
          # put the  self.all_individual_atom_names_list and self.all_residue_names_List in a list to match
          # the the atom name with a residue and find the unique matches
          if None in [
              unique_individual_atom_names_dict_iter,
              individual_atom_names_list_iter,
              missing_bead_to_atom_name_iter,
          ]:
              self.input_error = True
              print_error_message = (
                  "ERROR: The unique_atom_naming function failed while "
                  "running the charmm_writer function. Ensure the proper inputs are "
                  "in the bead_to_atom_name_dict."
              )
              raise ValueError(print_error_message)
  
          else:
              self.all_res_unique_atom_name_dict = {}
              for res_i in range(0, len(self.all_individual_atom_names_list)):
                  try:
                      if (
                          self.all_individual_atom_names_list[res_i]
                          not in self.all_res_unique_atom_name_dict[
                              self.all_residue_names_List[res_i]
                          ]
                      ):
                          self.all_res_unique_atom_name_dict.setdefault(
                              self.all_residue_names_List[res_i], []
                          ).append(self.all_individual_atom_names_list[res_i])
                  except:
                      self.all_res_unique_atom_name_dict.setdefault(
                          self.all_residue_names_List[res_i], []
                      ).append(self.all_individual_atom_names_list[res_i])
  
          print(
              "all_res_unique_atom_name_dict = {}".format(
                  self.all_res_unique_atom_name_dict
              )
          )
  
      def write_inp(self):
          """This write_inp function writes the Charmm style parameter (force field) file, which can be utilized
          in the GOMC and NAMD engines."""
          print("******************************")
          print("")
          print(
              "The charmm force field file writer (the write_inp function) is running"
          )
  
          if self.ff_filename is None:
              self.input_error = True
              print_error_message = (
                  "ERROR: The force field file name was not specified and in the "
                  "Charmm object. "
                  "Therefore, the force field file (.inp) can not be written. "
                  "Please use the force field file name when building the Charmm object, "
                  "then use the write_inp function."
              )
              raise TypeError(print_error_message)
          else:
  
              print("******************************")
              print("")
              print(
                  "The charmm force field file writer (the write_inp function) is running"
              )
              print("******************************")
              print("")
              print("writing the GOMC force field file ")
              date_time = datetime.datetime.today()
  
              unique_residue_data_dict = {}
              unique_residue_data_list = []
              residue_data_name_list = []
  
              # if self.structure_box_1 != None:
              if self.structure_box_1:
                  residue_iterate = 0
                  for m, residue in enumerate(self.structure_box_0_ff.residues):
                      residue_iterate = residue_iterate + 1
                      unique_residue_data_list.append(
                          str(self.structure_box_0_ff.residues[m])
                      )
                      unique_residue_data_dict.update(
                          {unique_residue_data_list[m]: m + 1}
                      )
                      residue_data_name_list.append(
                          self.structure_box_0_ff.residues[m].name
                      )
  
                  for m, residue in enumerate(self.structure_box_1_ff.residues):
                      unique_residue_data_list.append(
                          str(self.structure_box_1_ff.residues[m])
                      )
                      unique_residue_data_dict.update(
                          {unique_residue_data_list[m]: m + 1 + residue_iterate}
                      )
                      residue_data_name_list.append(
                          self.structure_box_1_ff.residues[m].name
                      )
  
              else:
                  for m, residue in enumerate(self.structure_box_0_ff.residues):
                      unique_residue_data_list.append(
                          str(self.structure_box_0_ff.residues[m])
                      )
                      unique_residue_data_dict.update(
                          {unique_residue_data_list[m]: m + 1}
                      )
                      residue_data_name_list.append(
                          self.structure_box_0_ff.residues[m].name
                      )
  
              for n in range(0, len(residue_data_name_list)):
                  if residue_data_name_list[n] not in self.residues:
                      print(
                          "residue_data_name_list = "
                          + str(residue_data_name_list)
                      )
                      self.input_error = True
                      print_error_message = "ERROR: Please specifiy all residues (residues) in a list"
                      raise ValueError(print_error_message)
  
              with open(self.ff_filename, "w") as data:
                  # Syntax which can change based on the functional form
                  # Infer functional form based on the properties of the residue_id_list or structure_box_0_ff
                  if self.detect_forcefield_style:
                      # Check angles
                      if len(self.structure_selection.urey_bradleys) > 0:
                          print("Urey bradley terms detected")
                          data.write(
                              "! Urey bradley terms detected, will use angle_style self"
                          )
                          data.write(
                              "! POTENTIAL ERROR: GOMC does no support the Urey bradley terms"
                          )
                          use_urey_bradleys = True
                      else:
                          print(
                              "No urey bradley terms detected, will use angle_style harmonic"
                          )
                          use_urey_bradleys = False
  
                      # Check dihedrals
                      if len(self.structure_selection.rb_torsions) > 0:
                          print(
                              "will use CHARMM_torsions  =  K0 + K1 * (1 + Cos[n1*(t) - (d1)] ) + "
                              + "K2 * (1 + Cos[n2*(t) - (d2)] ) + K3 * (1 + Cos[n3*(t) - (d3)] ) + "
                              + "K4 * (1 + Cos[n4*(t) - (d4)] ) + K5 * (1 + Cos[n5*(t) - (d5)] ) "
                          )
                          use_rb_torsions = True
  
                      else:
                          use_rb_torsions = False
  
                      if len(self.structure_selection.dihedrals) > 0:
                          print(
                              "Charmm dihedrals detected, will use dihedral_style charmm"
                          )
                          # this will need tested with a standard charmm input format before releasing it
                          use_dihedrals = True
                          self.input_error = True
                          print_error_message = (
                              "ERROR: use_dihedrals = {} "
                              "Charmm dihedrals not yet supported.".format(
                                  use_dihedrals
                              )
                          )
                          raise ValueError(print_error_message)
                      else:
                          use_dihedrals = False
  
                  if use_rb_torsions and use_dihedrals:
                      warn(
                          "Multiple dihedral styles detected, check your Forcefield XML and structure_selection"
                      )
  
                  # Check impropers
                  for dihedral in self.structure_selection.dihedrals:
                      if dihedral.improper:
                          warn(
                              "Amber-style impropers are currently not supported"
                          )
  
                  bonds = [
                      [bond.atom1.idx + 1, bond.atom2.idx + 1]
                      for bond in self.structure_selection.bonds
                  ]
                  angles = [
                      [
                          angle.atom1.idx + 1,
                          angle.atom2.idx + 1,
                          angle.atom3.idx + 1,
                      ]
                      for angle in self.structure_selection.angles
                  ]
                  if use_rb_torsions:
                      dihedrals = [
                          [
                              dihedral.atom1.idx + 1,
                              dihedral.atom2.idx + 1,
                              dihedral.atom3.idx + 1,
                              dihedral.atom4.idx + 1,
                          ]
                          for dihedral in self.structure_selection.rb_torsions
                      ]
                  elif use_dihedrals:
                      dihedrals = [
                          [
                              dihedral.atom1.idx + 1,
                              dihedral.atom2.idx + 1,
                              dihedral.atom3.idx + 1,
                              dihedral.atom4.idx + 1,
                          ]
                          for dihedral in self.structure_selection.dihedrals
                      ]
                  else:
                      dihedrals = []
                  impropers = [
                      [
                          improper.atom1.idx + 1,
                          improper.atom2.idx + 1,
                          improper.atom3.idx + 1,
                          improper.atom4.idx + 1,
                      ]
                      for improper in self.structure_selection.impropers
                  ]
  
                  if bonds:
                      if len(self.structure_selection.bond_types) == 0:
                          bond_types = np.ones(len(bonds), dtype=int)
                      else:
                          bond_types, unique_bond_types = _get_bond_types(
                              self.structure_selection,
                              self.sigma_conversion_factor,
                              self.epsilon_conversion_factor,
                          )
  
                  if angles:
                      angle_types, unique_angle_types = _get_angle_types(
                          self.structure_selection,
                          self.sigma_conversion_factor,
                          self.epsilon_conversion_factor,
                          use_urey_bradleys=use_urey_bradleys,
                      )
  
                  if dihedrals:
                      dihedral_types, unique_dihedral_types = _get_dihedral_types(
                          self.structure_selection,
                          use_rb_torsions,
                          use_dihedrals,
                          self.epsilon_conversion_factor,
                      )
  
                  if impropers:
                      improper_types, unique_improper_types = _get_impropers(
                          self.structure_selection, self.epsilon_conversion_factor
                      )
  
                  # if self.structure_box_1 != None:
                  if self.structure_box_1:
                      data.write(
                          "*  "
                          + self.filename_box_0
                          + " and "
                          + self.filename_box_1
                          + " - created by mBuild using the on "
                          + str(date_time)
                          + "\n"
                      )
                  else:
                      data.write(
                          "*  "
                          + self.filename_box_0
                          + " - created by mBuild using the on "
                          + str(date_time)
                          + "\n"
                      )
  
                  data.write(
                      "*  "
                      + "parameters from the "
                      + str(self.forcefield_selection)
                      + " force field(s) via MoSDef\n"
                  )
                  data.write(
                      "*  1-4 coulombic scaling = "
                      + str(self.combined_1_4_coul_dict_per_residue)
                      + ", and 1-4 LJ scaling = "
                      + str(self.combined_1_4_lj_dict_per_residue)
                      + "\n\n"
                  )
                  data.write(
                      "*  "
                      + "{:d} atoms\n".format(len(self.structure_selection.atoms))
                  )
  
                  if self.atom_style in ["full", "molecular"]:
                      data.write("*  " + "{:d} bonds\n".format(len(bonds)))
                      data.write("*  " + "{:d} angles\n".format(len(angles)))
                      data.write(
                          "*  " + "{:d} dihedrals\n".format(len(dihedrals))
                      )
                      data.write(
                          "*  " + "{:d} impropers\n\n".format(len(impropers))
                      )
  
                  data.write(
                      "*  " + "{:d} atom types\n".format(len(set(self.types)))
                  )
                  if self.atom_style in ["full", "molecular"]:
                      if bonds:
                          data.write(
                              "*  "
                              + "{:d} bond types\n".format(
                                  len(set(unique_bond_types))
                              )
                          )
                      if angles:
                          data.write(
                              "*  "
                              + "{:d} angle types\n".format(
                                  len(set(unique_angle_types))
                              )
                          )
                      if dihedrals:
                          data.write(
                              "*  "
                              + "{:d} dihedral types\n".format(
                                  len(set(unique_dihedral_types))
                              )
                          )
                      if impropers:
                          data.write(
                              "*  "
                              + "{:d} improper types\n".format(
                                  len(set(unique_improper_types))
                              )
                          )
  
                  data.write("\n")
  
                  data.write("\n* masses\n\n")
                  data.write(
                      "!atom_types \tmass \t\t  atomTypeForceFieldName_ResidueName "
                      + "(i.e., atoms_type_per_utilized_FF)  \n"
                  )
                  for atom_type, mass in self.mass_dict.items():
                      mass_format = "*  {}\t\t{:.6f}\t! {}\n"
                      data.write(
                          mass_format.format(
                              base10_to_base52_alph(
                                  atom_type - self.sub_1_for_base_52
                              ),
                              mass,
                              self.unique_types[atom_type - 1],
                          )
                      )
  
                  # Bond coefficients
                  if bonds:
                      data.write("\n")
                      data.write("BONDS * harmonic\n")
                      data.write("!\n")
                      data.write("!V(bond) = Kb(b - b0)**2\n")
                      data.write("!\n")
                      data.write("!Kb: kcal/mole/A**2\n")
                      data.write("!b0: A\n")
                      data.write(
                          "!Kb (kcal/mol) = Kb_K (K) * Boltz. const.; (9999999999 if no stretching)\n"
                      )
                      data.write("!\n")
  
                      if self.unit_style == "real":
                          data.write(
                              "!atom_types \t Kb\t\tb0 \t\t  atoms_types_per_utilized_FF\n"
                          )
                      for params, idx in unique_bond_types.items():
                          bond_format = "{}\t{}\t{}\t{}\t\t! {}\t{}\n"
                          if (
                              (self.gomc_fix_bonds_angles is not None)
                              and (
                                  (params[3] and params[4])
                                  in self.gomc_fix_bonds_angles
                              )
                          ) or (
                              (
                                  (self.gomc_fix_bonds is not None)
                                  and (
                                      (params[3] and params[4])
                                      in self.gomc_fix_bonds
                                  )
                              )
                          ):
                              fix_bond_k_value = "999999999999"
                              data.write(
                                  bond_format.format(
                                      base10_to_base52_alph(
                                          self.atom_types_to_index_value_dict[
                                              params[2][0] + "_" + str(params[3])
                                          ]
                                      ),
                                      base10_to_base52_alph(
                                          self.atom_types_to_index_value_dict[
                                              params[2][1] + "_" + str(params[4])
                                          ]
                                      ),
                                      fix_bond_k_value,
                                      params[1],
                                      params[2][0] + "_" + str(params[3]),
                                      params[2][1] + "_" + str(params[4]),
                                  )
                              )
  
                          else:
                              data.write(
                                  bond_format.format(
                                      base10_to_base52_alph(
                                          self.atom_types_to_index_value_dict[
                                              params[2][0] + "_" + str(params[3])
                                          ]
                                      ),
                                      base10_to_base52_alph(
                                          self.atom_types_to_index_value_dict[
                                              params[2][1] + "_" + str(params[4])
                                          ]
                                      ),
                                      params[0],
                                      params[1],
                                      params[2][0] + "_" + str(params[3]),
                                      params[2][1] + "_" + str(params[4]),
                                  )
                              )
  
                  # Angle coefficients
                  if angles:
                      if use_urey_bradleys:
                          data.write(
                              "\n!  Urey Bradley terms detected but not written,"
                              + "since they are currently not compatible with GOMC\n"
                          )
  
                      data.write("\nANGLES * harmonic\n")
                      data.write("!\n")
                      data.write("!V(angle) = Ktheta(Theta - Theta0)**2\n")
                      data.write("!\n")
                      data.write("!Ktheta: kcal/mole/rad**2\n")
                      data.write("!Theta0: degrees\n")
                      data.write("!\n")
                      data.write(
                          "! Ktheta (kcal/mol) = Ktheta_K (K) * Boltz. const.\t\t\n"
                      )
                      data.write("!\n")
                      data.write(
                          "!atom_types \t\tKtheta\t\tTheta0\t\t\t  atoms_types_per_utilized_FF\n"
                      )
                      for params, idx in unique_angle_types.items():
  
                          if (
                              (self.gomc_fix_bonds_angles is not None)
                              and ((params[4] and params[5] and params[6]))
                              in self.gomc_fix_bonds_angles
                          ) or (
                              (self.gomc_fix_angles is not None)
                              and ((params[4] and params[5] and params[6]))
                              in self.gomc_fix_angles
                          ):
                              fix_angle_k_value = "999999999999"
                              angle_format = (
                                  "{}\t{}\t{}\t{}\t\t{:.5f}\t\t! {}\t{}\t{}\n"
                              )
                              data.write(
                                  angle_format.format(
                                      base10_to_base52_alph(
                                          self.atom_types_to_index_value_dict[
                                              params[3][0] + "_" + params[4]
                                          ]
                                      ),
                                      base10_to_base52_alph(
                                          self.atom_types_to_index_value_dict[
                                              params[2] + "_" + params[5]
                                          ]
                                      ),
                                      base10_to_base52_alph(
                                          self.atom_types_to_index_value_dict[
                                              params[3][1] + "_" + params[6]
                                          ]
                                      ),
                                      fix_angle_k_value,
                                      params[1],
                                      params[3][0] + "_" + params[4],
                                      params[2] + "_" + params[5],
                                      params[3][1] + "_" + params[6],
                                  )
                              )
  
                          else:
                              data.write(
                                  "{}\t{}\t{}\t{}\t\t{:.5f}\t\t! {}\t{}\t{}\n".format(
                                      base10_to_base52_alph(
                                          self.atom_types_to_index_value_dict[
                                              params[3][0] + "_" + params[4]
                                          ]
                                      ),
                                      base10_to_base52_alph(
                                          self.atom_types_to_index_value_dict[
                                              params[2] + "_" + params[5]
                                          ]
                                      ),
                                      base10_to_base52_alph(
                                          self.atom_types_to_index_value_dict[
                                              params[3][1] + "_" + params[6]
                                          ]
                                      ),
                                      params[0],
                                      params[1],
                                      params[3][0] + "_" + params[4],
                                      params[2] + "_" + params[5],
                                      params[3][1] + "_" + params[6],
                                  )
                              )
  
                  # Dihedral coefficients
                  if dihedrals:
                      if use_rb_torsions:
                          list_if_large_error_dihedral_overall = []
  
                          list_if_largest_error_abs_values_for_dihedral_overall = (
                              []
                          )
                          list_dihedral_overall_error = []
  
                          list_if_abs_max_values_for_dihedral_overall = []
                          list_dihedral_atoms_all_dihedral_overall = []
  
                          data.write("\nDIHEDRALS * CHARMM\n")
                          data.write("!\n")
                          data.write(
                              "!V(dihedral) = Kchi(1 + cos(n(chi) - delta))\n"
                          )
                          data.write("!\n")
                          data.write("!Kchi: kcal/mole\n")
                          data.write("!n: multiplicity\n")
                          data.write("!delta: degrees\n")
                          data.write("!\n")
                          data.write(
                              "! Kchi (kcal/mol) = Kchi_K (K) * Boltz. const.\n"
                          )
                          data.write(
                              "! Boltzmann = 0.0019872041 kcal / (mol * K)\n"
                          )
                          data.write("!\n")
                          if self.unit_style == "real":
                              data.write(
                                  "!atom_types \t\t\tKchi\t\tn\tdelta\t\t  atoms_types_per_utilized_FF\n"
                              )
  
                          for params, idx in unique_dihedral_types.items():
                              charmm_coeffs = RB_to_CHARMM(
                                  params[0],
                                  params[1],
                                  params[2],
                                  params[3],
                                  params[4],
                                  params[5],
                              )
  
                              # check the error between the convertions of RB_tortions to CHARMM DIHEDRALS (start)
                              rb_to_charmm_abs_diff = []
                              no_pi = np.pi
                              dihedral_steps = 2 * 10 ** (-3)
                              dihedral_range = 4 * no_pi
                              dihedral_no_steps = (
                                  int(dihedral_range / dihedral_steps) + 1
                              )
  
                              for i in range(0, dihedral_no_steps + 1):
                                  t = i * dihedral_steps
  
                                  rb_dihedral_calc = (
                                      params[0]
                                      + params[1] * (np.cos(t - no_pi)) ** 1
                                      + params[2] * (np.cos(t - no_pi)) ** 2
                                      + params[3] * (np.cos(t - no_pi)) ** 3
                                      + params[4] * (np.cos(t - no_pi)) ** 4
                                      + params[5] * (np.cos(t - no_pi)) ** 5
                                  )
                                  """CHARMM_torsions
                                  = K0 * (1 + Cos[n0 * (t) - (d0)]) + K1 * (1 + Cos[n1 * (t) - (d1)]) + K2 * (
                                             # 1 + Cos[n2 * (t) - (d2)])
                                  + K3 * (1 + Cos[n3 * (t) - (d3)]) + K4 * (1 + Cos[n4 * (t) - (d4)]) + K5 * (
                                             # 1 + Cos[n5 * (t) - (d5)])
  
                                  = K0 + K1 * (1 + Cos[n1 * (t) - (d1)]) + K2 * (1 + Cos[n2 * (t) - (d2)])
                                  + K3 * (1 + Cos[n3 * (t) - (d3)]) + K4 * (1 + Cos[n4 * (t) - (d4)]) + K5 * (
                                             # 1 + Cos[n5 * (t) - (d5)]). """
  
                                  rb_to_charmm_calc = (
                                      charmm_coeffs[0, 0]
                                      * (
                                          1
                                          + np.cos(
                                              int(charmm_coeffs[0, 1]) * t
                                              - charmm_coeffs[0, 2] * no_pi / 180
                                          )
                                      )
                                      + charmm_coeffs[1, 0]
                                      * (
                                          1
                                          + np.cos(
                                              int(charmm_coeffs[1, 1]) * t
                                              - charmm_coeffs[1, 2] * no_pi / 180
                                          )
                                      )
                                      + charmm_coeffs[2, 0]
                                      * (
                                          1
                                          + np.cos(
                                              int(charmm_coeffs[2, 1]) * t
                                              - charmm_coeffs[2, 2] * no_pi / 180
                                          )
                                      )
                                      + charmm_coeffs[3, 0]
                                      * (
                                          1
                                          + np.cos(
                                              int(charmm_coeffs[3, 1]) * t
                                              - charmm_coeffs[3, 2] * no_pi / 180
                                          )
                                      )
                                      + charmm_coeffs[4, 0]
                                      * (
                                          1
                                          + np.cos(
                                              int(charmm_coeffs[4, 1]) * t
                                              - charmm_coeffs[4, 2] * no_pi / 180
                                          )
                                      )
                                      + charmm_coeffs[5, 0]
                                      * (
                                          1
                                          + np.cos(
                                              int(charmm_coeffs[5, 1]) * t
                                              - charmm_coeffs[5, 2] * no_pi / 180
                                          )
                                      )
                                  )
  
                                  rb_to_charmm_absolute_difference = np.absolute(
                                      rb_dihedral_calc - rb_to_charmm_calc
                                  )
                                  rb_to_charmm_abs_diff.append(
                                      rb_to_charmm_absolute_difference
                                  )
  
                              list_if_large_error_dihedral_iteration = []
                              list_abs_max_dihedral_iteration = []
  
                              if max(rb_to_charmm_abs_diff) > 10 ** (-10):
                                  list_if_large_error_dihedral_iteration.append(1)
                                  list_abs_max_dihedral_iteration.append(
                                      max(rb_to_charmm_abs_diff)
                                  )
  
                                  list_if_large_error_dihedral_overall.append(1)
                                  list_if_largest_error_abs_values_for_dihedral_overall.append(
                                      max(rb_to_charmm_abs_diff)
                                  )
                                  list_dihedral_overall_error.append(
                                      str(params[8])
                                      + ", "
                                      + str(params[9])
                                      + ", "
                                      + str(params[10])
                                      + ", "
                                      + str(params[11])
                                  )
  
                              else:
                                  list_if_large_error_dihedral_iteration.append(0)
  
                                  list_if_abs_max_values_for_dihedral_overall.append(
                                      max(rb_to_charmm_abs_diff)
                                  )
                                  list_dihedral_atoms_all_dihedral_overall.append(
                                      str(params[8])
                                      + ", "
                                      + str(params[9])
                                      + ", "
                                      + str(params[10])
                                      + ", "
                                      + str(params[11])
                                  )
  
                              # **************************************
                              # check the error between the convertions of RB_tortions to CHARMM DIHEDRALS (end)
                              # **************************************
                              dihedral_format = "{}\t{}\t{}\t{}\t{:.6f}\t{}\t{}\t\t! {}\t{}\t{}\t{}\n"
                              data.write(
                                  dihedral_format.format(
                                      base10_to_base52_alph(
                                          self.atom_types_to_index_value_dict[
                                              params[8] + "_" + params[12]
                                          ]
                                      ),
                                      base10_to_base52_alph(
                                          self.atom_types_to_index_value_dict[
                                              params[9] + "_" + params[13]
                                          ]
                                      ),
                                      base10_to_base52_alph(
                                          self.atom_types_to_index_value_dict[
                                              params[10] + "_" + params[14]
                                          ]
                                      ),
                                      base10_to_base52_alph(
                                          self.atom_types_to_index_value_dict[
                                              params[11] + "_" + params[15]
                                          ]
                                      ),
                                      charmm_coeffs[0, 0],
                                      int(charmm_coeffs[0, 1]),
                                      charmm_coeffs[0, 2],
                                      params[8] + "_" + params[12],
                                      params[9] + "_" + params[13],
                                      params[10] + "_" + params[14],
                                      params[11] + "_" + params[15],
                                  )
                              )
                              data.write(
                                  dihedral_format.format(
                                      base10_to_base52_alph(
                                          self.atom_types_to_index_value_dict[
                                              params[8] + "_" + params[12]
                                          ]
                                      ),
                                      base10_to_base52_alph(
                                          self.atom_types_to_index_value_dict[
                                              params[9] + "_" + params[13]
                                          ]
                                      ),
                                      base10_to_base52_alph(
                                          self.atom_types_to_index_value_dict[
                                              params[10] + "_" + params[14]
                                          ]
                                      ),
                                      base10_to_base52_alph(
                                          self.atom_types_to_index_value_dict[
                                              params[11] + "_" + params[15]
                                          ]
                                      ),
                                      charmm_coeffs[1, 0],
                                      int(charmm_coeffs[1, 1]),
                                      charmm_coeffs[1, 2],
                                      params[8] + "_" + params[12],
                                      params[9] + "_" + params[13],
                                      params[10] + "_" + params[14],
                                      params[11] + "_" + params[15],
                                  )
                              )
                              data.write(
                                  dihedral_format.format(
                                      base10_to_base52_alph(
                                          self.atom_types_to_index_value_dict[
                                              params[8] + "_" + params[12]
                                          ]
                                      ),
                                      base10_to_base52_alph(
                                          self.atom_types_to_index_value_dict[
                                              params[9] + "_" + params[13]
                                          ]
                                      ),
                                      base10_to_base52_alph(
                                          self.atom_types_to_index_value_dict[
                                              params[10] + "_" + params[14]
                                          ]
                                      ),
                                      base10_to_base52_alph(
                                          self.atom_types_to_index_value_dict[
                                              params[11] + "_" + params[15]
                                          ]
                                      ),
                                      charmm_coeffs[2, 0],
                                      int(charmm_coeffs[2, 1]),
                                      charmm_coeffs[2, 2],
                                      params[8] + "_" + params[12],
                                      params[9] + "_" + params[13],
                                      params[10] + "_" + params[14],
                                      params[11] + "_" + params[15],
                                  )
                              )
                              data.write(
                                  dihedral_format.format(
                                      base10_to_base52_alph(
                                          self.atom_types_to_index_value_dict[
                                              params[8] + "_" + params[12]
                                          ]
                                      ),
                                      base10_to_base52_alph(
                                          self.atom_types_to_index_value_dict[
                                              params[9] + "_" + params[13]
                                          ]
                                      ),
                                      base10_to_base52_alph(
                                          self.atom_types_to_index_value_dict[
                                              params[10] + "_" + params[14]
                                          ]
                                      ),
                                      base10_to_base52_alph(
                                          self.atom_types_to_index_value_dict[
                                              params[11] + "_" + params[15]
                                          ]
                                      ),
                                      charmm_coeffs[3, 0],
                                      int(charmm_coeffs[3, 1]),
                                      charmm_coeffs[3, 2],
                                      params[8] + "_" + params[12],
                                      params[9] + "_" + params[13],
                                      params[10] + "_" + params[14],
                                      params[11] + "_" + params[15],
                                  )
                              )
                              data.write(
                                  dihedral_format.format(
                                      base10_to_base52_alph(
                                          self.atom_types_to_index_value_dict[
                                              params[8] + "_" + params[12]
                                          ]
                                      ),
                                      base10_to_base52_alph(
                                          self.atom_types_to_index_value_dict[
                                              params[9] + "_" + params[13]
                                          ]
                                      ),
                                      base10_to_base52_alph(
                                          self.atom_types_to_index_value_dict[
                                              params[10] + "_" + params[14]
                                          ]
                                      ),
                                      base10_to_base52_alph(
                                          self.atom_types_to_index_value_dict[
                                              params[11] + "_" + params[15]
                                          ]
                                      ),
                                      charmm_coeffs[4, 0],
                                      int(charmm_coeffs[4, 1]),
                                      charmm_coeffs[4, 2],
                                      params[8] + "_" + params[12],
                                      params[9] + "_" + params[13],
                                      params[10] + "_" + params[14],
                                      params[11] + "_" + params[15],
                                  )
                              )
                              data.write(
                                  dihedral_format.format(
                                      base10_to_base52_alph(
                                          self.atom_types_to_index_value_dict[
                                              params[8] + "_" + params[12]
                                          ]
                                      ),
                                      base10_to_base52_alph(
                                          self.atom_types_to_index_value_dict[
                                              params[9] + "_" + params[13]
                                          ]
                                      ),
                                      base10_to_base52_alph(
                                          self.atom_types_to_index_value_dict[
                                              params[10] + "_" + params[14]
                                          ]
                                      ),
                                      base10_to_base52_alph(
                                          self.atom_types_to_index_value_dict[
                                              params[11] + "_" + params[15]
                                          ]
                                      ),
                                      charmm_coeffs[5, 0],
                                      int(charmm_coeffs[5, 1]),
                                      charmm_coeffs[5, 2],
                                      params[8] + "_" + params[12],
                                      params[9] + "_" + params[13],
                                      params[10] + "_" + params[14],
                                      params[11] + "_" + params[15],
                                  )
                              )
  
                          if sum(list_if_large_error_dihedral_overall) > 0:
                              list_max_error_abs_dihedral_overall = max(
                                  list_if_largest_error_abs_values_for_dihedral_overall
                              )
                              info_if_dihedral_error_too_large = (
                                  "! WARNING: RB-torsion to CHARMM "
                                  "dihedral conversion error"
                                  " is to large [error > 10^(-10)] \n"
                                  "! WARNING: Maximum( "
                                  "|(RB-torsion calc)-(CHARMM dihedral calc)| ) =  "
                                  + str(list_max_error_abs_dihedral_overall)
                                  + "\n"
                              )
                              warn(
                                  "! WARNING: RB-torsion to CHARMM dihedral conversion error"
                                  "is to large [error > 10^(-10)] \n"
                                  "! WARNING: Maximum( |(RB-torsion calc)-(CHARMM dihedral calc)| ) =  "
                                  + str(list_max_error_abs_dihedral_overall)
                                  + "\n"
                              )
                              data.write(info_if_dihedral_error_too_large)
                              print(info_if_dihedral_error_too_large)
                          else:
                              list_if_abs_max_values_for_dihedral_overall_max = (
                                  max(list_if_abs_max_values_for_dihedral_overall)
                              )
                              info_if_dihedral_error_ok = (
                                  "! RB-torsion to CHARMM dihedral conversion error is OK "
                                  "[error <= 10^(-10)]\n"
                                  + "! Maximum( |(RB-torsion calc)-(CHARMM dihedral calc)| ) =  "
                                  + str(
                                      list_if_abs_max_values_for_dihedral_overall_max
                                  )
                                  + "\n"
                              )
                              data.write(info_if_dihedral_error_ok)
                              print(info_if_dihedral_error_ok)
  
                      elif use_dihedrals:
                          data.write(
                              "ERROR: not set up to use to use_dihedrals form for data input from the xml file"
                          )
  
                  # Improper coefficients
                  if impropers:
                      data.write(
                          "ERROR: GOMC is not currently able to use improper in its calculations"
                      )
  
                  # Pair coefficients
                  print(
                      "NBFIX_Mixing not used or no mixing used for the non-bonded potentials out"
                  )
                  print("self.non_bonded_type = " + str(self.non_bonded_type))
                  if self.non_bonded_type == "LJ":
                      data.write("\n")
                      data.write("NONBONDED\n")
                      data.write("!\n")
                      data.write(
                          "!V(Lennard-Jones) = Eps,i,j[(Rmin,i,j/ri,j)**12 - 2(Rmin,i,j/ri,j)**6]\n"
                      )
                      data.write("!\n")
  
                      if self.unit_style == "real":
                          data.write(
                              "!atype \tignored\tepsilon \tRmin/2 \t\tignored\teps,1-4\t\tRmin/2,1-4\t\t"
                              + "  atom_type_per_utilized_FF\n"
                          )
  
                      print("forcefield_dict = " + str(self.forcefield_dict))
  
                      for idx, epsilon in self.epsilon_dict.items():
                          nb_format = "{}\t{:.2f}\t{:.9f}\t{:.11f}\t{:.2f}\t{:.9f}\t{:.11f}\t\t! {}\t{}\n"
                          data.write(
                              nb_format.format(
                                  base10_to_base52_alph(idx),
                                  0,
                                  -epsilon,
                                  self.sigma_dict[idx] * (2 ** (1 / 6)) / 2,
                                  0,
                                  float(self.LJ_1_4_dict[idx]) * (-epsilon),
                                  self.sigma_dict[idx] * (2 ** (1 / 6)) / 2,
                                  self.forcefield_dict[idx],
                                  self.forcefield_dict[idx],
                              )
                          )
  
                  elif self.non_bonded_type == "Mie":
                      data.write(
                          "ERROR: Currently the Mie potential (non_bonded_type) is not supported in this MoSDeF "
                          "GOMC parameter writer\n"
                      )
                      print_error_message = (
                          "ERROR: Currently the Mie potential (non_bonded_type) is not "
                          "supported in this MoSDeF GOMC parameter writer."
                      )
                      raise ValueError(print_error_message)
                  else:
                      data.write(
                          "ERROR: Currently this potential (non_bonded_type) is not supported in "
                          "this MoSDeF GOMC parameter writer\n"
                      )
                      print_error_message = (
                          "ERROR: Currently this potential (non_bonded_type) is not supported in "
                          "this MoSDeF GOMC parameter writer."
                      )
                      raise ValueError(print_error_message)
  
                  # writing end in file
                  data.write("\nEND\n")
  
          # **********************************
          # **********************************
          # FF writer (end)
          # **********************************
          # **********************************
  
      def write_psf(self):
          """This write_psf function writes the Charmm style PSF (topology) file, which can be utilized
          in the GOMC and NAMD engines."""
          # **********************************
          # **********************************
          # psf writer (start)
          # **********************************
          # **********************************
  
          print("******************************")
          print("")
          print(
              "The charmm X-plor format psf writer (the write_psf function) is running"
          )
  
          date_time = datetime.datetime.today()
  
          print(
              "write_psf: forcefield_selection = {}, residues = {}".format(
                  self.forcefield_selection, self.residues
              )
          )
  
          print("******************************")
          print("")
  
          if self.structure_box_1:
              list_of_structures = [
                  self.structure_box_0_ff,
                  self.structure_box_1_ff,
              ]
              list_of_file_names = [self.filename_box_0, self.filename_box_1]
              stuct_only = [self.structure_box_0_ff, self.structure_box_1_ff]
          else:
              list_of_structures = [self.structure_box_0_ff]
              list_of_file_names = [self.filename_box_0]
              stuct_only = [self.structure_box_0_ff]
  
          for q in range(0, len(list_of_structures)):
              stuct_iteration = list_of_structures[q]
              file_name_iteration = list_of_file_names[q]
              output = str(file_name_iteration) + ".psf"
              stuct_only_iteration = stuct_only[q]
              # Lammps syntax depends on the functional form
              # Infer functional form based on the properties of the stuct_iteration
              if self.detect_forcefield_style:
                  # Check  for angles
                  if len(stuct_iteration.urey_bradleys) > 0:
                      print(
                          "Warning: Urey bradley terms detected. GOMC does no support the Urey-Bradley terms"
                      )
                      warn(
                          "warning: Urey bradley terms detected. "
                          "GOMC does no support the Urey-Bradley terms"
                      )
                      use_urey_bradleys = True
                  else:
                      print("No urey bradley terms detected")
                      use_urey_bradleys = False
  
                  # Check for dihedrals
                  if len(stuct_iteration.rb_torsions) > 0:
                      print(
                          "RB Torsions detected, will converted to CHARMM Dihedrals"
                      )
                      use_rb_torsions = True
                      dihedrals_list = stuct_iteration.rb_torsions
                      dihedrals = [
                          [
                              dihedral.atom1.idx + 1,
                              dihedral.atom2.idx + 1,
                              dihedral.atom3.idx + 1,
                              dihedral.atom4.idx + 1,
                          ]
                          for dihedral in stuct_iteration.rb_torsions
                      ]
                  else:
                      use_rb_torsions = False
  
                  if len(stuct_iteration.dihedrals) > 0:
                      print(
                          "Charmm dihedrals detected, so CHARMM Dihedrals will remain"
                      )
                      use_dihedrals = True
                      dihedrals_list = stuct_iteration.dihedrals
                      dihedrals = [
                          [
                              dihedral.atom1.idx + 1,
                              dihedral.atom2.idx + 1,
                              dihedral.atom3.idx + 1,
                              dihedral.atom4.idx + 1,
                          ]
                          for dihedral in stuct_iteration.dihedrals
                      ]
                  else:
                      use_dihedrals = False
              if (use_rb_torsions is False) and (use_dihedrals is False):
                  dihedrals_list = []
                  dihedrals = []
              if use_rb_torsions and use_dihedrals:
                  warn(
                      "Multiple dihedral styles detected, check your "
                      "Forcefield XML and structure files"
                  )
  
              # Check for impropers
              for dihedral in stuct_iteration.dihedrals:
                  if dihedral.improper:
                      warn(
                          "ERROR: Amber-style impropers are currently not supported in GOMC"
                      )
  
              impropers_list = stuct_iteration.impropers
              impropers = [
                  [
                      improper.atom1.idx + 1,
                      improper.atom2.idx + 1,
                      improper.atom3.idx + 1,
                      improper.atom4.idx + 1,
                  ]
                  for improper in stuct_iteration.impropers
              ]
  
              no_atoms = len(stuct_iteration.atoms)
              no_bonds = len(stuct_iteration.bonds)
              no_angles = len(stuct_iteration.angles)
  
              no_dihedrals = len(dihedrals)
              no_impropers = len(impropers)
  
              no_donors = len(stuct_iteration.donors)
              no_acceptors = len(stuct_iteration.acceptors)
              no_groups = len(stuct_iteration.groups)
  
              # psf printing (start)
  
              residue_data_list = []
              residue_names_list = []
              for k, atom in enumerate(stuct_only_iteration.atoms):
                  residue_data_list.append(str(atom.residue))
                  residue_names_list.append(atom.residue.name)
  
              unique_residue_data_dict = {}
              unique_residue_data_list = []
              residue_data_name_list = []
  
              for m, residue in enumerate(stuct_only_iteration.residues):
                  unique_residue_data_list.append(
                      str(stuct_only_iteration.residues[m])
                  )
                  unique_residue_data_dict.update(
                      {unique_residue_data_list[m]: m + 1}
                  )
                  residue_data_name_list.append(
                      stuct_only_iteration.residues[m].name
                  )
  
              res_no_chain_iter_corrected = []
              residue_id_list = []
              residue_id_adder_fixed_struct_wo_bonds = 0
              for f, PSF_atom_iteration_0 in enumerate(
                  stuct_only_iteration.atoms
              ):
                  if f > 0:
                      if (
                          PSF_atom_iteration_0.residue.chain
                          == previous_residue_chain
                          and len(PSF_atom_iteration_0.bonds) == 0
                      ):
                          residue_id_adder_fixed_struct_wo_bonds += 1
  
                  previous_residue_chain = PSF_atom_iteration_0.residue.chain
  
                  residue_id_int = int(
                      unique_residue_data_dict[residue_data_list[f]]
                      + residue_id_adder_fixed_struct_wo_bonds
                  )
                  res_id_adder = int(
                      (residue_id_int % self.max_residue_no) % self.max_residue_no
                  )
                  if int(res_id_adder) == 0:
                      res_no_iteration_corrected = int(self.max_residue_no)
                  else:
                      res_no_iteration_corrected = res_id_adder
  
                  res_no_chain_iter_corrected.append(res_no_iteration_corrected)
                  residue_id_list.append(residue_id_int)
  
              output_write = genopen(output, "w")
  
              first_indent = "%8s"
              psf_formating = (
                  "%8s %-4s %-4s %-4s %-4s %4s %10.6f %13.4f" + 11 * " "
              )
  
              output_write.write("PSF ")
              output_write.write("\n\n")
  
              no_of_remarks = 3
              output_write.write(first_indent % no_of_remarks + " !NTITLE\n")
              output_write.write(
                  " REMARKS this file "
                  + file_name_iteration
                  + " - created by mBuild/foyer using the"
                  + "\n"
              )
              output_write.write(
                  " REMARKS parameters from the "
                  + str(self.forcefield_selection)
                  + " force field via MoSDef\n"
              )
              output_write.write(
                  " REMARKS created on " + str(date_time) + "\n\n\n"
              )
  
              # This converts the atom name in the GOMC psf and pdb files to unique atom names
              print(
                  "bead_to_atom_name_dict = {}".format(
                      self.bead_to_atom_name_dict
                  )
              )
              [
                  unique_individual_atom_names_dict,
                  individual_atom_names_list,
                  missing_bead_to_atom_name,
              ] = unique_atom_naming(
                  stuct_only_iteration,
                  residue_id_list,
                  residue_names_list,
                  bead_to_atom_name_dict=self.bead_to_atom_name_dict,
              )
  
              if None in [
                  unique_individual_atom_names_dict,
                  individual_atom_names_list,
                  missing_bead_to_atom_name,
              ]:
                  self.input_error = True
                  print_error_message = (
                      "ERROR: The unique_atom_naming function failed while "
                      "running the charmm_writer function. Ensure the proper inputs are "
                      "in the bead_to_atom_name_dict."
                  )
                  raise ValueError(print_error_message)
  
              # ATOMS: Calculate the atom data
              # psf_formating is conducted for the for CHARMM format (i.e., atom types are base 52, letters only)
              output_write.write(first_indent % no_atoms + " !NATOM\n")
              for i_atom, PSF_atom_iteration_1 in enumerate(
                  stuct_iteration.atoms
              ):
                  segment_id = PSF_atom_iteration_1.residue.segid or "SYS"
                  atom_type_iter = base10_to_base52_alph(
                      self.atom_types_to_index_value_dict[
                          PSF_atom_iteration_1.type
                          + "_"
                          + PSF_atom_iteration_1.residue.name
                      ]
                  )
  
                  atom_lines_iteration = psf_formating % (
                      i_atom + 1,
                      segment_id,
                      res_no_chain_iter_corrected[i_atom],
                      str(residue_names_list[i_atom])[: self.max_resname_char],
                      individual_atom_names_list[i_atom],
                      atom_type_iter,
                      PSF_atom_iteration_1.charge,
                      PSF_atom_iteration_1.mass,
                  )
  
                  output_write.write("%s\n" % atom_lines_iteration)
  
              output_write.write("\n")
  
              # BONDS: Calculate the bonding data
              output_write.write(first_indent % no_bonds + " !NBOND: bonds\n")
              for i_bond, PSF_bond_iteration_1 in enumerate(
                  stuct_iteration.bonds
              ):
                  output_write.write(
                      (first_indent * 2)
                      % (
                          PSF_bond_iteration_1.atom1.idx + 1,
                          PSF_bond_iteration_1.atom2.idx + 1,
                      )
                  )
  
                  if (i_bond + 1) % 4 == 0:
                      output_write.write("\n")
  
              if no_bonds % 4 == 0:
                  output_write.write("\n")
              else:
                  output_write.write("\n\n")
  
              if no_bonds == 0:
                  output_write.write("\n")
  
              # ANGLES: Calculate the angle data
              output_write.write(first_indent % no_angles + " !NTHETA: angles\n")
              for i_angle, angle_iteration in enumerate(stuct_iteration.angles):
  
                  output_write.write(
                      (first_indent * 3)
                      % (
                          angle_iteration.atom1.idx + 1,
                          angle_iteration.atom2.idx + 1,
                          angle_iteration.atom3.idx + 1,
                      )
                  )
  
                  if (i_angle + 1) % 3 == 0:
                      output_write.write("\n")
  
              if no_angles % 3 == 0:
                  output_write.write("\n")
              else:
                  output_write.write("\n\n")
  
              if no_angles == 0:
                  output_write.write("\n")
  
              # DIHEDRALS: Calculate the dihedral  data
              output_write.write(
                  first_indent % no_dihedrals + " !NPHI: dihedrals\n"
              )
              for i_dihedral, dihedral_iter in enumerate(dihedrals_list):
                  (
                      dihedral_atom_1,
                      dihedral_atom_2,
                      dihedral_atom_3,
                      dihedral_atom_4,
                  ) = (
                      dihedral_iter.atom1,
                      dihedral_iter.atom2,
                      dihedral_iter.atom3,
                      dihedral_iter.atom4,
                  )
  
                  output_write.write(
                      (first_indent * 4)
                      % (
                          dihedral_atom_1.idx + 1,
                          dihedral_atom_2.idx + 1,
                          dihedral_atom_3.idx + 1,
                          dihedral_atom_4.idx + 1,
                      )
                  )
  
                  if (i_dihedral + 1) % 2 == 0:
                      output_write.write("\n")
  
              if no_dihedrals % 2 == 0:
                  output_write.write("\n")
              else:
                  output_write.write("\n\n")
  
              if no_dihedrals == 0:
                  output_write.write("\n")
  
              # IMPROPERS: Calculate the improper data
              output_write.write(
                  first_indent % no_impropers + " !NIMPHI: impropers\n"
              )
              for i_improper, improper_iter in enumerate(impropers_list):
                  (
                      improper_atom_1,
                      improper_atom_2,
                      improper_atom_3,
                      improper_atom_4,
                  ) = (
                      improper_iter.atom1,
                      improper_iter.atom2,
                      improper_iter.atom3,
                      improper_iter.atom4,
                  )
  
                  output_write.write(
                      (first_indent * 4)
                      % (
                          improper_atom_1.idx + 1,
                          improper_atom_2.idx + 1,
                          improper_atom_3.idx + 1,
                          improper_atom_4.idx + 1,
                      )
                  )
  
                  if (i_improper + 1) % 2 == 0:
                      output_write.write("\n")
  
              if no_impropers % 2 == 0:
                  output_write.write("\n")
              else:
                  output_write.write("\n\n")
  
              if no_impropers == 0:
                  output_write.write("\n")
  
              # DONOR: calculate the donor data
              output_write.write(first_indent % no_donors + " !NDON: donors\n")
              for donor_i, donor_iter in enumerate(stuct_iteration.donors):
  
                  output_write.write(
                      (first_indent * 2)
                      % (donor_iter.atom1.idx + 1, donor_iter.atom2.idx + 1)
                  )
                  if (donor_i + 1) % 4 == 0:
                      output_write.write("\n")
  
              if no_donors % 4 == 0:
                  output_write.write("\n")
              else:
                  output_write.write("\n\n")
  
              if no_donors == 0:
                  output_write.write("\n")
  
              # ACCEPTOR: calculate the acceptor data
              output_write.write(
                  first_indent % no_acceptors + " !NACC: acceptors\n"
              )
              for acceptor_i, acceptor_iter in enumerate(
                  stuct_iteration.acceptors
              ):
  
                  output_write.write(
                      (first_indent * 2)
                      % (acceptor_iter.atom1.idx + 1, acceptor_iter.atom2.idx + 1)
                  )
                  if (acceptor_i + 1) % 4 == 0:
                      output_write.write("\n")
  
              if no_acceptors % 4 == 0:
                  output_write.write("\n")
              else:
                  output_write.write("\n\n")
  
              if no_acceptors == 0:
                  output_write.write("\n")
  
              # NNB: calculate the NNB data
              output_write.write(first_indent % 0 + " !NNB\n\n")
              for nbb_i, atoms_iter in enumerate(stuct_iteration.atoms):
  
                  output_write.write(first_indent % 0)
                  if (nbb_i + 1) % 8 == 0:
                      output_write.write("\n")
  
              if no_atoms % 8 == 0:
                  output_write.write("\n")
              else:
                  output_write.write("\n\n")
  
              if no_atoms == 0:
                  output_write.write("\n")
  
              # GROUP: calculate the group data
              try:
                  group_data = stuct_iteration.groups.nst2
              except AttributeError:
                  group_data = 0
              output_write.write(
                  (first_indent * 2) % (no_groups or 1, group_data) + " !NGRP \n"
              )
              if stuct_iteration.groups is True:
                  for group_i, group_iter in enumerate(stuct_iteration.groups):
  
                      output_write.write(
                          (first_indent * 3)
                          % (
                              group_iter.atom.idx,
                              group_iter.type,
                              group_iter.move,
                          )
                      )
                      if (group_i + 1) % 3 == 0:
                          output_write.write("\n")
  
                  if no_groups % 3 == 0:
                      output_write.write("\n")
                  else:
                      output_write.write("\n\n")
  
                  if no_groups == 0:
                      output_write.write("\n")
  
              else:
                  structure_abs_charge_value = abs(
                      sum(
                          atom_charge_iter.charge
                          for atom_charge_iter in stuct_iteration.atoms
                      )
                  )
                  if structure_abs_charge_value < 1.0e-4:
                      group_type = 1
                  else:
                      group_type = 2
                  output_write.write((first_indent * 3) % (0, group_type, 0))
                  output_write.write("\n")
  
              output_write.write("\n")
              output_write.close()
          # **********************************
          # **********************************
          # psf writer (end)
          # **********************************
          # **********************************
  
      def write_pdb(self):
          """This write_psf function writes the Charmm style PDB (coordinate file), which can be utilized
          in the GOMC and NAMD engines."""
          # **********************************
          # **********************************
          # pdb writer (start)
          # **********************************
          # **********************************
          date_time = datetime.datetime.today()
          print("******************************")
          print("")
          print("The charmm pdb writer (the write_pdb function) is running")
          print("write_charmm_pdb: residues == {}".format(self.residues))
          print("fix_residue = {}".format(self.fix_residue))
          print("fix_residue_in_box = {}".format(self.fix_residue_in_box))
          print("bead_to_atom_name_dict = {}".format(self.bead_to_atom_name_dict))
  
          if self.fix_residue is None and self.fix_residue_in_box is None:
              print(
                  "INFORMATION: No atoms are fixed in this pdb file for the GOMC simulation engine. "
              )
          else:
              warn(
                  "Some atoms are fixed in this pdb file for the GOMC simulation engine. "
              )
  
          print("******************************")
          print("")
  
          if self.structure_box_1:
              list_of_structures = [
                  self.structure_box_0_ff,
                  self.structure_box_1_ff,
              ]
              list_of_file_names = [self.filename_box_0, self.filename_box_1]
              stuct_only = [self.structure_box_0_ff, self.structure_box_1_ff]
          else:
              list_of_structures = [self.structure_box_0_ff]
              list_of_file_names = [self.filename_box_0]
              stuct_only = [self.structure_box_0_ff]
  
          for q in range(0, len(list_of_structures)):
              file_name_iteration = list_of_file_names[q]
              output = str(file_name_iteration) + ".pdb"
              stuct_only_iteration = stuct_only[q]
  
              output_write = genopen(output, "w")
              # output_write.write(
              #'REMARK this file ' + file_name_iteration + ' - created by mBuild/foyer using the' + '\n')
              # output_write.write(
              #'REMARK parameters from the ' + str(self.forcefield_selection) + ' force field via MoSDef\n')
              # output_write.write('REMARK created on ' + str(date_time) + '\n')
  
              unique_residue_data_dict = {}
              unique_residue_data_list = []
              residue_data_name_list = []
              for m, residue in enumerate(stuct_only_iteration.residues):
                  unique_residue_data_list.append(
                      str(stuct_only_iteration.residues[m])
                  )
                  unique_residue_data_dict.update(
                      {unique_residue_data_list[m]: m + 1}
                  )
                  residue_data_name_list.append(
                      stuct_only_iteration.residues[m].name
                  )
  
              for n in range(0, len(residue_data_name_list)):
                  if residue_data_name_list[n] not in self.residues:
                      self.input_error = True
                      print_error_message = "ERROR: Please specifiy all residues (residues) in a list"
                      raise ValueError(print_error_message)
  
              residue_data_list = []
              for k, atom in enumerate(stuct_only_iteration.atoms):
                  residue_data_list.append(str(atom.residue))
  
              if (self.fix_residue is not None) and (
                  self.fix_residue_in_box is not None
              ):
                  for n in range(0, len(self.fix_residue)):
                      if self.fix_residue[n] in self.fix_residue_in_box:
                          self.input_error = True
                          print_error_message = (
                              "ERROR: residue type can not be specified to both "
                              "fix_residue and fix_residue_in_box"
                          )
                          raise ValueError(print_error_message)
  
              residue_names_list = []
              fix_atoms_list = []
              for k, atom in enumerate(stuct_only_iteration.atoms):
                  residue_names_list.append(atom.residue.name)
                  if (self.fix_residue is not None) and (
                      atom.residue.name in self.fix_residue
                  ):
                      beta_iteration = 1.00
                  elif (self.fix_residue_in_box is not None) and (
                      atom.residue.name in self.fix_residue_in_box
                  ):
                      beta_iteration = 2.00
                  else:
                      beta_iteration = 0.00
                  fix_atoms_list.append(beta_iteration)
  
              if stuct_only_iteration.box is not None:
                  output_write.write(
                      "CRYST1%9.3f%9.3f%9.3f%7.2f%7."
                      "2f%7.2f %-11s%4s\n"
                      % (
                          stuct_only_iteration.box[0],
                          stuct_only_iteration.box[1],
                          stuct_only_iteration.box[2],
                          stuct_only_iteration.box[3],
                          stuct_only_iteration.box[4],
                          stuct_only_iteration.box[5],
                          stuct_only_iteration.space_group,
                          "",
                      )
                  )
  
              all_atom_coordinates = stuct_only_iteration.get_coordinates("all")
              if all_atom_coordinates is None:
                  self.input_error = True
                  print_error_message = (
                      "ERROR: the submitted structure has no PDB coordinates, "
                      "so the PDB writer has terminated. "
                  )
                  raise ValueError(print_error_message)
  
              pdb_atom_line_format = "ATOM  %5s %-4s%1s%-4s%1s%4d%1s   %8.3f%8.3f%8.3f%6.2f%6.2f      %-4s%2s%-2s\n"
  
              atom_alternate_location_list = []
              residue_code_insertion_list = []
              x_list = []
              y_list = []
              z_list = []
              atom_occupancy_list = []
              atom_bfactor_list = []
              element_list = []
  
              # lock occupany factor at 1 (instead of: atom.occupancy)
              locked_occupany_factor = 1.00
              max_no_atoms_in_base10 = 99999  # 99,999 for atoms in psf/pdb
  
              res_no_chain_iter_corrected = []
              res_chain_iteration_corrected_list = []
              residue_id_list = []
              residue_id_adder_fixed_struct_wo_bonds = (
                  0  # for example zeolite used as fixed atoms wo bonds
              )
              for i, atom_iter in enumerate(stuct_only_iteration.atoms):
                  if i > 0:
                      if (
                          atom_iter.residue.chain == previous_residue_chain
                          and len(atom_iter.bonds) == 0
                      ):
                          residue_id_adder_fixed_struct_wo_bonds += 1
  
                  previous_residue_chain = atom_iter.residue.chain
                  residue_id_int = int(
                      unique_residue_data_dict[residue_data_list[i]]
                      + residue_id_adder_fixed_struct_wo_bonds
                  )
                  res_chain_iteration_corrected_list.append(
                      base10_to_base26_alph(
                          int(residue_id_int / (self.max_residue_no + 1))
                      )[-1:]
                  )
                  res_id_adder = int(
                      (residue_id_int % self.max_residue_no) % self.max_residue_no
                  )
                  if int(res_id_adder) == 0:
                      res_no_chain_iter_corrected.append(int(self.max_residue_no))
                  else:
                      res_no_chain_iter_corrected.append(res_id_adder)
  
                  residue_id_list.append(residue_id_int)
  
              # This converts the atom name in the CHARMM psf and pdb files to unique atom names
              [
                  unique_individual_atom_names_dict,
                  individual_atom_names_list,
                  missing_bead_to_atom_name,
              ] = unique_atom_naming(
                  stuct_only_iteration,
                  residue_id_list,
                  residue_names_list,
                  bead_to_atom_name_dict=self.bead_to_atom_name_dict,
              )
  
              if None in [
                  unique_individual_atom_names_dict,
                  individual_atom_names_list,
                  missing_bead_to_atom_name,
              ]:
                  self.input_error = True
                  print_error_message = (
                      "ERROR: The unique_atom_naming function failed while "
                      "running the charmm_writer function. Ensure the proper inputs are "
                      "in the bead_to_atom_name_dict."
                  )
  
                  raise ValueError(print_error_message)
  
              for coord_iter, atom_coordinates in enumerate(all_atom_coordinates):
  
                  for PDB_residue_count in stuct_only_iteration.residues:
                      segment_id = ""
                      atom_iteration = sorted(
                          PDB_residue_count.atoms, key=lambda atom: atom.number
                      )
                      for atom_iteration_2 in atom_iteration:
                          x, y, z = atom_coordinates[atom_iteration_2.idx]
                          atom_alternate_location_list.append(
                              atom_iteration_2.altloc
                          )
                          residue_code_insertion_list.append(
                              PDB_residue_count.insertion_code[:1]
                          )
                          x_list.append(x)
                          y_list.append(y)
                          z_list.append(z)
                          atom_occupancy_list.append(atom_iteration_2.occupancy)
                          atom_bfactor_list.append(atom_iteration_2.bfactor)
                          element_list.append(
                              Element[atom_iteration_2.atomic_number].upper()
                          )
  
                  for v, atom_iter_1 in enumerate(stuct_only_iteration.atoms):
  
                      if v + 1 > max_no_atoms_in_base10:
                          atom_number = base10_to_base16_alph_num(v + 1)
  
                      else:
                          atom_number = v + 1
  
                      output_write.write(
                          pdb_atom_line_format
                          % (
                              atom_number,
                              individual_atom_names_list[v],
                              atom_alternate_location_list[v],
                              str(residue_names_list[v])[: self.max_resname_char],
                              res_chain_iteration_corrected_list[v],
                              res_no_chain_iter_corrected[v],
                              residue_code_insertion_list[v],
                              x_list[v],
                              y_list[v],
                              z_list[v],
                              locked_occupany_factor,
                              fix_atoms_list[v],
                              segment_id,
                              element_list[v],
                              "",
                          )
                      )
  
                  output_write.write("%-80s\n" % "END")
  
              output_write.close()
  
              # **********************************
              # **********************************
              # pdb writer (end)
              # **********************************
              # **********************************
  