diff --git a/mbuild/conversion.py b/mbuild/conversion.py index 065460446..7e23ef54b 100644 --- a/mbuild/conversion.py +++ b/mbuild/conversion.py @@ -21,8 +21,6 @@ from mbuild.box import Box from mbuild.exceptions import MBuildError from mbuild.formats.json_formats import compound_from_json, compound_to_json -from mbuild.formats.par_writer import write_par -from mbuild.formats.xyz import read_xyz, write_xyz from mbuild.utils.io import ( has_gmso, has_mdtraj, @@ -401,26 +399,6 @@ def load_file( compound = compound_from_json(filename) return compound # Handle xyz file - # TODO 1.0: Get rid of this conditional and everything under it? xyz will be GMSO backend now - if extension == ".xyz" and not "top" in kwargs: - if coords_only: - tmp = read_xyz(filename) - if tmp.n_particles != compound.n_particles: - raise ValueError( - f"Number of atoms in {filename}" - f"does not match {compound}" - ) - ref_and_compound = zip( - tmp._particles(include_ports=False), - compound.particles(include_ports=False), - ) - for ref_particle, particle in ref_and_compound: - particle.pos = ref_particle.pos - else: - compound = read_xyz(filename, compound=compound) - elif extension == ".xyz" and "top" in kwargs: - backend = "mdtraj" - # Then gmso reader if backend == "gmso": top = gmso.Topology.load(filename=filename, **kwargs) @@ -985,9 +963,10 @@ def save( residues : str of list of str Labels of residues in the Compound. Residues are assigned by checking against Compound.name. - **kwargs #TODO 1.0: Update description here and the link to GMSO + **kwargs : dict, optional Depending on the file extension these will be passed to either - `write_gsd`, , `write_mcf`, or `parmed.Structure.save`. + GMSO's intenral writers or `parmed.Structure.save`. + See https://github.com/mosdef-hub/gmso/tree/main/gmso/formats See https://parmed.github.io/ParmEd/html/structobj/parmed.structure. Structure.html#parmed.structure.Structure.save @@ -1031,8 +1010,6 @@ def save( return # Savers supported by mbuild.formats - # TODO 1.0: Do we have a pdb writer anywhere? Right now, we use parmed - # TODO 1.0: GMSO can't save mol2 files, do we prioritize a mol2 writer, or continue using parmed backend here? savers = { ".gro": save_in_gmso, ".gsd": save_in_gmso, @@ -1047,17 +1024,7 @@ def save( saver = savers[extension] except KeyError: saver = None - # Provide a warning if rigid_ids are not sequential from 0 - if compound.contains_rigid: - unique_rigid_ids = sorted( - set([p.rigid_id for p in compound.rigid_particles()]) - ) - if max(unique_rigid_ids) != len(unique_rigid_ids) - 1: - warn("Unique rigid body IDs are not sequential starting from zero.") - if saver: # mBuild supported saver. - if extension == ".gsd": - kwargs["rigid_bodies"] = [p.rigid_id for p in compound.particles()] # Calling save_in_gmso saver( filename=filename, @@ -1079,9 +1046,21 @@ def save( structure.save(filename, overwrite=overwrite, **kwargs) -# TODO 1.0: Add doc strings, links, etc.. def save_in_gmso(compound, filename, box, overwrite, **kwargs): - """Convert to GMSO, call gmso writers.""" + """Convert to GMSO, call GMSO internal writers. + + Parameters + ---------- + compound : mbuild.compound.Compound, required. + The mBuild compound to convert to a GSMO topology. + box : mbuild.box.Box, required + The mBuild box to be converted to a gmso box, if different that compound.box + overwrite : bool, required + If `True` overwrite the file if it already exists. + **kwargs + Keyword arguments used in GMSO's savers. + See https://github.com/mosdef-hub/gmso/tree/main/gmso/formats + """ gmso_top = to_gmso(compound=compound, box=box) gmso_top.save(filename=filename, overwrite=overwrite, **kwargs) @@ -1930,7 +1909,6 @@ def to_gmso( """ from gmso.external.convert_mbuild import from_mbuild - # TODO: Pass in rigid body IDs here once added to GMSO return from_mbuild( compound=compound, box=box, diff --git a/mbuild/formats/par_writer.py b/mbuild/formats/par_writer.py deleted file mode 100644 index 26410efa9..000000000 --- a/mbuild/formats/par_writer.py +++ /dev/null @@ -1,201 +0,0 @@ -"""CHARMM Par format.""" - -import warnings - -__all__ = ["write_par"] - - -def write_par(structure, filename): - """Write CHARMM Par file given a parametrized structure. - - Notes - ----- - Follows format according to - https://www.ks.uiuc.edu/Training/Tutorials/namd/namd-tutorial-unix-html/ - node25.html - Furthermore, ParmEd should support writing CHARMM par, rtf, str files - by converting the parmed.Structure into parmed.CharmmParameterSet - - Parmed stores rmin/2 in "rmin" - """ - # ATOMS - with open(filename, "w") as f: - f.write("ATOMS\n") - unique_atoms = set() - for atom in structure.atoms: - unique_atoms.add((atom.atom_type.name, atom.atom_type.mass)) - for atom in unique_atoms: - f.write("MASS -1 {:8s} {:8.4f}\n".format(atom[0], atom[1])) - - f.write("\nBONDS\n") - unique_bonds = set() - for bond in structure.bonds: - unique_bonds.add( - ( - bond.atom1.atom_type.name, - bond.atom2.atom_type.name, - bond.type, - ) - ) - - for bond in unique_bonds: - f.write( - "{:8s} {:8s} {:.5f} {:.5f}\n".format( - bond[0], bond[1], bond[2].k, bond[2].req - ) - ) - - f.write("\nANGLES\n") - unique_angles = set() - unique_ubs = set() - for angle in structure.angles: - associated_ub = False - for ub in structure.urey_bradleys: - if ((angle.atom1, angle.atom3) == (ub.atom1, ub.atom2)) or ( - angle.atom3, - angle.atom1, - ) == (ub.atom1, ub.atom2): - unique_ubs.add( - ( - angle.atom1.atom_type.name, - angle.atom2.atom_type.name, - angle.atom3.atom_type.name, - angle.type, - ub.type, - ) - ) - associated_ub = True - - if not associated_ub: - unique_angles.add( - ( - angle.atom1.atom_type.name, - angle.atom2.atom_type.name, - angle.atom3.atom_type.name, - angle.type, - ) - ) - - for ub in unique_ubs: - f.write( - "{:8s} {:8s} {:8s} {:.5f} {:.5f} {:.5f} {:.5f}\n".format( - ub[0], - ub[1], - ub[2], - ub[3].k, - ub[3].theteq, - ub[4].k, - ub[4].req, - ) - ) - for angle in unique_angles: - f.write( - "{:8s} {:8s} {:8s} {:.5f} {:.5f}\n".format( - angle[0], angle[1], angle[2], angle[3].k, angle[3].theteq - ) - ) - - # These dihedrals need to be PeriodicTorsion Style (Charmm style) - if len(structure.rb_torsions) > 0: - warnings.warn("RB Torsions detected, but unsupported in par writer") - f.write("\nDIHEDRALS\n") - unique_dihedrals = set() - scnb = set() - for dihedral in structure.dihedrals: - if not dihedral.improper: - unique_dihedrals.add( - ( - dihedral.atom1.atom_type.name, - dihedral.atom2.atom_type.name, - dihedral.atom3.atom_type.name, - dihedral.atom4.atom_type.name, - dihedral.type, - ) - ) - scnb.add(dihedral.type.scnb) - else: - msg = ( - "AMBER-style improper detected between " - + "{} {} {} {}".format( - dihedral.atom1, - dihedral.atom2, - dihedral.atom3, - dihedral.atom4, - ) - + ", but unsupported in par writer" - ) - warnings.warn(msg) - - for dihedral in unique_dihedrals: - f.write( - "{:8s} {:8s} {:8s} {:8s} {:.5f} {:5d} {:.5f}\n".format( - dihedral[0], - dihedral[1], - dihedral[2], - dihedral[3], - dihedral[4].phi_k, - dihedral[4].per, - dihedral[4].phase, - ) - ) - - f.write("\nIMPROPER\n") - unique_impropers = set() - for improper in structure.impropers: - unique_impropers.add( - ( - improper.atom1.atom_type.name, - improper.atom2.atom_type.name, - improper.atom3.atom_type.name, - improper.atom4.atom_type.name, - improper.type, - ) - ) - for improper in unique_impropers: - f.write( - "{:8s} {:8s} {:8s} {:8s} {:.5f} {:5d} {:.5f}\n".format( - improper[2], - improper[0], - improper[1], - improper[3], - improper[4].psi_k, - 0, - improper[4].psi_eq, - ) - ) - - sc_nb = [a for a in scnb] - if len(sc_nb) > 1: - warnings.warn( - "Multiple 1-4 LJ scalings were detected, " - "defaulting to first LJ scaling detected, {}".format(sc_nb[0]) - ) - sc_nb = sc_nb[0] - elif len(sc_nb) == 1: - sc_nb = sc_nb[0] - elif len(sc_nb) == 0: - warnings.warn("No 1-4 LJ scaling was detected, defaulting 1") - sc_nb = 1.0 - - f.write("\nNONBONDED\n") - unique_atypes = set() - for atom in structure.atoms: - unique_atypes.add(atom.atom_type) - for atype in unique_atypes: - # atype, 0.0, epsilon, rmin/2, 0.0, epsilon(1-4), rmin/2 (1-4) - f.write( - "{:8s} {:8.3f} {:8.3f} {:8.3f} {:8.3f} {:8.3f} {:8.3f}\n".format( - atype.name, - 0.0, - -1 * atype.epsilon, - atype.rmin, - 0.0, - -1 * sc_nb * atype.epsilon, - atype.rmin, - ) - ) - - if structure.has_NBFIX(): - warnings.warn("NBFixes detected but unsupported in par writer") - - f.write("\nEND") diff --git a/mbuild/formats/xyz.py b/mbuild/formats/xyz.py deleted file mode 100644 index f4def4843..000000000 --- a/mbuild/formats/xyz.py +++ /dev/null @@ -1,134 +0,0 @@ -"""XYZ format.""" - -from warnings import warn - -import numpy as np -from ele import element_from_symbol -from ele.exceptions import ElementError - -import mbuild as mb -from mbuild.exceptions import MBuildError - -__all__ = ["read_xyz", "write_xyz"] - - -def read_xyz(filename, compound=None): - """Read an XYZ file. - - The expected format is as follows: - The first line contains the number of atoms in the file. The second line - contains a comment, which is not read. Remaining lines, one for each - atom in the file, include an elemental symbol followed by X, Y, and Z - coordinates in Angstroms. Columns are expected tbe separated by - whitespace. See https://openbabel.org/wiki/XYZ_(format). - - Parameters - ---------- - filename : str - Path of the input file - - Returns - ------- - compound : Compound - - Notes - ----- - The XYZ file format neglects many important details, including bonds, - residues, and box information. - - There are some other flavors of the XYZ file format and not all are - guaranteed to be compatible with this reader. For example, the TINKER - XYZ format is not expected to be properly read. - """ - if compound is None: - compound = mb.Compound() - - guessed_elements = set() - - compound_list = [] - with open(filename, "r") as xyz_file: - n_atoms = int(xyz_file.readline()) - xyz_file.readline() - coords = np.zeros(shape=(n_atoms, 3), dtype=np.float64) - for row, _ in enumerate(coords): - line = xyz_file.readline().split() - if not line: - msg = ( - "Incorrect number of lines in input file. Based on the " - "number in the first line of the file, {} rows of atoms " - "were expected, but at least one fewer was found." - ) - raise MBuildError(msg.format(n_atoms)) - coords[row] = line[1:4] - coords[row] *= 0.1 - name = line[0] - try: - element = name.capitalize() - element = element_from_symbol(element) - except ElementError: - if name not in guessed_elements: - warn( - "No matching element found for {}; the particle will " - "be added to the compound without an element " - "attribute.".format(name) - ) - guessed_elements.add(name) - element = None - - particle = mb.Compound(pos=coords[row], name=name, element=element) - compound_list.append(particle) - - # Verify we have read the last line by ensuring the next line in blank - line = xyz_file.readline().split() - if line: - msg = ( - "Incorrect number of lines in input file. Based on the " - "number in the first line of the file, {} rows of atoms " - "were expected, but at least one more was found." - ) - raise MBuildError(msg.format(n_atoms)) - - compound.add(compound_list) - - return compound - - -def write_xyz(structure, filename, write_atomnames=False, prec=6): - """Output an XYZ file. - - Parameters - ---------- - structure : parmed.Structure - ParmEd structure object - filename : str - Path of the output file - write_atomnames : bool - Write the `atom.name` attribute of the parmed structure - to the first column of the xyz file rather than the element - prec : int - The number of decimal places to write to the output file - - Notes - ----- - Coordinates are written in Angstroms. By default, the first column is the - element. These choices follow the conventions for the XYZ file format. - - The XYZ file format neglects many important details, notably as bonds, - residues, and box information. - """ - if isinstance(structure, mb.Compound): - raise ValueError("Expected a ParmEd structure, got an mbuild.Compound") - - xyz = np.array([[atom.xx, atom.xy, atom.xz] for atom in structure.atoms]) - if write_atomnames: - names = [atom.name for atom in structure.atoms] - else: - names = [atom.element_name for atom in structure.atoms] - - with open(filename, "w") as xyz_file: - xyz_file.write(str(len(structure.atoms))) - xyz_file.write("\n" + filename + " - created by mBuild\n") - for name, (x, y, z) in zip(names, xyz): - xyz_file.write( - f"{name} {x:11.{prec}f} {y:11.{prec}f} {z:11.{prec}f}\n" - )