Data Structures

CHARMM-style PDB, PSF, and Force Field File Writers

class mosdef_gomc.formats.gmso_charmm_writer.Charmm(structure_box_0, filename_box_0, structure_box_1=None, filename_box_1=None, forcefield_selection=None, residues=None, gomc_fix_bonds_angles=None, gomc_fix_bonds=None, gomc_fix_angles=None, bead_to_atom_name_dict=None, atom_type_naming_style='all_unique', fix_residue=None, fix_residue_in_box=None, set_residue_pdb_occupancy_to_1=None, ff_filename=None, gmso_match_ff_by='molecule')[source]
Generates a Charmm object via foyer and gmso 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 and NAMD output data files.

Mw = g/mol

charge = e

Energy units = Lennard-Jones (LJ) -> kcal/mol, Mie and Exp-6 (Buckingham) -> Kelvin (K)

Harmonic bonds: Kb = Energy units, b0 = Angstroms

Harmonic angles: Ktheta = (Energy units)/rad**2 , Theta0 = degrees

Dihedral angles Harmonic (not unavailable): Ktheta = Energy units, n (LJ) = 0 (unitless integer), delta = degrees

Dihedral angles: Periodic: Ktheta = Energy units, n (unitless integer) -> n (LJ) = 1-5 and n (Mie or Buckingham) = 0-5, delta = degrees

Improper angles Harmonic (not unavailable): Ktheta = Energy units, n (LJ) = 0 (unitless integer), delta = degrees

Improper angles: Periodic: Ktheta = Energy units, n = 1-5 (unitless integer), delta = degrees

Lennard-Jones (LJ)-NONBONDED: epsilon = Energy units, Rmin/2 = Angstroms

Mie-NONBONDED: epsilon = Energy units, sigma = Angstroms, n = integer (unitless)

Exp6-NONBONDED: epsilon = Energy units, sigma = Angstroms (FF XMLs in Rmin), alpha = integer (unitless)

Lennard-Jones (LJ)-NBFIX (not unavailable): epsilon = Energy units, Rmin = Angstroms

Mie-NBFIX (not unavailable): same as Mie-NONBONDED

Exp6-NBFIX (not unavailable): same as Exp6-NONBONDED

Note: The Charmm object is only compatible with molecules with a single residue (single residue name).

Note: The units are only the same for GOMC and NAMD when using the Lennard-Jones (LJ) non-bonded type.

Note: NAMD is only compatible with the Lennard-Jones (LJ) non-bonded type.

Note: There are CHARMM style atom type naming conventions based on the force field xml provided and the GOMC simulation input parameters, outlined in the atom_type_naming_style input variable.
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 multiple 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.

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’}

gmso_match_ff_by: str (“group” or “molecule”), default = “molecule”
How the GMSO force field is applied, using the molecules name/residue name (mbuild.Compound().name) for GOMC and NAMD. This is regardless number of levels in the (mbuild.Compound().name or mbuild.Box()).
“molecule” applies the force field using the molecule’s name
WARNING: This is the atom’s name for a single atom molecule (1 atom/bead = molecule).

Molecule > 1 atom —-> uses the “mbuild.Compound().name” (1 level above the atoms/beads) is the molecule’s name. This “mb.Compound().name” (1 level above the atoms/beads) needs to be used in the Charmm object’s residue_list and forcefield_selection (if >1 force field), and will be the residue name in the PSF, PDB, and FF files.

Molecule = 1 atom/bead —-> uses the “atom/bead’s name” is the molecule’s name.

This “atom/bead’s name” needs to be used in the Charmm object’s residue_list and forcefield_selection (if >1 force field), and will be the residue name in the PSF, PDB, and FF files.

NOTE: Non-bonded zeolites or other fixed structures without bonds will use the “atom/bead’s name” as the molecule’s name, if they are single non-bonded atoms. However, the user may want to use the “group” option instead for this type of system, if applicable (see the “group” option).

Example (Charmm_writer selects the user changeable “ETH”, in the Charmm object residue_list and forcefield_selection (if >1 force field), which sets the residue “ETH” in the PSF, PDB, and FF files).
>>> ethane = mbuild.load("CC", smiles=True)
>>> ethane.name = "ETH"
>>> ethane_box = mbuild.fill_box(
>>>     compound=[ethane],
>>>     n_compounds=[100],
>>>     box=[4, 4, 4]
>>> )
Example (Charmm_writer must to select the non-user changeable “_CH4” (per the foyer TraPPE force field), in the Charmm object residue_list and forcefield_selection (if >1 force field), which sets the residue “_CH4” in the PSF, PDB, and FF files).
>>> methane_ua_bead_name = "_CH4"
>>> methane_child_bead = mbuild.Compound(name=methane_ua_bead_name)
>>> methane_box = mbuild.fill_box(
>>>    compound=methane_child_bead,
>>>    n_compounds=4,
>>>    box=[1, 2, 3]
>>> )
methane_box.name = "MET"
Example (Charmm_writer must to select the non-user changeable “Na” and “Cl” in the Charmm object residue_list and forcefield_selection (if >1 force field), which sets the residues “Na” and “Cl” in the PSF, PDB, and FF files).
>>> sodium_atom_name = "Na"
>>> sodium_child_atom = mbuild.Compound(name=sodium_atom_name)
>>> sodium = mb.Compound(name="SOD")
>>> sodium.add(sodium_child_atom, inherit_periodicity=False)
>>> chloride_atom_name = "Cl"
>>> chloride_child_bead = mbuild.Compound(name=chloride_atom_name)
>>> chloride = mb.Compound(name="CHL")
>>> chloride.add(chloride_child_atom, inherit_periodicity=False)
>>> sodium_chloride_box = mbuild.fill_box(
>>>     compound=[sodium, chloride],
>>>     n_compounds=[4, 4],
>>>     box=[1, 2, 3]
>>> )
Example zeolite (Charmm_writer must to select the non-user changeable “Si” and “O” in the Charmm object residue_list and forcefield_selection (if >1 force field), which sets the residues “Si” and “O” in the PSF, PDB, and FF files):
>>> lattice_cif_ETV_triclinic = load_cif(file_or_path=get_mosdef_gomc_fn("ETV_triclinic.cif"))
>>> ETV_triclinic = lattice_cif_ETV_triclinic.populate(x=1, y=1, z=1)
>>> ETV_triclinic.name = "ETV"
“group” applies the force field to any molecules under the mbuild.Compound().name (children).
The “group” search starts at the top level if there is no mbuild.Compound().name found on that level, then the next level down (children) are searched, keeping moving down every time.

WARNING: This “group” option will take all the molecules below it, regardless if they are selected to be separate from the group via the residue_list and forcefield_selection (if >1 force field).

This is ideal to use when you are building simulation box(es) using mbuild.fill_box(), with molecules, and it allows you to add another level to single single atom molecules (1 atom/bead = molecule) to rename the mbuild.Compound().name, changing the residue’s name and allowing keeping the atom/bead name so the force field is applied properly.

NOTE: This “group” option may be best for non-bonded zeolites or other fixed structures without bonds, if they are single non-bonded atoms. Using this “group” option, the user can select the residue name for the Charmm_writer’s residue_list and forcefield_selection (if >1 force field) to force field all the atoms with a single residue name, and output this residue name in the PSF, PDB, and FF files.

Example (Charmm_writer select the user changeable “MET” in the Charmm object residue_list and forcefield_selection (if >1 force field), which sets the residue “MET” in the PSF, PDB, and FF files):
>>> methane_ua_bead_name = "_CH4"
>>> methane_child_bead = mbuild.Compound(name=methane_ua_bead_name)
>>> methane_box = mbuild.fill_box(compound=methane_child_bead, n_compounds=4, box=[1, 2, 3])
>>> methane_box.name = "MET"
Example (Charmm_writer select the user changeable “MET” in the Charmm object residue_list and forcefield_selection (if >1 force field), which sets the residue “MET” in the PSF, PDB, and FF files):
>>> methane_child_bead = mb.Compound(name="_CH4")
>>> methane = mb.Compound(name="MET")
>>> methane.add(methane_child_bead, inherit_periodicity=False)
>>> box_liq = mb.fill_box(
>>>     compound=methane,
>>>     n_compounds=1230,
>>>     box=[4.5, 4.5, 4.5]
>>> )
Example zeolite (Charmm_writer select the user changeable “ETV” in the Charmm object residue_list and forcefield_selection (if >1 force field), which sets the residue “ETV” in the PSF, PDB, and FF files):
>>> lattice_cif_ETV_triclinic = load_cif(file_or_path=get_mosdef_gomc_fn("ETV_triclinic.cif"))
>>> ETV_triclinic = lattice_cif_ETV_triclinic.populate(x=1, y=1, z=1)
>>> ETV_triclinic.name = "ETV"
ff_filename: str, default =None
If a string, it will write the force field files that work in GOMC and NAMD structures.

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}

atom_type_naming_style: str, optional, default=’all_unique’, (‘general’ or ‘all_unique’)
‘general’: WARNING: The ‘general’ convention is UNSAFE, and the EXPERT user SHOULD USE AT THEIR OWN RISK, making SURE ALL THE BONDED PARAMETERS HAVE THE SAME VALUES IN THE UTILIZED FORCE FIELD XMLs. Also, this DOES NOT ENSURE that THERE ARE NO specific Foyer XML ATOM TYPE BONDED CONNECTIONS in the Foyer FORCE FIELD XMLs, instead of the Foyer atom class type bonded connections, which could RESULT IN AN INCORRECT FORCE FIELD PARAMETERIZATION. This is UNSAFE to use even with the same force field XML file, so the EXPERT user SHOULD USE AT THEIR OWN RISK.

The ‘general’ convention only tests if the sigma, epsilons, mass, and Mie-n values are identical between the different molecules (residues in this context) and their applied force fields and DOES NOT check that any or all of the bonded parameters have the same or conflicting values.

The ‘general’ convention is where all the atom classes in the Foyer force field XML files are converted to the CHARMM-style atom types (FOYER ATOM CLASSES). The ‘general’ convention ONLY auto-checks that the sigma, epsilon, mass, and Mie-n values are the same and does not currently ensure all the bonded parameters are the same or conflicting between different force field XML files. If the sigma, epsilons, mass, and Mie-n values are the same between force fields the general method can be applied; if not, it defaults to the ‘all_unique’ method.

Example of CHARMM style atom types in an all-atom ethane and ethanol system: * Ethane: alkane carbon = CT, alkane hydrogen = HC * Ethanol: alkane carbon = CT, alkane hydrogen = HC , oxygen in alcohol = OH, hydrogen in alcohol = OH

This is only permitted when the following is true; otherwise it will default to the the ‘all_unique’: * All the MoSDeF force field XML’s atom classes’ non-bonded parameters (sigma, epsilon, mass, and Mie-n power constant) values are THE SAME. * If the general CHARMM style atom type in any residue/molecule’s gomc_fix_bonds_angles, gomc_fix_bonds, or gomc_fix_angles NOT IN any other residue/molecule, the ‘all_unique’ type will be used.

‘all_unique’: The ‘all_unique’ convention is the SAFE way to parameterize the system. The MoSDeF force field XML atom names within residue/molecule are all unique, where each one adds an alpha numberic value after the MoSDeF force field XML atom classes to ensure uniqueness within the within residue/molecule. The OPLS atom types/name do not require all the sigma, epsilon, mass values to be the same, but have less bonded class parameters.

Example of CHARMM style atom types in an all-atom ethane and ethanol system: * Ethane: alkane carbon type 0 = CT0, alkane hydrogen type 0 = HC0 * Ethanol: alkane carbon type 1 = CT1, alkane carbon type 2 = CT2, alkane hydrogen type 1 = HC1 , oxygen in alcohol type 0 = OH0, hydrogen in alcohol type 0 = OH0

This is selected when auto-selected when: * All the MoSDeF force field XML’s atom classes’ non-bonded parameters (sigma, epsilon, mass, and Mie-n power constant) values are NOT THE SAME. * If the general CHARMM style atom type in any residue/molecule’s gomc_fix_bonds_angles, gomc_fix_bonds, or gomc_fix_angles are IN any other residue/molecule.

fix_residue: list or None, default=None
Changes occcur in the pdb file only. All the atoms in the residue are have their Beta values in the PDB file set to 1.00; Otherwise, they will be 0.00. NOTE: In GOMC, this only fixes atoms automatically as listed below. NOTE: In NAMD, these all fixes need to be manually set in the control file (please see NAMD manual). 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 (Temperature factor) 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. All the atoms in the residue are have their Beta values in the PDB file set to 2.00; Otherwise, they will be 0.00. NOTE: In GOMC, this only fixes atoms automatically as listed below. NOTE: In NAMD, these all fixes need to be manually set in the control file (please see NAMD manual). 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 (Temperature factor) 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. NOTE that this is mainly for GOMC but also applies for NAMD (please see NAMD manual).

set_residue_pdb_occupancy_to_1: list or None, default=None
Changes occcur in the pdb file only. All the atoms in the residue are have their occupancy values in the PDB file set to 1.00; Otherwise, they will be 0.00. NOTE: In GOMC, This defines which atoms belong to which box for the GCMC and GEMC ensembles. NOTE: In NAMD, This can be used for fixes which are manually set in the control file (please see NAMD manual).
Notes

Urey-Bradleys, and NBFIX are not currenly supported. Currently the NBFIX is not available but will be in the near future. OPLS, AMBER, and CHARMM forcefield styles are supported (without harmonic dihedrals),

The atom typing is currently provided via a base 62 numbering (0 to 9, then capital and lowercase lettering A to Z). This base 62 numbering allows for (62)^2 unique atom types. There are 4 characters allowed for the atom type name and 2 for the base 62 addition making them unique if needed. Otherwise, if the GMSO general class is used for the atom typing, then the 6 characters is allowed, because additional numbers are not required to make them unique. This specifically avoids X because this represents any atom type in the CHARMM-style force field files.

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.

GOMC standard LJ form = epsilon * ( (Rmin/r)**12 - 2*(Rmin/r)**6 ) ———> = 4*epsilon * ( (sigma/r)**12 - (sigma/r)**6 ), when converted to sigmas Both forms above are accepted and compared automatically, but all input FFs have to be of the above input forms, aside from the whole potential energy scaling factor.

GOMC standard Mie form = (n/(n-m)) * (n/m)**(m/(n-m)) * epsilon * ((sigma/r)**n - (sigma/r)**m) where m = 6 –> (n/(n-6)) * (n/6)**(6/(n-6)) * epsilon * ((sigma/r)**n - (sigma/r)**6) The above form is accepted but only if all input FFs have the same form, aside from the whole potential energy scaling factor.

GOMC standard Exp-6 form = alpha*epsilon/(alpha -6) * Exp( alpha*(1-r/Rmin) - (Rmin/r)**6 ), where r >= Rmax infinity , where r < Rmax The above form is accepted but only if all input FFs have the same form, aside from the whole potential energy scaling factor.

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).

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’}

ff_filename: str, default =None
If a string, it will write the force field files that work in GOMC and NAMD structures.

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. All the atoms in the residue are have their Beta values in the PDB file set to 1.00; Otherwise, they will be 0.00. NOTE: In GOMC, this only fixes atoms automatically as listed below. NOTE: In NAMD, these all fixes need to be manually set in the control file (please see NAMD manual). 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 (Temperature factor) 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. All the atoms in the residue are have their Beta values in the PDB file set to 2.00; Otherwise, they will be 0.00. NOTE: In GOMC, this only fixes atoms automatically as listed below. NOTE: In NAMD, these all fixes need to be manually set in the control file (please see NAMD manual). 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 (Temperature factor) 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. NOTE that this is mainly for GOMC but also applies for NAMD (please see NAMD manual).

set_residue_pdb_occupancy_to_1: list or None, default=None
Changes occcur in the pdb file only. All the atoms in the residue are have their occupancy values in the PDB file set to 1.00; Otherwise, they will be 0.00. NOTE: In GOMC, This defines which atoms belong to which box for the GCMC and GEMC ensembles. NOTE: In NAMD, This can be used for fixes which are manually set in the control file (please see NAMD manual).

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))

topology_box_0_ff: gmso.Topology
The box 0 topology (from structure_box_0) after all the provided force fields are applied.

topology_box_1_ff: gmso.Topology
The box 1 topology (from structure_box_1) after all the provided force fields are applied. 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_kcal_per_mol_atom_type_dict: dict {str: float or int}
The uniquely numbered atom type (key) and it’s non-bonded epsilon coefficient in units of kcal/mol (value). The atom type is defined by the AtomClass_ResidueName (Example of a carbon atom in an ethane molecule, AtomClass_ResidueName –> CT_ETH).

sigma_angstrom_atom_type_dict: dict {str: float or int}
The uniquely numbered atom type (key) and it’s non-bonded sigma coefficient in angstroms (value). The atom type is defined by the AtomClass_ResidueName (Example of a carbon atom in an ethane molecule, AtomClass_ResidueName –> CT_ETH).

mie_n_atom_type_dict: dict {str: float, int, or None}
The uniquely numbered atom type (key) and it’s non-bonded unitless n coefficient (value). The atom type is defined by the AtomClass_ResidueName (Example of a carbon atom in an ethane molecule, AtomClass_ResidueName –> CT_ETH). NOTE: The value is None if a Mie FF is not used.

exp6_alpha_atom_type_dict: dict {str: float, int, or None}
The uniquely numbered atom type (key) and it’s non-bonded unitless alpha coefficient (value). The atom type is defined by the AtomClass_ResidueName (Example of a carbon atom in an ethane molecule, AtomClass_ResidueName –> CT_ETH). NOTE: The value is None if a Exp6 FF is not used.

nonbonded_1_4_dict: dict {str: float or int}
The uniquely numbered atom type (key) and it’s non-bonded 1-4 scaling factor (value). The atom type is defined by the AtomClass_ResidueName (Example of a carbon atom in an ethane molecule, AtomClass_ResidueName –> CT_ETH).

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.

combined_1_4_nonbonded_dict_per_residue
The residue name/molecule (key) and it’s non-bonded 1-4 coulombic scaling factor (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.

electrostatic_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.

NOTE: NAMD and GOMC can not have mulitple electrostatic 1-4 scalers, since it is provided as a single input in their control files.

combined_1_4_electrostatic_dict_per_residue: dict, {str: float or int}
The residue name/molecule (key) and it’s non-bonded 1-4 coulombic scaling factor (value).

NOTE: NAMD and GOMC can not have mulitple electrostatic 1-4 scalers, since it is provided as a single input in their control files.

combining_rule: str (‘geometric’ or ‘lorentz’”’),
The possible mixing/combining rules are ‘geometric’ or ‘lorentz’, which provide the geometric and arithmetic mixing rule, respectively. NOTE: Arithmetic means the ‘lorentz’ combining or mixing rule. NOTE: GMSO default to the ‘lorentz’ mixing rule if none is provided, and this writers default is the GMSO default. NOTE: NAMD and GOMC can not have multiple values for the combining_rule.

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’}

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: set(str, …, str)}
A dictionary that provides the residue names (keys) and a set of the unique atom names in the residue (value), for the combined structures (box 0 and box 1 (if supplied)).

__init__(structure_box_0, filename_box_0, structure_box_1=None, filename_box_1=None, forcefield_selection=None, residues=None, gomc_fix_bonds_angles=None, gomc_fix_bonds=None, gomc_fix_angles=None, bead_to_atom_name_dict=None, atom_type_naming_style='all_unique', fix_residue=None, fix_residue_in_box=None, set_residue_pdb_occupancy_to_1=None, ff_filename=None, gmso_match_ff_by='molecule')[source]

write_inp()[source]

This write_inp function writes the Charmm style parameter (force field) file, which can be utilized in the GOMC and NAMD engines.

write_pdb(space_group='P 1')[source]

This write_pdb function writes the Charmm style PDB (coordinate file), which can be utilized in the GOMC and NAMD engines.

Parameters:

space_group: str (default=”P 1”)
The space group of the structure

write_psf()[source]

This write_psf function writes the Charmm style PSF (topology) file, which can be utilized in the GOMC and NAMD engines.

GOMC Control File Writer

mosdef_gomc.formats.gmso_gomc_conf_writer.print_required_input(description=False)[source]

Prints the required ensemble arguments with an optional description based on the ensemble type

Parameters:

description: bool, default = False.
If True, it prints the descriptions of the required ensemble inputs (i.e. dict), If False, it only prints the required ensemble inputs without the descriptions (i.e. list)

Returns:

Prints the required ensemble arguments with an optional description based on the ensemble type
mosdef_gomc.formats.gmso_gomc_conf_writer.write_gomc_control_file(charmm_object, conf_filename, ensemble_type, RunSteps, Temperature, ff_psf_pdb_file_directory=None, check_input_files_exist=True, Restart=False, Checkpoint=False, ExpertMode=False, Parameters=None, Coordinates_box_0=None, Structure_box_0=None, Coordinates_box_1=None, Structure_box_1=None, binCoordinates_box_0=None, extendedSystem_box_0=None, binVelocities_box_0=None, binCoordinates_box_1=None, extendedSystem_box_1=None, binVelocities_box_1=None, input_variables_dict=None)[source]
The usable command that creates the GOMCControl object and writes the GOMC control file via the GOMCControl.write_conf_file function.

Constructs the GOMCControl object with the defaults, or adding additional data in the input_variable section. Default setting for the GOMC configuraion files are based upon an educated guess, which should result in reasonable sampling for a given ensemble/simulation type. However, there is no guarantee that the default setting will provide the best or adequate sampling for the selected system. The user has the option to modify the configuration/control files based on the simulation specifics or to optimize the system beyond the standard settings. These override options are available via the keyword arguments in input_variable_dict.
Parameters:
charmm_object: Charmm object
Charmm object is has been parameterized from the selected force field.

ensemble_typ: str, [‘NVT’, ‘NPT’, ‘GEMC_NPT’, ‘GCMC-NVT’, ‘GCMC’]
The ensemble type of the simulation.

RunSteps: int (>0), must be an integer greater than zero.
Sets the total number of simulation steps.

Temperature: unyt.unyt_quantity in temperature units (> 1 K)
Temperature of system in unyt temperature units (written in Kelvin).

ff_psf_pdb_file_directory: str (optional), default=None (i.e., the current directory).
The full or relative directory added to the force field, psf, and pdb file names, created via the Charmm object.

check_input_files_exist: bool, (default=True)
Check if the force field, psf, and pdb files exist. If input, it also checks if the Checkpoint, binCoordinates_box_0, extendedSystem_box_0, binVelocities_box_0, binCoordinates_box_1, extendedSystem_box_1, and binVelocities_box_1 files. If the files are checked and do not exist, the writer will throw a ValueError. True, check if the all the input files exist. False, do not check if the all the input files exist.

Restart: boolean, default = False
Determines whether to restart the simulation from restart file (*_restart.pdb and *_restart.psf) or not.

Checkpoint: list or bool (see below), default = False
If True: Checkpoint = [True, str] —> [True, checkpoint_filename] If False: Checkpoint = False or [False, str] —> [False, checkpoint_filename] Checkpoint files guarantees the trajectory files can be concatenated, along with ensuring no deviation from a single simulation’s results in a simulation which was interrupted and restarted. Determines whether to restart the simulation from checkpoint file or not. Restarting the simulation with would result in an identical outcome, as if previous simulation was continued. This is required for hybrid Monte-Carlo Molecular Dyanamics in open-ensembles (GCMC/GEMC) to concatenate trajectory files since the molecular transfers rearranges the order of the molecules. Checkpointing will ensure the molecules are loaded in the same order each cycle. NOTE: The checkpoint_filename with or without the path can be entered with the .chk extension, or without an extension; if it does not have the file extension, the .chk extension will then be added to the filename.

ExpertMode: boolean, default = False
This allows the move ratios to be any value, regardless of the ensemble, provided all the move ratios sum to 1. For example, this mode is utilized to easily equilibrate a GCMC or GEMC ensemble in a pseudo NVT mode by removing the requirement that the volume and swap moves must be non-zero. In other words, when the volume and swap moves are zero, the GCMC and GEMC ensembles will run pseudo NVT simulations in 1 and 2 simulation boxes, respectively. The simulation’s output and restart files will keep their original output structure for the selected ensemble, which is advantageous when automating a workflow.

Parameters: str, (default=None)
Override all other force field directory and filename input with the correct extension (.inp or .par). Note: the default directory is the current directory with the Charmm object file name.

Coordinates_box_0: str, (default=None)
Override all other box 0 pdb directory and filename inputs with the correct extension. Note: the default directory is the current directory with the Charmm object file name.

Structure_box_0: str, (default=None)
Override all other box 0 psf directory and filename inputs with the correct extension. Note: the default directory is the current directory with the Charmm object file name.

Coordinates_box_1: str, (default=None)
Override all other box 1 pdb directory and filename inputs with the correct extension. Note: the default directory is the current directory with the Charmm object file name.

Structure_box_1: str, (default=None)
Override all other box 1 psf directory and filename inputs with the correct extension. Note: the default directory is the current directory with the Charmm object file name.

binCoordinates_box_0: str, (default=None)
The box 0 binary coordinate file is used only for restarting a GOMC simulation, which provides increased numerical accuracy.

extendedSystem_box_0: str, (default=None)
The box 0 vectors and origin file is used only for restarting a GOMC simulation.

binVelocities_box_0: str, (default=None)
The box 0 binary velocity file is used only for restarting a GOMC simulation, which provides increased numerical accuracy. These velocities are only passed thru GOMC since Monte Carlo simulations do not utilize any velocity information.

binCoordinates_box_1: str, (default=None)
The box 1 binary coordinate file is used only for restarting a GOMC simulation, which provides increased numerical accuracy.

extendedSystem_box_1: str, (default=None)
The box 1 vectors and origin file is used only for restarting a GOMC simulation.

binVelocities_box_1: str, (default=None)
The box 1 binary velocity file is used only for restarting a GOMC simulation, which provides increased numerical accuracy. These velocities are only passed thru GOMC since Monte Carlo simulations do not utilize any velocity information.

input_variables_dict: dict, default=None
These input variables are optional and override the default settings. Changing these variables likely required for more advanced systems. The details of the acceptable input variables for the selected ensembles can be found by running the code below in python,
>>> print_valid_ensemble_input_variables('GCMC', description = True)
which prints the input_variables with their subsection description for the selected ‘GCMC’ ensemble (other ensembles can be set as well).

Example: input_variables_dict = {‘PRNG’: 123, ‘ParaTypeCHARMM’: True }
# *******************************************************************

# input_variables_dict options (keys and values) - (start)

# Note: the input_variables_dict keys are also attributes

# *******************************************************************

PRNG: string or int (>= 0) (“RANDOM” or int), default = “RANDOM”
PRNG = Pseudo-Random Number Generator (PRNG). There are two (2) options, entering the string, “RANDOM”, or a integer.

— “RANDOM”: selects a random seed number. This will enter the line “PRNG RANDOM” in the gomc configuration file.

— integer: which defines the integer seed number for the simulation. This is equivalent to entering the following two lines in the configuration file

line 1 = PRNG INTSEED

line 2 = Random_Seed user_selected_integer.

Example 1: for a random seed enter the string “RANDOM”.

Example 2: for a specific seed number enter a integer of your choosing.

ParaTypeCHARMM: boolean, default = True
True if a CHARMM forcefield, False otherwise. Note: This is changed by the MoSDeF-GOMC Charmm object if the ParaTypeCHARMM or ParaTypeMie is detected.

ParaTypeMie: boolean, default = False
True if a Mie or Exp6 forcefield types, False otherwise. Note: This is changed by the MoSDeF-GOMC Charmm object if the ” ParaTypeCHARMM or ParaTypeMie is detected.

ParaTypeMARTINI: boolean, default = False
True if a MARTINI forcefield, False otherwise.

RcutCoulomb_box_0: unyt.unyt_quantity in length units (>= 0), default=None
Sets a specific radius in box 0 (written in Angstroms) where the short-range electrostatic energy will be calculated (i.e., The distance to truncate the short-range electrostatic energy in box 0.) Note: if None, GOMC will default to the Rcut value

RcutCoulomb_box_1: unyt.unyt_quantity in length units (>= 0), default=None
Sets a specific radius in box 1 (written in Angstroms) where the short-range electrostatic energy will be calculated (i.e., The distance to truncate the short-range electrostatic energy in box 0.) Note: if None, GOMC will default to the Rcut value

Pressure: unyt.unyt_quantity in pressure units (>= 0), default = 1.01325 bar
The pressure (written in bar) is utilized for the NPT and GEMC_NPT simulations.

Rcut: unyt.unyt_quantity in length units (>= 0 and RcutLow < Rswitch < Rcut),
default = 10 angstroms Sets a specific radius (written in Angstroms) that non-bonded interaction energy and force will be considered and calculated using defined potential function. The distance in Angstoms to truncate the LJ, Mie, or other VDW type potential at. Note: Rswitch is only used when the “Potential” = SWITCH.

RcutLow: unyt.unyt_quantity in length units (>= 0 and RcutLow < Rswitch < Rcut),
default = 0 angstroms Sets a specific minimum possible distance (written in Angstroms) that reject any move that places any atom closer than specified distance. The minimum possible distance between any atoms. Sets a specific radius in Angstroms that non-bonded interaction Note: Rswitch is only used when the “Potential” = SWITCH. WARNING: When using a molecule that has charge atoms with non-bonded epsilon values of zero (i.e., water), the RcutLow need to be greater than zero, typically 1 angstrom. WARNING: When using the free energy calculations, RcutLow needs to be set to zero (RcutLow=0); otherwise, the free energy calculations can produce results that are slightly off or wrong.

LRC: boolean, default = True
If True, the simulation considers the long range tail corrections for the non-bonded VDW or dispersion interactions. Note: In case of using SHIFT or SWITCH potential functions, LRC will be ignored.

IPC: boolean, default = False
If True, the simulation adds the impulse correction term to the pressure, which considers to correct for the discontinuous Rcut potential (i.e., a hard cutoff potential, meaning a potential without tail corrections) the long range tail corrections for the non-bonded VDW or dispersion interactions. If False, the impulse correction term to the pressure is not applied. Note: This can not be used if LRC is True or the Potential is set to SWITCH, or SHIFT.

Exclude: str [“1-2”, “1-3”, or “1-4”], default = “1-3”
Note: In CHARMM force field, the 1-4 interaction needs to be considered. Choosing “Excude 1-3”, will modify 1-4 interaction based on 1-4 parameters in parameter file. If a kind force field is used, where 1-4 interaction needs to be ignored, such as TraPPE, either Exclude “1-4” needs to be chosen or 1-4 parameter needs to be assigned to zero in the parameter file.

— “1-2”: All interaction pairs of bonded atoms, except the ones that separated with one bond, will be considered and modified using 1-4 parameters defined in parameter file.

— “1-3”: All interaction pairs of bonded atoms, except the ones that separated with one or two bonds, will be considered and modified using 1-4 parameters defined in parameter file.

— “1-4”: All interaction pairs of bonded atoms, except the ones that separated with one, two or three bonds, will be considered using non-bonded parameters defined in parameter file.

Potential: str, [“VDW”, “EXP6”, “SHIFT” or “SWITCH”], default = “VDW”
Defines the potential function type to calculate non-bonded dispersion interaction energy and force between atoms.

— “VDW”: Non-bonded dispersion interaction energy and force calculated based on n-6 (Lennard - Jones) equation. This function will be discussed further in the Intermolecular energy and Virial calculation section.

— “EXP6”: Non-bonded dispersion interaction energy and force calculated based on exp-6 equation.

— “SHIFT”: This option forces the potential energy to be zero at Rcut distance.

— “SWITCH”: This option smoothly forces the potential energy to be zero at Rcut distance and starts modifying the potential at Rswitch distance. Depending on force field type, specific potential function will be applied.

Rswitch: unyt.unyt_quantity in length units (>= 0 and RcutLow < Rswitch < Rcut),
default = 9 angstroms Note: Rswitch is only used when the SWITCH function is used (i.e., “Potential” = SWITCH). The Rswitch distance is written in Angstrom. If the “SWITCH” function is chosen, Rswitch needs to be defined, otherwise, the program will be terminated. When using choosing “SWITCH” as potential function, the Rswitch distance defines where the non-bonded interaction energy modification is started, which is eventually truncated smoothly at Rcut distance.

ElectroStatic: boolean, default = True
Considers the coulomb interactions or not. If True, coulomb interactions are considered and false if not. Note: To simulate the polar molecule in MARTINI force field, ElectroStatic needs to be turn on (i.e., True). The MARTINI force field uses short-range coulomb interaction with constant Dielectric of 15.0.

Ewald: boolean, default = True
Considers the standard Ewald summation method for electrostatic calculations. If True, Ewald summation calculation needs to be considered and false if not. Note: By default, GOMC will set ElectroStatic to True if Ewald summation method was used to calculate coulomb interaction.

CachedFourier: boolean, default = False
Considers storing the reciprocal terms for Ewald summation calculation in order to improve the code performance. This option would increase the code performance with the cost of memory usage. If True, to store reciprocal terms of Ewald summation calculation and False if not. Warning: Monte Carlo moves, such as MEMC-1, MEMC-2, MEMC-3, IntraMEMC-1, IntraMEMC-2, and IntraMEMC-3 are not support with CachedFourier.

Tolerance: float (0.0 < float < 1.0), default = 1e-05
Sets the accuracy in Ewald summation calculation. Ewald separation parameter and number of reciprocal vectors for the Ewald summation are determined based on the accuracy parameter.

Dielectric: int or float (>= 0.0), default = 15
Sets dielectric value used in coulomb interaction when the Martini force field is used. Note: In MARTINI force field, Dielectric needs to be set to 15.0.

PressureCalc: list [bool , int (> 0)] or [bool , step_frequency],
default = [True, 10k] or [True , set via formula based on the number of RunSteps or 10k max] Calculate the system pressure or not. bool = True, enables the pressure calculation during the simulation, false disables the calculation. The int/step frequency sets the frequency of calculating the pressure.

EqSteps: int (> 0), default = set via formula based on the number of RunSteps or 1M max
Sets the number of steps necessary to equilibrate the system. Averaging will begin at this step. Note: In GCMC simulations, the Histogram files will be outputed at EqSteps.

AdjSteps: int (> 0), default = set via formula based on the number of RunSteps or 1k max
Sets the number of steps per adjustment of the parameter associated with each move (e.g. maximum translate distance, maximum rotation, maximum volume exchange, etc.).

VDWGeometricSigma: boolean, default = True
Use geometric mean, as required by OPLS force field, to combining Lennard-Jones sigma parameters for different atom types. If set to True, GOMC uses geometric mean to combine Lennard-Jones or VDW sigmas. Note: The default setting is pulled from the force field XML, if not present it defaults to geometric via MoSDeF’s default setting of VDWGeometricSigma is True, which uses the arithmetic mean when combining Lennard-Jones or ” VDW sigma parameters for different atom types. NOTE: In GOMC, for Mie FFs the following is always true –> n_ij = (n_ii + n_jj)/2. NOTE: In GOMC, for Exp FFs the following is always true –> alpha_ij = (alpha_ii * alpha_jj)**0.5.

useConstantArea: boolean, default = False
Changes the volume of the simulation box by fixing the cross-sectional area (x-y plane). If True, the volume will change only in z axis, If False, the volume of the box will change in a way to maintain the constant axis ratio.

FixVolBox0: boolean, default = False
Changing the volume of fluid phase (Box 1) to maintain the constant imposed pressure and Temperature, while keeping the volume of adsorbed phase (Box 0) fixed. Note: By default, GOMC will set useConstantArea to False if no value was set. It means, the volume of the box will change in a way to maintain the constant axis ratio.

ChemPot:dict {str (4 dig limit) , unyt.unyt_quantity in energy units (kcal/mol, kJ/mol, K)},
default=None The chemical potentials are written in GOMC units of energy, K. There is a 4 character limit for the string/residue name since the PDB/PSF files have a 4 character limitation and require and exact match in the conf file. Note: These strings must match the residue in the psf and psb files or it will fail. The name of the residues and their corresponding chemical potential must specified for every residue in the system (i.e., {“residue_name”: chemical_potential}). Note: IF 2 KEYS WITH THE SAME STRING/RESIDUE ARE PROVIDED, ONE WILL BE AUTOMATICALLY OVERWRITTEN AND NO ERROR WILL BE THROWN IN THIS CONTROL FILE WRITER.

Example 1 (system with only water): {“H2O”: -4000 * unyt.K}

Example 2 (system with water and ethanol): {“H2O”: -4000 * unyt.K, “ETH”: -8000 * unyt.K}

Fugacity: dict {str, unyt.unyt_quantity in pressure units (>= 0)}, default=None
The fugacity is written in GOMC units of pressure, bar. There is a 4 character limit for the string/residue name since the PDB/PSF files have a 4 character limitation and require and exact match in the conf file. Note: These strings must match the residue in the psf and psb files or it will fail. The name of the residues and their corresponding fugacity must specified for every residue in the system (i.e., {“residue_name”: fugacity}). Note: IF 2 KEYS WITH THE SAME STRING/RESIDUE ARE PROVIDED, ONE WILL BE AUTOMATICALLY OVERWRITTEN AND NO ERROR WILL BE THROWN IN THIS CONTROL FILE WRITER.

Example 1 (system with only water): {“H2O”: 1 * unyt.bar}

Example 2 (system with water and ethanol): {“H2O”: 0.5 * unyt.bar, “ETH”: 10 * unyt.atm}

CBMC_First: int (>= 0), default = 12
The number of CD-CBMC trials to choose the first atom position (Lennard-Jones trials for first seed growth).

CBMC_Nth: int (>= 0), default = 10
The Number of CD-CBMC trials to choose the later atom positions (Lennard-Jones trials for first seed growth).

CBMC_Ang: int (>= 0), default = 50
The Number of CD-CBMC bending angle trials to perform for geometry (per the coupled-decoupled CBMC scheme).

CBMC_Dih: int (>= 0), default = 50
The Number of CD-CBMC dihedral angle trials to perform for geometry (per the coupled-decoupled CBMC scheme).

OutputName: str (NO SPACES), , default = “Output_data”, default = [True, 1M] or
[True , set via formula based on the number of RunSteps or 1M max] The UNIQUE STRING NAME, WITH NO SPACES, which is used for the output block average, PDB, and PSF file names.

CoordinatesFreq: list [bool , int (> 0)] or [Generate_data_bool , steps_per_data_output_int],
default = [False, 1M] or [False , set via formula based on the number of RunSteps or M max] Controls output of PDB file (coordinates). If bool is True, this enables outputting the coordinate files at the integer frequency (set steps_per_data_output_int), while “False” disables outputting the coordinates.

DCDFreq: list [bool , int (> 0)] or [Generate_data_bool , steps_per_data_output_int],
default = [True, 1M] or [True , set via formula based on the number of RunSteps or M max] Controls output of DCD file (coordinates). If bool is True, this enables outputting the coordinate files at the integer frequency (set steps_per_data_output_int), while “False” disables outputting the coordinates.

RestartFreq: list [bool , int (> 0)] or [Generate_data_bool , steps_per_data_output_int],
default = [True, 1M] or [True , set via formula based on the number of RunSteps or 1M max] This creates the PDB and PSF (coordinate and topology) files for restarting the system at the set steps_per_data_output_int (frequency) If bool is True, this enables outputting the PDB/PSF restart files at the integer frequency (set steps_per_data_output_int), while “false” disables outputting the PDB/PSF restart files.

CheckpointFreq: list [bool , int (> 0)] or [Generate_data_bool , steps_per_data_output_int],
default = [True, 1M] or [True , set via formula based on the number of RunSteps or 1M max] Controls the output of the last state of simulation at a specified step, in a binary file format (checkpoint.dat). Checkpoint file contains the following information in full precision:

Last simulation step that saved into checkpoint file

Simulation cell dimensions and angles

(3) Maximum amount of displacement (Å), rotation (δ), and volume (Å^3) that is used in the Displacement, Rotation, MultiParticle, and Volume moves

Number of Monte Carlo move trial and acceptance

All molecule’s coordinates

Random number sequence

If bool is True, this enables outputing the checkpoint file at the integer frequency (set steps_per_data_ouput_int), while “False” disables outputting the checkpoint file.

ConsoleFreq: list [bool , int (> 0)] or [Generate_data_bool , steps_per_data_output_int],
default = [True, 10k] or [True , set via formula based on the number of RunSteps or 10k max] Controls the output to the “console” or log file, which prints the acceptance statistics, and run timing info. In addition, instantaneously-selected thermodynamic properties will be output to this file. If bool is True, this enables outputting the console data at the integer frequency (set steps_per_data_output_int), while “False” disables outputting the console data file.

BlockAverageFreq: list [bool , int (> 0)] or [Generate_data_bool , steps_per_data_output_int],
default = [True, 10k] or [True , set via formula based on the number of RunSteps or 10k max] Controls the block averages output of selected thermodynamic properties. Block averages are averages of thermodynamic values of interest for chunks of the simulation (for post-processing of averages or std. dev. in those values). If bool is True, this enables outputting the block averaging data/file at the integer frequency (set steps_per_data_output_int), while “False” disables outputting the block averaging data/file.

HistogramFreq: list [bool , int (> 0)] or [Generate_data_bool , steps_per_data_output_int],
default = [True, 10k] or [True , set via formula based on the number of RunSteps or 10k max] Controls the histograms. Histograms are a binned listing of observation frequency for a specific thermodynamic variable. In the GOMC code, they also control the output of a file containing energy/molecule samples, which is only used for the “GCMC” ensemble simulations for histogram reweighting purposes. If bool is True, this enables outputting the data to the histogram data at the integer frequency (set steps_per_data_output_int), while “False” disables outputting the histogram data.

DistName: str (NO SPACES), default = “dis”
Short phrase which will be combined with RunNumber and RunLetter to use in the name of the binned histogram for molecule distribution.

HistName: str (NO SPACES), default = “his”
Short phrase, which will be combined with RunNumber and RunLetter, to use in the name of the energy/molecule count sample file.

RunNumber: int ( > 0 ), default = 1
Sets a number, which is a part of DistName and HistName file name.

RunLetter: str (1 alphabetic character only), default = “a”
Sets a letter, which is a part of DistName and HistName file name.

SampleFreq: int ( > 0 ), default = 500
The number of steps per histogram sample or frequency.

OutEnergy: [bool, bool], default = [True, True]
The list provides the booleans to [block_averages_bool, console_output_bool]. This outputs the energy data into the block averages and console output/log

OutPressure: [bool, bool], default = [True, True]
The list provides the booleans to [block_averages_bool, console_output_bool]. This outputs the pressure data into the block averages and console output/log files.

OutMolNum: [bool, bool], default = [True, True]
The list provides the booleans to [block_averages_bool, console_output_bool]. This outputs the number of molecules data into the block averages and console output/log files.

OutDensity: [bool, bool], default = [True, True]
The list provides the booleans to [block_averages_bool, console_output_bool]. This outputs the density data into the block averages and console output/log files.

OutVolume: [bool, bool], default = [True, True]
The list provides the booleans to [block_averages_bool, console_output_bool]. This outputs the volume data into the block averages and console output/log files.

OutSurfaceTension: [bool, bool], default = [False, False]
The list provides the booleans to [block_averages_bool, console_output_bool]. This outputs the surface tension data into the block averages and console output/log files.

FreeEnergyCalc: list [bool , int (> 0)] or [Generate_data_bool , steps_per_data_output_int],
default=None bool = True enabling free energy calculation during the simulation, false disables the calculation. The int/step frequency sets the frequency of calculating the free energy. WARNING: When using the free energy calculations, RcutLow needs to be set to zero (RcutLow=0); otherwise, the free energy calculations can produce results that are slightly off or wrong.

MoleculeType: list [str , int (> 0)] or [“residue_name” , residue_ID], default=None
The user must set this variable as there is no working default. Note: ONLY 4 characters can be used for the string (i.e., “residue_name”). Sets the solute molecule kind (residue name) and molecule number (residue ID), which absolute solvation free will be calculated for.

InitialState: int (>= 0), default=None
The user must set this variable as there is no working default. The index of LambdaCoulomb and LambdaVDW vectors. Sets the index of the LambdaCoulomb and LambdaVDW vectors, to determine the simulation lambda value for VDW and Coulomb interactions. WARNING: This must an integer within the vector count of the LambdaVDW and LambdaCoulomb, in which the counting starts at 0.

LambdaVDW: list of floats (0 <= floats <= 1), default=None
The user must set this variable as there is no working default (default = {}). Lambda values for VDW interaction in ascending order. Sets the intermediate lambda states to which solute-solvent VDW interactions are scaled.

WARNING: This list must be the same length as the “LambdaCoulomb” list length.

WARNING: All lambda values must be stated in the ascending order, starting with 0.0 and ending with 1.0; otherwise, the program will terminate.

Example of ascending order 1: [0.0, 0.1, 1.0]

Example of ascending order 2: [0.0, 0.1, 0.2, 0.4, 0.9, 1.0]

LambdaCoulomb: list of floats (0 <= floats <= 1), default=None
Lambda values for Coulombic interaction in ascending order. Sets the intermediate lambda states to which solute-solvent Coulombic interactions are scaled. GOMC defauts to the “LambdaVDW” values for the Coulombic interaction if no “LambdaCoulomb” variable is set.

WARNING: This list must be the same length as the “LambdaVDW” list length.

WARNING: All lambda values must be stated in the ascending order, starting with 0.0 and ending with 1.0; otherwise, the program will terminate.

Example of ascending order 1: [0.0, 0.1, 1.0]

Example of ascending order 2: [0.0, 0.1, 0.2, 0.4, 0.9, 1.0]

ScaleCoulomb: bool, default = False
Determines to scale the Coulombic interaction non-linearly (soft-core scheme) or not. True if the Coulombic interaction needs to be scaled non-linearly. False if the Coulombic interaction needs to be scaled linearly.

ScalePower: int (>= 0), default = 2
The p value in the soft-core scaling scheme, where the distance between solute and solvent is scaled non-linearly.

ScaleAlpha: int or float (>= 0), default = 0.5
The alpha value in the soft-core scaling scheme, where the distance between solute and solvent is scaled non-linearly.

MinSigma: unyt.unyt_quantity in length units (>= 0), default = 3 angstroms
The minimum sigma value in the soft-core scaling scheme (written in Angstroms), where the distance between solute and solvent is scaled non-linearly.

DisFreq: int or float (0 <= value <= 1), default are specific for each ensemble
{‘NVT’: 0.15, ‘NPT’: 0.15, ‘GEMC_NVT’: 0.19, ‘GEMC_NPT’: 0.19, ‘GCMC’: 0.15} Fractional percentage at which the displacement move will occur (i.e., fraction of displacement moves).

RotFreq: int or float (0 <= value <= 1), default are specific for each ensemble
{‘NVT’: 0.15, ‘NPT’: 0.15, ‘GEMC_NVT’: 0.2, ‘GEMC_NPT’: 0.2, ‘GCMC’: 0.15} Fractional percentage at which the rotation move will occur. (i.e., fraction of rotation moves).

IntraSwapFreq: int or float (0 <= value <= 1), default are specific for each ensemble
{‘NVT’: 0.3, ‘NPT’: 0.29, ‘GEMC_NVT’: 0.1, ‘GEMC_NPT’: 0.1, ‘GCMC’: 0.1} Fractional percentage at which the molecule will be removed from a box and inserted into the same box using coupled-decoupled configurational-bias algorithm. (i.e., fraction of intra-molecule swap moves).

SwapFreq: int or float (0 <= value <= 1), default are specific for each ensemble
{‘NVT’: 0.0, ‘NPT’: 0.0, ‘GEMC_NVT’: 0.2, ‘GEMC_NPT’: 0.2, ‘GCMC’: 0.35} For Gibbs and Grand Canonical (GC) ensemble runs only: Fractional percentage at which molecule swap move will occur using coupled-decoupled configurational-bias. (i.e., fraction of molecule swaps moves).

RegrowthFreq: int or float (0 <= value <= 1), default are specific for each ensemble
{‘NVT’: 0.3, ‘NPT’: 0.3, ‘GEMC_NVT’: 0.2, ‘GEMC_NPT’: 0.2, ‘GCMC’: 0.15} Fractional percentage at which part of the molecule will be deleted and then regrown using coupled- decoupled configurational-bias algorithm (i.e., fraction of molecular growth moves).

CrankShaftFreq: int or float (0 <= value <= 1), default are specific for each ensemble
{‘NVT’: 0.1, ‘NPT’: 0.1, ‘GEMC_NVT’: 0.1, ‘GEMC_NPT’: 0.1, ‘GCMC’: 0.1} Fractional percentage at which crankshaft move will occur. In this move, two atoms that are forming angle or dihedral are selected randomly and form a shaft. Then any atoms or group that are within these two selected atoms, will rotate around the shaft to sample intra-molecular degree of freedom (i.e., fraction of crankshaft moves).

VolFreq: int or float (0 <= value <= 1), default are specific for each ensemble
{‘NVT’: 0.0, ‘NPT’: 0.01, ‘GEMC_NVT’: 0.01, ‘GEMC_NPT’: 0.01, ‘GCMC’: 0.0} For isobaric-isothermal (NPT) ensemble and Gibbs ensemble (GEMC_NPT and GEMC_NVT) runs only: Fractional percentage at which a volume move will occur (i.e., fraction of Volume moves).

MultiParticleFreq: int or float (0 <= value <= 1), default are specific for each ensemble
{‘NVT’: 0.0, ‘NPT’: 0.0, ‘GEMC_NVT’: 0.0, ‘GEMC_NPT’: 0.0, ‘GCMC’: 0.0} Fractional percentage at which multi-particle move will occur. In this move, all molecules in the selected simulation box will be rigidly rotated or displaced simultaneously, along the calculated torque or force respectively (i.e., fraction of multi-particle moves). Note: Both MultiParticle (MultiParticleFreq) and MultiParticleBrownian (MultiParticleBrownianFreq) cannot be used at the same time.

MultiParticleBrownianFreq: int or float (0 <= value <= 1), default are specific for each ensemble
{‘NVT’: 0.0, ‘NPT’: 0.0, ‘GEMC_NVT’: 0.0, ‘GEMC_NPT’: 0.0, ‘GCMC’: 0.0} Fractional percentage at which multi-particle brownian move will occur. In this move, all molecules in the selected simulation box will be rigidly rotated or displaced simultaneously, along the calculated torque or force, respectively (i.e., fraction of Multi-Particle Brownian moves). Note: all of the move types are not available in for every ensemble. Note: all of the move fractions must sum to 1, or the control file writer will fail. Note: Both MultiParticle (MultiParticleFreq) and MultiParticleBrownian (MultiParticleBrownianFreq) cannot be used at the same time.

MultiParticleLiquid: bool, default=True (GEMC-only)
Use the multi-particle move in the liquid phase. Note: GOMC determines the boxes are liquid or gas/vapor before running each move, based on the calculated density of the GEMC boxes.

MultiParticleGas: bool, default=False (GEMC-only)
Use the multi-particle move in the gas/vapor phase. Note: GOMC determines the boxes are liquid or gas/vapor before

IntraMEMC-1Freq: int or float (0 <= value <= 1), default are specific for each ensemble
{‘NVT’: 0.0, ‘NPT’: 0.0, ‘GEMC_NVT’: 0.0, ‘GEMC_NPT’: 0.0, ‘GCMC’: 0.0} Fractional percentage at which specified number of small molecule kind will be exchanged with a specified large molecule kind in defined sub-volume within same simulation box. This move need additional information such as ExchangeVolumeDim, ExchangeRatio, ExchangeSmallKind, and ExchangeLargeKind.

MEMC-1Freq: int or float (0 <= value <= 1), default are specific for each ensemble
{‘NVT’: 0.0, ‘NPT’: 0.0, ‘GEMC_NVT’: 0.0, ‘GEMC_NPT’: 0.0, ‘GCMC’: 0.0} Fractional percentage at which specified number of small molecule kind will be exchanged with a specified large molecule kind in defined sub-volume, between simulation boxes. This move need additional information such as ExchangeVolumeDim, ExchangeRatio, ExchangeSmallKind, and ExchangeLargeKind.

IntraMEMC-2Freq: int or float (0 <= value <= 1), default are specific for each ensemble
{‘NVT’: 0.0, ‘NPT’: 0.0, ‘GEMC_NVT’: 0.0, ‘GEMC_NPT’: 0.0, ‘GCMC’: 0.0} Fractional percentage at which specified number of small molecule kind will be exchanged with a specified large molecule kind in defined sub-volume within same simulation box. Backbone of small and large molecule kind will be used to insert the large molecule more efficiently. This move need additional information such as ExchangeVolumeDim, ExchangeRatio, ExchangeSmallKind, ExchangeLargeKind, SmallKindBackBone, and LargeKindBackBone.

MEMC-2Freq: int or float (0 <= value <= 1), default are specific for each ensemble
{‘NVT’: 0.0, ‘NPT’: 0.0, ‘GEMC_NVT’: 0.0, ‘GEMC_NPT’: 0.0, ‘GCMC’: 0.0} Fractional percentage at which specified number of small molecule kind will be exchanged with a specified large molecule kind in defined sub-volume, between simulation boxes. Backbone of small and large molecule kind will be used to insert the large molecule more efficiently. This move need additional information such as ExchangeVolumeDim, ExchangeRatio, ExchangeSmallKind, ExchangeLargeKind, SmallKindBackBone, and LargeKindBackBone.

IntraMEMC-3Freq: int or float (0 <= value <= 1), default are specific for each ensemble
{‘NVT’: 0.0, ‘NPT’: 0.0, ‘GEMC_NVT’: 0.0, ‘GEMC_NPT’: 0.0, ‘GCMC’: 0.0} Fractional percentage at which specified number of small molecule kind will be exchanged with a specified large molecule kind in defined sub-volume within same simulation box. Specified atom of the large molecule kind will be used to insert the large molecule using coupled-decoupled configurational-bias. This move need additional information such as ExchangeVolumeDim, ExchangeRatio, ExchangeSmallKind, ExchangeLargeKind, and LargeKindBackBone.

MEMC-3Freq: int or float (0 <= value <= 1), default are specific for each ensemble
{‘NVT’: 0.0, ‘NPT’: 0.0, ‘GEMC_NVT’: 0.0, ‘GEMC_NPT’: 0.0, ‘GCMC’: 0.0} Fractional percentage at which specified number of small molecule kind will be exchanged with a specified large molecule kind in defined sub-volume, between simulation boxes. Specified atom of the large molecule kind will be used to insert the large molecule using coupled-decoupled configurational-bias. This move need additional information such as ExchangeVolumeDim, ExchangeRatio, ExchangeSmallKind, ExchangeLargeKind, and LargeKindBackBone.

TargetedSwapFreq: int or float (0 <= value <= 1), default are specific for each ensemble
{‘NVT’: 0.0, ‘NPT’: 0.0, ‘GEMC_NVT’: 0.0, ‘GEMC_NPT’: 0.0, ‘GCMC’: 0} Fractional percentage at which targeted swap move will occur. Note: This is only usable with the ‘GCMC’, ‘GEMC_NVT’, and ‘GEMC_NPT’ ensembles. Note: This is used in conjunction with the “TargetedSwap_DataInput” variables.

IntraTargetedSwapFreq: int or float (0 <= value <= 1), default are specific for each ensemble
{‘NVT’: 0.0, ‘NPT’: 0.0, ‘GEMC_NVT’: 0.0, ‘GEMC_NPT’: 0.0, ‘GCMC’: 0} Note: This is used in conjunction with the “TargetedSwap_DataInput” variables.

TargetedSwap_DataInput: dict, default=None
A dictionary which can contain one or several targeted swap regions, each designated with their own tag ID number (aka, subvolume number). A few examples for TargetedSwap_DataInput input is provided below. NOTE: MULTIPLE SIMULATION BOXES CAN BE UTILIZED BY SETTING MULTIPLE tag_ID_integer VALUES (INTEGERS VALUES), NOTE: THIS IS REQUIRED WHEN USING EITHER THE “TargetedSwapFreq” OR “IntraTargetedSwapFreq” MC MOVES. WARNING: THE tag_ID_integer VALUES MUST BE UNIQUE FOR EACH SUBVOLUME, OR THE DICTIONARY WILL OVERWRITE THE PREVIOUS SUBVOLUME (tag_ID_integer) SECTION WITH THE CURRENT tag_ID_integer AND ITS RESPECTIVE VALUES.

Example 1 input_variables_dict={“TargetedSwap_DataInput”: {tag_ID_integer: {“SubVolumeType”: “dynamic”, “SubVolumeBox”: 0, “SubVolumeCenterList”: [‘0-10’, 12, 15, ‘22-40’], “SubVolumeDim”: [1 * unyt.angstrom, 2 * unyt.angstrom, 3 * unyt.angstrom], “SubVolumeResidueKind”: “ALL”, “SubVolumeRigidSwap”: False, “SubVolumePBC”: “XY”, “SubVolumeChemPot”: {“MET”: -21 * unyt.kcal / unyt.mol, “met”: -31 * unyt.K}}}

Example 2 input_variables_dict={“TargetedSwap_DataInput”: {tag_ID_integer: {“SubVolumeType”: “static”, “SubVolumeBox”: 0, “SubVolumeCenter”: [1 , 12, 15] * unyt.angstrom, “SubVolumeDim”: [3, 3, 3] * 10**(-10) * unyt.m, “SubVolumeResidueKind”: [“MET”, “met”], “SubVolumeRigidSwap”: False, “SubVolumePBC”: “XYZ”, “SubVolumeFugacity”: {“MET”: 0.1 * unyt.bar, “met”: 1 * unyt.bar}}}

The details of each key and value for the “TargetedSwap_DataInput” are provided below.

— “SubVolumeType”: str (“static” or “dynamic”), No default is provided. The type of targeted swap box (subvolume) that will be created. The “static” type will maintain the box (subvolume) in a fixed location during the whole simulation, with the center of the box determined by the coordinates set in the “SubvolumeCenter” parameter. The “dynamic” type will allow for dynamic movement of the box (subvolume) based atom indices provided in the SubvolumeCenterList variable. For the “dynamic” type, the user must define a list of atom indices using “SubVolumeCenterList” keyword; the geometric center of the provided atom indices will be used as the center of subVolume. User must ensure that the atoms defined in the atom list remain in the simulation box (by setting the Beta value to 2 in PDB file).

— “SubVolumeBox”: int (0 or 1), No default is provided. The simulation box in which the targeted swap subvolume will be applied. NOTE: Only box zero (0) can be used for the GCMC, NPT, and NVT ensembles.

— “SubVolumeCenter”: list of three (3) unyt.unyt_quantity in length units, [x-axis, y-axis, z-axis], No default is provided. The simulation box is centered on this x, y, and z-axis points (written in Angstroms), which is only utilized when “SubVolumeType” is set to “static”.

— “SubVolumeCenterList”: list of int and/or str (>=0), [atom_index, …, atom_index], No default is provided. The simulation box subVolume is centered on the geometric center of the provided atom indices, which is only used when the “SubVolumeType” is set to “dynamic”. For example, [0-10’, 12, 15] means that atom indices 0 to 10, 12 and 15 are used as the geometric center of the simulation box subVolume. NOTE: THE ATOM INDICES RANGES MUST BE A STRING IN THE FORM ‘2-20’, WITH THE FIRST ATOM INDICES BEING SMALLER THAN THE SECOND (i.e, ‘a-b’, where a < b). ALL SINGULAR ATOM INDICES MUST BE INTEGERS. NOTE: THE SAME ATOM INDICES CAN BE USED 2, 3 OR X TIMES TO WEIGHT that atom 2, 3, OR X TIMES MORE IN THE GEOMETRIC CENTERING CALCULATION. NOTE: THE ATOM INDICES START AT ZERO (0), WHILE THE PDB AND PSF FILES START AT ONE (1). THEREFORE, YOU NEED TO BE CAREFUL WHEN SETTING THE INDICES FROM THE PDB OR PSF VALUES AS THEY ARE ONE (1) NUMBER OFF.

— “SubVolumeDim”: list of three (3) unyt.unyt_quantity in length units (>= 0), [x-axis, y-axis, z-axis], No default is provided. This sets the size of the simulation box (subVolume) in the x, y, and z-axis (written in Angstroms).

— “SubVolumeResidueKind”: str or list of str, “ALL” or “residue” or [“ALL”] or [residue_str, …, residue_str], No default is provided. The residues that will be used in the “TargetedSwap_DataInput” subvolume. Alternatively, the user can just set the value to [“ALL”] or “ALL”, which covers all the residues.

— “SubVolumeRigidSwap”: bool, default = True Choose whether to use a rigid or flexible molecule insertion using CD-CBMC for the subVolume. True uses a rigid molecule insertion, while False uses a flexible molecule insertion

— “SubVolumePBC”: str (‘X’, ‘XY’, ‘XZ’, ‘XYZ’, ‘Y’, ‘YZ’, or ‘Z’), default = ‘XYZ’. Apply periodic boundary condition (PBC) in selected axes for the subVolume. Example 1, ‘X’ applies PBC only in the X axis. Example 2, ‘XY’ applies PBC only in the X and Y axes. Example 3, ‘XYZ’ applies PBC in the X, Y, and Z axes.

— “SubVolumeChemPot”: dict {str (4 dig limit) , unyt.unyt_quantity in energy units (kcal/mol, kJ/mol, or K)}}, No default is provided. The chemical potentials are written in GOMC units of energy, K. If no SubVolumeChemPot is provided the default system chemical potential values are used. There is a 4 character limit for the string/residue name since the PDB/PSF files have a 4 character limitation and require an exact match in the conf file. Note: These strings must match the residue in the psf and psb files or it will fail. The name of the residues and their corresponding chemical potential must be specified for every residue in the system (i.e., {“residue_name”: chemical_potential}). Note: THIS IS ONLY REQUIRED FOR THE GCMC ENSEMBLE. Note: IF 2 KEYS WITH THE SAME STRING/RESIDUE ARE PROVIDED, ONE WILL BE AUTOMATICALLY OVERWRITTEN AND NO ERROR WILL BE THROWN IN THIS CONTROL FILE WRITER. Note: ONLY THE “SubVolumeChemPot” OR THE “SubVolumeFugacity” CAN BE USED FOR ALL THE TARGET SWAP BOXES (SUBVOLUMES). IF MIX OF “SubVolumeChemPot” AND “SubVolumeFugacity” ARE USED THE CONTROL FILE WRITER WILL THROW AN ERROR.

— “SubVolumeFugacity”: dict {str , iunyt.unyt_quantity in pressure units (>= 0)}, No default is provided. The fugacity is written in GOMC units of pressure, bar. If no “SubVolumeFugacity” is provided the default system fugacity values are used. There is a 4 character limit for the string/residue name since the PDB/PSF files have a 4 character limitation and require an exact match in the conf file. Note: These strings must match the residue in the psf and psb files or it will fail. The name of the residues and their corresponding fugacity must be specified for every residue in the system (i.e., {“residue_name”: fugacity}). Note: THIS IS ONLY REQUIRED FOR THE GCMC ENSEMBLE. Note: IF 2 KEYS WITH THE SAME STRING/RESIDUE ARE PROVIDED, ONE WILL BE AUTOMATICALLY OVERWRITTEN AND NO ERROR WILL BE THROWN IN THIS CONTROL FILE WRITER. Note: ONLY THE “SubVolumeChemPot” OR THE “SubVolumeFugacity” CAN BE USED FOR ALL THE TARGET SWAP BOXES (SUBVOLUMES). IF MIX OF “SubVolumeChemPot” AND “SubVolumeFugacity” ARE USED THE CONTROL FILE WRITER WILL THROW AN ERROR.

ExchangeVolumeDim: list of 3 floats or integers or [X-dimension, Y-dimension, Z-dimension)],
default = [1.0, 1.0, 1.0] To use all variations of MEMC and Intra-MEMC Monte Carlo moves, the exchange subvolume must be defined. The exchange sub-volume is defined as an orthogonal box with x, y, and z-dimensions, where small molecule/molecules kind will be selected from to be exchanged with a large molecule kind. Note: Currently, the X and Y dimension cannot be set independently (X = Y = max(X, Y)). Note: A heuristic for setting good values of the x, y, and z-dimensions is to use the geometric size of the large molecule plus 1-2 Å in each dimension. Note: In case of exchanging 1 small molecule kind with 1 large molecule kind in IntraMEMC-2, IntraMEMC-3, MEMC-2, MEMC-3 Monte Carlo moves, the sub-volume dimension has no effect on acceptance rate.

MEMC_DataInput: nested lists, default=None
Enter data as a list with some sub-lists as follows: [[ExchangeRatio_int (> 0), ExchangeLargeKind_str, [LargeKindBackBone_atom_1_str_or_NONE, LargeKindBackBone_atom_2_str_or_NONE ], ExchangeSmallKind_str, [SmallKindBackBone_atom_1_str_or_NONE, SmallKindBackBone_atom_2_str_or_NONE ]], …, [ExchangeRatio_int (> 0), ExchangeLargeKind_str, [LargeKindBackBone_atom_1_str_or_NONE, LargeKindBackBone_atom_2_str_or_NONE ], ExchangeSmallKind_str, [SmallKindBackBone_atom_1_str_or_NONE, SmallKindBackBone_atom_2_str_or_NONE ]. NOTE: CURRENTLY ALL THESE INPUTS NEED TO BE SPECIFIED, REGARDLESS OF THE MEMC TYPE SELECTION. IF THE SmallKindBackBone or LargeKindBackBone IS NOT REQUIRED FOR THE MEMC TYPE, None CAN BE USED IN PLACE OF A STRING.

Note: These strings must match the residue in the psf and psb files or it will fail. It is recommended that the user print the Charmm object psf and pdb files and review the residue names that match the atom name before using the in the MEMC_DataInput variable input.

Note: see the below data explanations for the ExchangeRatio, ExchangeSmallKind, ExchangeLargeKind, LargeKindBackBone, SmallKindBackBone.

Example 1 (MEMC-1): [ [1, ‘WAT’, [None, None], ‘wat’, [None, None]] , [1, ‘WAT’, [None, None], ‘wat’, [None, None]]

Example 2 (MEMC-2): [ [1, ‘WAT’, [‘O1’, ‘H1’], ‘wat’, [‘O1’, ‘H1’ ]] , [1, ‘WAT’, [‘H1’, ‘H2’], ‘wat’, [‘H1’, ‘H2’ ]]

Example 3 (MEMC-3): [ [2, ‘WAT’, ‘O1’, ‘H1’], ‘wat’, [None, None]] , [2, ‘WAT’, [‘H1’, ‘H2’], ‘wat’, [None, None]]

— ExchangeRatio = MEMC parameters (all ensembles): int (> 0), default=None The Ratio of exchanging small molecule/molecules with 1 large molecule. To use all variation of MEMC and Intra-MEMC Monte Carlo moves, the exchange ratio must be defined. The exchange ratio defines how many small molecule will be exchanged with 1 large molecule. For each large-small molecule pairs, one exchange ratio must be defined.

— ExchangeSmallKind = MEMC parameters (all ensembles): str, default=None The small molecule kind (resname) to be exchanged. Note: ONLY 4 characters can be used for the strings. To use all variation of MEMC and Intra-MEMC Monte Carlo moves, the small molecule kind to be exchanged with a large molecule kind must be defined. Multiple small molecule kind can be specified.

— ExchangeLargeKind = MEMC parameters (all ensembles): str, default=None The large molecule kind (resname) to be exchanged. Note: ONLY 4 characters can be used for the strings. To use all variation of MEMC and Intra-MEMC Monte Carlo moves, the large molecule kind to be exchanged with small molecule kind must be defined. Multiple large molecule kind can be specified.

— LargeKindBackBone = MEMC parameters (all ensembles): list [str, str] or [None, None], default=None Note: ONLY 4 characters can be used for the strings. The [None, None] values can only be used if that MEMC type does not require them. The strings for the the atom name 1 and atom name 2 that belong to the large molecule’s backbone (i.e., [str_for_atom_name_1, str_for_atom_name_2]) To use MEMC-2, MEMC-3, IntraMEMC-2, and IntraMEMC-3 Monte Carlo moves, the large molecule backbone must be defined. The backbone of the molecule is defined as a vector that connects two atoms belong to the large molecule. The large molecule backbone will be used to align the sub-volume in MEMC-2 and IntraMEMC-2 moves, while in MEMC-3 and IntraMEMC-3 moves, it uses the atom name to start growing the large molecule using coupled-decoupled configurational-bias. For each large-small molecule pairs, two atom names must be defined. Note: all atom names in the molecule must be unique. Note: In MEMC-3 and IntraMEMC-3 Monte Carlo moves, both atom names must be same, otherwise program will be terminated. Note: If the large molecule has only one atom (mono atomic molecules), same atom name must be used for str_for_atom_name_1 and str_for_atom_name_2 of the LargeKindBackBone.

— SmallKindBackBone = MEMC parameters (all ensembles): list [str, str] or [None, None], default=None Note: ONLY 4 characters can be used for the strings. The [None, None] values can only be used if that MEMC type does not require them. The strings for the the atom name 1 and atom name 2 that belong to the small molecule’s backbone (i.e., [str_for_atom_name_1, str_for_atom_name_2]) To use MEMC-2, and IntraMEMC-2 Monte Carlo moves, the small molecule backbone must be defined. The backbone of the molecule is defined as a vector that connects two atoms belong to the small molecule and will be used to align the sub-volume. For each large-small molecule pairs, two atom names must be defined. Note: all atom names in the molecule must be unique. Note: If the small molecule has only one atom (mono atomic molecules), same atom name must be used str_for_atom_name_1 and str_for_atom_name_2 of the SmallKindBackBone.

# *******************************************************************

# input_variables_dict options (keys and values) - (end)

# Note: the input_variables_dict keys are also attributes

# *******************************************************************
Returns:

If completed without errors: str
returns “GOMC_CONTROL_FILE_WRITTEN” when the GOMC input control file is writen

If completed with errors: None
Notes

The user input variables (input_variables_dict) and the specific ensembles.

The details of the required inputs for the selected ensembles can be found by running this python workbook,
>>> print_valid_required_input_variables('NVT', description = True)
which prints the required inputs with their subsection description for the selected ‘NVT’ ensemble (other ensembles can be set as well).

The details of the input variables for the selected ensembles can be found by running this python workbook,
>>> print_valid_ensemble_input_variables('NPT', description = True)
which prints the input variables with their subsection description for the selected ‘NPT’ ensemble (other ensembles can be set as well).

Note: The box units imported are in nm (standard MoSDeF units). The units for this writer are auto-scaled to Angstroms, so they can be directly used in the GOMC or NAMD engines.

Note: all of the move types are not available in for every ensemble.

Note: all of the move fractions must sum to 1, or the control file writer will fail.

Note: In GOMC, for Mie FFs the following is always true –> n_ij = (n_ii + n_jj)/2. Note: In GOMC, for Exp FFs the following is always true –> alpha_ij = (alpha_ii * alpha_jj)**0.5.

The input variables (input_variables_dict) and text extracted with permission from the GOMC manual version 2.60. Some of the text was modified from its original version. Cite: Potoff, Jeffrey; Schwiebert, Loren; et. al. GOMC Documentation. https://raw.githubusercontent.com/GOMC-WSU/GOMC/master/GOMC_Manual.pdf, 2021.
Attributes:
input_error: bool
This error is typically incurred from an error in the user input values. However, it could also be due to a bug, provided the user is inputting the data as this Class intends.

all_failed_input_List: list
A list of all the inputs that failed, but there may be some inputs that

ensemble_typ: str, [‘NVT’, ‘NPT’, ‘GEMC_NPT’, ‘GCMC-NVT’, ‘GCMC’]
The ensemble type of the simulation.

RunSteps: int (>0), must be an integer greater than zero.
Sets the total number of simulation steps.

Temperature: unyt.unyt_quantity in temperature units (> 1 K)
Temperature of system in unyt temperature units (written in Kelvin).

ff_psf_pdb_file_directory: str (optional), default=None (i.e., the current directory).
The full or relative directory added to the force field, psf, and pdb file names, created via the Charmm object.

check_input_files_exist: bool (default=True)
Check if the force field, psf, and pdb files exist. If input, it also checks if the Checkpoint, binCoordinates_box_0, extendedSystem_box_0, binVelocities_box_0, binCoordinates_box_1, extendedSystem_box_1, and binVelocities_box_1 files. If the files are checked and do not exist, the writer will throw a ValueError. True, check if the all the input files exist. False, do not check if the all the input files exist.

Restart: boolean, default = False
Determines whether to restart the simulation from restart file (*_restart.pdb and *_restart.psf) or not.

Checkpoint: list (see below), default = [False, ‘’]
If True: Checkpoint = [True, str] —> [True, checkpoint_filename] If False: Checkpoint = [False, str] —> [False, checkpoint_filename] Checkpoint files guarantees the trajectory files can be concatenated, along with ensuring no deviation from a single simulation’s results in a simulation which was interrupted and restarted. Determines whether to restart the simulation from checkpoint file or not. Restarting the simulation with would result in an identical outcome, as if previous simulation was continued. This is required for hybrid Monte-Carlo Molecular Dyanamics in open-ensembles (GCMC/GEMC) to concatenate trajectory files since the molecular transfers rearranges the order of the molecules. Checkpointing will ensure the molecules are loaded in the same order each cycle. NOTE: the checkpoint_filename with or without the path can be entered with the .chk extension, or without an extension; if it does not have the file extension, the .chk extension will then be added to the filename.

Parameters: str, (default=None)
Override all other force field directory and filename input with the correct extension (.inp or .par). Note: the default directory is the current directory with the Charmm object file name.

Coordinates_box_0: str, (default=None)
Override all other box 0 pdb directory and filename inputs with the correct extension. Note: the default directory is the current directory with the Charmm object file name.

Structure_box_0: str, (default=None)
Override all other box 0 psf directory and filename inputs with the correct extension. Note: the default directory is the current directory with the Charmm object file name.

Coordinates_box_1: str, (default=None)
Override all other box 1 pdb directory and filename inputs with the correct extension. Note: the default directory is the current directory with the Charmm object file name.

Structure_box_1: str, (default=None)
Override all other box 1 psf directory and filename inputs with the correct extension. Note: the default directory is the current directory with the Charmm object file name.

binCoordinates_box_0: str, (default=None)
The box 0 binary coordinate file is used only for restarting a GOMC simulation, which provides increased numerical accuracy.

extendedSystem_box_0: str, (default=None)
The box 0 vectors and origin file is used only for restarting a GOMC simulation.

binVelocities_box_0: str, (default=None)
The box 0 binary velocity file is used only for restarting a GOMC simulation, which provides increased numerical accuracy. These velocities are only passed thru GOMC since Monte Carlo simulations do not utilize any velocity information.

binCoordinates_box_1: str, (default=None)
The box 1 binary coordinate file is used only for restarting a GOMC simulation, which provides increased numerical accuracy.

extendedSystem_box_1: str, (default=None)
The box 1 vectors and origin file is used only for restarting a GOMC simulation.

binVelocities_box_1: str, (default=None)
The box 1 binary velocity file is used only for restarting a GOMC simulation, which provides increased numerical accuracy. These velocities are only passed thru GOMC since Monte Carlo simulations do not utilize any velocity information.

input_variables_dict: dict, default=None
These input variables are optional and override the default settings. Changing these variables likely required for more advanced systems. The details of the acceptable input variables for the selected ensembles can be found by running the code below in python,
>>> print_valid_ensemble_input_variables('GCMC', description = True)
which prints the input_variables with their subsection description for the selected ‘GCMC’ ensemble (other ensembles can be set as well). Example: input_variables_dict = {‘PRNG’: 123, ‘ParaTypeCHARMM’: True }
conf_filename: str
The name of the GOMC contol file, which will be created. The extension of the GOMC control file can be .conf, or no extension can be provided. If no extension is provided, this writer will automatically add the .conf extension to the provided string.

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))

electrostatic_1_4: float or int
The 1-4 electrostatic/coulombic scaling factor, which is the same for all the residues/molecules, regardless if differenct force fields are utilized.

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.

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)).

any input_variables_dict key: varies (see each input_variables_dict key and value)
Any of the input variables keys is also an Attribute and can be called the same way. Please see the input_variables_dict keys in the Parameters section above for all the available attributes.

NAMD Control File Writer

The NAMD control file writer is not currently available.

MoSDeF software functions

The MosDeF-GOMC was built on the MosDeF ecosystem. Therefore, please see the MoSDeF software documentation for details on the respective packages:

mBuild – A hierarchical, component based molecule builder
foyer – A package for atom-typing as well as applying and disseminating forcefields
GMSO – Flexible storage of chemical topology for molecular simulation