Source code for mbuild.formats.charmm_writer

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
[docs]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. """
[docs] 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 ) )
[docs] 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) # ********************************** # **********************************
[docs] 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) # ********************************** # **********************************
[docs] 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) # ********************************** # **********************************