From aabbab68d54d1772fa5ff25d42f1e938abd56452 Mon Sep 17 00:00:00 2001 From: Chris Jones <50423140+chrisjonesBSU@users.noreply.github.com> Date: Sat, 26 Oct 2024 15:56:41 -0600 Subject: [PATCH] Remove deprecated formats, update reader and writer backends (#1191) * remove hoomdxml * remove protobuf * remove hoomd3 force writers * change backend for mol2 and xyz to gmso * remove protobuf from dev environments * add some TODO notes * change behavior of xyz test * comment out test for elements * fix test based on expected behavior for xyz files, remove wrong n atoms tests, that is being tested for in GMSO test_xyz.py * remove lammpsdata and tests * remove duplicate benzene fixture that reads benzene from a mol2 file * update benzene naming to fix failing rigid tests * add hoomd_writer file that has GSD and snapshot writer * fix typo * add else statement * move to_hoomdsnapshot to conversion, add method in Compound class. Delete old snapshot and gsd writers * fix import * remove hoomd writer file, add TODO flags to conversion * add TODO notes, pass in overwrite to Topology.save(), comment out file formats that require typed systems * remove tests for gsd, hoomd, cassandra and par writers, move tests for loading xys into test_compound * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * add check for charges back in save method * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * add if statement before checking charge is neutral * remove .top from unit test * remove lammps and lammpsdata files from tests and support in conversion.py * Don't use gmso mol2 backend * remove gmso as backend for mol2 saver * remove Hoomd V3 test from CI * remove par writer and xyz reader/write * pin gmso to latest ver in window env file * add has_hoomd back to io.py, handle tests that use hoomd via gmso * fix syntax * remove unused import * update doc strings with new links, correct file types and methods * add saver mapping to compound.save doc strings --------- Co-authored-by: CalCraven Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> --- .github/workflows/CI.yaml | 6 - environment-dev-win.yml | 3 +- environment-dev.yml | 1 - mbuild/compound.py | 111 +- mbuild/conversion.py | 321 +++--- mbuild/formats/gsdwriter.py | 301 ------ mbuild/formats/hoomd_forcefield.py | 466 --------- mbuild/formats/hoomd_simulation.py | 467 --------- mbuild/formats/hoomd_snapshot.py | 622 ------------ mbuild/formats/hoomdxml.py | 433 -------- mbuild/formats/lammpsdata.py | 1522 ---------------------------- mbuild/formats/par_writer.py | 201 ---- mbuild/formats/protobuf.py | 195 ---- mbuild/formats/xyz.py | 134 --- mbuild/tests/base_test.py | 7 +- mbuild/tests/test_cassandramcf.py | 467 --------- mbuild/tests/test_compound.py | 198 +--- mbuild/tests/test_gsd.py | 320 ------ mbuild/tests/test_hoomd.py | 445 -------- mbuild/tests/test_lammpsdata.py | 742 -------------- mbuild/tests/test_par.py | 40 - mbuild/tests/test_protobuf.py | 18 - mbuild/tests/test_rigid.py | 5 +- mbuild/tests/test_xyz.py | 72 -- mbuild/utils/io.py | 30 +- 25 files changed, 241 insertions(+), 6886 deletions(-) delete mode 100644 mbuild/formats/gsdwriter.py delete mode 100644 mbuild/formats/hoomd_forcefield.py delete mode 100644 mbuild/formats/hoomd_simulation.py delete mode 100644 mbuild/formats/hoomd_snapshot.py delete mode 100644 mbuild/formats/hoomdxml.py delete mode 100644 mbuild/formats/lammpsdata.py delete mode 100644 mbuild/formats/par_writer.py delete mode 100644 mbuild/formats/protobuf.py delete mode 100644 mbuild/formats/xyz.py delete mode 100644 mbuild/tests/test_cassandramcf.py delete mode 100644 mbuild/tests/test_gsd.py delete mode 100644 mbuild/tests/test_hoomd.py delete mode 100644 mbuild/tests/test_lammpsdata.py delete mode 100644 mbuild/tests/test_par.py delete mode 100644 mbuild/tests/test_protobuf.py delete mode 100644 mbuild/tests/test_xyz.py diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml index 864831ffd..d2ce820de 100644 --- a/.github/workflows/CI.yaml +++ b/.github/workflows/CI.yaml @@ -54,12 +54,6 @@ jobs: - name: Test (OS -> ${{ matrix.os }} / Python -> ${{ matrix.python-version }}) run: python -m pytest -v --cov=mbuild --cov-report=xml --cov-append --cov-config=setup.cfg --color yes --pyargs mbuild - - name: Tests for HOOMD-blue (V3) - run: | - micromamba install -c conda-forge "hoomd>3" - python -m pytest -v --cov=mbuild --cov-report=xml --cov-append --cov-config=setup.cfg --color yes mbuild/tests/test_hoomd.py - if: runner.os == 'Linux' - - name: Upload Coverage Report uses: codecov/codecov-action@v4 with: diff --git a/environment-dev-win.yml b/environment-dev-win.yml index a95ac86a9..851d3e5b1 100644 --- a/environment-dev-win.yml +++ b/environment-dev-win.yml @@ -6,7 +6,7 @@ dependencies: - bump2version - codecov - ele - - gmso>=0.9.0 + - gmso>=0.12.4 - foyer>=0.11.0 - garnett>=0.7.1 - gsd>=2.9 @@ -20,7 +20,6 @@ dependencies: - parmed>=3.4.3 - pip - pre-commit - - protobuf - py3Dmol - pycifrw - pytest diff --git a/environment-dev.yml b/environment-dev.yml index 0f151bc9e..42591309f 100644 --- a/environment-dev.yml +++ b/environment-dev.yml @@ -21,7 +21,6 @@ dependencies: - parmed>=3.4.3 - pip - pre-commit - - protobuf - py3Dmol - pycifrw - pytest diff --git a/mbuild/compound.py b/mbuild/compound.py index 02eab970b..261f3b8ea 100644 --- a/mbuild/compound.py +++ b/mbuild/compound.py @@ -2007,7 +2007,6 @@ def _visualize_py3dmol( os.path.join(tmp_dir, "tmp.mol2"), include_ports=show_ports, overwrite=True, - parmed_kwargs={"infer_residues": False}, ) view = py3Dmol.view() @@ -2979,15 +2978,9 @@ def save( self, filename, include_ports=False, - forcefield_name=None, - forcefield_files=None, - forcefield_debug=False, box=None, overwrite=False, residues=None, - combining_rule="lorentz", - foyer_kwargs=None, - parmed_kwargs=None, **kwargs, ): """Save the Compound to a file. @@ -2997,22 +2990,11 @@ def save( filename : str Filesystem path in which to save the trajectory. The extension or prefix will be parsed and control the format. Supported extensions: - 'hoomdxml', 'gsd', 'gro', 'top', 'lammps', 'lmp', 'mcf', 'pdb', 'xyz', - 'json', 'mol2', 'sdf', 'psf'. See parmed/structure.py for more - information on savers. + 'gsd', 'gro', 'top', 'mcf', 'pdb', 'xyz', + 'json', 'mol2', 'sdf', 'psf'. See `mbuild.conversion.save()` + for more information about writer methods. include_ports : bool, optional, default=False Save ports contained within the compound. - forcefield_files : str, optional, default=None - Apply a forcefield to the output file using a forcefield provided - by the `foyer` package. - forcefield_name : str, optional, default=None - Apply a named forcefield to the output file using the `foyer` - package, e.g. 'oplsaa'. `Foyer forcefields - `_ - forcefield_debug : bool, optional, default=False - Choose verbosity level when applying a forcefield through `foyer`. - Specifically, when missing atom types in the forcefield xml file, - determine if the warning is condensed or verbose. box : mb.Box, optional, default=self.boundingbox (with buffer) Box information to be written to the output file. If 'None', a bounding box is used with 0.25nm buffers at each face to avoid @@ -3022,50 +3004,25 @@ def save( residues : str of list of str Labels of residues in the Compound. Residues are assigned by checking against Compound.name. - combining_rule : str, optional, default='lorentz' - Specify the combining rule for nonbonded interactions. Only relevant - when the `foyer` package is used to apply a forcefield. Valid - options are 'lorentz' and 'geometric', specifying Lorentz-Berthelot - and geometric combining rules respectively. - foyer_kwargs : dict, optional, default=None - Keyword arguments to provide to `foyer.Forcefield.apply`. - Depending on the file extension these will be passed to either - `write_gsd`, `write_hoomdxml`, `write_lammpsdata`, - `write_mcf`, or `parmed.Structure.save`. - See `parmed structure documentation - `_ - parmed_kwargs : dict, optional, default=None - Keyword arguments to provide to :meth:`mbuild.Compound.to_parmed` **kwargs + See `mbuild.conversion.save()`. Depending on the file extension these will be passed to either - `write_gsd`, `write_hoomdxml`, `write_lammpsdata`, `write_mcf`, or - `parmed.Structure.save`. + Parmed or GMSO backend writers See https://parmed.github.io/ParmEd/html/structobj/parmed.structure. - Structure.html#parmed.structure.Structure.save + Structure.html#parmed.structure.Structure.save and + https://github.com/mosdef-hub/gmso/tree/main/gmso/formats Other Parameters ---------------- ref_distance : float, optional, default=1.0 - Normalization factor used when saving to .gsd and .hoomdxml formats + Normalization factor used when saving to the .gsd format for converting distance values to reduced units. ref_energy : float, optional, default=1.0 - Normalization factor used when saving to .gsd and .hoomdxml formats + Normalization factor used when saving to the .gsd format for converting energy values to reduced units. ref_mass : float, optional, default=1.0 - Normalization factor used when saving to .gsd and .hoomdxml formats + Normalization factor used when saving to the .gsd format for converting mass values to reduced units. - atom_style: str, default='full' - Defines the style of atoms to be saved in a LAMMPS data file. The - following atom styles are currently supported: - 'full', 'atomic', 'charge', 'molecular' - See `LAMMPS atom style documentation - `_ for more - information. - unit_style: str, default='real' - Defines to unit style to be save in a LAMMPS data file. Defaults - to 'real' units. Current styles are supported: 'real', 'lj'. See - `LAMMPS unit style documentation_ - `_ for more information. Notes ----- @@ -3074,29 +3031,24 @@ def save( * filename * include_ports + The savers used for each supported file type are: + GMSO: .gro, .gsd, .data, .xyz, .mcf, .top + Parmed: .mol2, .pdb, .prmtop, .cif, .crd + PyBel: .sdf + See Also -------- - conversion.save : Main saver logic - formats.gsdwrite.write_gsd : Write to GSD format - formats.hoomdxml.write_hoomdxml : Write to Hoomd XML format - formats.xyzwriter.write_xyz : Write to XYZ format - formats.lammpsdata.write_lammpsdata : Write to LAMMPS data format - formats.cassandramcf.write_mcf : Write to Cassandra MCF format - formats.json_formats.compound_to_json : Write to a json file + mbuild.conversion.save : Main saver logic + mbuild.formats.cassandramcf.write_mcf : Write to Cassandra MCF format + mbuild.formats.json_formats.compound_to_json : Write to a json file """ conversion.save( - self, - filename, - include_ports, - forcefield_name, - forcefield_files, - forcefield_debug, - box, - overwrite, - residues, - combining_rule, - foyer_kwargs, - parmed_kwargs, + compound=self, + filename=filename, + include_ports=include_ports, + box=box, + overwrite=overwrite, + residues=residues, **kwargs, ) @@ -3246,6 +3198,21 @@ def to_gmso(self, **kwargs): """ return conversion.to_gmso(self, **kwargs) + def to_hoomdsnapshot(self, **kwargs): + """Create a HOOMD-Blue snapshot from an mBuild Compound. + + Parameters + ---------- + compound : mb.Compound + The mb.Compound to be converted. + + Returns + ------- + snapshot : gsd.hoomd.Frame + HOOMD-Blue compatible topology. + """ + return conversion.to_hoomdsnapshot(self, **kwargs) + # Interface to Trajectory for reading/writing .pdb and .mol2 files. # ----------------------------------------------------------------- def from_trajectory( diff --git a/mbuild/conversion.py b/mbuild/conversion.py index eb4361712..6d42091ab 100644 --- a/mbuild/conversion.py +++ b/mbuild/conversion.py @@ -7,6 +7,7 @@ from pathlib import Path from warnings import warn +import gmso import numpy as np import parmed as pmd from ele import ( @@ -19,12 +20,7 @@ import mbuild as mb from mbuild.box import Box from mbuild.exceptions import MBuildError -from mbuild.formats.gsdwriter import write_gsd -from mbuild.formats.hoomdxml import write_hoomdxml from mbuild.formats.json_formats import compound_from_json, compound_to_json -from mbuild.formats.lammpsdata import write_lammpsdata -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, @@ -59,9 +55,9 @@ def load( Parameters ---------- filename_or_object : str, mdtraj.Trajectory, parmed.Structure, - mbuild.Compound, pybel.Molecule, + mbuild.Compound, pybel.Molecule. Name of the file or topology from which to load atom and bond - information. + information, or a valid SMILES string. relative_to_module : str, optional, default=None Instead of looking in the current working directory, look for the file where this module is defined. This is typically used in Compound @@ -80,7 +76,7 @@ def load( smiles: bool, optional, default=False Use RDKit or OpenBabel to parse filename as a SMILES string or file containing a SMILES string. If this is set to True, `rdkit` is the - default backend. + default backend and `filename_or_object` should be the SMILES string. infer_hierarchy : bool, optional, default=True If True, infer hierarchy from chains and residues ignore_box_warn : bool, optional, default=False @@ -172,21 +168,20 @@ def load_object( mb.Compound """ # Create type_dict type -> loading method - # Will need to add a gmso method soon - type_dict = {pmd.Structure: from_parmed} - if has_gmso: - gmso = import_("gmso") - type_dict.update({gmso.Topology: from_gmso}) + type_dict = { + pmd.Structure: from_parmed, + gmso.Topology: from_gmso, + } if has_openbabel: pybel = import_("pybel") type_dict.update({pybel.Molecule: from_pybel}) - if has_mdtraj: md = import_("mdtraj") type_dict.update({md.Trajectory: from_trajectory}) # Check if the given object is an mb.Compound + # TODO 1.0: What is the use case for this? if isinstance(obj, mb.Compound): if not compound: warn("Given object is already an mb.Compound, doing nothing.") @@ -210,7 +205,7 @@ def load_object( return compound # If nothing is return raise an error - raise ValueError(f"Object of type {type(obj).__name__} is not supported") + raise ValueError(f"Object of type {type(obj).__name__} is not supported.") def load_smiles( @@ -374,13 +369,11 @@ def load_file( # Need to come up with a different dict structure default_backends = { ".json": "internal", - ".xyz": "internal", + ".xyz": "gmso", ".sdf": "pybel", - ".hoomdxml": "mdtraj", - ".mol2": "mdtraj", + ".mol2": "gmso", ".pdb": "mdtraj", } - # Handle mbuild *.py files containing a class that wraps a structure file # in its own folder. E.g., you build a system from ~/foo.py and it imports # from ~/bar/baz.py where baz.py loads ~/bar/baz.pdb. @@ -391,11 +384,9 @@ def load_file( extension = Path(filename).suffix if not backend: - try: - # Try matching backend based on extension + try: # Try matching backend based on extension backend = default_backends[extension] - except KeyError: - # Else use default backend + except KeyError: # Else use default backend backend = "mdtraj" if has_mdtraj else "parmed" # First check internal readers @@ -406,29 +397,8 @@ def load_file( compound = compound_from_json(filename) return compound # Handle xyz file - 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": - gmso = import_("gmso") - top = gmso.Topology.load(filename=filename, **kwargs) compound = from_gmso( topology=top, @@ -443,8 +413,7 @@ def load_file( if extension == ".sdf": pybel_mol = pybel.readfile("sdf", filename) # pybel returns a generator, so we grab the first molecule of a - # list of len 1. - # Raise ValueError if there are more molecules + # list of len 1. Raise ValueError if there are more molecules pybel_mol = [i for i in pybel_mol] if len(pybel_mol) == 1: compound = from_pybel( @@ -495,7 +464,6 @@ def load_file( # by the corresponding backend if rigid: compound.label_rigid_bodies() - return compound @@ -967,15 +935,9 @@ def save( compound, filename, include_ports=False, - forcefield_name=None, - forcefield_files=None, - forcefield_debug=False, box=None, overwrite=False, residues=None, - combining_rule=None, - foyer_kwargs=None, - parmed_kwargs=None, **kwargs, ): """Save the Compound to a file. @@ -987,22 +949,10 @@ def save( filename : str Filesystem path in which to save the trajectory. The extension or prefix will be parsed and control the format. Supported extensions are: - 'hoomdxml', 'gsd', 'gro', 'top', 'lammps', 'lmp', 'mcf', 'xyz', 'pdb', - 'sdf', 'mol2', 'psf'. See parmed/structure.py for more information on - savers. + 'gsd', 'gro', 'top', 'mcf', 'xyz', 'pdb', 'sdf', 'mol2', 'psf'. + See parmed/structure.py for more information on savers. include_ports : bool, optional, default=False Save ports contained within the compound. - forcefield_files : str, optional, default=None - Apply a forcefield to the output file using a forcefield provided by the - `foyer` package. - forcefield_name : str, optional, default=None - Apply a named forcefield to the output file using the `foyer` package, - e.g. 'oplsaa'. Forcefields listed here: - https://github.com/mosdef-hub/foyer/tree/master/foyer/forcefields - forcefield_debug : bool, optional, default=False - Choose level of verbosity when applying a forcefield through `foyer`. - Specifically, when missing atom types in the forcefield xml file, - determine if the warning is condensed or verbose. box : mb.Box, optional, default=compound.boundingbox (with buffer) Box information to be written to the output file. If 'None', a bounding box is used with 0.25nm buffers at each face to avoid overlapping atoms. @@ -1011,44 +961,24 @@ def save( residues : str of list of str Labels of residues in the Compound. Residues are assigned by checking against Compound.name. - combining_rule : str or None, optional, default=None - Specify the combining rule for nonbonded interactions. Only relevant - when the `foyer` package is used to apply a forcefield. Valid options - are 'lorentz' and 'geometric' and None, specifying Lorentz-Berthelot and - geometric combining rules respectively. If None is provided, the combining_rule - specified from the forcefield will be used. - foyer_kwargs : dict, optional, default=None - Keyword arguments to provide to `foyer.Forcefield.apply`. - parmed_kwargs : dict, optional, default=None - Keyword arguments to provide to :meth:`mbuild.Compound.to_parmed` - **kwargs + **kwargs : dict, optional Depending on the file extension these will be passed to either - `write_gsd`, `write_hoomdxml`, `write_lammpsdata`, `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 Other Parameters ---------------- ref_distance : float, optional, default=1.0 - Normalization factor used when saving to .gsd and .hoomdxml formats for + Normalization factor used when saving to .gsd formats for converting distance values to reduced units. ref_energy : float, optional, default=1.0 - Normalization factor used when saving to .gsd and .hoomdxml formats for + Normalization factor used when saving to .gsd formats for converting energy values to reduced units. ref_mass : float, optional, default=1.0 - Normalization factor used when saving to .gsd and .hoomdxml formats for + Normalization factor used when saving to .gsd formats for converting mass values to reduced units. - atom_style: str, default='full' - Defines the style of atoms to be saved in a LAMMPS data file. The - following atom styles are currently supported: 'full', 'atomic', - 'charge', 'molecular'. See http://lammps.sandia.gov/doc/atom_style.html - for more information on atom styles. - unit_style: str, default='real' - Defines to unit style to be save in a LAMMPS data file. Defaults to - 'real' units. Current styles are supported: 'real', 'lj'. See - https://lammps.sandia.gov/doc/99/units.html for more information on - unit styles Notes ----- @@ -1058,15 +988,19 @@ def save( See Also -------- - formats.gsdwriter.write_gsd : Write to GSD format - formats.hoomdxml.write_hoomdxml : Write to Hoomd XML format - formats.xyzwriter.write_xyz : Write to XYZ format - formats.lammpsdata.write_lammpsdata : Write to LAMMPS data format formats.cassandramcf.write_mcf : Write to Cassandra MCF format formats.json_formats.compound_to_json : Write to a json file """ - extension = os.path.splitext(filename)[-1] + if os.path.exists(filename) and not overwrite: + raise IOError(f"{filename} exists; not overwriting") + if compound.charge: + if round(compound.charge, 4) != 0.0: + warn( + f"System is not charge neutral. Total charge is {compound.charge}." + ) + extension = os.path.splitext(filename)[-1] + # Keep json stuff with internal mbuild method if extension == ".json": compound_to_json( compound, file_path=filename, include_ports=include_ports @@ -1075,98 +1009,60 @@ def save( # Savers supported by mbuild.formats savers = { - ".hoomdxml": write_hoomdxml, - ".gsd": write_gsd, - ".xyz": write_xyz, - ".lammps": write_lammpsdata, - ".lmp": write_lammpsdata, - ".par": write_par, + ".gro": save_in_gmso, + ".gsd": save_in_gmso, + ".data": save_in_gmso, + ".xyz": save_in_gmso, + # ".mol2": save_in_gmso, + ".mcf": save_in_gmso, + ".top": save_in_gmso, } - if has_networkx: - from mbuild.formats.cassandramcf import write_mcf - - savers.update({".mcf": write_mcf}) try: saver = savers[extension] except KeyError: saver = None - - if os.path.exists(filename) and not overwrite: - raise IOError(f"{filename} exists; not overwriting") - - if not parmed_kwargs: - parmed_kwargs = {} - structure = compound.to_parmed( - box=box, - residues=residues, - include_ports=include_ports, - **parmed_kwargs, - ) - # Apply a force field with foyer if specified - if forcefield_name or forcefield_files: - foyer = import_("foyer") - ff = foyer.Forcefield( - forcefield_files=forcefield_files, - name=forcefield_name, - debug=forcefield_debug, - ) - - if not foyer_kwargs: - foyer_kwargs = {} - structure = ff.apply(structure, **foyer_kwargs) - - if combining_rule is None: - combining_rule = structure.combining_rule - - if structure.combining_rule != combining_rule: - warn( - f"Overwriting forcefield-specified combining rule ({combining_rule})" - f"to new combining rule ({combining_rule})." - "This can cause inconsistent between the 1-4 pair interactions," - "calculated in foyer, and the new combining rule." - "Consider directly changing the metadata of the Forcefield." - ) - - structure.combining_rule = combining_rule - if structure.__dict__.get("defaults"): - structure.defaults.comb_rule = ( - 2 if combining_rule == "lorentz" else 3 - ) - - total_charge = sum([atom.charge for atom in structure]) - if round(total_charge, 4) != 0.0: - warn(f"System is not charge neutral. Total charge is {total_charge}.") - - # 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 in [".gsd", ".hoomdxml"]: - kwargs["rigid_bodies"] = [p.rigid_id for p in compound.particles()] - saver(filename=filename, structure=structure, **kwargs) + # Calling save_in_gmso + saver( + filename=filename, + compound=compound, + box=box, + overwrite=overwrite, + **kwargs, + ) elif extension == ".sdf": pybel = import_("pybel") - new_compound = mb.Compound() - # Convert pmd.Structure to mb.Compound - new_compound.from_parmed(structure) - # Convert mb.Compound to pybel molecule - pybel_molecule = new_compound.to_pybel() + pybel_molecule = compound.to_pybel() # Write out pybel molecule to SDF file output_sdf = pybel.Outputfile("sdf", filename, overwrite=overwrite) output_sdf.write(pybel_molecule) output_sdf.close() - else: # ParmEd supported saver. + structure = compound.to_parmed() structure.save(filename, overwrite=overwrite, **kwargs) +def save_in_gmso(compound, filename, box, overwrite, **kwargs): + """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) + + def catalog_bondgraph_type(compound, bond_graph=None): """Identify type of subgraph found at this stage of the compound. @@ -1298,6 +1194,69 @@ def pull_residues( return residuesList +def to_hoomdsnapshot( + compound, + identify_connections=True, + ref_distance=1.0, + ref_mass=1.0, + rigid_bodies=None, + shift_coords=True, + write_special_pairs=True, + **kwargs, +): + """Output a gsd.hoomd.Frame (HOOMD-Blue topology format). + + Parameters + ---------- + compound : mb.Compound + mBuild compound to save to the GSD format. + identify_connections : bool + If `True`, then infer angles and dihedrals in the topology. + ref_distance : float, optional, default=1.0 + Reference distance for conversion to reduced units + ref_mass : float, optional, default=1.0 + Reference mass for conversion to reduced units + rigid_bodies : list of int, optional, default=None + List of rigid body information. An integer value is required for each + atom corresponding to the index of the rigid body the particle is to be + associated with. A value of None indicates the atom is not part of a + rigid body. + shift_coords : bool, optional, default=True + Shift coordinates from (0, L) to (-L/2, L/2) if necessary. + write_special_pairs : bool, optional, default=True + Writes out special pair information necessary to correctly use the OPLS + fudged 1,4 interactions in HOOMD-Blue. + + Notes + ----- + See gmso.external.convert_hoomd on how to make sure the GSD file contains + forcefield information (e.g. atom types, angle types, etc.). + """ + import unyt as u + from gmso.external import from_mbuild, to_gsd_snapshot + + gmso_top = from_mbuild(compound=compound) + + if identify_connections: + gmso_top.identify_connections() + + base_units = { + "length": ref_distance * u.Unit("nm"), + "mass": ref_mass * u.Unit("amu"), + "energy": 1 * u.Unit("kJ/mol"), + } + + snapshot, refs = to_gsd_snapshot( + top=gmso_top, + base_units=base_units, + rigid_bodies=rigid_bodies, + shift_coords=shift_coords, + parse_special_pairs=write_special_pairs, + auto_scale=False, + ) + return snapshot + + def to_parmed( compound, box=None, @@ -1811,7 +1770,6 @@ def to_rdkit(compound): for particle in compound.particles(): temp_atom = Chem.Atom(particle.element.atomic_number) - # this next line is necessary to prevent rdkit from adding hydrogens # this will also set the label to be the element with particle index temp_atom.SetProp( @@ -1924,7 +1882,14 @@ def _iterate_children(compound, nodes, edges, names_only=False): return nodes, edges -def to_gmso(compound, box=None, **kwargs): +def to_gmso( + compound, + box=None, + parse_label=True, + custom_groups=None, + infer_elements=False, + **kwargs, +): """Create a GMSO Topology from a mBuild Compound. Parameters @@ -1941,7 +1906,13 @@ def to_gmso(compound, box=None, **kwargs): """ from gmso.external.convert_mbuild import from_mbuild - return from_mbuild(compound, box=None, **kwargs) + return from_mbuild( + compound=compound, + box=box, + parse_label=parse_label, + custom_groups=custom_groups, + infer_elements=infer_elements, + ) def to_intermol(compound, molecule_types=None): # pragma: no cover diff --git a/mbuild/formats/gsdwriter.py b/mbuild/formats/gsdwriter.py deleted file mode 100644 index 2de3d5cd3..000000000 --- a/mbuild/formats/gsdwriter.py +++ /dev/null @@ -1,301 +0,0 @@ -"""GSD format. - -https://gsd.readthedocs.io/en/stable/ -""" - -import numpy as np - -from mbuild.utils.geometry import coord_shift -from mbuild.utils.io import import_ -from mbuild.utils.sorting import natural_sort - -__all__ = ["write_gsd"] - - -def write_gsd( - structure, - filename, - ref_distance=1.0, - ref_mass=1.0, - ref_energy=1.0, - rigid_bodies=None, - shift_coords=True, - write_special_pairs=True, - **kwargs -): - """Output a GSD file (HOOMD v2 default data format). - - Parameters - ---------- - structure : parmed.Structure - ParmEd Structure object - filename : str - Path of the output file. - ref_distance : float, optional, default=1.0 - Reference distance for conversion to reduced units - ref_mass : float, optional, default=1.0 - Reference mass for conversion to reduced units - ref_energy : float, optional, default=1.0 - Reference energy for conversion to reduced units - rigid_bodies : list of int, optional, default=None - List of rigid body information. An integer value is required for each - atom corresponding to the index of the rigid body the particle is to be - associated with. A value of None indicates the atom is not part of a - rigid body. - shift_coords : bool, optional, default=True - Shift coordinates from (0, L) to (-L/2, L/2) if necessary. - write_special_pairs : bool, optional, default=True - Writes out special pair information necessary to correctly use the OPLS - fudged 1,4 interactions - in HOOMD. - - Notes - ----- - Force field parameters are not written to the GSD file and must be included - manually into a HOOMD input script. - """ - import_("gsd") - import gsd.hoomd - - xyz = np.array([[atom.xx, atom.xy, atom.xz] for atom in structure.atoms]) - if shift_coords: - xyz = coord_shift(xyz, structure.box[:3]) - - gsd_snapshot = gsd.hoomd.Frame() - - gsd_snapshot.configuration.step = 0 - gsd_snapshot.configuration.dimensions = 3 - - # Write box information - if np.allclose(structure.box[3:6], np.array([90, 90, 90])): - gsd_snapshot.configuration.box = np.hstack( - (structure.box[:3] / ref_distance, np.zeros(3)) - ) - else: - a, b, c = structure.box[0:3] / ref_distance - alpha, beta, gamma = np.radians(structure.box[3:6]) - - lx = a - xy = b * np.cos(gamma) - xz = c * np.cos(beta) - ly = np.sqrt(b**2 - xy**2) - yz = (b * c * np.cos(alpha) - xy * xz) / ly - lz = np.sqrt(c**2 - xz**2 - yz**2) - - gsd_snapshot.configuration.box = np.array([lx, ly, lz, xy, xz, yz]) - - _write_particle_information( - gsd_snapshot, - structure, - xyz, - ref_distance, - ref_mass, - ref_energy, - rigid_bodies, - ) - if write_special_pairs: - _write_pair_information(gsd_snapshot, structure) - if structure.bonds: - _write_bond_information(gsd_snapshot, structure) - if structure.angles: - _write_angle_information(gsd_snapshot, structure) - if structure.rb_torsions: - _write_dihedral_information(gsd_snapshot, structure) - - with gsd.hoomd.open(filename, mode="w") as gsd_file: - gsd_file.append(gsd_snapshot) - - -def _write_particle_information( - gsd_snapshot, - structure, - xyz, - ref_distance, - ref_mass, - ref_energy, - rigid_bodies, -): - """Write out the particle information.""" - gsd_snapshot.particles.N = len(structure.atoms) - gsd_snapshot.particles.position = xyz / ref_distance - - types = [a.name if a.type == "" else a.type for a in structure.atoms] - - unique_types = list(set(types)) - unique_types.sort(key=natural_sort) - gsd_snapshot.particles.types = unique_types - - typeids = np.array([unique_types.index(t) for t in types]) - gsd_snapshot.particles.typeid = typeids - - masses = np.array([atom.mass for atom in structure.atoms]) - masses[masses == 0] = 1.0 - gsd_snapshot.particles.mass = masses / ref_mass - - charges = np.array([atom.charge for atom in structure.atoms]) - e0 = 2.396452e-04 - """ - Permittivity of free space = 2.396452e-04 e^2/((kcal/mol) Angstrom), - where e is the elementary charge - """ - charge_factor = (4.0 * np.pi * e0 * ref_distance * ref_energy) ** 0.5 - gsd_snapshot.particles.charge = charges / charge_factor - - if rigid_bodies: - rigid_bodies = [-1 if body is None else body for body in rigid_bodies] - gsd_snapshot.particles.body = rigid_bodies - - -def _write_pair_information(gsd_snapshot, structure): - """Write the special pairs in the system. - - Parameters - ---------- - gsd_snapshot : gsd.snapshot, - The file object of the GSD file being written - structure : parmed.Structure - Parmed structure object holding system information - """ - pair_types = [] - pair_typeid = [] - pairs = [] - for ai in structure.atoms: - for aj in ai.dihedral_partners: - # make sure we don't double add - if ai.idx > aj.idx: - ps = "-".join(sorted([ai.type, aj.type], key=natural_sort)) - if ps not in pair_types: - pair_types.append(ps) - pair_typeid.append(pair_types.index(ps)) - pairs.append((ai.idx, aj.idx)) - gsd_snapshot.pairs.types = pair_types - gsd_snapshot.pairs.typeid = pair_typeid - gsd_snapshot.pairs.group = pairs - gsd_snapshot.pairs.N = len(pairs) - - -def _write_bond_information(gsd_snapshot, structure): - """Write the bonds in the system. - - Parameters - ---------- - gsd_snapshot : - The file object of the GSD file being written - structure : parmed.Structure - Parmed structure object holding system information - """ - gsd_snapshot.bonds.N = len(structure.bonds) - - unique_bond_types = set() - for bond in structure.bonds: - t1, t2 = bond.atom1.type, bond.atom2.type - if t1 == "" or t2 == "": - t1, t2 = bond.atom1.name, bond.atom2.name - t1, t2 = sorted([t1, t2], key=natural_sort) - try: - bond_type = "-".join((t1, t2)) - except AttributeError: # no forcefield applied, bond.type is None - bond_type = ("-".join((t1, t2)), 0.0, 0.0) - unique_bond_types.add(bond_type) - unique_bond_types = sorted(list(unique_bond_types), key=natural_sort) - gsd_snapshot.bonds.types = unique_bond_types - - bond_typeids = [] - bond_groups = [] - for bond in structure.bonds: - t1, t2 = bond.atom1.type, bond.atom2.type - if t1 == "" or t2 == "": - t1, t2 = bond.atom1.name, bond.atom2.name - t1, t2 = sorted([t1, t2], key=natural_sort) - try: - bond_type = "-".join((t1, t2)) - except AttributeError: # no forcefield applied, bond.type is None - bond_type = ("-".join((t1, t2)), 0.0, 0.0) - bond_typeids.append(unique_bond_types.index(bond_type)) - bond_groups.append((bond.atom1.idx, bond.atom2.idx)) - - gsd_snapshot.bonds.typeid = bond_typeids - gsd_snapshot.bonds.group = bond_groups - - -def _write_angle_information(gsd_snapshot, structure): - """Write the angles in the system. - - Parameters - ---------- - gsd_snapshot : - The file object of the GSD file being written - structure : parmed.Structure - Parmed structure object holding system information - """ - gsd_snapshot.angles.N = len(structure.angles) - - unique_angle_types = set() - for angle in structure.angles: - t1, t2, t3 = angle.atom1.type, angle.atom2.type, angle.atom3.type - t1, t3 = sorted([t1, t3], key=natural_sort) - angle_type = "-".join((t1, t2, t3)) - unique_angle_types.add(angle_type) - unique_angle_types = sorted(list(unique_angle_types), key=natural_sort) - gsd_snapshot.angles.types = unique_angle_types - - angle_typeids = [] - angle_groups = [] - for angle in structure.angles: - t1, t2, t3 = angle.atom1.type, angle.atom2.type, angle.atom3.type - t1, t3 = sorted([t1, t3], key=natural_sort) - angle_type = "-".join((t1, t2, t3)) - angle_typeids.append(unique_angle_types.index(angle_type)) - angle_groups.append((angle.atom1.idx, angle.atom2.idx, angle.atom3.idx)) - - gsd_snapshot.angles.typeid = angle_typeids - gsd_snapshot.angles.group = angle_groups - - -def _write_dihedral_information(gsd_snapshot, structure): - """Write the dihedrals in the system. - - Parameters - ---------- - gsd_snapshot : - The file object of the GSD file being written - structure : parmed.Structure - Parmed structure object holding system information - """ - gsd_snapshot.dihedrals.N = len(structure.rb_torsions) - - unique_dihedral_types = set() - for dihedral in structure.rb_torsions: - t1, t2 = dihedral.atom1.type, dihedral.atom2.type - t3, t4 = dihedral.atom3.type, dihedral.atom4.type - if [t2, t3] == sorted([t2, t3], key=natural_sort): - dihedral_type = "-".join((t1, t2, t3, t4)) - else: - dihedral_type = "-".join((t4, t3, t2, t1)) - unique_dihedral_types.add(dihedral_type) - unique_dihedral_types = sorted( - list(unique_dihedral_types), key=natural_sort - ) - gsd_snapshot.dihedrals.types = unique_dihedral_types - - dihedral_typeids = [] - dihedral_groups = [] - for dihedral in structure.rb_torsions: - t1, t2 = dihedral.atom1.type, dihedral.atom2.type - t3, t4 = dihedral.atom3.type, dihedral.atom4.type - if [t2, t3] == sorted([t2, t3], key=natural_sort): - dihedral_type = "-".join((t1, t2, t3, t4)) - else: - dihedral_type = "-".join((t4, t3, t2, t1)) - dihedral_typeids.append(unique_dihedral_types.index(dihedral_type)) - dihedral_groups.append( - ( - dihedral.atom1.idx, - dihedral.atom2.idx, - dihedral.atom3.idx, - dihedral.atom4.idx, - ) - ) - - gsd_snapshot.dihedrals.typeid = dihedral_typeids - gsd_snapshot.dihedrals.group = dihedral_groups diff --git a/mbuild/formats/hoomd_forcefield.py b/mbuild/formats/hoomd_forcefield.py deleted file mode 100644 index 265947785..000000000 --- a/mbuild/formats/hoomd_forcefield.py +++ /dev/null @@ -1,466 +0,0 @@ -"""HOOMD v3 forcefield format.""" - -import itertools -import operator -import warnings -from collections import namedtuple - -import numpy as np -import parmed as pmd - -import mbuild as mb -from mbuild.utils.conversion import RB_to_OPLS -from mbuild.utils.decorators import deprecated -from mbuild.utils.io import import_ -from mbuild.utils.sorting import natural_sort - -from .hoomd_snapshot import _get_hoomd_version, to_hoomdsnapshot - -hoomd = import_("hoomd") - - -dep_msg = """ -Support for Hoomd-Blue writers will be removed in mBuild 1.0. -See GMSO (https://github.com/mosdef-hub/gmso/tree/main/gmso/external/convert_hoomd) for -Hoomd-Blue 3.x and 4.x format support. -To convert to GMSO, use the method `Compound.to_gmso()`. -""" -print(dep_msg) - - -@deprecated(dep_msg) -def create_hoomd_forcefield( - structure, - r_cut, - ref_distance=1.0, - ref_mass=1.0, - ref_energy=1.0, - auto_scale=False, - nlist_buffer=0.4, - snapshot_kwargs={}, - pppm_kwargs={"Nx": 8, "Ny": 8, "Nz": 8, "order": 4}, - init_snap=None, -): - """Convert a parametrized pmd.Structure to a HOOMD snapshot and forces. - - Parameters - ---------- - structure : parmed.Structure - ParmEd Structure object - r_cut : float - Cutoff radius in simulation units - ref_distance : float, optional, default=1.0 - Reference distance for unit conversion ((Angstrom) / (desired units)) - ref_mass : float, optional, default=1.0 - Reference mass for unit conversion ((Dalton) / (desired units)) - ref_energy : float, optional, default=1.0 - Reference energy for unit conversion ((kcal/mol) / (desired units)) - auto_scale : bool, optional, default=False - Scale to reduced units by automatically using the largest sigma value - as ref_distance, largest mass value as ref_mass, and largest epsilon - value as ref_energy - nlist_buffer : float, optional, default=True - buffer argument to pass to hoomd.md.nlist.Cell - snapshot_kwargs : dict - Keyword arguments to pass to to_hoomdsnapshot - pppm_kwargs : dict - Keyword arguments to pass to hoomd.md.long_range.pppm.make_pppm_coulomb_forces - init_snap : hoomd.Snapshot, optional, default=None - Initial snapshot to which to add the ParmEd structure object - (useful for rigid bodies) - - Returns - ------- - hoomd_snapshot : hoomd.Snapshot - HOOMD snapshot object to initialize the simulation - hoomd_forcefield : list[hoomd.md.force.Force] - List of hoomd force computes created during conversion - ReferenceValues : namedtuple - Values used in scaling - - Notes - ----- - If you pass a non-parametrized pmd.Structure, you will not have - angle, dihedral, or force field information. You may be better off - creating a hoomd.Snapshot. - See mbuild.formats.hoomd_snapshot.to_hoomdsnapshot() - - About units: This method operates on a Parmed.Structure object - where the units used differ from those used in mBuild and Foyer - which may have been used when creating the typed Parmed.Structure. - - The default units used when writing out the HOOMD Snapshot are: - Distance (Angstrom) - Mass (Dalton) - Energy (kcal/mol) - - If you wish to convert this unit system to another, you can use the - reference parameters (ref_distance, ref_mass, ref_energy). - The values used here should be expected to convert from the Parmed - Structure units (above) to your desired units. - The Parmed.Structure values are divided by the reference values. - - If you wish to used a reduced unit system, set auto_scale = True. - When auto_scale is True, the reference parameters won't be used. - - Examples - -------- - To convert the energy units from kcal/mol to kj/mol: - use ref_energy = 0.2390057 (kcal/kj) - - To convert the distance units from Angstrom to nm: - use ref_distance = 10 (angstroms/nm) - - To use a reduced unit system, where mass, sigma, and epsilon are - scaled by the largest value of each: - use auto_scale = True, ref_distance = ref_energy = ref_mass = 1 - - """ - if isinstance(structure, mb.Compound): - raise ValueError( - "You passed mb.Compound to create_hoomd_simulation, there will be " - "no angles, dihedrals, or force field parameters. Please use " - "hoomd_snapshot.to_hoomdsnapshot to create a hoomd.Snapshot, then " - "create your own hoomd context and pass your hoomd.Snapshot to " - "hoomd.init.read_snapshot()" - ) - elif not isinstance(structure, pmd.Structure): - raise ValueError( - "Please pass a parmed.Structure to create_hoomd_simulation" - ) - - hoomd_version = _get_hoomd_version() - if hoomd_version.major < 3: - raise RuntimeError( - "Unsupported HOOMD-blue version:", str(hoomd_version) - ) - - hoomd_forcefield = [] - - if auto_scale: - if not all([i == 1 for i in (ref_distance, ref_energy, ref_mass)]): - warnings.warn( - "Autoscale option selected--provided reference values will not " - "be used." - ) - - pair_coeffs = list( - set((a.type, a.epsilon, a.sigma) for a in structure.atoms) - ) - - ref_mass = max([atom.mass for atom in structure.atoms]) - ref_energy = max(pair_coeffs, key=operator.itemgetter(1))[1] - ref_distance = max(pair_coeffs, key=operator.itemgetter(2))[2] - - ReferenceValues = namedtuple("ref_values", ["distance", "mass", "energy"]) - ref_values = ReferenceValues(ref_distance, ref_mass, ref_energy) - - snapshot, _ = to_hoomdsnapshot( - structure, - ref_distance=ref_distance, - ref_mass=ref_mass, - ref_energy=ref_energy, - **snapshot_kwargs, - hoomd_snapshot=init_snap, - ) - - nl = hoomd.md.nlist.Cell(exclusions=["bond", "1-3"], buffer=nlist_buffer) - - if structure.atoms[0].type != "": - print("Processing LJ and QQ") - lj = _init_hoomd_lj( - structure, - nl, - r_cut, - ref_distance=ref_distance, - ref_energy=ref_energy, - ) - qq = _init_hoomd_qq(structure, nl, snapshot, r_cut, **pppm_kwargs) - hoomd_forcefield.append(lj) - if qq is not None: - hoomd_forcefield.extend(qq) - if structure.adjusts: - print("Processing 1-4 interactions, adjusting neighborlist exclusions") - lj_14, qq_14 = _init_hoomd_14_pairs( - structure, - nl, - snapshot, - r_cut, - ref_distance=ref_distance, - ref_energy=ref_energy, - ) - hoomd_forcefield.append(lj_14) - hoomd_forcefield.append(qq_14) - if structure.bond_types: - print("Processing harmonic bonds") - harmonic_bond = _init_hoomd_bonds( - structure, ref_distance=ref_distance, ref_energy=ref_energy - ) - hoomd_forcefield.append(harmonic_bond) - if structure.angle_types: - print("Processing harmonic angles") - harmonic_angle = _init_hoomd_angles(structure, ref_energy=ref_energy) - hoomd_forcefield.append(harmonic_angle) - if structure.dihedral_types: - print("Processing periodic torsions") - periodic_torsions = _init_hoomd_dihedrals( - structure, ref_energy=ref_energy - ) - hoomd_forcefield.append(periodic_torsions) - if structure.rb_torsion_types: - print("Processing RB torsions") - rb_torsions = _init_hoomd_rb_torsions(structure, ref_energy=ref_energy) - hoomd_forcefield.append(rb_torsions) - - return snapshot, hoomd_forcefield, ref_values - - -def _init_hoomd_lj(structure, nl, r_cut, ref_distance=1.0, ref_energy=1.0): - """LJ parameters.""" - # Identify the unique atom types before setting - atom_type_params = {} - for atom in structure.atoms: - if atom.type not in atom_type_params: - atom_type_params[atom.type] = atom.atom_type - - # Set the hoomd parameters for self-interactions - lj = hoomd.md.pair.LJ(nlist=nl) - for name, atom_type in atom_type_params.items(): - lj.params[(name, name)] = dict( - sigma=atom_type.sigma / ref_distance, - epsilon=atom_type.epsilon / ref_energy, - ) - if atom_type.epsilon / ref_energy == 0: - lj.r_cut[(name, name)] = 0 - else: - lj.r_cut[(name, name)] = r_cut - - # Cross interactions, mixing rules, NBfixes - all_atomtypes = sorted(atom_type_params.keys()) - for a1, a2 in itertools.combinations_with_replacement(all_atomtypes, 2): - nb_fix_info = atom_type_params[a1].nbfix.get(a2, None) - # nb_fix_info = (rmin, eps, rmin14, eps14) - if nb_fix_info is None: - # No nbfix means use mixing rule to find cross-interaction - if structure.combining_rule == "lorentz": - sigma = ( - atom_type_params[a1].sigma + atom_type_params[a2].sigma - ) / (2 * ref_distance) - epsilon = ( - ( - atom_type_params[a1].epsilon - * atom_type_params[a2].epsilon - ) - / ref_energy**2 - ) ** 0.5 - elif structure.combining_rule == "geometric": - sigma = ( - (atom_type_params[a1].sigma * atom_type_params[a2].sigma) - / ref_distance**2 - ) ** 0.5 - epsilon = ( - ( - atom_type_params[a1].epsilon - * atom_type_params[a2].epsilon - ) - / ref_energy**2 - ) ** 0.5 - else: - raise ValueError( - f"Mixing rule {structure.combining_rule} not supported, " - 'use "lorentz" or "geometric"' - ) - else: - # If we have nbfix info, use it - sigma = nb_fix_info[0] / (ref_distance * (2 ** (1 / 6))) - epsilon = nb_fix_info[1] / ref_energy - lj.params[(a1, a2)] = dict(sigma=sigma, epsilon=epsilon) - if epsilon == 0: - lj.r_cut[(a1, a2)] = 0 - else: - lj.r_cut[(a1, a2)] = r_cut - - return lj - - -def _init_hoomd_qq(structure, nl, snapshot, r_cut, Nx=1, Ny=1, Nz=1, order=4): - """Charge interactions.""" - num_charged = np.sum(snapshot.particles.charge[:] != 0) - if num_charged == 0: - print("No charged groups found, ignoring electrostatics") - return None - else: - qq = hoomd.md.long_range.pppm.make_pppm_coulomb_forces( - nlist=nl, resolution=(Nx, Ny, Nz), order=order, r_cut=r_cut - ) - return qq - - -def _init_hoomd_14_pairs( - structure, nl, snapshot, r_cut, ref_distance=1.0, ref_energy=1.0 -): - """Special_pairs to handle 14 scaling. - - See discussion: https://groups.google.com/forum/ - #!topic/hoomd-users/iZ9WCpHczg0 - """ - # Update neighborlist to exclude 1-4 interactions, - # but impose a special_pair force to handle these pairs - nl.exclusions = nl.exclusions + [ - "1-4", - ] - - if snapshot.pairs.N == 0: - print("No 1,4 pairs found in hoomd snapshot") - return None, None - - lj_14 = hoomd.md.special_pair.LJ() - qq_14 = hoomd.md.special_pair.Coulomb() - params_14 = {} - # Identify unique 14 scalings - for adjust in structure.adjusts: - t1 = adjust.atom1.type - t2 = adjust.atom2.type - ps = "-".join(sorted([t1, t2])) - if ps not in params_14: - params_14[ps] = adjust.type - for name, adjust_type in params_14.items(): - lj_14.params[name] = dict( - sigma=adjust_type.sigma / ref_distance, - # The adjust epsilon already carries the scaling - epsilon=adjust_type.epsilon / ref_energy, - ) - if adjust_type.epsilon / ref_energy == 0: - lj_14.r_cut[name] = 0 - else: - lj_14.r_cut[name] = r_cut - qq_14.params[name] = dict(alpha=adjust_type.chgscale) - qq_14.r_cut[name] = r_cut - - return lj_14, qq_14 - - -def _init_hoomd_bonds(structure, ref_distance=1.0, ref_energy=1.0): - """Harmonic bonds.""" - # Identify the unique bond types before setting - bond_type_params = {} - for bond in structure.bonds: - t1, t2 = bond.atom1.type, bond.atom2.type - t1, t2 = sorted([t1, t2], key=natural_sort) - if t1 != "" and t2 != "": - bond_type = "-".join((t1, t2)) - if bond_type not in bond_type_params: - bond_type_params[bond_type] = bond.type - - # Set the hoomd parameters - harmonic_bond = hoomd.md.bond.Harmonic() - for name, bond_type in bond_type_params.items(): - # A (paramerized) parmed structure with no bondtype - # is because of constraints - if bond_type is None: - print("Bond with no bondtype detected, setting coefficients to 0") - harmonic_bond.params[name] = dict(k=0, r0=0) - else: - harmonic_bond.params[name] = dict( - k=2 * bond_type.k * ref_distance**2 / ref_energy, - r0=bond_type.req / ref_distance, - ) - - return harmonic_bond - - -def _init_hoomd_angles(structure, ref_energy=1.0): - """Harmonic angles.""" - # Identify the unique angle types before setting - angle_type_params = {} - for angle in structure.angles: - t1, t2, t3 = angle.atom1.type, angle.atom2.type, angle.atom3.type - t1, t3 = sorted([t1, t3], key=natural_sort) - angle_type = "-".join((t1, t2, t3)) - if angle_type not in angle_type_params: - angle_type_params[angle_type] = angle.type - - # set the hoomd parameters - harmonic_angle = hoomd.md.angle.Harmonic() - for name, angle_type in angle_type_params.items(): - harmonic_angle.params[name] = dict( - t0=np.deg2rad(angle_type.theteq), - k=2 * angle_type.k / ref_energy, - ) - - return harmonic_angle - - -def _init_hoomd_dihedrals(structure, ref_energy=1.0): - """Periodic dihedrals (dubbed harmonic dihedrals in HOOMD).""" - dihedral_type_params = {} - for dihedral in structure.dihedrals: - t1, t2 = dihedral.atom1.type, dihedral.atom2.type - t3, t4 = dihedral.atom3.type, dihedral.atom4.type - if [t2, t3] == sorted([t2, t3], key=natural_sort): - dihedral_type = "-".join((t1, t2, t3, t4)) - else: - dihedral_type = "-".join((t4, t3, t2, t1)) - if dihedral_type not in dihedral_type_params: - if isinstance(dihedral.type, pmd.DihedralType): - dihedral_type_params[dihedral_type] = dihedral.type - elif isinstance(dihedral.type, pmd.DihedralTypeList): - if len(dihedral.type) > 1: - warnings.warn( - "Multiple dihedral types detected" - + " for single dihedral, will ignore all except " - + " first dihedral type." - + "First dihedral type: {}".format(dihedral.type[0]) - ) - dihedral_type_params[dihedral_type] = dihedral.type[0] - - # Set the hoomd parameters - # These are periodic torsions - hoomd_version = _get_hoomd_version() - if ( - hoomd_version.major == 3 and hoomd_version.minor >= 7 - ) or hoomd_version.major == 4: - periodic_torsion = hoomd.md.dihedral.Periodic() - else: - periodic_torsion = hoomd.md.dihedral.Harmonic() - for name, dihedral_type in dihedral_type_params.items(): - periodic_torsion.params[name] = dict( - k=2 * dihedral_type.phi_k / ref_energy, - d=1, - n=dihedral_type.per, - phi0=np.deg2rad(dihedral_type.phase), - ) - - return periodic_torsion - - -def _init_hoomd_rb_torsions(structure, ref_energy=1.0): - """RB dihedrals (implemented as OPLS dihedrals in HOOMD).""" - # Identify the unique dihedral types before setting - dihedral_type_params = {} - for dihedral in structure.rb_torsions: - t1, t2 = dihedral.atom1.type, dihedral.atom2.type - t3, t4 = dihedral.atom3.type, dihedral.atom4.type - if [t2, t3] == sorted([t2, t3], key=natural_sort): - dihedral_type = "-".join((t1, t2, t3, t4)) - else: - dihedral_type = "-".join((t4, t3, t2, t1)) - if dihedral_type not in dihedral_type_params: - dihedral_type_params[dihedral_type] = dihedral.type - - # Set the hoomd parameter - rb_torsion = hoomd.md.dihedral.OPLS() - for name, dihedral_type in dihedral_type_params.items(): - F_coeffs = RB_to_OPLS( - dihedral_type.c0 / ref_energy, - dihedral_type.c1 / ref_energy, - dihedral_type.c2 / ref_energy, - dihedral_type.c3 / ref_energy, - dihedral_type.c4 / ref_energy, - dihedral_type.c5 / ref_energy, - ) - rb_torsion.params[name] = dict( - k1=F_coeffs[1], k2=F_coeffs[2], k3=F_coeffs[3], k4=F_coeffs[4] - ) - - return rb_torsion diff --git a/mbuild/formats/hoomd_simulation.py b/mbuild/formats/hoomd_simulation.py deleted file mode 100644 index 12da75818..000000000 --- a/mbuild/formats/hoomd_simulation.py +++ /dev/null @@ -1,467 +0,0 @@ -"""HOOMD simulation format.""" - -import itertools -import operator -import warnings -from collections import namedtuple - -import numpy as np -import parmed as pmd - -import mbuild as mb -from mbuild.utils.conversion import RB_to_OPLS -from mbuild.utils.decorators import deprecated -from mbuild.utils.io import import_ -from mbuild.utils.sorting import natural_sort - -from .hoomd_snapshot import to_hoomdsnapshot - -gsd = import_("gsd") -gsd.hoomd = import_("gsd.hoomd") -hoomd = import_("hoomd") -hoomd.md = import_("hoomd.md") - - -dep_msg = """ -Support for Hoomd-Blue 2.0 will be removed in mBuild 1.0. -See GMSO (https://github.com/mosdef-hub/gmso/tree/main/gmso/external/convert_hoomd) for -Hoomd-Blue 3.x and 4.x format support. -To convert to GMSO, use the method `Compound.to_gmso()`. -""" -print(dep_msg) - - -@deprecated(dep_msg) -def create_hoomd_simulation( - structure, - r_cut, - ref_distance=1.0, - ref_mass=1.0, - ref_energy=1.0, - auto_scale=False, - snapshot_kwargs={}, - pppm_kwargs={"Nx": 8, "Ny": 8, "Nz": 8, "order": 4}, - init_snap=None, - restart=None, - nlist=hoomd.md.nlist.cell, -): - """Convert a parametrized pmd.Structure to hoomd.SimulationContext. - - Parameters - ---------- - structure : parmed.Structure - ParmEd Structure object - r_cut : float - Cutoff radius in simulation units - ref_distance : float, optional, default=1.0 - Reference distance for unit conversion (from Angstrom) - ref_mass : float, optional, default=1.0 - Reference mass for unit conversion (from Dalton) - ref_energy : float, optional, default=1.0 - Reference energy for unit conversion (from kcal/mol) - auto_scale : bool, optional, default=False - Scale to reduced units by automatically using the largest sigma value - as ref_distance, largest mass value as ref_mass, and largest epsilon - value as ref_energy - snapshot_kwargs : dict - Kwargs to pass to to_hoomdsnapshot - pppm_kwargs : dict - Kwargs to pass to hoomd's pppm function - init_snap : hoomd.data.SnapshotParticleData, optional, default=None - Initial snapshot to which to add the ParmEd structure object - (useful for rigid bodies) - restart : str, optional, default=None - Path to the gsd file from which to restart the simulation. - https://hoomd-blue.readthedocs.io/en/v2.9.4/restartable-jobs.html - Note: It is assumed that the ParmEd structure and the system in - restart.gsd contain the same types. The ParmEd structure is still used - to initialize the forces, but restart.gsd is used to initialize the - system state (e.g., particle positions, momenta, etc). - nlist : hoomd.md.nlist, default=hoomd.md.nlist.cell - Type of neighborlist to use, see - https://hoomd-blue.readthedocs.io/en/stable/module-md-nlist.html - for more information. - - Returns - ------- - hoomd_objects : list - List of hoomd objects created during conversion - ReferenceValues : namedtuple - Values used in scaling - - Note - ---- - While the hoomd objects are returned, the hoomd.SimulationContext is - accessible via `hoomd.context.current`. If you pass a non-parametrized - pmd.Structure, you will not have angle, dihedral, or force field - information. You may be better off creating a hoomd.Snapshot. - Reference units should be expected to convert parmed Structure units: - - --- angstroms, kcal/mol, and daltons - """ - if isinstance(structure, mb.Compound): - raise ValueError( - "You passed mb.Compound to create_hoomd_simulation, there will be " - "no angles, dihedrals, or force field parameters. Please use " - "hoomd_snapshot.to_hoomdsnapshot to create a hoomd.Snapshot, then " - "create your own hoomd context and pass your hoomd.Snapshot to " - "hoomd.init.read_snapshot()" - ) - elif not isinstance(structure, pmd.Structure): - raise ValueError( - "Please pass a parmed.Structure to create_hoomd_simulation" - ) - - _check_hoomd_version() - version_numbers = _check_hoomd_version() - if float(version_numbers[0]) >= 3: - raise RuntimeError("This method does not support HOOMD-blue v3.x.") - - hoomd_objects = [] # Potential adaptation for Hoomd v3 API - - if auto_scale: - if not all([i == 1 for i in (ref_distance, ref_energy, ref_mass)]): - warnings.warn( - "Autoscale option selected--provided reference values will not " - "be used." - ) - pair_coeffs = list( - set((a.type, a.epsilon, a.sigma) for a in structure.atoms) - ) - - ref_mass = max([atom.mass for atom in structure.atoms]) - ref_energy = max(pair_coeffs, key=operator.itemgetter(1))[1] - ref_distance = max(pair_coeffs, key=operator.itemgetter(2))[2] - - ReferenceValues = namedtuple("ref_values", ["distance", "mass", "energy"]) - ref_values = ReferenceValues(ref_distance, ref_mass, ref_energy) - - if not hoomd.context.current: - hoomd.context.initialize("") - if restart is None: - snapshot, _ = to_hoomdsnapshot( - structure, - ref_distance=ref_distance, - ref_mass=ref_mass, - ref_energy=ref_energy, - **snapshot_kwargs, - hoomd_snapshot=init_snap, - ) - hoomd_objects.append(snapshot) - hoomd_system = hoomd.init.read_snapshot(snapshot) - hoomd_objects.append(hoomd_system) - else: - with gsd.hoomd.open(restart) as f: - snapshot = f[-1] - hoomd_objects.append(snapshot) - hoomd_system = hoomd.init.read_gsd(restart, restart=restart) - hoomd_objects.append(hoomd_system) - print("Simulation initialized from restart file") - - nl = nlist() - nl.reset_exclusions(exclusions=["1-2", "1-3"]) - hoomd_objects.append(nl) - - if structure.atoms[0].type != "": - print("Processing LJ and QQ") - lj = _init_hoomd_lj( - structure, - nl, - r_cut, - ref_distance=ref_distance, - ref_energy=ref_energy, - ) - qq = _init_hoomd_qq(structure, nl, r_cut, **pppm_kwargs) - hoomd_objects.append(lj) - hoomd_objects.append(qq) - if structure.adjusts: - print("Processing 1-4 interactions, adjusting neighborlist exclusions") - lj_14, qq_14 = _init_hoomd_14_pairs( - structure, - nl, - r_cut, - ref_distance=ref_distance, - ref_energy=ref_energy, - ) - hoomd_objects.append(lj_14) - hoomd_objects.append(qq_14) - if structure.bond_types: - print("Processing harmonic bonds") - harmonic_bond = _init_hoomd_bonds( - structure, ref_distance=ref_distance, ref_energy=ref_energy - ) - hoomd_objects.append(harmonic_bond) - if structure.angle_types: - print("Processing harmonic angles") - harmonic_angle = _init_hoomd_angles(structure, ref_energy=ref_energy) - hoomd_objects.append(harmonic_angle) - if structure.dihedral_types: - print("Processing periodic torsions") - periodic_torsions = _init_hoomd_dihedrals( - structure, ref_energy=ref_energy - ) - hoomd_objects.append(periodic_torsions) - if structure.rb_torsion_types: - print("Processing RB torsions") - rb_torsions = _init_hoomd_rb_torsions(structure, ref_energy=ref_energy) - hoomd_objects.append(rb_torsions) - print("HOOMD SimulationContext updated from ParmEd Structure") - - return hoomd_objects, ref_values - - -def _init_hoomd_lj(structure, nl, r_cut, ref_distance=1.0, ref_energy=1.0): - """LJ parameters.""" - # Identify the unique atom types before setting - atom_type_params = {} - for atom in structure.atoms: - if atom.type not in atom_type_params: - atom_type_params[atom.type] = atom.atom_type - - # Set the hoomd parameters for self-interactions - lj = hoomd.md.pair.lj(r_cut, nl) - for name, atom_type in atom_type_params.items(): - lj.pair_coeff.set( - name, - name, - sigma=atom_type.sigma / ref_distance, - epsilon=atom_type.epsilon / ref_energy, - ) - if atom_type.epsilon / ref_energy == 0: - lj.pair_coeff.set(name, name, r_cut=0) - - # Cross interactions, mixing rules, NBfixes - all_atomtypes = sorted(atom_type_params.keys()) - for a1, a2 in itertools.combinations_with_replacement(all_atomtypes, 2): - nb_fix_info = atom_type_params[a1].nbfix.get(a2, None) - # nb_fix_info = (rmin, eps, rmin14, eps14) - if nb_fix_info is None: - # No nbfix means use mixing rule to find cross-interaction - if structure.combining_rule == "lorentz": - sigma = ( - atom_type_params[a1].sigma + atom_type_params[a2].sigma - ) / (2 * ref_distance) - epsilon = ( - ( - atom_type_params[a1].epsilon - * atom_type_params[a2].epsilon - ) - / ref_energy**2 - ) ** 0.5 - elif structure.combining_rule == "geometric": - sigma = ( - (atom_type_params[a1].sigma * atom_type_params[a2].sigma) - / ref_distance**2 - ) ** 0.5 - epsilon = ( - ( - atom_type_params[a1].epsilon - * atom_type_params[a2].epsilon - ) - / ref_energy**2 - ) ** 0.5 - else: - raise ValueError( - f"Mixing rule {structure.combining_rule} not supported, " - "use lorentz" - ) - else: - # If we have nbfix info, use it - sigma = nb_fix_info[0] / (ref_distance * (2 ** (1 / 6))) - epsilon = nb_fix_info[1] / ref_energy - lj.pair_coeff.set(a1, a2, sigma=sigma, epsilon=epsilon) - if epsilon == 0: - lj.pair_coeff.set(a1, a2, r_cut=0) - - return lj - - -def _init_hoomd_qq(structure, nl, r_cut, Nx=1, Ny=1, Nz=1, order=4): - """Charge interactions.""" - charged = hoomd.group.charged() - if len(charged) == 0: - print("No charged groups found, ignoring electrostatics") - return None - else: - qq = hoomd.md.charge.pppm(charged, nl) - qq.set_params(Nx, Ny, Nz, order, r_cut) - return qq - - -def _init_hoomd_14_pairs( - structure, nl, r_cut, ref_distance=1.0, ref_energy=1.0 -): - """Special_pairs to handle 14 scalings. - - See discussion: https://groups.google.com/forum/#!topic/hoomd-users/ - iZ9WCpHczg0 - """ - # Update neighborlist to exclude 1-4 interactions, - # but impose a special_pair force to handle these pairs - nl.exclusions.append("1-4") - - if hoomd.context.current.system_definition.getPairData().getN() == 0: - print("No 1,4 pairs found in hoomd snapshot") - return None, None - - lj_14 = hoomd.md.special_pair.lj() - qq_14 = hoomd.md.special_pair.coulomb() - params_14 = {} - # Identify unique 14 scalings - for adjust in structure.adjusts: - t1 = adjust.atom1.type - t2 = adjust.atom2.type - ps = "-".join(sorted([t1, t2])) - if ps not in params_14: - params_14[ps] = adjust.type - for name, adjust_type in params_14.items(): - lj_14.pair_coeff.set( - name, - sigma=adjust_type.sigma / ref_distance, - # The adjust epsilon already carries the scaling - epsilon=adjust_type.epsilon / ref_energy, - # Do NOT use hoomd's alpha to modify any LJ terms - alpha=1, - r_cut=r_cut, - ) - qq_14.pair_coeff.set(name, alpha=adjust_type.chgscale, r_cut=r_cut) - - return lj_14, qq_14 - - -def _init_hoomd_bonds(structure, ref_distance=1.0, ref_energy=1.0): - """Harmonic bonds.""" - # Identify the unique bond types before setting - bond_type_params = {} - for bond in structure.bonds: - t1, t2 = bond.atom1.type, bond.atom2.type - t1, t2 = sorted([t1, t2], key=natural_sort) - if t1 != "" and t2 != "": - bond_type = "-".join((t1, t2)) - if bond_type not in bond_type_params: - bond_type_params[bond_type] = bond.type - - # Set the hoomd parameters - harmonic_bond = hoomd.md.bond.harmonic() - for name, bond_type in bond_type_params.items(): - # A (paramerized) parmed structure with no bondtype - # is because of constraints - if bond_type is None: - print("Bond with no bondtype detected, setting coefficients to 0") - harmonic_bond.bond_coeff.set(name, k=0, r0=0) - else: - harmonic_bond.bond_coeff.set( - name, - k=2 * bond_type.k * ref_distance**2 / ref_energy, - r0=bond_type.req / ref_distance, - ) - - return harmonic_bond - - -def _init_hoomd_angles(structure, ref_energy=1.0): - """Harmonic angles.""" - # Identify the unique angle types before setting - angle_type_params = {} - for angle in structure.angles: - t1, t2, t3 = angle.atom1.type, angle.atom2.type, angle.atom3.type - t1, t3 = sorted([t1, t3], key=natural_sort) - angle_type = "-".join((t1, t2, t3)) - if angle_type not in angle_type_params: - angle_type_params[angle_type] = angle.type - - # set the hoomd parameters - harmonic_angle = hoomd.md.angle.harmonic() - for name, angle_type in angle_type_params.items(): - harmonic_angle.angle_coeff.set( - name, - t0=np.deg2rad(angle_type.theteq), - k=2 * angle_type.k / ref_energy, - ) - - return harmonic_angle - - -def _init_hoomd_dihedrals(structure, ref_energy=1.0): - """Periodic dihedrals (dubbed harmonic dihedrals in HOOMD).""" - # Identify the unique dihedral types before setting - # need Hoomd 2.8.0 to use proper dihedral implemtnation - # from this PR https://github.com/glotzerlab/hoomd-blue/pull/492 - version_numbers = _check_hoomd_version() - if float(version_numbers[0]) < 2 or float(version_numbers[1]) < 8: - from mbuild.exceptions import MBuildError - - raise MBuildError("Please upgrade Hoomd to at least 2.8.0") - - dihedral_type_params = {} - for dihedral in structure.dihedrals: - t1, t2 = dihedral.atom1.type, dihedral.atom2.type - t3, t4 = dihedral.atom3.type, dihedral.atom4.type - if [t2, t3] == sorted([t2, t3], key=natural_sort): - dihedral_type = "-".join((t1, t2, t3, t4)) - else: - dihedral_type = "-".join((t4, t3, t2, t1)) - if dihedral_type not in dihedral_type_params: - if isinstance(dihedral.type, pmd.DihedralType): - dihedral_type_params[dihedral_type] = dihedral.type - elif isinstance(dihedral.type, pmd.DihedralTypeList): - if len(dihedral.type) > 1: - warnings.warn( - "Multiple dihedral types detected for single dihedral, " - "will ignore all except first dihedral type. First " - "dihedral type: {}".format(dihedral.type[0]) - ) - dihedral_type_params[dihedral_type] = dihedral.type[0] - - # Set the hoomd parameters - # These are periodic torsions - periodic_torsion = hoomd.md.dihedral.harmonic() - for name, dihedral_type in dihedral_type_params.items(): - periodic_torsion.dihedral_coeff.set( - name, - k=2 * dihedral_type.phi_k / ref_energy, - d=1, - n=dihedral_type.per, - phi_0=np.deg2rad(dihedral_type.phase), - ) - - return periodic_torsion - - -def _init_hoomd_rb_torsions(structure, ref_energy=1.0): - """RB dihedrals (implemented as OPLS dihedrals in HOOMD).""" - # Identify the unique dihedral types before setting - dihedral_type_params = {} - for dihedral in structure.rb_torsions: - t1, t2 = dihedral.atom1.type, dihedral.atom2.type - t3, t4 = dihedral.atom3.type, dihedral.atom4.type - if [t2, t3] == sorted([t2, t3], key=natural_sort): - dihedral_type = "-".join((t1, t2, t3, t4)) - else: - dihedral_type = "-".join((t4, t3, t2, t1)) - if dihedral_type not in dihedral_type_params: - dihedral_type_params[dihedral_type] = dihedral.type - - # Set the hoomd parameter - rb_torsion = hoomd.md.dihedral.opls() - for name, dihedral_type in dihedral_type_params.items(): - F_coeffs = RB_to_OPLS( - dihedral_type.c0 / ref_energy, - dihedral_type.c1 / ref_energy, - dihedral_type.c2 / ref_energy, - dihedral_type.c3 / ref_energy, - dihedral_type.c4 / ref_energy, - dihedral_type.c5 / ref_energy, - error_if_outside_tolerance=False, - ) - rb_torsion.dihedral_coeff.set( - name, k1=F_coeffs[1], k2=F_coeffs[2], k3=F_coeffs[3], k4=F_coeffs[4] - ) - - return rb_torsion - - -def _check_hoomd_version(): - version = hoomd.__version__ - version_numbers = version.split(".") - return version_numbers diff --git a/mbuild/formats/hoomd_snapshot.py b/mbuild/formats/hoomd_snapshot.py deleted file mode 100644 index b52b9db72..000000000 --- a/mbuild/formats/hoomd_snapshot.py +++ /dev/null @@ -1,622 +0,0 @@ -"""HOOMD snapshot format.""" - -import operator -from collections import namedtuple - -import numpy as np -import packaging.version -import parmed as pmd - -from mbuild.box import Box -from mbuild.compound import Compound, Particle -from mbuild.utils.geometry import coord_shift -from mbuild.utils.io import import_ -from mbuild.utils.sorting import natural_sort - -hoomd = import_("hoomd") - -__all__ = ["to_hoomdsnapshot", "from_snapshot"] - - -def _get_hoomd_version(): - if "version" in dir(hoomd): - return packaging.version.parse(hoomd.version.version) - else: - return packaging.version.parse(hoomd.__version__) - - -def from_snapshot(snapshot, scale=1.0): - """Convert a Snapshot to a Compound. - - Snapshot can be a hoomd.Snapshot or a gsd.hoomd.Frame. - - Parameters - ---------- - snapshot : hoomd.Snapshot or gsd.hoomd.Frame - Snapshot from which to build the mbuild Compound. - scale : float, optional, default 1.0 - Value by which to scale the length values - - Returns - ------- - comp : Compound - - Note - ---- - GSD and HOOMD snapshots center their boxes on the origin (0,0,0), so the - compound is shifted by half the box lengths - """ - comp = Compound() - bond_array = snapshot.bonds.group - n_atoms = snapshot.particles.N - - if "SnapshotSystemData_float" in dir(hoomd._hoomd) and isinstance( - snapshot, hoomd._hoomd.SnapshotSystemData_float - ): - # hoomd v2 - box = snapshot.box - comp.box = Box.from_lengths_tilt_factors( - lengths=np.array([box.Lx, box.Ly, box.Lz]) * scale, - tilt_factors=np.array([box.xy, box.xz, box.yz]), - ) - else: - # gsd / hoomd v3 - box = np.asarray(snapshot.configuration.box) - comp.box = Box.from_lengths_tilt_factors( - lengths=box[:3] * scale, tilt_factors=box[3:] - ) - - # GSD and HOOMD snapshots center their boxes on the origin (0,0,0) - shift = np.array(comp.box.lengths) / 2 - # Add particles - for i in range(n_atoms): - name = snapshot.particles.types[snapshot.particles.typeid[i]] - xyz = snapshot.particles.position[i] * scale + shift - charge = snapshot.particles.charge[i] - - atom = Particle(name=name, pos=xyz, charge=charge) - comp.add(atom, label=str(i)) - - # Add bonds - particle_dict = {idx: p for idx, p in enumerate(comp.particles())} - for i in range(bond_array.shape[0]): - atom1 = int(bond_array[i][0]) - atom2 = int(bond_array[i][1]) - comp.add_bond([particle_dict[atom1], particle_dict[atom2]]) - return comp - - -def to_hoomdsnapshot( - structure, - ref_distance=1.0, - ref_mass=1.0, - ref_energy=1.0, - rigid_bodies=None, - shift_coords=True, - write_special_pairs=True, - auto_scale=False, - parmed_kwargs={}, - hoomd_snapshot=None, -): - """Convert a Compound or parmed.Structure to hoomd.Snapshot. - - Parameters - ---------- - structure : parmed.Structure - ParmEd Structure object - Reference distance for unit conversion ((Angstrom) / (desired units)) - ref_mass : float, optional, default=1.0 - Reference mass for unit conversion ((Dalton) / (desired units)) - ref_energy : float, optional, default=1.0 - Reference energy for unit conversion ((kcal/mol) / (desired units)) - rigid_bodies : list of int, optional, default=None - List of rigid body information. An integer value is required for - each atom corresponding to the index of the rigid body the particle - is to be associated with. A value of None indicates the atom is not - part of a rigid body. - shift_coords : bool, optional, default=True - Shift coordinates from (0, L) to (-L/2, L/2) if necessary. - auto_scale : bool, optional, default=False - Automatically use largest sigma value as ref_distance, - largest mass value as ref_mass - and largest epsilon value as ref_energy - write_special_pairs : bool, optional, default=True - Writes out special pair information necessary to correctly use - the OPLS fudged 1,4 interactions in HOOMD. - hoomd_snapshot : hoomd.Snapshot, optional, default=None - Initial snapshot to which to add the ParmEd structure object. - The box information of the initial snapshot will be overwritten. - (useful for rigid bodies) - - Returns - ------- - hoomd_snapshot : hoomd.Snapshot - ReferenceValues : namedtuple - Values used in scaling - - Notes - ----- - This method does not create hoomd forcefield objects and - the snapshot returned does not store the forcefield parameters. - See mbuild.formats.hoomd_forcefield.create_hoomd_forcefield() - - About units: This method operates on a Parmed.Structure object - where the units used differ from those used in mBuild and Foyer - which may have been used when creating the typed Parmed.Structure. - - The default units used when writing out the HOOMD Snapshot are: - Distance (Angstrom) - Mass (Dalton) - Energy (kcal/mol) - - If you wish to convert this unit system to another, you can use the - reference parameters (ref_distance, ref_mass, ref_energy). - The values used here should be expected to convert from the Parmed - Structure units (above) to your desired units. - The Parmed.Structure values are divided by the reference values. - - If you wish to used a reduced unit system, set auto_scale = True. - When auto_scale is True, the reference parameters won't be used. - - Examples - -------- - To convert the energy units from kcal/mol to kj/mol: - use ref_energy = 0.2390057 (kcal/kj) - - To convert the distance units from Angstrom to nm: - use ref_distance = 10 (angstroms/nm) - - To use a reduced unit system, where mass, sigma, and epsilon are - scaled by the largest value of each: - use auto_scale = True, ref_distance = ref_energy = ref_mass = 1 - - """ - if not isinstance(structure, (Compound, pmd.Structure)): - raise ValueError( - "You are trying to create a hoomd.Snapshot from a " - f"{type(structure)} please pass a Compound or pmd.Structure" - ) - elif isinstance(structure, Compound): - structure = structure.to_parmed(**parmed_kwargs) - - hoomd_version = _get_hoomd_version() - if hoomd_version.major == 2 and not hoomd.context.current: - hoomd.context.initialize("") - - if auto_scale: - ref_mass = max([atom.mass for atom in structure.atoms]) - pair_coeffs = list( - set( - (atom.type, atom.epsilon, atom.sigma) - for atom in structure.atoms - ) - ) - ref_energy = max(pair_coeffs, key=operator.itemgetter(1))[1] - ref_distance = max(pair_coeffs, key=operator.itemgetter(2))[2] - - ReferenceValues = namedtuple("ref_values", ["distance", "mass", "energy"]) - ref_values = ReferenceValues(ref_distance, ref_mass, ref_energy) - - xyz = np.array([[atom.xx, atom.xy, atom.xz] for atom in structure.atoms]) - if shift_coords: - xyz = coord_shift(xyz, structure.box[:3]) - - # Get box information - if np.allclose(structure.box[3:6], np.array([90, 90, 90])): - lx, ly, lz = structure.box[:3] / ref_distance - xy, xz, yz = 0, 0, 0 - else: - a, b, c = structure.box[0:3] / ref_distance - alpha, beta, gamma = np.radians(structure.box[3:6]) - - lx = a - xy = b * np.cos(gamma) - xz = c * np.cos(beta) - ly = np.sqrt(b**2 - xy**2) - yz = (b * c * np.cos(alpha) - xy * xz) / ly - lz = np.sqrt(c**2 - xz**2 - yz**2) - - ( - n_particles, - scaled_positions, - unique_types, - typeids, - scaled_mass, - scaled_charges, - rigid_bodies, - ) = _parse_particle_information( - structure, xyz, ref_distance, ref_mass, ref_energy, rigid_bodies - ) - ( - n_bonds, - unique_bond_types, - bond_typeids, - bond_groups, - ) = _parse_bond_information(structure) - ( - n_angles, - unique_angle_types, - angle_typeids, - angle_groups, - ) = _parse_angle_information(structure) - ( - n_dihedrals, - unique_dihedral_types, - dihedral_typeids, - dihedral_groups, - ) = _parse_dihedral_information(structure) - ( - n_impropers, - unique_improper_types, - improper_typeids, - improper_groups, - ) = _parse_improper_information(structure) - pair_types, pair_typeid, pairs, n_pairs = _parse_pair_information(structure) - - if hoomd_snapshot is not None: - n_init = hoomd_snapshot.particles.N - - if n_init > 0: - n_particles += n_init - # shift the typeids - typeids += len(set(hoomd_snapshot.particles.types)) - hoomd_snapshot.particles.types += unique_types - # shift bond/angle/dihedral indices - if n_bonds > 0: - bond_groups = [ - tuple(row) for row in np.array(bond_groups) + n_init - ] - init_bondtypes = len(hoomd_snapshot.bonds.types) - hoomd_snapshot.bonds.types += unique_bond_types - if n_angles > 0: - angle_groups = [ - tuple(row) for row in np.array(angle_groups) + n_init - ] - init_angletypes = len(hoomd_snapshot.angles.types) - hoomd_snapshot.angles.types += unique_angle_types - if n_dihedrals > 0: - dihedral_groups = [ - tuple(row) for row in np.array(dihedral_groups) + n_init - ] - init_dihedraltypes = len(hoomd_snapshot.dihedrals.types) - hoomd_snapshot.dihedrals.types += unique_dihedral_types - if n_impropers > 0: - improper_groups = [ - tuple(row) for row in np.array(improper_groups) + n_init - ] - init_impropertypes = len(hoomd_snapshot.impropers.types) - hoomd_snapshot.impropers.types += unique_improper_types - if n_pairs > 0: - pairs = [tuple(row) for row in np.array(pairs) + n_init] - init_pairtypes = len(hoomd_snapshot.pairs.types) - hoomd_snapshot.pairs.types += pair_types - else: - raise RuntimeError( - "Initial snapshot provided, but it contains no particles" - ) - - if hoomd_version.major == 2: - hoomd_snapshot.box = hoomd.data.boxdim( - Lx=lx, Ly=ly, Lz=lz, xy=xy, xz=xz, yz=yz - ) - - # save the box for later use when wrapping coordinates - box = hoomd_snapshot.box - elif hoomd_version.major == 3 or hoomd_version.major == 4: - hoomd_snapshot.configuration.box = [lx, ly, lz, xy, xz, yz] - else: - raise RuntimeError("Unsupported HOOMD version:", str(hoomd_version)) - - init_bonds = hoomd_snapshot.bonds.N - if init_bonds > 0: - n_bonds += init_bonds - bond_typeids = list(np.array(bond_typeids) + init_bondtypes) - - init_angles = hoomd_snapshot.angles.N - if init_angles > 0: - n_angles += init_angles - angle_typeids = list(np.array(angle_typeids) + init_angletypes) - - init_dihedrals = hoomd_snapshot.dihedrals.N - if init_dihedrals > 0: - n_dihedrals += init_dihedrals - dihedral_typeids = list( - np.array(dihedral_typeids) + init_dihedraltypes - ) - - init_impropers = hoomd_snapshot.impropers.N - if init_impropers > 0: - n_impropers += init_impropers - improper_typeids = list( - np.array(improper_typeids) + init_impropertypes - ) - - init_pairs = hoomd_snapshot.pairs.N - if init_pairs > 0: - n_pairs += init_pairs - pair_typeid = list(np.array(pair_typeid) + init_pairtypes) - - else: - n_init = 0 - init_bonds = 0 - init_angles = 0 - init_dihedrals = 0 - init_impropers = 0 - init_pairs = 0 - - if hoomd_version.major == 2: - hoomd_snapshot = hoomd.data.make_snapshot( - N=n_particles, - box=hoomd.data.boxdim(Lx=lx, Ly=ly, Lz=lz, xy=xy, xz=xz, yz=yz), - particle_types=unique_types, - bond_types=unique_bond_types, - angle_types=unique_angle_types, - dihedral_types=unique_dihedral_types, - improper_types=unique_improper_types, - pair_types=pair_types, - ) - box = hoomd.data.boxdim(Lx=lx, Ly=ly, Lz=lz, xy=xy, xz=xz, yz=yz) - elif hoomd_version.major == 3 or hoomd_version.major == 4: - hoomd_snapshot = hoomd.Snapshot() - hoomd_snapshot.configuration.box = [lx, ly, lz, xy, xz, yz] - hoomd_snapshot.particles.types = unique_types - hoomd_snapshot.bonds.types = unique_bond_types - hoomd_snapshot.angles.types = unique_angle_types - hoomd_snapshot.dihedrals.types = unique_dihedral_types - hoomd_snapshot.impropers.types = unique_improper_types - hoomd_snapshot.pairs.types = pair_types - box = hoomd.Box(Lx=lx, Ly=ly, Lz=lz, xy=xy, xz=xz, yz=yz) - else: - raise RuntimeError("Unsupported HOOMD version:", str(hoomd_version)) - - # wrap particles into the box manually for v2 - if hoomd_version.major == 2: - scaled_positions = np.stack( - [box.wrap(xyz)[0] for xyz in scaled_positions] - ) - - def set_size(obj, n): - if hoomd_version.major == 2: - obj.resize(n) - elif hoomd_version.major == 3 or hoomd_version.major == 4: - obj.N = n - else: - raise RuntimeError("Unsupported HOOMD version:", str(hoomd_version)) - - set_size(hoomd_snapshot.particles, n_particles) - hoomd_snapshot.particles.position[n_init:] = scaled_positions - hoomd_snapshot.particles.typeid[n_init:] = typeids - hoomd_snapshot.particles.mass[n_init:] = scaled_mass - hoomd_snapshot.particles.charge[n_init:] = scaled_charges - hoomd_snapshot.particles.body[n_init:] = rigid_bodies - - # wrap the particles into the box using the v3 API - if hoomd_version.major == 3 or hoomd_version.major == 4: - hoomd_snapshot.wrap() - - if n_bonds > 0: - set_size(hoomd_snapshot.bonds, n_bonds) - hoomd_snapshot.bonds.typeid[init_bonds:] = bond_typeids - hoomd_snapshot.bonds.group[init_bonds:] = bond_groups - - if n_angles > 0: - set_size(hoomd_snapshot.angles, n_angles) - hoomd_snapshot.angles.typeid[init_angles:] = angle_typeids - hoomd_snapshot.angles.group[init_angles:] = np.reshape( - angle_groups, (-1, 3) - ) - - if n_dihedrals > 0: - set_size(hoomd_snapshot.dihedrals, n_dihedrals) - hoomd_snapshot.dihedrals.typeid[init_dihedrals:] = dihedral_typeids - hoomd_snapshot.dihedrals.group[init_dihedrals:] = np.reshape( - dihedral_groups, (-1, 4) - ) - - if n_impropers > 0: - set_size(hoomd_snapshot.impropers, n_impropers) - hoomd_snapshot.impropers.typeid[init_impropers:] = improper_typeids - hoomd_snapshot.impropers.group[init_impropers:] = np.reshape( - improper_groups, (-1, 4) - ) - - if n_pairs > 0: - set_size(hoomd_snapshot.pairs, n_pairs) - hoomd_snapshot.pairs.typeid[init_pairs:] = pair_typeid - hoomd_snapshot.pairs.group[init_pairs:] = np.reshape(pairs, (-1, 2)) - - return hoomd_snapshot, ref_values - - -def _parse_particle_information( - structure, xyz, ref_distance, ref_mass, ref_energy, rigid_bodies -): - """Write out the particle information.""" - n_particles = len(structure.atoms) - scaled_positions = xyz / ref_distance - - types = [ - atom.name if atom.type == "" else atom.type for atom in structure.atoms - ] - - unique_types = list(set(types)) - unique_types.sort(key=natural_sort) - - typeids = np.array([unique_types.index(t) for t in types]) - - masses = np.array([atom.mass for atom in structure.atoms]) - masses[masses == 0] = 1.0 - scaled_mass = masses / ref_mass - - charges = np.array([atom.charge for atom in structure.atoms]) - e0 = 2.396452e-04 - """ - Permittivity of free space = 2.396452e-04 e^2/((kcal/mol) Angstrom), - where e is the elementary charge - """ - charge_factor = (4.0 * np.pi * e0 * ref_distance * ref_energy) ** 0.5 - scaled_charges = charges / charge_factor - - if rigid_bodies: - rigid_bodies = [-1 if body is None else body for body in rigid_bodies] - else: - rigid_bodies = [-1 for _ in structure.atoms] - - return ( - n_particles, - scaled_positions, - unique_types, - typeids, - scaled_mass, - scaled_charges, - rigid_bodies, - ) - - -def _parse_pair_information(structure): - pair_types = [] - pair_typeid = [] - pairs = [] - for ai in structure.atoms: - for aj in ai.dihedral_partners: - # make sure we don't double add - if ai.idx > aj.idx: - ps = "-".join(sorted([ai.type, aj.type])) - if ps not in pair_types: - pair_types.append(ps) - pair_typeid.append(pair_types.index(ps)) - pairs.append((ai.idx, aj.idx)) - n_pairs = len(pairs) - - return pair_types, pair_typeid, pairs, n_pairs - - -def _parse_bond_information(structure): - n_bonds = len(structure.bonds) - unique_bond_types = set() - for bond in structure.bonds: - t1, t2 = bond.atom1.type, bond.atom2.type - if t1 == "" or t2 == "": - t1, t2 = bond.atom1.name, bond.atom2.name - t1, t2 = sorted([t1, t2], key=natural_sort) - try: - bond_type = "-".join((t1, t2)) - except AttributeError: # no forcefield applied, bond.type is None - bond_type = ("-".join((t1, t2)), 0.0, 0.0) - unique_bond_types.add(bond_type) - unique_bond_types = sorted(list(unique_bond_types), key=natural_sort) - - bond_typeids = [] - bond_groups = [] - for bond in structure.bonds: - t1, t2 = bond.atom1.type, bond.atom2.type - if t1 == "" or t2 == "": - t1, t2 = bond.atom1.name, bond.atom2.name - t1, t2 = sorted([t1, t2], key=natural_sort) - try: - bond_type = "-".join((t1, t2)) - except AttributeError: # no forcefield applied, bond.type is None - bond_type = ("-".join((t1, t2)), 0.0, 0.0) - bond_typeids.append(unique_bond_types.index(bond_type)) - bond_groups.append((bond.atom1.idx, bond.atom2.idx)) - - return n_bonds, unique_bond_types, bond_typeids, bond_groups - - -def _parse_angle_information(structure): - n_angles = len(structure.angles) - - unique_angle_types = set() - for angle in structure.angles: - t1, t2, t3 = angle.atom1.type, angle.atom2.type, angle.atom3.type - t1, t3 = sorted([t1, t3], key=natural_sort) - angle_type = "-".join((t1, t2, t3)) - unique_angle_types.add(angle_type) - unique_angle_types = sorted(list(unique_angle_types), key=natural_sort) - - angle_typeids = [] - angle_groups = [] - for angle in structure.angles: - t1, t2, t3 = angle.atom1.type, angle.atom2.type, angle.atom3.type - t1, t3 = sorted([t1, t3], key=natural_sort) - angle_type = "-".join((t1, t2, t3)) - angle_typeids.append(unique_angle_types.index(angle_type)) - angle_groups.append((angle.atom1.idx, angle.atom2.idx, angle.atom3.idx)) - - return n_angles, unique_angle_types, angle_typeids, angle_groups - - -def _parse_dihedral_information(structure): - n_dihedrals = len(structure.rb_torsions + structure.dihedrals) - - unique_dihedral_types = set() - for dihedral in structure.rb_torsions + structure.dihedrals: - t1, t2 = dihedral.atom1.type, dihedral.atom2.type - t3, t4 = dihedral.atom3.type, dihedral.atom4.type - if [t2, t3] == sorted([t2, t3], key=natural_sort): - dihedral_type = "-".join((t1, t2, t3, t4)) - else: - dihedral_type = "-".join((t4, t3, t2, t1)) - unique_dihedral_types.add(dihedral_type) - unique_dihedral_types = sorted( - list(unique_dihedral_types), key=natural_sort - ) - - dihedral_typeids = [] - dihedral_groups = [] - for dihedral in structure.rb_torsions + structure.dihedrals: - t1, t2 = dihedral.atom1.type, dihedral.atom2.type - t3, t4 = dihedral.atom3.type, dihedral.atom4.type - if [t2, t3] == sorted([t2, t3], key=natural_sort): - dihedral_type = "-".join((t1, t2, t3, t4)) - else: - dihedral_type = "-".join((t4, t3, t2, t1)) - dihedral_typeids.append(unique_dihedral_types.index(dihedral_type)) - dihedral_groups.append( - ( - dihedral.atom1.idx, - dihedral.atom2.idx, - dihedral.atom3.idx, - dihedral.atom4.idx, - ) - ) - - return n_dihedrals, unique_dihedral_types, dihedral_typeids, dihedral_groups - - -def _parse_improper_information(structure): - n_dihedrals = len(structure.impropers) - - unique_dihedral_types = set() - for dihedral in structure.impropers: - t1, t2 = dihedral.atom1.type, dihedral.atom2.type - t3, t4 = dihedral.atom3.type, dihedral.atom4.type - if [t2, t3] == sorted([t2, t3], key=natural_sort): - dihedral_type = "-".join((t1, t2, t3, t4)) - else: - dihedral_type = "-".join((t4, t3, t2, t1)) - unique_dihedral_types.add(dihedral_type) - unique_dihedral_types = sorted( - list(unique_dihedral_types), key=natural_sort - ) - - dihedral_typeids = [] - dihedral_groups = [] - for dihedral in structure.impropers: - t1, t2 = dihedral.atom1.type, dihedral.atom2.type - t3, t4 = dihedral.atom3.type, dihedral.atom4.type - if [t2, t3] == sorted([t2, t3], key=natural_sort): - dihedral_type = "-".join((t1, t2, t3, t4)) - else: - dihedral_type = "-".join((t4, t3, t2, t1)) - dihedral_typeids.append(unique_dihedral_types.index(dihedral_type)) - dihedral_groups.append( - ( - dihedral.atom1.idx, - dihedral.atom2.idx, - dihedral.atom3.idx, - dihedral.atom4.idx, - ) - ) - - return n_dihedrals, unique_dihedral_types, dihedral_typeids, dihedral_groups diff --git a/mbuild/formats/hoomdxml.py b/mbuild/formats/hoomdxml.py deleted file mode 100644 index 1d026b653..000000000 --- a/mbuild/formats/hoomdxml.py +++ /dev/null @@ -1,433 +0,0 @@ -"""HOOMD xml format.""" - -import operator -from collections import namedtuple -from math import radians - -import numpy as np - -from mbuild.utils.conversion import RB_to_OPLS -from mbuild.utils.decorators import deprecated -from mbuild.utils.geometry import coord_shift - -__all__ = ["write_hoomdxml"] - - -@deprecated("The .hoomdxml file format will be removed in mBuild 1.0.") -def write_hoomdxml( - structure, - filename, - ref_distance=1.0, - ref_mass=1.0, - ref_energy=1.0, - rigid_bodies=None, - shift_coords=True, - auto_scale=False, -): - """Output a HOOMD XML file. - - Parameters - ---------- - structure : parmed.Structure - ParmEd structure object - filename : str - Path of the output file. - ref_distance : float, optional, default=1.0, units=nanometers - Reference distance for conversion to reduced units - ref_mass : float, optional, default=1.0, units=amu - Reference mass for conversion to reduced units - ref_energy : float, optional, default=1.0, units=kJ/mol - Reference energy for conversion to reduced units - rigid_bodies : list - List of rigid body information. An integer value is required - for each particle corresponding to the number of the rigid body with - which the particle should be included. A value of None indicates the - particle is not part of any rigid body. - shift_coords : bool, optional, default=True - Shift coordinates from (0, L) to (-L/2, L/2) if necessary. - auto_scale : bool, optional, default=False - Automatically use largest sigma value as ref_distance, largest mass value - as ref_mass and largest epsilon value as ref_energy. - - Returns - ------- - ReferenceValues : namedtuple - Values used in scaling - - Example - ------- - - >>> ref_values = ethane.save( - ... filename='ethane-opls.hoomdxml', - ... forcefield_name='oplsaa', - ... auto_scale=True - ... ) - >>> print(ref_values.mass, ref_values.distance, ref_values.energy) - - Notes - ----- - The following elements are always written: - - * **position** : particle positions - * **type** : particle types - * **mass** : particle masses (default 1.0) - * **charge** : particle charges - - The following elements may be written if applicable: - - * **pair_coeffs** : Pair coefficients for each particle type (assumes a - 12-6 LJ pair style). The following information is written for each - particle type: - - * type : particle type - * epsilon : LJ epsilon - * sigma : LJ sigma - - * **bond_coeffs** : Coefficients for each bond type (assumes a harmonic bond - style). The following information is written for each bond type: - - * type : bond type - * k : force constant (units of energy/distance^2) - * r0 : bond rest length (units of distance) - - * **bond** : system bonds - * **angle_coeffs** : Coefficients for each angle type (assumes a harmonic - angle style). The following information is written for each angle type: - - * type : angle type - * k : force constant (units of energy/radians^2) - * theta : rest angle (units of radians) - - * **angle** : system angles - * **dihedral_coeffs** : Coefficients for each dihedral type (assumes an - OPLS dihedral style). The following information is written for each - dihedral type: - - * type : dihedral type - * k1, k2, k3, k4 : force coefficients (units of energy) - - * **dihedral** : system dihedrals - - * **body** : ID of the rigid body to which each particle belongs - """ - ref_distance *= 10 # Parmed unit hack - ref_energy /= 4.184 # Parmed unit hack - forcefield = True - if structure[0].type == "": - forcefield = False - if auto_scale and forcefield: - ref_mass = max([atom.mass for atom in structure.atoms]) - pair_coeffs = list( - set( - (atom.type, atom.epsilon, atom.sigma) - for atom in structure.atoms - ) - ) - ref_energy = max(pair_coeffs, key=operator.itemgetter(1))[1] - ref_distance = max(pair_coeffs, key=operator.itemgetter(2))[2] - - xyz = np.array([[atom.xx, atom.xy, atom.xz] for atom in structure.atoms]) - if shift_coords: - xyz = coord_shift(xyz, structure.box[:3]) - - with open(filename, "w") as xml_file: - xml_file.write('\n') - xml_file.write('\n') - xml_file.write( - "\n" - ) - # We need to convert back into nm and kJ - xml_file.write( - "\n".format( - ref_distance / 10, ref_mass, ref_energy * 4.184 - ) - ) - xml_file.write('\n') - _write_box_information(xml_file, structure, ref_distance) - _write_particle_information( - xml_file, - structure, - xyz, - forcefield, - ref_distance, - ref_mass, - ref_energy, - ) - _write_bond_information(xml_file, structure, ref_distance, ref_energy) - _write_angle_information(xml_file, structure, ref_energy) - _write_dihedral_information(xml_file, structure, ref_energy) - _write_rigid_information(xml_file, rigid_bodies) - xml_file.write("\n") - xml_file.write("") - - ReferenceValues = namedtuple("ref_values", ["distance", "mass", "energy"]) - - return ReferenceValues(ref_distance, ref_mass, ref_energy) - - -def _write_particle_information( - xml_file, structure, xyz, forcefield, ref_distance, ref_mass, ref_energy -): - """Write out the particle information. - - Parameters - ---------- - xml_file : file object - The file object of the hoomdxml file being written - structure : parmed.Structure - Parmed structure object - xyz : np.ndarray, shape=(n,3), dtype=float - The particle positions to be written. - forcefield : bool - If True, write the particle "type". Write the particle "name" otherwise. - ref_distance : float, default=1.0 - Reference distance for conversion to reduced units - ref_mass : float, default=1.0 - Reference mass for conversion to reduced units - ref_energy : float, default=1.0 - Reference energy for conversion to reduced units - """ - xml_file.write('\n'.format(xyz.shape[0])) - for pos in xyz: - xml_file.write("{}\t{}\t{}\n".format(*pos / ref_distance)) - xml_file.write("\n") - if forcefield: - types = [atom.type for atom in structure.atoms] - else: - types = [atom.name for atom in structure.atoms] - - xml_file.write("\n") - for atom_type in types: - xml_file.write("{}\n".format(atom_type)) - xml_file.write("\n") - - masses = [atom.mass for atom in structure.atoms] - xml_file.write("\n") - for mass in masses: - if mass == 0: - mass = 1.0 - xml_file.write("{}\n".format(mass / ref_mass)) - xml_file.write("\n") - - charges = [atom.charge for atom in structure.atoms] - xml_file.write("\n") - e0 = 2.396452e-04 # e^2 mol/(kcal A), permittivity of free space - charge_factor = (4.0 * np.pi * e0 * ref_distance * ref_energy) ** 0.5 - for charge in charges: - xml_file.write("{}\n".format(charge / charge_factor)) - xml_file.write("\n") - if forcefield: - pair_coeffs = list( - set( - (atom.type, atom.epsilon, atom.sigma) - for atom in structure.atoms - ) - ) - pair_coeffs.sort(key=lambda pair_type: pair_type[0]) - xml_file.write("\n") - for param_set in pair_coeffs: - xml_file.write( - "{}\t{:.4f}\t{:.4f}\n".format( - param_set[0], - param_set[1] / ref_energy, - param_set[2] / ref_distance, - ) - ) - xml_file.write("\n") - - -def _write_bond_information(xml_file, structure, ref_distance, ref_energy): - """Write the bonds in the system. - - Parameters - ---------- - xml_file : file object - The file object of the hoomdxml file being written - structure : parmed.Structure - Parmed structure object - ref_distance : float, default=1.0 - Reference distance for conversion to reduced units - ref_energy : float, default=1.0 - Reference energy for conversion to reduced units - """ - unique_bond_types = set() - xml_file.write("\n") - for bond in structure.bonds: - t1, t2 = bond.atom1.type, bond.atom2.type - if t1 == "" or t2 == "": - t1, t2 = bond.atom1.name, bond.atom2.name - t1, t2 = sorted([t1, t2]) - try: - bond_type = ("-".join((t1, t2)), bond.type.k, bond.type.req) - except AttributeError: # no forcefield applied, bond.type is None - bond_type = ("-".join((t1, t2)), 0.0, 0.0) - unique_bond_types.add(bond_type) - xml_file.write( - "{} {} {}\n".format(bond_type[0], bond.atom1.idx, bond.atom2.idx) - ) - xml_file.write("\n") - xml_file.write("\n") - xml_file.write("\n") - for bond_type, k, req in unique_bond_types: - xml_file.write( - "{} {} {}\n".format( - bond_type, - k * 2.0 / ref_energy * ref_distance**2.0, - req / ref_distance, - ) - ) - xml_file.write("\n") - - -def _write_angle_information(xml_file, structure, ref_energy): - """Write the angles in the system. - - Parameters - ---------- - xml_file : file object - The file object of the hoomdxml file being written - structure : parmed.Structure - Parmed structure object - ref_energy : float, default=1.0 - Reference energy for conversion to reduced units - """ - unique_angle_types = set() - xml_file.write("\n") - for angle in structure.angles: - t1, t2, t3 = angle.atom1.type, angle.atom2.type, angle.atom3.type - t1, t3 = sorted([t1, t3]) - angle_type = ("-".join((t1, t2, t3)), angle.type.k, angle.type.theteq) - unique_angle_types.add(angle_type) - xml_file.write( - "{} {} {} {}\n".format( - angle_type[0], angle.atom1.idx, angle.atom2.idx, angle.atom3.idx - ) - ) - xml_file.write("\n") - xml_file.write("\n") - xml_file.write("\n") - for angle_type, k, teq in unique_angle_types: - xml_file.write( - "{} {} {}\n".format(angle_type, k * 2.0 / ref_energy, radians(teq)) - ) - xml_file.write("\n") - - -def _write_dihedral_information(xml_file, structure, ref_energy): - """Write dihedrals in the system. - - Parameters - ---------- - xml_file : file object - The file object of the hoomdxml file being written - structure : parmed.Structure - Parmed structure object - ref_energy : float, default=1.0 - Reference energy for conversion to reduced units - """ - unique_dihedral_types = set() - xml_file.write("\n") - for dihedral in structure.rb_torsions: - t1, t2 = (dihedral.atom1.type, dihedral.atom2.type) - t3, t4 = dihedral.atom3.type, dihedral.atom4.type - if [t2, t3] == sorted([t2, t3]): - types_in_dihedral = "-".join((t1, t2, t3, t4)) - else: - types_in_dihedral = "-".join((t4, t3, t2, t1)) - dihedral_type = ( - types_in_dihedral, - dihedral.type.c0, - dihedral.type.c1, - dihedral.type.c2, - dihedral.type.c3, - dihedral.type.c4, - dihedral.type.c5, - dihedral.type.scee, - dihedral.type.scnb, - ) - unique_dihedral_types.add(dihedral_type) - xml_file.write( - "{} {} {} {} {}\n".format( - dihedral_type[0], - dihedral.atom1.idx, - dihedral.atom2.idx, - dihedral.atom3.idx, - dihedral.atom4.idx, - ) - ) - xml_file.write("\n") - xml_file.write("\n") - xml_file.write("\n") - for ( - dihedral_type, - c0, - c1, - c2, - c3, - c4, - c5, - scee, - scnb, - ) in unique_dihedral_types: - opls_coeffs = RB_to_OPLS( - c0, c1, c2, c3, c4, c5, error_if_outside_tolerance=False - ) - opls_coeffs /= ref_energy - xml_file.write( - "{} {:.5f} {:.5f} {:.5f} {:.5f}\n".format( - dihedral_type, - opls_coeffs[1], - opls_coeffs[2], - opls_coeffs[3], - opls_coeffs[4], - ) - ) - xml_file.write("\n") - - -def _write_rigid_information(xml_file, rigid_bodies): - """Write rigid body information. - - Parameters - ---------- - xml_file : file object - The file object of the hoomdxml file being written - rigid_bodies : list, len=n_particles - The rigid body that each particle belongs to (-1 for none) - """ - if not all(body is None for body in rigid_bodies): - xml_file.write("\n") - for body in rigid_bodies: - if body is None: - body = -1 - xml_file.write("{}\n".format(int(body))) - xml_file.write("\n") - - -def _write_box_information(xml_file, structure, ref_distance): - """Write box information. - - Parameters - ---------- - xml_file : file object - The file object of the hoomdxml file being written - structure : parmed.Structure - Parmed structure object - ref_energy : float, default=1.0 - Reference energy for conversion to reduced units - """ - if np.allclose(structure.box[3:6], np.array([90, 90, 90])): - box_str = '\n' - xml_file.write(box_str.format(*structure.box[:3] / ref_distance)) - else: - a, b, c = structure.box[0:3] / ref_distance - alpha, beta, gamma = np.radians(structure.box[3:6]) - - lx = a - xy = b * np.cos(gamma) - xz = c * np.cos(beta) - ly = np.sqrt(b**2 - xy**2) - yz = (b * c * np.cos(alpha) - xy * xz) / ly - lz = np.sqrt(c**2 - xz**2 - yz**2) - box_str = '\n' - xml_file.write(box_str.format(lx, ly, lz, xy, xz, yz)) diff --git a/mbuild/formats/lammpsdata.py b/mbuild/formats/lammpsdata.py deleted file mode 100644 index 5cd5846fb..000000000 --- a/mbuild/formats/lammpsdata.py +++ /dev/null @@ -1,1522 +0,0 @@ -"""LAMMPS data format.""" - -import itertools as it -from collections import OrderedDict -from warnings import warn - -import numpy as np -from parmed import Structure -from parmed.parameters import ParameterSet -from scipy.constants import epsilon_0 - -from mbuild import Box -from mbuild.utils.conversion import RB_to_OPLS -from mbuild.utils.decorators import deprecated -from mbuild.utils.orderedset import OrderedSet -from mbuild.utils.sorting import natural_sort - -__all__ = ["write_lammpsdata"] -# Define constants for conversions -KCAL_TO_KJ = 4.184 -NM2_TO_ANG2 = 100.0 -NM_TO_ANG = 10.0 -ELEM_TO_COUL = 1.602176634e-19 - - -dep_msg = """ -Support for writing out LAMMPS data files will be removed -in mbuild 1.0. -See GMSO (https://github.com/mosdef-hub/gmso/tree/main/gmso/formats/lammpsdata) for -continued support for LAMMPS. -""" -print(dep_msg) - - -@deprecated(dep_msg) -def write_lammpsdata( - structure, - filename, - atom_style="full", - unit_style="real", - mins=None, - maxs=None, - pair_coeff_label=None, - detect_forcefield_style=True, - nbfix_in_data_file=True, - use_urey_bradleys=False, - use_rb_torsions=True, - use_dihedrals=False, - sigma_conversion_factor=None, - epsilon_conversion_factor=None, - mass_conversion_factor=None, - charge_conversion_factor=True, - zero_dihedral_weighting_factor=False, - moleculeID_offset=1, -): - """Output a LAMMPS data file. - - Outputs a LAMMPS data file in the 'full' atom style format. Default units - are 'real' units. See http://lammps.sandia.gov/doc/atom_style.html for - more information on atom styles. - - Parameters - ---------- - structure : parmed.Structure - ParmEd structure object - filename : str - Path of the output file - atom_style: str, optional, default='full' - Defines the style of atoms to be saved in a LAMMPS data file. The - following atom styles are currently supported: - 'full', 'atomic', 'charge', 'molecular' - see http://lammps.sandia.gov/doc/atom_style.html for more information - on atom styles. - unit_style: str, optional, default='real' - Defines to unit style to be save in a LAMMPS data file. Defaults to - 'real' units. Current styles are supported: 'real', 'lj', 'metal' - see https://lammps.sandia.gov/doc/99/units.html for more information - on unit styles - mins : list, optional, default=None - Minimum box dimension in x, y, z directions, nm - maxs : list, optional, default=None - Maximum box dimension in x, y, z directions, nm - pair_coeff_label : str, optional, default=None - Provide a custom label to the pair_coeffs section in the lammps data - file. A value of None means a suitable default will be chosen. - detect_forcefield_style : bool, optional, default=True - If True, format lammpsdata parameters based on the contents of - the parmed Structure - use_urey_bradleys : bool, optional, default=False - If True, will treat angles as CHARMM-style angles with urey bradley - terms while looking for `structure.urey_bradleys` - use_rb_torsions : bool, optional, default=True - If True, will treat dihedrals OPLS-style torsions while looking for - `structure.rb_torsions` - use_dihedrals : bool, optional, default=False - If True, will treat dihedrals as CHARMM-style dihedrals while looking - for `structure.dihedrals` - zero_dihedral_weighting_factor : bool, optional, default=False - If True, will set weighting parameter to zero in CHARMM-style dihedrals. - This should be True if the CHARMM dihedral style is used in non-CHARMM forcefields. - sigma_conversion_factor : float, optional, default=None - If unit_style is set to 'lj', then sigma conversion factor is used to non-dimensionalize. - Assume to be in units of nm. If None, will take the largest sigma value in - the structure.atoms.sigma values. - epsilon_conversion_factor : float, optional, default=None - If unit_style is set to 'lj', then epsilon conversion factor is used to non-dimensionalize. - Assume to be in units of kCal/mol. If None, will take the largest epsilon value in - the structure.atoms.epsilon values. - mass_conversion_factor : float, optional, default=None - If unit_style is set to 'lj', then mass conversion factor is used to non-dimensionalize. - Assume to be in units of amu. If None, will take the largest mass value in - the structure.atoms.mass values. - charge_conversion_factor : bool, optional, default=True - If unit_style is set to 'lj', then charge conversion factor may or may not be used to - non-dimensionalize. Assume to be in elementary charge units. If ``True``, the charges - are scaled by ``np.sqrt(4*np.pi()*eps_0*sigma_conversion_factor*epsilon_conversion_factor)``. - If ``False``, the charges are not scaled and the user must be wary to choose the dielectric - constant properly, which may be more convenient to implement an implicit solvent. - moleculeID_offset : int , optional, default=1 - Since LAMMPS treats the MoleculeID as an additional set of information - to identify what molecule an atom belongs to, this currently - behaves as a residue id. This value needs to start at 1 to be - considered a real molecule. - - Notes - ----- - See http://lammps.sandia.gov/doc/2001/data_format.html for a full - description of the LAMMPS data format. Currently the following sections are - supported (in addition to the header): *Masses*, *Nonbond Coeffs*, - *Bond Coeffs*, *Angle Coeffs*, *Dihedral Coeffs*, *Atoms*, *Bonds*, - *Angles*, *Dihedrals*, *Impropers* - OPLS and CHARMM forcefield styles are supported, AMBER forcefield styles - are NOT - - Some of this function has beed adopted from `mdtraj`'s support of the - LAMMPSTRJ trajectory format. See - https://github.com/mdtraj/mdtraj/blob/master/mdtraj/formats/lammpstrj.py - for details. - - unique_types : a sorted list of unique atomtypes for all atoms in the structure. - Defined by: - atomtype : atom.type - unique_bond_types: an enumarated OrderedDict of unique bond types for all bonds in the structure. - Defined by bond parameters and component atomtypes, in order: - --- k : bond.type.k - --- req : bond.type.req - --- atomtypes : sorted((bond.atom1.type, bond.atom2.type)) - unique_angle_types: an enumerated OrderedDict of unique angle types for all - angles in the structure. - Defined by angle parameters and component atomtypes, in order: - --- k : angle.type.k - --- theteq : angle.type.theteq - --- vertex atomtype: angle.atom2.type - --- atomtypes: sorted((bond.atom1.type, bond.atom3.type)) - - unique_dihedral_types: an enumerated OrderedDict of unique dihedrals type - for all dihedrals in the structure. - Defined by dihedral parameters and component atomtypes, in order: - --- c0 : dihedral.type.c0 - --- c1 : dihedral.type.c1 - --- c2 : dihedral.type.c2 - --- c3 : dihedral.type.c3 - --- c4 : dihedral.type.c4 - --- c5 : dihedral.type.c5 - --- scee : dihedral.type.scee - --- scnb : dihedral.type.scnb - --- atomtype 1 : dihedral.atom1.type - --- atomtype 2 : dihedral.atom2.type - --- atomtype 3 : dihedral.atom3.type - --- atomtype 4 : dihedral.atom4.type - """ - if atom_style not in ["atomic", "charge", "molecular", "full"]: - raise ValueError( - 'Atom style "{}" is invalid or is not currently supported'.format( - atom_style - ) - ) - if unit_style not in ["real", "lj", "metal"]: - raise ValueError( - 'Unit style "{}" is invalid or is not currently supported'.format( - unit_style - ) - ) - - forcefield = True - if structure[0].type == "": - forcefield = False - # copy structure so the input structure isn't modified in-place - structure = structure.copy(cls=Structure, split_dihedrals=True) - - # units of angstroms - xyz = np.array([[atom.xx, atom.xy, atom.xz] for atom in structure.atoms]) - - if forcefield: - types = [atom.type for atom in structure.atoms] - else: - types = [atom.name for atom in structure.atoms] - - unique_types = list(set(types)) - unique_types.sort(key=natural_sort) - - charges = np.array([atom.charge for atom in structure.atoms]) - - # lammps does not require the box to be centered at any a specific origin - # min and max dimensions are therefore needed to write the file in a - # consistent way the parmed structure only stores the box length. It is - # not rigorous to assume bounds are 0 to L or -L/2 to L/2 - - # NOTE: 0 to L is current default, mins and maxs should be passed by user - - if _check_minsmaxs(mins, maxs): - mins = np.array(mins) * NM_TO_ANG - maxs = np.array(maxs) * NM_TO_ANG - box = Box.from_mins_maxs_angles( - mins=mins, maxs=maxs, angles=structure.box[3:6] - ) # box lengths input in nm, so convert to angstrom - else: # Parmed internally converts box to angstroms - box = Box( - lengths=np.array([val for val in structure.box[0:3]]), - angles=structure.box[3:6], - ) - - warn( - "Explicit box bounds (i.e., mins and maxs) were not provided. Box " - "bounds are assumed to be min = 0 and max = length in each " - "direction. This may not produce a system with the expected " - "spatial location and may cause non-periodic systems to fail. " - "Bounds can be defined explicitly by passing the them to the " - "write_lammpsdata function or by passing box info to the save " - "function." - ) - - # Lammps syntax depends on the functional form - # Infer functional form based on the properties of the structure - if detect_forcefield_style: - # Check angles - if len(structure.urey_bradleys) > 0: - print("Urey bradley terms detected, will use angle_style charmm") - use_urey_bradleys = True - else: - print( - "No urey bradley terms detected, will use angle_style harmonic" - ) - use_urey_bradleys = False - - # Check dihedrals - if len(structure.rb_torsions) > 0: - print("RB Torsions detected, will use dihedral_style opls") - use_rb_torsions = True - else: - use_rb_torsions = False - if len([d for d in structure.dihedrals if not d.improper]) > 0: - print("Charmm dihedrals detected, will use dihedral_style charmm") - use_dihedrals = True - else: - use_dihedrals = False - if use_rb_torsions and use_dihedrals: - raise ValueError( - "Multiple dihedral styles detected, check your " - "Forcefield XML and structure" - ) - - # save atom index information for all bonded params in structure - bonds = [[b.atom1.idx + 1, b.atom2.idx + 1] for b in structure.bonds] - angles = [ - [angle.atom1.idx + 1, angle.atom2.idx + 1, angle.atom3.idx + 1] - for angle in structure.angles - ] - if use_rb_torsions: - dihedrals = [ - [d.atom1.idx + 1, d.atom2.idx + 1, d.atom3.idx + 1, d.atom4.idx + 1] - for d in structure.rb_torsions - ] - elif use_dihedrals: - dihedrals = [ - [d.atom1.idx + 1, d.atom2.idx + 1, d.atom3.idx + 1, d.atom4.idx + 1] - for d in structure.dihedrals - if not d.improper - ] - else: - dihedrals = [] - impropers = [ - [i.atom1.idx + 1, i.atom2.idx + 1, i.atom3.idx + 1, i.atom4.idx + 1] - for i in structure.impropers - ] - imp_dihedrals = [ - [d.atom1.idx + 1, d.atom2.idx + 1, d.atom3.idx + 1, d.atom4.idx + 1] - for d in structure.dihedrals - if d.improper - ] - - if impropers and imp_dihedrals: - raise ValueError("Use of multiple improper styles is not supported") - - # Get lj conversion factors if they exist and apply lj params to charges and box - # Params are in read in mbuild units of angstrom, kcal/mol, amu and converted to lammps units - if unit_style == "lj": - sigma_conversion_factor = _evaluate_lj_conversion_factors( - structure, "sigma", sigma_conversion_factor - ) - epsilon_conversion_factor = _evaluate_lj_conversion_factors( - structure, "epsilon", epsilon_conversion_factor - ) - mass_conversion_factor = _evaluate_lj_conversion_factors( - structure, "mass", mass_conversion_factor - ) - # Convert coordinates and charges to LJ units - xyz = xyz / sigma_conversion_factor - if charge_conversion_factor: - charges = (charges * ELEM_TO_COUL) / np.sqrt( - 4 - * np.pi - * sigma_conversion_factor - * NM_TO_ANG**-1 - * epsilon_conversion_factor - * KCAL_TO_KJ - * epsilon_0 - * 10**-6 - ) - charges[np.isinf(charges)] = 0 - eng_unit_str = " " - - elif unit_style == "real": - sigma_conversion_factor = 1 - epsilon_conversion_factor = 1 - mass_conversion_factor = 1 - eng_unit_str = "kcal/mol" - - elif unit_style == "metal": - sigma_conversion_factor = 1 # unit in Angstrom - epsilon_conversion_factor = 1 / 0.043364115 # kcal/mol to eV - mass_conversion_factor = 1 - eng_unit_str = "eV" - else: - raise ValueError( - "Currently only support 'real', 'lj', and 'metal' unit_styles." - ) - - # Divide by conversion factor - Lx = box.Lx * (1 / sigma_conversion_factor) - Ly = box.Ly * (1 / sigma_conversion_factor) - Lz = box.Lz * (1 / sigma_conversion_factor) - box = Box(lengths=(Lx, Ly, Lz), angles=box.angles) - - # Get bonded parameter information from structure - if bonds: - if len(structure.bond_types) == 0: - bond_types = np.ones(len(bonds), dtype=int) - else: - bond_types, unique_bond_types = _get_bond_types( - structure, - bonds, - unique_types, - sigma_conversion_factor, - epsilon_conversion_factor, - ) - - if angles: - angle_types, unique_angle_types = _get_angle_types( - structure, - unique_types, - use_urey_bradleys, - sigma_conversion_factor, - epsilon_conversion_factor, - ) - - if imp_dihedrals: - ( - imp_dihedral_types, - unique_imp_dihedral_types, - ) = _get_improper_dihedral_types( - structure, - unique_types, - epsilon_conversion_factor, - ) - - if dihedrals: - dihedral_types, unique_dihedral_types = _get_dihedral_types( - structure, - unique_types, - use_rb_torsions, - use_dihedrals, - epsilon_conversion_factor, - zero_dihedral_weighting_factor, - ) - - if impropers: - improper_types, unique_improper_types = _get_impropers( - structure, - unique_types, - epsilon_conversion_factor, - ) - - # Write lammps data file https://docs.lammps.org/2001/data_format.html - with open(filename, "w") as data: - data.write(f"{filename} - created by mBuild; units = {unit_style}\n") - if unit_style == "lj": - data.write("#Normalization factors: ") - data.write( - "sigma - {:.3E} (angstrom), ".format(sigma_conversion_factor) - ) - data.write( - "epsilon - {:.3E} (kcal/mol), ".format( - epsilon_conversion_factor - ) - ) - data.write("mass - {:.3E} (amu)".format(mass_conversion_factor)) - data.write("\n") - - # Write counts of bonded interactions - data.write("{:d} atoms\n".format(len(structure.atoms))) - if atom_style in ["full", "molecular"]: - data.write("{:d} bonds\n".format(len(bonds))) - data.write("{:d} angles\n".format(len(angles))) - data.write("{:d} dihedrals\n".format(len(dihedrals))) - data.write( - "{:d} impropers\n\n".format(len(impropers) + len(imp_dihedrals)) - ) - - # Write counts of unique bonded types - data.write("{:d} atom types\n".format(len(set(types)))) - if atom_style in ["full", "molecular"]: - if bonds: - data.write("{:d} bond types\n".format(len(set(bond_types)))) - if angles: - data.write("{:d} angle types\n".format(len(set(angle_types)))) - if dihedrals: - data.write( - "{:d} dihedral types\n".format(len(set(dihedral_types))) - ) - if impropers: - data.write( - "{:d} improper types\n".format(len(set(improper_types))) - ) - elif imp_dihedrals: - data.write( - "{:d} improper types\n".format(len(set(imp_dihedral_types))) - ) - - data.write("\n") - # Write box data - # NOTE: Needs better logic handling maxs and mins of a bounding box - # NOTE: CCC, "could be an option to grab mins and maxs of the structure boundingbox" - # Lammps supports any box origin, so it is okay if it matches the structure positions - _write_box_information(box, data, mins) - - # Write mass data - _write_mass_information( - structure, - data, - mass_conversion_factor, - forcefield, - unique_types, - types, - ) - - if forcefield: - # Write pair coefficients data - _write_pair_information( - structure, - data, - forcefield, - sigma_conversion_factor, - epsilon_conversion_factor, - unique_types, - types, - pair_coeff_label, - unit_style, - nbfix_in_data_file, - eng_unit_str, - ) - - # Write bond coefficients - if bonds: - _write_bond_information( - structure, data, unique_bond_types, unit_style, eng_unit_str - ) - - # Write angle coefficients - if angles: - _write_angle_information( - structure, - data, - unique_angle_types, - use_urey_bradleys, - unit_style, - eng_unit_str, - ) - - # Write dihedral coefficients - if dihedrals: - _write_dihedral_information( - structure, - data, - unique_dihedral_types, - unit_style, - use_rb_torsions, - use_dihedrals, - eng_unit_str, - ) - - # Write improper coefficients - if impropers: - _write_improper_information( - structure, data, unique_improper_types, unit_style - ) - # Write improper dihedrals - elif imp_dihedrals: - _write_imp_dihedral_information( - structure, data, unique_imp_dihedral_types, unit_style - ) - - # Write Atom data - data.write("\nAtoms # {}\n\n".format(atom_style)) - if atom_style == "atomic": - atom_line = "{index:d}\t{type_index:d}\t{x:.6f}\t{y:.6f}\t{z:.6f}\n" - elif atom_style == "charge": - if unit_style in ["real", "metal"]: - atom_line = "{index:d}\t{type_index:d}\t{charge:.6f}\t{x:.6f}\t{y:.6f}\t{z:.6f}\n" - elif unit_style == "lj": - atom_line = "{index:d}\t{type_index:d}\t{charge:.4ef}\t{x:.6f}\t{y:.6f}\t{z:.6f}\n" - elif atom_style == "molecular": - atom_line = "{index:d}\t{zero:d}\t{type_index:d}\t{x:.6f}\t{y:.6f}\t{z:.6f}\n" - elif atom_style == "full": - if unit_style in ["real", "metal"]: - atom_line = "{index:d}\t{zero:d}\t{type_index:d}\t{charge:.6f}\t{x:.6f}\t{y:.6f}\t{z:.6f}\n" - elif unit_style == "lj": - atom_line = "{index:d}\t{zero:d}\t{type_index:d}\t{charge:.4e}\t{x:.6f}\t{y:.6f}\t{z:.6f}\n" - - for i, coords in enumerate(xyz): - data.write( - atom_line.format( - index=i + 1, - type_index=unique_types.index(types[i]) + 1, - zero=structure.atoms[i].residue.idx + moleculeID_offset, - charge=charges[i], - x=coords[0], - y=coords[1], - z=coords[2], - ) - ) - - if atom_style in ["full", "molecular"]: - # Bond data - if bonds: - data.write("\nBonds\n\n") - for i, bond in enumerate(bonds): - data.write( - "{:d}\t{:d}\t{:d}\t{:d}\n".format( - i + 1, bond_types[i], bond[0], bond[1] - ) - ) - - # Angle data - if angles: - data.write("\nAngles\n\n") - for i, angle in enumerate(angles): - data.write( - "{:d}\t{:d}\t{:d}\t{:d}\t{:d}\n".format( - i + 1, angle_types[i], angle[0], angle[1], angle[2] - ) - ) - - # Dihedral data - if dihedrals: - data.write("\nDihedrals\n\n") - for i, dihedral in enumerate(dihedrals): - data.write( - "{:d}\t{:d}\t{:d}\t{:d}\t{:d}\t{:d}\n".format( - i + 1, - dihedral_types[i], - dihedral[0], - dihedral[1], - dihedral[2], - dihedral[3], - ) - ) - # Improper data for harmonic impropers - if impropers: - data.write("\nImpropers\n\n") - for i, improper in enumerate(impropers): - # The atoms are written central-atom third in LAMMPS data file. - # The first atom in structure.impropers is the central atom - data.write( - "{:d}\t{:d}\t{:d}\t{:d}\t{:d}\t{:d}\n".format( - i + 1, - improper_types[i], - improper[2], - improper[1], - improper[0], - improper[3], - ) - ) - elif imp_dihedrals: # impropers in structure.dihedrals are cvff - data.write("\nImpropers\n\n") - for i, improper in enumerate(imp_dihedrals): - # The atoms are written central-atom third in LAMMPS data file. - # The third atom in structure.dihedral is the central atom - # if dihedral.improper is True - # This is correct for AMBER impropers even though - # LAMMPS documentation implies central-atom-first. - data.write( - "{:d}\t{:d}\t{:d}\t{:d}\t{:d}\t{:d}\n".format( - i + 1, - imp_dihedral_types[i], - improper[0], - improper[1], - improper[2], - improper[3], - ) - ) - - -def _evaluate_lj_conversion_factors( - structure, conversion_name, conversion_factor -): - """Get Lennard Jones style conversion factors. `conversion_name`` can be sigma, epsilon, or mass.""" - if conversion_factor is None: - # Check if structure is parametrized - if any([atom.sigma for atom in structure.atoms]) is None: - raise ValueError( - "LJ units specified but one or more atoms has undefined LJ " - "parameters." - ) - else: - conversion_factor = np.max( - [getattr(atom, conversion_name) for atom in structure.atoms] - ) - warn( - f"Assuming {conversion_name} conversion factor of " - + str(conversion_factor) - ) - if conversion_factor == 0: - conversion_factor = 1 - warn( - f"{conversion_name} conversion factor cannot be inferred from the maximum {conversion_name} value in the ParmEd Structure. " - "Setting the {conversion_name} conversion factor to 1" - ) - elif conversion_factor <= 0: - raise ValueError( - f"The {conversion_name} conversion factor to convert to LJ units should be greater than 0." - ) - else: - # assume conversion factor passed in mbuild units, convert to lammps real units - if conversion_name == "sigma": - conversion_factor *= NM_TO_ANG - elif conversion_name == "epsilon": - conversion_factor *= 1 / KCAL_TO_KJ - elif conversion_name == "mass": - pass - return conversion_factor - - -def _check_minsmaxs(mins, maxs): - """Return True if both mins and maxs have been defined, and each have length 3 otherwise returns False.""" - if mins and maxs: - if len(mins) == 3 and len(maxs) == 3: - return True - else: - warn( - "mins and maxs passed to write_lammpsdata, but list size is " - "incorrect. mins and maxs will be ignored." - ) - return False - else: - return False - - -def _get_bond_types( - structure, - bonds, - atom_types, - sigma_conversion_factor, - epsilon_conversion_factor, - bond_precision=3, -): - """Will get the bond types from a parmed structure and convert them to lammps real units.""" - unique_bond_types = OrderedDict( - enumerate( - OrderedSet( - *[ - ( - round( - bond.type.k - * ( - sigma_conversion_factor**2 - / epsilon_conversion_factor - ), - bond_precision, - ), - round( - bond.type.req / sigma_conversion_factor, - bond_precision, - ), - tuple(sorted((bond.atom1.type, bond.atom2.type))), - ) - for bond in structure.bonds - ] - ) - ) - ) - - magnitude = np.ceil(np.log10(len(atom_types))) - bond_tuples = [x[2] for x in unique_bond_types.values()] - ordered_bond_tuples = bond_tuples.copy() - ordered_bond_tuples.sort( - key=lambda x: atom_types.index(x[0]) * 10**magnitude - + atom_types.index(x[1]) - ) - - unique_bond_types = OrderedDict( - [ - (unique_bond_types[bond_tuples.index(x)], i + 1) - for i, x in enumerate(ordered_bond_tuples) - ] - ) - - bond_types = [ - unique_bond_types[ - ( - round( - bond.type.k - * (sigma_conversion_factor**2 / epsilon_conversion_factor), - bond_precision, - ), - round(bond.type.req / sigma_conversion_factor, bond_precision), - tuple(sorted((bond.atom1.type, bond.atom2.type))), - ) - ] - for bond in structure.bonds - ] - return bond_types, unique_bond_types - - -def _get_angle_types( - structure, - atom_types, - use_urey_bradleys, - sigma_conversion_factor, - epsilon_conversion_factor, - angle_precision=3, -): - """ - Will get the angle types from a parmed structure and convert them to lammps real units. - - Can get the parameters if urey bradleys or harmonic angles. - """ - if use_urey_bradleys: - charmm_angle_types = [] - for angle in structure.angles: - ub_k = 0 - ub_req = 0 - for ub in structure.urey_bradleys: - if (angle.atom1, angle.atom3) == (ub.atom1, ub.atom2): - ub_k = ub.type.k - ub_req = ub.type.req - charmm_angle_types.append( - ( - round( - angle.type.k - * ( - sigma_conversion_factor**2 - / epsilon_conversion_factor - ), - angle_precision, - ), - round(angle.type.theteq, angle_precision), - round(ub_k / epsilon_conversion_factor, angle_precision), - round(ub_req, angle_precision), - tuple(sorted((angle.atom1.type, angle.atom3.type))), - ) - ) - - unique_angle_types = OrderedDict( - enumerate(OrderedSet(*charmm_angle_types)) - ) - - magnitude = np.ceil(np.log10(len(atom_types))) - angle_tuples = [x[-1] for x in unique_angle_types.values()] - ordered_angle_tuples = angle_tuples.copy() - ordered_angle_tuples.sort( - key=lambda x: atom_types.index(x[0]) * 10**magnitude - + atom_types.index(x[1]) - ) - - else: - unique_angle_types = OrderedDict( - enumerate( - OrderedSet( - *[ - ( - round( - angle.type.k * (1 / epsilon_conversion_factor), - angle_precision, - ), - round(angle.type.theteq, angle_precision), - angle.atom2.type, - tuple(sorted((angle.atom1.type, angle.atom3.type))), - ) - for angle in structure.angles - ] - ) - ) - ) - - magnitude = np.ceil(np.log10(len(atom_types))) - angle_tuples = [ - (x[2], x[-1][0], x[-1][1]) for x in unique_angle_types.values() - ] - ordered_angle_tuples = angle_tuples.copy() - ordered_angle_tuples.sort( - key=lambda x: atom_types.index(x[0]) * 10 ** (2 * magnitude) - + atom_types.index(x[1]) * 10**magnitude - + atom_types.index(x[2]) - ) - - unique_angle_types = OrderedDict( - [ - (unique_angle_types[angle_tuples.index(x)], i + 1) - for i, x in enumerate(ordered_angle_tuples) - ] - ) - - if use_urey_bradleys: - angle_types = [ - unique_angle_types[ub_info] for ub_info in charmm_angle_types - ] - else: - angle_types = [ - unique_angle_types[ - ( - round( - angle.type.k * (1 / epsilon_conversion_factor), - angle_precision, - ), - round(angle.type.theteq, angle_precision), - angle.atom2.type, - tuple(sorted((angle.atom1.type, angle.atom3.type))), - ) - ] - for angle in structure.angles - ] - - return angle_types, unique_angle_types - - -typeLister = lambda x: (x.atom1.type, x.atom2.type, x.atom3.type, x.atom4.type) - - -def _dihedral_sorter(dihedral): - if dihedral.atom2.type < dihedral.atom3.type or ( - dihedral.atom2.type == dihedral.atom3.type - and dihedral.atom1.type < dihedral.atom4.type - ): - return typeLister(dihedral) - else: - return tuple(reversed(typeLister(dihedral))) - - -def _get_dihedral_types( - structure, - atom_types, - use_rb_torsions, - use_dihedrals, - epsilon_conversion_factor, - zero_dihedral_weighting_factor, - dihedral_precision=5, -): - """ - Will get the dihedral types from a parmed structure and convert them to lammps real units. - - Can be in the form of rb_torsions or charmm dihedrals. - """ - lj_unit = 1.0 / epsilon_conversion_factor - if use_rb_torsions: - unique_dihedral_types = set( - [ - ( - round(dihedral.type.c0 * lj_unit, dihedral_precision), - round(dihedral.type.c1 * lj_unit, dihedral_precision), - round(dihedral.type.c2 * lj_unit, dihedral_precision), - round(dihedral.type.c3 * lj_unit, dihedral_precision), - round(dihedral.type.c4 * lj_unit, dihedral_precision), - round(dihedral.type.c5 * lj_unit, dihedral_precision), - round(dihedral.type.scee, 1), - round(dihedral.type.scnb, 1), - *_dihedral_sorter(dihedral), - ) - for dihedral in structure.rb_torsions - ] - ) - unique_dihedral_typesDict = { - v: k + 1 - for k, v in enumerate( - sorted( - unique_dihedral_types, - key=lambda item: tuple(item[i] for i in [-3, -2, -4, -1]), - ) - ) - } - - elif use_dihedrals: - structure.join_dihedrals() # creates a list of dihedrals for dihedral.type - weight = 0.0 - if not zero_dihedral_weighting_factor: - weight = 1.0 - charmm_dihedralsList = [] - for dihedral in structure.dihedrals: - if not dihedral.improper: - for dih_type in dihedral.type: - charmm_dihedralsList.append( - ( - dih_type.phi_k * lj_unit, - int(round(dih_type.per, 0)), - round(dih_type.phase, 3), - round(weight / len(dihedral.type), 4), - round(dih_type.scee, 1), - round(dih_type.scnb, 1), - *_dihedral_sorter(dihedral), - ) - ) - unique_dihedral_typesDict = { - v: k + 1 - for k, v in enumerate( - sorted( - charmm_dihedralsList, - key=lambda item: tuple(item[i] for i in [-3, -2, -4, -1]), - ) - ) - } - - if use_rb_torsions: - dihedral_types = [ - unique_dihedral_typesDict[ - ( - round(dihedral.type.c0 * lj_unit, dihedral_precision), - round(dihedral.type.c1 * lj_unit, dihedral_precision), - round(dihedral.type.c2 * lj_unit, dihedral_precision), - round(dihedral.type.c3 * lj_unit, dihedral_precision), - round(dihedral.type.c4 * lj_unit, dihedral_precision), - round(dihedral.type.c5 * lj_unit, dihedral_precision), - round(dihedral.type.scee, 1), - round(dihedral.type.scnb, 1), - *_dihedral_sorter(dihedral), - ) - ] - for dihedral in structure.rb_torsions - ] - elif use_dihedrals: - dihedral_types = [] - for dihedral in structure.dihedrals: - if not dihedral.improper: - for dih_type in dihedral.type: - dihedral_types.append( - unique_dihedral_typesDict[ - ( - dih_type.phi_k * lj_unit, - int(round(dih_type.per, 0)), - round(dih_type.phase, 3), - round(weight / len(dihedral.type), 4), - round(dih_type.scee, 1), - round(dih_type.scnb, 1), - *_dihedral_sorter(dihedral), - ) - ] - ) - - return dihedral_types, unique_dihedral_typesDict - - -def _get_improper_dihedral_types( - structure, atom_types, epsilon_conversion_factor, imp_dih_precision=3 -): - """ - Will get the improper types from a parmed structure and convert them to lammps real units. - - Type cvff https://docs.lammps.org/improper_cvff.html - """ - lj_unit = 1 / epsilon_conversion_factor - improper_dihedrals = [] - for dihedral in structure.dihedrals: - if dihedral.improper: - dih_type = dihedral.type - phase = abs(int(round(dih_type.phase, 0))) - if not (phase == 0 or phase == 180): - raise ValueError("Improper dihedral phase must be 0 or 180") - if phase: - d = -1 - else: - d = 1 - improper_dihedrals.append( - ( - round(dih_type.phi_k * lj_unit, imp_dih_precision), - d, - int(round(dih_type.per, 0)), - round(dih_type.scee, 1), - round(dih_type.scnb, 1), - dihedral.atom1.type, - dihedral.atom2.type, - dihedral.atom3.type, - dihedral.atom4.type, - ) - ) - unique_imp_dihedral_types = dict(enumerate(OrderedSet(*improper_dihedrals))) - magnitude = np.ceil(np.log10(len(atom_types))) - dihedral_tuples = [ - (x[1], x[5:]) for x in unique_imp_dihedral_types.values() - ] - ordered_dihedral_tuples = dihedral_tuples.copy() - ordered_dihedral_tuples.sort( - key=lambda x: atom_types.index(x[1][0]) * 10 ** (4 * magnitude) - + atom_types.index(x[1][1]) * 10 ** (3 * magnitude) - + atom_types.index(x[1][2]) * 10 ** (2 * magnitude) - + atom_types.index(x[1][3]) * 10**magnitude - + x[0] - ) - - unique_imp_dihedral_types = OrderedDict( - [ - (unique_imp_dihedral_types[dihedral_tuples.index(x)], i + 1) - for i, x in enumerate(ordered_dihedral_tuples) - ] - ) - imp_dihedral_types = [ - unique_imp_dihedral_types[dihedral_info] - for dihedral_info in improper_dihedrals - ] - - return imp_dihedral_types, unique_imp_dihedral_types - - -def _get_impropers( - structure, atom_types, epsilon_conversion_factor, improper_precision=3 -): - """ - Will get the improper types from a parmed structure and convert them to lammps real units. - - Type harmonic https://docs.lammps.org/improper_harmonic.html. - """ - lj_unit = 1 / epsilon_conversion_factor - unique_improper_types = dict( - enumerate( - OrderedSet( - *[ - ( - round( - improper.type.psi_k * lj_unit, improper_precision - ), - round(improper.type.psi_eq, improper_precision), - improper.atom3.type, - improper.atom2.type, - improper.atom1.type, - improper.atom4.type, - ) - for improper in structure.impropers - ] - ) - ) - ) - - magnitude = np.ceil(np.log10(len(atom_types))) - improper_tuples = [x[2:] for x in unique_improper_types.values()] - ordered_improper_tuples = improper_tuples.copy() - ordered_improper_tuples.sort( - key=lambda x: atom_types.index(x[0]) * 10 ** (4 * magnitude) - + atom_types.index(x[1]) * 10 ** (3 * magnitude) - + atom_types.index(x[2]) * 10 ** (2 * magnitude) - + atom_types.index(x[3]) * 10**magnitude - ) - - unique_improper_types = OrderedDict( - [ - (unique_improper_types[improper_tuples.index(x)], i + 1) - for i, x in enumerate(ordered_improper_tuples) - ] - ) - - improper_types = [ - unique_improper_types[ - ( - round(improper.type.psi_k * lj_unit, improper_precision), - round(improper.type.psi_eq, improper_precision), - improper.atom3.type, - improper.atom2.type, - improper.atom1.type, - improper.atom4.type, - ) - ] - for improper in structure.impropers - ] - - return improper_types, unique_improper_types - - -def _write_box_information(box, data, mins): - """Write box information to lammps data file, can handle non-orthogonal boxes.""" - if np.allclose(box.angles, 90.0, atol=1e-5) and (mins is None): - for i, dim in enumerate(["x", "y", "z"]): - data.write( - "{0:.6f} {1:.6f} {2}lo {2}hi\n".format( - 0.0, - box.lengths[i], - dim, - ) - ) - # NOTE: - # currently non-orthogonal bounding box translates - # Compound such that mins are new origin - else: - a = box.Lx - b = box.Ly - c = box.Lz - # alpha, beta, gamma = np.radians(box.angles) - - xy = box.xy - xz = box.xz - yz = box.yz - - # NOTE: using (0,0,0) as origin - xlo, ylo, zlo = (0.0, 0.0, 0.0) - xhi = xlo + a - yhi = ylo + b - zhi = zlo + c - - xlo_bound = xlo + np.min([0.0, xy, xz, xy + xz]) - xhi_bound = xhi + np.max([0.0, xy, xz, xy + xz]) - ylo_bound = ylo + np.min([0.0, yz]) - yhi_bound = yhi + np.max([0.0, yz]) - zlo_bound = zlo - zhi_bound = zhi - - data.write("{0:.6f} {1:.6f} xlo xhi\n".format(xlo_bound, xhi_bound)) - data.write("{0:.6f} {1:.6f} ylo yhi\n".format(ylo_bound, yhi_bound)) - data.write("{0:.6f} {1:.6f} zlo zhi\n".format(zlo_bound, zhi_bound)) - data.write("{0:.6f} {1:.6f} {2:6f} xy xz yz\n".format(xy, xz, yz)) - - -def _write_mass_information( - structure, data, mass_conversion_factor, forcefield, unique_types, types -): - """Write mass information from structure to lammps data file.""" - if not forcefield: - masses = ( - np.array([atom.mass for atom in structure.atoms]) - / mass_conversion_factor - ) - else: - tmp_masses = list() - for atom in structure.atoms: - # handle case where atomtype does not contain a mass - try: - tmp_masses.append(atom.atom_type.mass.in_units("g/mol").value) - except AttributeError: - warn( - f"No mass or defined atomtype for atom: {atom}. Using atom mass of {atom.mass / mass_conversion_factor}" - ) - tmp_masses.append(atom.mass) - masses = np.asarray(tmp_masses) / mass_conversion_factor - - mass_dict = OrderedDict( - [ - (unique_types.index(atom_type) + 1, mass) - for atom_type, mass in zip(types, masses) - ] - ) - data.write("\nMasses\n\n") - for atom_type, mass in sorted(mass_dict.items()): - data.write( - "{:d}\t{:.6f}\t# {}\n".format( - atom_type, mass, unique_types[atom_type - 1] - ) - ) - - -def _write_pair_information( - structure, - data, - forcefield, - sigma_conversion_factor, - epsilon_conversion_factor, - unique_types, - types, - pair_coeff_label, - unit_style, - nbfix_in_data_file, - eng_unit_str, -): - """Write nonbonded pair information to lammps data file.""" - epsilons = ( - np.array([atom.epsilon for atom in structure.atoms]) - / epsilon_conversion_factor - ) - sigmas = ( - np.array([atom.sigma for atom in structure.atoms]) - / sigma_conversion_factor - ) - forcefields = [atom.type for atom in structure.atoms] - epsilon_dict = dict( - [ - (unique_types.index(atom_type) + 1, epsilon) - for atom_type, epsilon in zip(types, epsilons) - ] - ) - sigma_dict = dict( - [ - (unique_types.index(atom_type) + 1, sigma) - for atom_type, sigma in zip(types, sigmas) - ] - ) - forcefield_dict = dict( - [ - (unique_types.index(atom_type) + 1, forcefield) - for atom_type, forcefield in zip(types, forcefields) - ] - ) - - # Modified cross-interactions - if structure.has_NBFIX(): - params = ParameterSet.from_structure(structure) - # Sort keys (maybe they should be sorted in ParmEd) - new_nbfix_types = OrderedDict() - for key in params.nbfix_types.keys(): - sorted_key = tuple(sorted(key)) - if sorted_key in new_nbfix_types: - warn("Sorted key matches an existing key") - if new_nbfix_types[sorted_key]: - warn("nbfixes are not symmetric, overwriting old " "nbfix") - new_nbfix_types[sorted_key] = params.nbfix_types[key] - params.nbfix_types = new_nbfix_types - warn( - "Explicitly writing cross interactions using mixing rule: " - "{}".format(structure.combining_rule) - ) - coeffs = OrderedDict() - for combo in it.combinations_with_replacement(unique_types, 2): - # Attempt to find pair coeffis in nbfixes - if combo in params.nbfix_types: - type1 = unique_types.index(combo[0]) + 1 - type2 = unique_types.index(combo[1]) + 1 - epsilon = params.nbfix_types[combo][0] # kcal OR lj units - rmin = params.nbfix_types[combo][1] # Angstrom OR lj units - sigma = rmin / 2 ** (1 / 6) - coeffs[(type1, type2)] = ( - round(sigma, 8), - round(epsilon, 8), - ) - else: - type1 = unique_types.index(combo[0]) + 1 - type2 = unique_types.index(combo[1]) + 1 - # Might not be necessary to be this explicit - if type1 == type2: - sigma = sigma_dict[type1] - epsilon = epsilon_dict[type1] - else: - if structure.combining_rule == "lorentz": - sigma = (sigma_dict[type1] + sigma_dict[type2]) * 0.5 - elif structure.combining_rule == "geometric": - sigma = (sigma_dict[type1] * sigma_dict[type2]) ** 0.5 - else: - raise ValueError( - "Only lorentz and geometric combining " - "rules are supported" - ) - epsilon = (epsilon_dict[type1] * epsilon_dict[type2]) ** 0.5 - coeffs[(type1, type2)] = ( - round(sigma, 8), - round(epsilon, 8), - ) - if nbfix_in_data_file: - if pair_coeff_label: - data.write("\nPairIJ Coeffs # {}\n".format(pair_coeff_label)) - else: - data.write("\nPairIJ Coeffs # modified lj\n") - - data.write( - "# type1 type2\tepsilon (%s)\tsigma (Angstrom)\n" % eng_unit_str - ) - - for (type1, type2), (sigma, epsilon) in coeffs.items(): - data.write( - "{0} \t{1} \t{2} \t\t{3}\t\t# {4}\t{5}\n".format( - type1, - type2, - epsilon, - sigma, - forcefield_dict[type1], - forcefield_dict[type2], - ) - ) - else: - if pair_coeff_label: - data.write("\nPair Coeffs # {}\n".format(pair_coeff_label)) - else: - data.write("\nPair Coeffs # lj\n") - - for idx, epsilon in sorted(epsilon_dict.items()): - data.write( - "{}\t{:.5f}\t{:.5f}\n".format(idx, epsilon, sigma_dict[idx]) - ) - print("Copy these commands into your input script:\n") - print( - "# type1 type2\tepsilon (%s)\tsigma (Angstrom)\n" % eng_unit_str - ) - for (type1, type2), (sigma, epsilon) in coeffs.items(): - print( - "pair_coeff\t{0} \t{1} \t{2} \t\t{3} \t\t# {4} \t{5}".format( - type1, - type2, - epsilon, - sigma, - forcefield_dict[type1], - forcefield_dict[type2], - ) - ) - - # Pair coefficients - else: - if pair_coeff_label: - data.write("\nPair Coeffs # {}\n".format(pair_coeff_label)) - else: - data.write("\nPair Coeffs # lj\n") - - if unit_style in ["real", "metal"]: - data.write("#\tepsilon (%s)\t\tsigma (Angstrom)\n" % eng_unit_str) - elif unit_style == "lj": - data.write("#\treduced_epsilon \t\treduced_sigma \n") - for idx, epsilon in sorted(epsilon_dict.items()): - data.write( - "{}\t{:.5f}\t\t{:.5f}\t\t# {}\n".format( - idx, epsilon, sigma_dict[idx], forcefield_dict[idx] - ) - ) - - -def _write_bond_information( - structure, data, unique_bond_types, unit_style, eng_unit_str -): - """Write Bond Coeffs section of lammps data file.""" - data.write("\nBond Coeffs # harmonic\n") - if unit_style == "real" or unit_style == "metal": - data.write("#\tk(%s/angstrom^2)\t\treq(angstrom)\n" % eng_unit_str) - elif unit_style == "lj": - data.write("#\treduced_k\t\treduced_req\n") - sorted_bond_types = { - k: v - for k, v in sorted(unique_bond_types.items(), key=lambda item: item[1]) - } - for params, idx in sorted_bond_types.items(): - data.write( - "{}\t{}\t\t{}\t\t# {}\t{}\n".format( - idx, - params[0], - params[1], - params[2][0], - params[2][1], - ) - ) - - -def _write_angle_information( - structure, - data, - unique_angle_types, - use_urey_bradleys, - unit_style, - eng_unit_str, -): - """Write Angle Coeffs section of lammps data file.""" - sorted_angle_types = { - k: v - for k, v in sorted(unique_angle_types.items(), key=lambda item: item[1]) - } - if use_urey_bradleys: - data.write("\nAngle Coeffs # charmm\n") - data.write( - "#\tk(%s/rad^2)\t\ttheteq(deg)\tk(%s/angstrom^2)\treq(angstrom)\n" - % (eng_unit_str, eng_unit_str) - ) - for params, idx in sorted_angle_types.items(): - data.write("{}\t{}\t{:.5f}\t{:.5f}\t{:.5f}\n".format(idx, *params)) - - else: - data.write("\nAngle Coeffs # harmonic\n") - - if unit_style == "lj": - data.write("#\treduced_k\t\ttheteq(deg)\n") - else: - data.write("#\tk(%s/rad^2)\t\ttheteq(deg)\n" % eng_unit_str) - - for params, idx in sorted_angle_types.items(): - data.write( - "{}\t{}\t\t{:.5f}\t# {}\t{}\t{}\n".format( - idx, - params[0], - params[1], - params[3][0], - params[2], - params[3][1], - ) - ) - - -def _write_dihedral_information( - structure, - data, - unique_dihedral_types, - unit_style, - use_rb_torsions, - use_dihedrals, - eng_unit_str, -): - """Write Dihedral Coeffs section of lammps data file.""" - sorted_dihedral_types = { - k: v - for k, v in sorted( - unique_dihedral_types.items(), - key=lambda item: tuple(item[0][i] for i in [-3, -2, -4, -1]), - ) - } - if use_rb_torsions: - data.write("\nDihedral Coeffs # opls\n") - if unit_style == "real" or unit_style == "metal": - data.write( - "#\tf1(%s)\tf2(%s)\tf3(%s)\tf4(%s)\n" - % (eng_unit_str, eng_unit_str, eng_unit_str, eng_unit_str) - ) - elif unit_style == "lj": - data.write("#\tf1\tf2\tf3\tf4 (all lj reduced units)\n") - for idx, (params, _) in enumerate(sorted_dihedral_types.items()): - opls_coeffs = RB_to_OPLS( - params[0], - params[1], - params[2], - params[3], - params[4], - params[5], - error_if_outside_tolerance=False, - ) - data.write( - "{}\t{:.5f}\t{:.5f}\t\t{:.5f}\t\t{:.5f}\t# {}\t{}\t{}\t{}\n".format( - idx + 1, - opls_coeffs[1], - opls_coeffs[2], - opls_coeffs[3], - opls_coeffs[4], - params[8], - params[9], - params[10], - params[11], - ) - ) - elif use_dihedrals: - data.write("\nDihedral Coeffs # charmm\n") - data.write("#k, n, phi, weight\n") - for idx, (params, _) in enumerate(sorted_dihedral_types.items()): - data.write( - "{}\t{:.5f}\t{:d}\t{:d}\t{:.5f}\t# {}\t{}\t{}\t{}\n".format( - idx + 1, - params[0], - params[1], - int(params[2]), - params[3], - params[6], - params[7], - params[8], - params[9], - ) - ) - - -def _write_improper_information( - structure, data, unique_improper_types, unit_style -): - """Write Impropers Coeffs section of lammps data file.""" - sorted_improper_types = { - k: v - for k, v in sorted( - unique_improper_types.items(), key=lambda item: item[1] - ) - } - data.write("\nImproper Coeffs # harmonic\n") - data.write("#k, phi\n") - for params, idx in sorted_improper_types.items(): - data.write( - "{}\t{:.5f}\t{:.5f}\t# {}\t{}\t{}\t{}\n".format( - idx, - params[0], - params[1], - params[2], - params[3], - params[4], - params[5], - ) - ) - - -def _write_imp_dihedral_information( - structure, data, unique_imp_dihedral_types, unit_style -): - """Write Impropers Coeffs section of lammps data file.""" - sorted_imp_dihedral_types = { - k: v - for k, v in sorted( - unique_imp_dihedral_types.items(), - key=lambda item: item[1], - ) - } - data.write("\nImproper Coeffs # cvff\n") - data.write("#K, d, n\n") - for params, idx in sorted_imp_dihedral_types.items(): - data.write( - "{}\t{:.5f}\t{:d}\t{:d}\t# {}\t{}\t{}\t{}\n".format( - idx, - params[0], - params[1], - params[2], - params[5], - params[6], - params[7], - params[8], - ) - ) 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/protobuf.py b/mbuild/formats/protobuf.py deleted file mode 100644 index bc88198b7..000000000 --- a/mbuild/formats/protobuf.py +++ /dev/null @@ -1,195 +0,0 @@ -"""Protocol buffers (Protobuf) format. - -A language-agnostic data serialization format developed by Google -https://developers.google.com/protocol-buffers -""" - -import ele -import numpy as np -from google.protobuf.text_format import Merge, PrintMessage - -from mbuild import Box, Compound -from mbuild.formats import compound_pb2 -from mbuild.utils.decorators import deprecated - -__all__ = ["write_pb2", "read_pb2"] - - -dep_msg = """ -Support for the Protobuf file format will be removed in mBuild 1.0. -""" -print(dep_msg) - - -@deprecated(dep_msg) -def write_pb2(cmpd, filename, binary=True): - """Convert mb.Compound to Protobuf Message file. - - Parameters - ---------- - cmpd : mb.Compound - filename : str - binary: bool, default True - If True, will print a binary file - If False, will print to a text file - """ - cmpd_to_proto = {} - - root_proto = compound_pb2.Compound() - root_proto = _mb_to_proto(cmpd, root_proto) - cmpd_to_proto[cmpd] = root_proto - - for sub_cmpd in cmpd.successors(): - parent_cmpd = sub_cmpd.parent - sub_proto = cmpd_to_proto[parent_cmpd].children.add() - sub_proto = _mb_to_proto(sub_cmpd, sub_proto) - cmpd_to_proto[sub_cmpd] = sub_proto - - _add_proto_bonds(cmpd, root_proto) - - if binary: - with open(filename, "wb") as f: - f.write(root_proto.SerializeToString()) - else: - with open(filename, "w") as f: - PrintMessage(root_proto, f) - - -def read_pb2(filename, binary=True): - """Convert a Protobuf Message file into mb.Compound. - - Parameters - ---------- - filename : str - binary: bool, default True - If True, will print a binary file - If False, will print to a text file - - Returns - ------- - root_compound : mb.Compound - """ - root_proto = compound_pb2.Compound() - if binary: - with open(filename, "rb") as f: - root_proto.ParseFromString(f.read()) - else: - with open(filename, "r") as f: - Merge(f.read(), root_proto) - - proto_to_cmpd = {} - root_compound = _proto_to_mb(root_proto) - proto_to_cmpd[root_proto.id] = root_compound - - for sub_proto, parent_proto in _proto_successors(root_proto): - if parent_proto.id not in proto_to_cmpd: - parent_cmpd = _proto_to_mb(parent_proto) - proto_to_cmpd[parent_proto.id] = parent_cmpd - parent_cmpd = proto_to_cmpd[parent_proto.id] - - if sub_proto.id not in proto_to_cmpd: - sub_cmpd = _proto_to_mb(sub_proto) - proto_to_cmpd[sub_proto.id] = sub_cmpd - sub_cmpd = proto_to_cmpd[sub_proto.id] - - parent_cmpd.add(sub_cmpd) - - _add_mb_bonds(root_proto, root_compound, proto_to_cmpd) - return root_compound - - -def _mb_to_proto(cmpd, proto): - """Given mb.Compound, parse propertes into compound_pb2.Compound.""" - proto.name = cmpd.name - proto.pos.x, proto.pos.y, proto.pos.z = cmpd.pos - proto.charge = cmpd.charge if cmpd.charge else 0.0 - proto.id = id(cmpd) - ( - proto.periodicity.x, - proto.periodicity.y, - proto.periodicity.z, - ) = cmpd.periodicity - if cmpd.element: - proto.element.name = cmpd.element.name - proto.element.symbol = cmpd.element.symbol - proto.element.atomic_number = cmpd.element.atomic_number - proto.element.mass = cmpd.element.mass - - return proto - - -def _add_proto_bonds(cmpd, proto): - """Parse the mb.Compound bonds, add to the proto bonds. - - Parameters - ---------- - cmpd : mb.Compound - proto : compound_pb2.Compound - """ - for b in cmpd.bonds(): - proto_bond = proto.bonds.add() - proto_bond.id1 = id(b[0]) - proto_bond.id2 = id(b[1]) - - -def _proto_successors(proto): - """Recurisve method to look for a compound_pb2's children. - - Parameters - ---------- - proto : compound_pb2 - - Notes - ----- - Base Case: there are no children to the proto, just return - Recursion: First look at proto's children and return these children - (sub_proto) Then make the recursive call to look at all the sub_proto's - successors. This is similar to mb.Compound().successors(). Unlike - mb.Compound(), we need to also keep track of parents in this recursion - """ - if len(proto.children) == 0: - return - for sub_proto in proto.children: - yield (sub_proto, proto) - for sub_sub_proto, parent_proto in _proto_successors(sub_proto): - yield (sub_sub_proto, parent_proto) - - -def _proto_to_mb(proto): - """Given compound_pb2.Compound, create mb.Compound. - - Parameters - ---------- - proto: compound_pb2.Compound() - """ - if proto.element.symbol == "": - elem = None - else: - elem = ele.element_from_symbol(proto.element.symbol) - lengths = [proto.periodicity.x, proto.periodicity.y, proto.periodicity.z] - if np.all(np.array(lengths) != 0): - box = Box(lengths) - else: - box = None - return Compound( - name=proto.name, - pos=[proto.pos.x, proto.pos.y, proto.pos.z], - charge=proto.charge, - box=box, - element=elem, - ) - - -def _add_mb_bonds(proto, cmpd, proto_to_cmpd): - """Parse the compound_pb2.Compound bonds, add to mb.Compound. - - Parameters - ---------- - proto : compound_pb2.Compound - cmpd : mb.Compound - proto_to_cmpd : dict - keys : compound_pb2.Compound.id - value : mb.Compound - """ - for bond in proto.bonds: - cmpd.add_bond([proto_to_cmpd[bond.id1], proto_to_cmpd[bond.id2]]) 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" - ) diff --git a/mbuild/tests/base_test.py b/mbuild/tests/base_test.py index 3437917c5..e58394dd1 100644 --- a/mbuild/tests/base_test.py +++ b/mbuild/tests/base_test.py @@ -26,6 +26,7 @@ def methane(self): @pytest.fixture def benzene(self): benzene = mb.load("C1=CC=CC=C1", smiles=True) + benzene.name = "Benzene" return benzene @@ -121,12 +122,6 @@ def sixpoints(self): molecule.generate_bonds("C", "C", 0.9, 1.1) return molecule - @pytest.fixture - def benzene(self): - compound = mb.load(get_fn("benzene.mol2")) - compound.name = "Benzene" - return compound - @pytest.fixture def rigid_benzene(self): compound = mb.load(get_fn("benzene.mol2")) diff --git a/mbuild/tests/test_cassandramcf.py b/mbuild/tests/test_cassandramcf.py deleted file mode 100644 index f54f4027a..000000000 --- a/mbuild/tests/test_cassandramcf.py +++ /dev/null @@ -1,467 +0,0 @@ -import pytest -from numpy import isclose - -import mbuild as mb -from mbuild.tests.base_test import BaseTest -from mbuild.utils.io import has_foyer - - -@pytest.mark.skipif(not has_foyer, reason="Foyer package not installed") -class TestCassandraMCF(BaseTest): - def test_not_parameterized(self, ethane): - with pytest.raises( - ValueError, match=r"MCF writing not supported without" - ): - ethane.save( - filename="ethane.mcf", - angle_style="harmonic", - dihedral_style="opls", - ) - - def test_invalid_structure(self, ethane): - with pytest.raises(ValueError, match=r"requires parmed structure"): - mb.formats.cassandramcf.write_mcf( - ethane, - "ethane.mcf", - angle_style="harmonic", - dihedral_style="opls", - ) - - def test_invalid_angle_style(self, ethane): - with pytest.raises( - ValueError, match=r"Invalid selection for angle_style" - ): - ethane.save( - filename="ethane-opls.mcf", - forcefield_name="oplsaa", - angle_style="harm", - dihedral_style="opls", - ) - - def test_invalid_dihedral_style(self, ethane): - with pytest.raises( - ValueError, match=r"Invalid selection for dihedral_style" - ): - ethane.save( - filename="ethane-opls.mcf", - forcefield_name="oplsaa", - angle_style="harmonic", - dihedral_style="op", - ) - - def test_dihedral_style_none(self, ethane): - ethane.save( - filename="ethane-opls.mcf", - forcefield_name="oplsaa", - angle_style="harmonic", - dihedral_style="none", - ) - - mcf_data = [] - with open("ethane-opls.mcf") as f: - for line in f: - mcf_data.append(line.strip().split()) - - for idx, line in enumerate(mcf_data): - if len(line) > 1: - if line[1] == "Dihedral_Info": - dihedral_section_start = idx - - assert mcf_data[dihedral_section_start + 2][5] == "none" - - def test_no_dihedrals(self, methane): - methane.save( - filename="methane-opls.mcf", - forcefield_name="oplsaa", - angle_style="harmonic", - dihedral_style="none", - ) - - mcf_data = [] - with open("methane-opls.mcf") as f: - for line in f: - mcf_data.append(line.strip().split()) - - for idx, line in enumerate(mcf_data): - if len(line) > 1: - if line[1] == "Dihedral_Info": - dihedral_section_start = idx - - assert mcf_data[dihedral_section_start + 1][0] == "0" - - def test_no_14(self, methane): - methane.save( - filename="methane-opls.mcf", - forcefield_name="oplsaa", - angle_style="harmonic", - dihedral_style="none", - ) - - mcf_data = [] - with open("methane-opls.mcf") as f: - for line in f: - mcf_data.append(line.strip().split()) - - for idx, line in enumerate(mcf_data): - if len(line) > 1: - if line[1] == "Intra_Scaling": - intrascaling_section_start = idx - - assert isclose(float(mcf_data[intrascaling_section_start + 1][2]), 0.0) - assert isclose(float(mcf_data[intrascaling_section_start + 2][2]), 0.0) - - def test_infer_14(self, ethane): - ethane.save( - filename="ethane-opls.mcf", - forcefield_name="oplsaa", - angle_style="harmonic", - dihedral_style="none", - ) - - mcf_data = [] - with open("ethane-opls.mcf") as f: - for line in f: - mcf_data.append(line.strip().split()) - - for idx, line in enumerate(mcf_data): - if len(line) > 1: - if line[1] == "Intra_Scaling": - intrascaling_section_start = idx - - assert isclose(float(mcf_data[intrascaling_section_start + 1][2]), 0.5) - assert isclose(float(mcf_data[intrascaling_section_start + 2][2]), 0.5) - - def test_unmatched_dihedral_style(self, ethane): - with pytest.raises(ValueError, match=r"but RB torsions found"): - ethane.save( - filename="ethane-opls.mcf", - forcefield_name="oplsaa", - angle_style="harmonic", - dihedral_style="charmm", - ) - - def test_unreasonable_lj14(self, ethane): - import foyer - - oplsaa = foyer.Forcefield(name="oplsaa") - ethane = oplsaa.apply(ethane) - with pytest.raises(ValueError, match=r"Unreasonable value"): - mb.formats.cassandramcf.write_mcf( - ethane, - "ethane.mcf", - angle_style="harmonic", - dihedral_style="opls", - lj14=2.0, - ) - - def test_unreasonable_coul14(self, ethane): - import foyer - - oplsaa = foyer.Forcefield(name="oplsaa") - ethane = oplsaa.apply(ethane) - with pytest.raises(ValueError, match=r"Unreasonable value"): - mb.formats.cassandramcf.write_mcf( - ethane, - "ethane.mcf", - angle_style="harmonic", - dihedral_style="opls", - coul14=-1.0, - ) - - def test_multiple_molecules(self, ethane): - n_ethane = 2 - ethane.name = "Ethane" - filled = mb.fill_box( - ethane, n_compounds=n_ethane, box=[0, 0, 0, 4, 4, 4] - ) - with pytest.raises( - ValueError, match=r"Not all components of the molecule" - ): - filled.save( - filename="box-ethane-opls.mcf", - forcefield_name="oplsaa", - angle_style="harmonic", - dihedral_style="opls", - ) - - def test_save_forcefield(self, ethane): - ethane.save( - filename="ethane-opls.mcf", - forcefield_name="oplsaa", - angle_style="harmonic", - dihedral_style="opls", - ) - - mcf_data = [] - with open("ethane-opls.mcf") as f: - for line in f: - mcf_data.append(line.strip().split()) - - for idx, line in enumerate(mcf_data): - if len(line) > 1: - if line[1] == "Atom_Info": - atom_section_start = idx - elif line[1] == "Bond_Info": - bond_section_start = idx - elif line[1] == "Angle_Info": - angle_section_start = idx - elif line[1] == "Dihedral_Info": - dihedral_section_start = idx - elif line[1] == "Improper_Info": - improper_section_start = idx - elif line[1] == "Fragment_Info": - fragment_section_start = idx - elif line[1] == "Fragment_Connectivity": - fragment_conn_start = idx - - # Check a some atom info - assert mcf_data[atom_section_start + 1][0] == "8" - assert mcf_data[atom_section_start + 2][1] == "opls_135" - assert isclose(float(mcf_data[atom_section_start + 2][3]), 12.011) - assert isclose(float(mcf_data[atom_section_start + 2][4]), -0.18) - assert mcf_data[atom_section_start + 2][5] == "LJ" - assert isclose(float(mcf_data[atom_section_start + 2][6]), 33.21249) - assert isclose(float(mcf_data[atom_section_start + 2][7]), 3.500) - - # Bond info - assert mcf_data[bond_section_start + 1][0] == "7" - passed_test = False - for line in mcf_data[bond_section_start + 2 : angle_section_start - 6]: - a1 = line[1] - a2 = line[2] - if (a1 == "1" and a2 == "2") or (a2 == "1" and a1 == "2"): - assert line[3] == "fixed" - assert isclose(float(line[4]), 1.090) - passed_test = True - break - assert passed_test - - # Angle info - assert mcf_data[angle_section_start + 1][0] == "12" - passed_test = False - for line in mcf_data[ - angle_section_start + 2 : dihedral_section_start - 8 - ]: - if line[2] == "5": - a1 = line[1] - a3 = line[3] - if (a1 == "1" and a3 == "6") or (a3 == "1" and a1 == "6"): - assert line[4] == "harmonic" - assert isclose(float(line[5]), 18870.7) - assert isclose(float(line[6]), 110.7) - passed_test = True - break - assert passed_test - - # Dihedral info - assert mcf_data[dihedral_section_start + 1][0] == "9" - passed_test = False - for line in mcf_data[ - dihedral_section_start + 2 : improper_section_start - 5 - ]: - a1 = line[1] - a2 = line[2] - a3 = line[3] - a4 = line[4] - if (a1 == "2" and a2 == "1" and a3 == "5" and a4 == "6") or ( - a4 == "2" and a3 == "1" and a2 == "5" and a1 == "6" - ): - assert line[5] == "OPLS" - assert isclose(float(line[6]), 0.000) - assert isclose(float(line[7]), 0.000) - assert isclose(float(line[8]), 0.000) - assert isclose(float(line[9]), 0.628) - - assert mcf_data[improper_section_start + 1][0] == "0" - assert mcf_data[fragment_section_start + 1][0] == "2" - - # Check fragment connectivity - assert mcf_data[fragment_conn_start + 1][0] == "1" - assert mcf_data[fragment_conn_start + 2][0] == "1" - assert mcf_data[fragment_conn_start + 2][1] == "1" - assert mcf_data[fragment_conn_start + 2][2] == "2" - - def test_save_ring_forcefield(self, benzene): - benzene.save( - filename="benzene-opls.mcf", - forcefield_name="oplsaa", - angle_style="fixed", - dihedral_style="opls", - ) - - mcf_data = [] - with open("benzene-opls.mcf") as f: - for line in f: - mcf_data.append(line.strip().split()) - - for idx, line in enumerate(mcf_data): - if len(line) > 1: - if line[1] == "Atom_Info": - atom_section_start = idx - elif line[1] == "Bond_Info": - bond_section_start = idx - elif line[1] == "Angle_Info": - angle_section_start = idx - elif line[1] == "Dihedral_Info": - dihedral_section_start = idx - elif line[1] == "Improper_Info": - improper_section_start = idx - elif line[1] == "Fragment_Info": - fragment_section_start = idx - elif line[1] == "Fragment_Connectivity": - fragment_conn_start = idx - - # Check a some atom info - assert mcf_data[atom_section_start + 1][0] == "12" - assert mcf_data[atom_section_start + 2][1] == "opls_145" - assert isclose(float(mcf_data[atom_section_start + 2][3]), 12.011) - assert isclose(float(mcf_data[atom_section_start + 2][4]), -0.115) - assert mcf_data[atom_section_start + 2][5] == "LJ" - assert isclose(float(mcf_data[atom_section_start + 2][6]), 35.22537) - assert isclose(float(mcf_data[atom_section_start + 2][7]), 3.550) - assert mcf_data[atom_section_start + 2][8] == "ring" - - # Bond info - assert mcf_data[bond_section_start + 1][0] == "12" - passed_test = False - for line in mcf_data[bond_section_start + 2 : angle_section_start - 6]: - a1 = line[1] - a2 = line[2] - if (a1 == "1" and a2 == "2") or (a2 == "1" and a1 == "2"): - assert line[3] == "fixed" - assert isclose(float(line[4]), 1.400) - passed_test = True - break - assert passed_test - - # Angle info - assert mcf_data[angle_section_start + 1][0] == "18" - passed_test = False - for line in mcf_data[ - angle_section_start + 2 : dihedral_section_start - 8 - ]: - if line[2] == "2": - a1 = line[1] - a3 = line[3] - if (a1 == "1" and a3 == "3") or (a3 == "1" and a1 == "3"): - assert line[4] == "fixed" - assert isclose(float(line[5]), 120.00) - passed_test = True - break - assert passed_test - - # Dihedral info - assert mcf_data[dihedral_section_start + 1][0] == "24" - passed_test = False - for line in mcf_data[ - dihedral_section_start + 2 : improper_section_start - 5 - ]: - a1 = line[1] - a2 = line[2] - a3 = line[3] - a4 = line[4] - if (a1 == "1" and a2 == "2" and a3 == "3" and a4 == "4") or ( - a4 == "1" and a3 == "2" and a2 == "3" and a1 == "4" - ): - assert line[5] == "OPLS" - assert isclose(float(line[6]), 0.000) - assert isclose(float(line[7]), -0.000) - assert isclose(float(line[8]), 15.167) - assert isclose(float(line[9]), -0.000) - - assert mcf_data[improper_section_start + 1][0] == "0" - assert mcf_data[fragment_section_start + 1][0] == "1" - - # Check fragment connectivity - assert mcf_data[fragment_conn_start + 1][0] == "0" - - def test_shorten_atomname(self, ethane): - import foyer - - from mbuild.formats.cassandramcf import write_mcf - - typed_ethane = foyer.forcefields.load_OPLSAA().apply(ethane) - typed_ethane[0].type = "C_very_very_very_extended" - write_mcf( - typed_ethane, - "ethane-opls.mcf", - angle_style="harmonic", - dihedral_style="opls", - ) - - mcf_data = [] - with open("ethane-opls.mcf") as f: - for line in f: - mcf_data.append(line.strip().split()) - - for idx, line in enumerate(mcf_data): - if len(line) > 1: - if line[1] == "Atom_Info": - atom_section_start = idx - - assert mcf_data[atom_section_start + 2][1] == "y_very_very_extended" - - def test_fused_rings(self): - import foyer - - import mbuild - from mbuild.formats.cassandramcf import write_mcf - - naph = mbuild.load("C1=CC=C2C=CC=CC2=C1", smiles=True) - # Note the atomtyping is wrong -- doesn't matter for test though - naph_ff = foyer.forcefields.load_OPLSAA().apply(naph) - write_mcf( - naph_ff, "naph.mcf", angle_style="harmonic", dihedral_style="opls" - ) - - mcf_data = [] - with open("naph.mcf") as f: - for line in f: - mcf_data.append(line.strip().split()) - - for idx, line in enumerate(mcf_data): - if len(line) > 1: - if line[1] == "Fragment_Info": - frag_section_start = idx - - assert int(mcf_data[frag_section_start + 1][0]) == 1 - assert int(mcf_data[frag_section_start + 2][1]) == 18 - - tmada = mbuild.load("C[C](C)(C)C12CC3CC(C1)CC(C3)C2", smiles=True) - tmada_ff = foyer.forcefields.load_OPLSAA().apply(tmada) - write_mcf( - tmada_ff, "tmada.mcf", angle_style="harmonic", dihedral_style="opls" - ) - - mcf_data = [] - with open("tmada.mcf") as f: - for line in f: - mcf_data.append(line.strip().split()) - - for idx, line in enumerate(mcf_data): - if len(line) > 1: - if line[1] == "Fragment_Info": - frag_section_start = idx - - assert int(mcf_data[frag_section_start + 1][0]) == 5 - assert int(mcf_data[frag_section_start + 2][1]) == 26 - assert int(mcf_data[frag_section_start + 3][1]) == 5 - assert int(mcf_data[frag_section_start + 4][1]) == 5 - assert int(mcf_data[frag_section_start + 5][1]) == 5 - assert int(mcf_data[frag_section_start + 6][1]) == 5 - - def test_infer_14_scaling_zero_eps(self): - import foyer - - import mbuild - from mbuild.formats.cassandramcf import write_mcf - - mol = mbuild.load("CO", smiles=True) - mol_ff = foyer.forcefields.load_OPLSAA().apply(mol) - with pytest.warns(UserWarning): - write_mcf( - mol_ff, - "test.mcf", - angle_style="harmonic", - dihedral_style="opls", - ) diff --git a/mbuild/tests/test_compound.py b/mbuild/tests/test_compound.py index 6d162c605..376f1eb6f 100644 --- a/mbuild/tests/test_compound.py +++ b/mbuild/tests/test_compound.py @@ -16,6 +16,7 @@ get_fn, has_foyer, has_freud, + has_hoomd, has_intermol, has_mdtraj, has_networkx, @@ -89,14 +90,11 @@ def test_load_conversion(self, ethane, h2o): assert np.allclose(test_converted1.xyz, test_converted2.xyz) def test_load_xyz(self): - class MyCompound(Compound): - def __init__(self): - super(MyCompound, self).__init__() - - mb.load(get_fn("ethane.xyz"), compound=self) - - myethane = MyCompound() + myethane = mb.load(get_fn("ethane.xyz")) assert myethane.n_particles == 8 + assert myethane.n_bonds == 0 + assert len(myethane.children) == 1 + assert set([p.name for p in myethane.particles()]) == {"C", "H"} def test_update_from_file(self, ch3): ch3.update_coordinates(get_fn("methyl.pdb")) @@ -373,12 +371,39 @@ def test_load_protein(self): # Chains info is lossed assert len(protein.children) == 393 - def test_save_simple(self, ch3): - extensions = [".xyz", ".pdb", ".mol2", ".json", ".sdf"] - for ext in extensions: - outfile = "methyl_out" + ext - ch3.save(filename=outfile) - assert os.path.exists(outfile) + @pytest.mark.parametrize( + "extension", + [ + ".xyz", + ".pdb", + ".mol2", + ".gro", + ".gsd", + ".json", + ".sdf", + ], + ) + def test_save_simple(self, ch3, extension): + # Can't save gsd files with Windows + if extension == ".gsd" and not has_hoomd: + return True + outfile = "methyl_out" + extension + ch3.save(filename=outfile) + assert os.path.exists(outfile) + + @pytest.mark.parametrize( + "extension", + [".xyz", ".pdb", ".mol2", ".json"], + ) + def test_save_load_simple(self, ch3, extension): + outfile = "methyl_out" + extension + ch3.save(filename=outfile) + ch3_new = mb.load(outfile) + assert ch3_new.n_particles == ch3.n_particles + assert np.allclose(ch3_new.xyz, ch3.xyz) + ch3_names = [p.name for p in ch3.particles()] + ch3_new_names = [p.name for p in ch3_new.particles()] + assert set(ch3_names) == set(ch3_new_names) def test_save_json_loop(self, ethane): ethane.save("ethane.json", show_ports=True) @@ -388,7 +413,7 @@ def test_save_json_loop(self, ethane): assert len(ethane.children) == len(ethane_copy.children) def test_save_box(self, ch3): - extensions = [".mol2", ".pdb", ".hoomdxml", ".gro", ".sdf"] + extensions = [".mol2", ".pdb", ".gro", ".sdf"] box_attributes = ["lengths"] custom_box = Box(lengths=[0.8, 0.8, 0.8], angles=[90, 90, 90]) for ext in extensions: @@ -404,7 +429,7 @@ def test_save_box(self, ch3): assert np.array_equal(pad_attr, custom_attr) def test_save_overwrite(self, ch3): - extensions = [".gsd", ".hoomdxml", ".lammps", ".lmp", ".top", ".gro"] + extensions = [".mol2", ".xyz", ".gro"] for ext in extensions: outfile = "lyhtem" + ext ch3.save(filename=outfile) @@ -412,86 +437,6 @@ def test_save_overwrite(self, ch3): with pytest.raises(IOError): ch3.save(filename=outfile, overwrite=False) - @pytest.mark.skipif(not has_foyer, reason="Foyer is not installed") - def test_save_forcefield(self, methane): - exts = [ - ".gsd", - ".hoomdxml", - ".lammps", - ".lmp", - ".top", - ".gro", - ".mol2", - ".pdb", - ".xyz", - ".sdf", - ] - for ext in exts: - methane.save( - "lythem" + ext, forcefield_name="oplsaa", overwrite=True - ) - - @pytest.mark.skipif(not has_foyer, reason="Foyer is not installed") - def test_save_forcefield_with_file(self, methane): - exts = [ - ".gsd", - ".hoomdxml", - ".lammps", - ".lmp", - ".top", - ".gro", - ".mol2", - ".pdb", - ".xyz", - ".sdf", - ] - for ext in exts: - methane.save( - "lythem" + ext, - forcefield_files=get_fn("methane_oplssaa.xml"), - overwrite=True, - ) - - @pytest.mark.skipif(not has_foyer, reason="Foyer is not installed") - @pytest.mark.parametrize( - "ff_filename,foyer_kwargs", - [ - ("ethane-angle-typo.xml", {"assert_angle_params": False}), - ("ethane-dihedral-typo.xml", {"assert_dihedral_params": False}), - ], - ) - def test_save_missing_topo_params(self, ff_filename, foyer_kwargs): - """Test that the user is notified if not all topology parameters are found.""" - from foyer.tests.utils import get_fn - - ethane = mb.load(get_fn("ethane.mol2")) - with pytest.raises(Exception): - ethane.save("ethane.mol2", forcefield_files=get_fn(ff_filename)) - with pytest.warns(UserWarning): - ethane.save( - "ethane.mol2", - forcefield_files=get_fn(ff_filename), - overwrite=True, - foyer_kwargs=foyer_kwargs, - ) - - @pytest.mark.skipif(not has_foyer, reason="Foyer is not installed") - def test_save_forcefield_with_file_foyer_kwargs(self, methane): - foyer_kwargs = {"assert_improper_params": True} - with pytest.raises(Exception): - methane.save( - "lythem.hoomdxml", - forcefield_files=get_fn("methane_oplssaa.xml"), - overwrite=True, - foyer_kwargs=foyer_kwargs, - ) - methane.save( - "lythem.hoomdxml", - forcefield_files=get_fn("methane_oplssaa.xml"), - overwrite=True, - foyer_kwargs={}, - ) - def test_save_resnames(self, ch3, h2o): system = Compound([ch3, h2o]) system.save("resnames.gro", residues=["CH3", "H2O"]) @@ -507,63 +452,6 @@ def test_save_resnames_single(self, c3, n4): assert struct.residues[0].number == 1 assert struct.residues[1].number == 2 - def test_save_residue_map(self, ethane): - filled = mb.fill_box(ethane, n_compounds=100, box=[0, 0, 0, 4, 4, 4]) - t0 = time.time() - foyer_kwargs = {"use_residue_map": True} - filled.save( - "filled.mol2", - forcefield_name="oplsaa", - residues="Ethane", - foyer_kwargs=foyer_kwargs, - ) - t1 = time.time() - - foyer_kwargs = {"use_residue_map": False} - filled.save( - "filled.mol2", - forcefield_name="oplsaa", - overwrite=True, - residues="Ethane", - foyer_kwargs=foyer_kwargs, - ) - t2 = time.time() - assert (t2 - t1) > (t1 - t0) - - @pytest.mark.skipif(not has_foyer, reason="Foyer is not installed") - def test_save_references(self, methane): - foyer_kwargs = {"references_file": "methane.bib"} - methane.save( - "methyl.mol2", forcefield_name="oplsaa", foyer_kwargs=foyer_kwargs - ) - assert os.path.isfile("methane.bib") - - @pytest.mark.skipif(not has_foyer, reason="Foyer is not installed") - def test_save_combining_rule(self, methane): - combining_rules = ["lorentz", "geometric"] - gmx_rules = {"lorentz": 2, "geometric": 3} - for combining_rule in combining_rules: - if combining_rule == "geometric": - methane.save( - "methane.top", - forcefield_name="oplsaa", - combining_rule=combining_rule, - overwrite=True, - ) - else: - with pytest.warns(UserWarning): - methane.save( - "methane.top", - forcefield_name="oplsaa", - combining_rule=combining_rule, - overwrite=True, - ) - with open("methane.top") as fp: - for i, line in enumerate(fp): - if i == 18: - gmx_rule = int(line.split()[1]) - assert gmx_rule == gmx_rules[combining_rule] - def test_clone_with_box(self, ethane): ethane.box = ethane.get_boundingbox() ethane.periodicity = (True, True, False) @@ -2672,9 +2560,9 @@ def test_elements_from_smiles(self, backend): for particle in mol.particles(): assert particle.element is not None - def test_mins_maxs(self, benzene): - assert np.allclose(benzene.mins, [-0.2267, -0.15422, 0.0]) - assert np.allclose(benzene.maxs, [0.20318, 0.34207, 0.0]) + def test_mins_maxs(self, rigid_benzene): + assert np.allclose(rigid_benzene.mins, [-0.2267, -0.15422, 0.0]) + assert np.allclose(rigid_benzene.maxs, [0.20318, 0.34207, 0.0]) def test_periodicity_raises(self): with pytest.raises(ValueError): diff --git a/mbuild/tests/test_gsd.py b/mbuild/tests/test_gsd.py deleted file mode 100644 index 704a934a1..000000000 --- a/mbuild/tests/test_gsd.py +++ /dev/null @@ -1,320 +0,0 @@ -import numpy as np -import pytest - -import mbuild as mb -from mbuild import Box -from mbuild.tests.base_test import BaseTest -from mbuild.utils.io import has_foyer, has_gsd - - -class TestGSD(BaseTest): - @pytest.mark.skipif(not has_gsd, reason="GSD package not installed") - def test_save(self, ethane): - ethane.save(filename="ethane.gsd") - - @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed") - @pytest.mark.skipif(not has_gsd, reason="GSD package not installed") - def test_save_forcefield(self, ethane): - ethane.save(filename="ethane-opls.gsd", forcefield_name="oplsaa") - - @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed") - @pytest.mark.skipif(not has_gsd, reason="GSD package not installed") - def test_save_box(self, ethane): - box = Box(lengths=[2.0, 2.0, 2.0], angles=[90.0, 90.0, 90.0]) - ethane.save( - filename="ethane-box.gsd", forcefield_name="oplsaa", box=box - ) - - @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed") - def test_save_triclinic_box_(self, ethane): - box = Box(lengths=np.array([2.0, 2.0, 2.0]), angles=[60, 70, 80]) - ethane.save( - filename="triclinic-box.gsd", forcefield_name="oplsaa", box=box - ) - - @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed") - @pytest.mark.skipif(not has_gsd, reason="GSD package not installed") - def test_particles(self, ethane): - from collections import OrderedDict - - import gsd - import gsd.hoomd - - from mbuild.utils.sorting import natural_sort - - ethane.save(filename="ethane.gsd", forcefield_name="oplsaa") - with gsd.hoomd.open("ethane.gsd", mode="r") as f: - frame = f[0] - - assert frame.configuration.step == 0 - assert frame.configuration.dimensions == 3 - - mass_dict = {"C": 12.011, "H": 1.008} - masses = frame.particles.mass.astype(float) - for mass, particle in zip(masses, ethane.particles()): - assert round(mass, 3) == mass_dict[particle.name] - - n_particles = frame.particles.N - assert n_particles == ethane.n_particles - - positions = frame.particles.position.astype(float) - shift = positions[0] - (ethane[0].pos * 10) - shifted_xyz = (ethane.xyz * 10) + shift - assert np.array_equal( - np.round(positions, decimals=4), np.round(shifted_xyz, decimals=4) - ) - - opls_type_dict = OrderedDict([("C", "opls_135"), ("H", "opls_140")]) - types_from_gsd = frame.particles.types - typeids_from_gsd = frame.particles.typeid.astype(int) - expected_types = [ - opls_type_dict[particle.name] for particle in ethane.particles() - ] - unique_types = list(set(expected_types)) - unique_types.sort(key=natural_sort) - expected_typeids = [ - unique_types.index(atype) for atype in expected_types - ] - assert np.array_equal(types_from_gsd, unique_types) - assert np.array_equal(typeids_from_gsd, expected_typeids) - - @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed") - @pytest.mark.skipif(not has_gsd, reason="GSD package not installed") - def test_box(self, ethane): - import gsd - import gsd.hoomd - - lengths = [2.0, 3.0, 4.0] - ethane.box = Box(lengths=[2.0, 3.0, 4.0]) - - ethane.save(filename="ethane.gsd", forcefield_name="oplsaa") - with gsd.hoomd.open("ethane.gsd", mode="r") as f: - frame = f[0] - - box_from_gsd = frame.configuration.box.astype(float) - (lx, ly, lz) = ethane.box.lengths - lx *= 10 - ly *= 10 - lz *= 10 - assert np.array_equal(box_from_gsd[:3], [lx, ly, lz]) - assert not np.any(box_from_gsd[3:]) - - ethane.periodicity = (True, True, True) - ethane.save(filename="ethane-periodicity.gsd", forcefield_name="oplsaa") - with gsd.hoomd.open("ethane-periodicity.gsd", mode="r") as f: - frame = f[0] - box_from_gsd_periodic = frame.configuration.box.astype(float) - assert np.array_equal(box_from_gsd, box_from_gsd_periodic) - - box = Box(lengths=np.array([2.0, 2.0, 2.0]), angles=[92, 104, 119]) - # check that providing a box to save overwrites compound.box - ethane.save( - filename="triclinic-box.gsd", forcefield_name="oplsaa", box=box - ) - with gsd.hoomd.open("triclinic-box.gsd", mode="r") as f: - frame = f[0] - lx, ly, lz, xy, xz, yz = frame.configuration.box - - a = lx - b = np.sqrt(ly**2 + xy**2) - c = np.sqrt(lz**2 + xz**2 + yz**2) - - assert np.isclose(np.cos(np.radians(92)), (xy * xz + ly * yz) / (b * c)) - assert np.isclose(np.cos(np.radians(104)), xz / c) - assert np.isclose(np.cos(np.radians(119)), xy / b) - - @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed") - @pytest.mark.skipif(not has_gsd, reason="GSD package not installed") - def test_rigid(self, benzene): - import gsd - import gsd.hoomd - - benzene.label_rigid_bodies(rigid_particles="C") - benzene.save(filename="benzene.gsd", forcefield_name="oplsaa") - with gsd.hoomd.open("benzene.gsd", mode="r") as f: - frame = f[0] - - rigid_bodies = frame.particles.body - expected_bodies = [ - -1 if p.rigid_id is None else p.rigid_id - for p in benzene.particles() - ] - for gsd_body, expected_body in zip(rigid_bodies, expected_bodies): - assert gsd_body == expected_body - - @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed") - @pytest.mark.skipif(not has_gsd, reason="GSD package not installed") - def test_bonded(self, ethane): - import gsd - import gsd.hoomd - from foyer import Forcefield - - ethane.save(filename="ethane.gsd", forcefield_name="oplsaa") - with gsd.hoomd.open("ethane.gsd", mode="r") as f: - frame = f[0] - - structure = ethane.to_parmed() - forcefield = Forcefield(name="oplsaa") - structure = forcefield.apply(structure) - - # Bonds - n_bonds = frame.bonds.N - assert n_bonds == len(structure.bonds) - - expected_unique_bond_types = ["opls_135-opls_135", "opls_135-opls_140"] - bond_types = frame.bonds.types - assert np.array_equal(bond_types, expected_unique_bond_types) - - bond_typeids = frame.bonds.typeid - bond_atoms = frame.bonds.group - expected_bond_atoms = [ - [bond.atom1.idx, bond.atom2.idx] for bond in structure.bonds - ] - assert np.array_equal(bond_atoms, expected_bond_atoms) - - bond_type_dict = {("C", "C"): 0, ("C", "H"): 1, ("H", "C"): 1} - expected_bond_typeids = [] - for bond in structure.bonds: - expected_bond_typeids.append( - bond_type_dict[(bond.atom1.name, bond.atom2.name)] - ) - assert np.array_equal(bond_typeids, expected_bond_typeids) - - # Angles - n_angles = frame.angles.N - assert n_angles == len(structure.angles) - - expected_unique_angle_types = [ - "opls_135-opls_135-opls_140", - "opls_140-opls_135-opls_140", - ] - angle_types = frame.angles.types - assert np.array_equal(angle_types, expected_unique_angle_types) - - angle_typeids = frame.angles.typeid - angle_atoms = frame.angles.group - expected_angle_atoms = [ - [angle.atom1.idx, angle.atom2.idx, angle.atom3.idx] - for angle in structure.angles - ] - assert np.array_equal(angle_atoms, expected_angle_atoms) - - angle_type_dict = { - ("C", "C", "H"): 0, - ("H", "C", "C"): 0, - ("H", "C", "H"): 1, - } - expected_angle_typeids = [] - for angle in structure.angles: - expected_angle_typeids.append( - angle_type_dict[ - (angle.atom1.name, angle.atom2.name, angle.atom3.name) - ] - ) - assert np.array_equal(angle_typeids, expected_angle_typeids) - - # Dihedrals - n_dihedrals = frame.dihedrals.N - assert n_dihedrals == len(structure.rb_torsions) - - expected_unique_dihedral_types = ["opls_140-opls_135-opls_135-opls_140"] - dihedral_types = frame.dihedrals.types - assert np.array_equal(dihedral_types, expected_unique_dihedral_types) - - dihedral_typeids = frame.dihedrals.typeid - dihedral_atoms = frame.dihedrals.group - expected_dihedral_atoms = [] - for dihedral in structure.rb_torsions: - expected_dihedral_atoms.append( - [ - dihedral.atom1.idx, - dihedral.atom2.idx, - dihedral.atom3.idx, - dihedral.atom4.idx, - ] - ) - assert np.array_equal(dihedral_atoms, expected_dihedral_atoms) - assert np.array_equal(dihedral_typeids, np.zeros(n_dihedrals)) - - @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed") - @pytest.mark.skipif(not has_gsd, reason="GSD package not installed") - def test_pairs(self, benzene): - import gsd - import gsd.hoomd - from foyer import Forcefield - - benzene.save(filename="benzene.gsd", forcefield_name="oplsaa") - with gsd.hoomd.open("benzene.gsd", mode="r") as f: - frame = f[0] - - structure = benzene.to_parmed() - forcefield = Forcefield(name="oplsaa") - structure = forcefield.apply(structure) - - # Pairs - assert len(frame.pairs.types) == 3 - assert frame.pairs.N == 21 - - @pytest.mark.skipif(not has_gsd, reason="GSD package not installed") - def test_units(self, ethane): - import gsd - import gsd.hoomd - - ref_distance = 3.5 - ref_energy = 0.066 - ref_mass = 12.011 - - box = Box(lengths=[2.0, 3.0, 4.0], angles=[90.0, 90.0, 90.0]) - ethane.save( - filename="ethane.gsd", - forcefield_name="oplsaa", - ref_distance=ref_distance, - ref_energy=ref_energy, - ref_mass=ref_mass, - box=box, - ) - with gsd.hoomd.open("ethane.gsd", mode="r") as f: - frame = f[0] - - box_from_gsd = frame.configuration.box.astype(float) - assert np.array_equal( - np.round(box_from_gsd[:3], decimals=5), - np.round(np.asarray(box.lengths) * 10 / ref_distance, 5), - ) - - mass_dict = {"C": 12.011, "H": 1.008} - masses = frame.particles.mass.astype(float) - for mass, p in zip(masses, ethane.particles()): - assert round(mass, 3) == round(mass_dict[p.name] / ref_mass, 3) - - charge_dict = {"C": -0.18, "H": 0.06} - charges = frame.particles.charge.astype(float) - e0 = 2.396452e-4 - charge_factor = (4.0 * np.pi * e0 * ref_distance * ref_energy) ** 0.5 - for charge, particle in zip(charges, ethane.particles()): - reduced_charge = charge_dict[particle.name] / charge_factor - assert round(charge, 3) == round(reduced_charge, 3) - - positions = frame.particles.position.astype(float) - shift = positions[0] - (ethane[0].pos * 10 / ref_distance) - shifted_xyz = (ethane.xyz * 10 / ref_distance) + shift - assert np.array_equal( - np.round(positions, decimals=4), np.round(shifted_xyz, decimals=4) - ) - - @pytest.mark.skipif(not has_gsd, reason="GSD package not installed") - def test_box_dimensions(self, benzene): - import gsd - import gsd.hoomd - - n_benzenes = 10 - filled = mb.fill_box( - benzene, n_compounds=n_benzenes, box=[0, 0, 0, 4, 4, 4] - ) - filled.save(filename="benzene.gsd") - with gsd.hoomd.open("benzene.gsd", mode="r") as f: - frame = f[0] - positions = frame.particles.position.astype(float) - for coords in positions: - assert coords.max() < 20 - assert coords.min() > -20 diff --git a/mbuild/tests/test_hoomd.py b/mbuild/tests/test_hoomd.py deleted file mode 100644 index 5dc243fc9..000000000 --- a/mbuild/tests/test_hoomd.py +++ /dev/null @@ -1,445 +0,0 @@ -import xml.etree.ElementTree - -import numpy as np -import packaging.version -import pytest - -import mbuild as mb -from mbuild.tests.base_test import BaseTest -from mbuild.utils.io import get_fn, has_foyer, has_gsd, has_hoomd, import_ - -if has_hoomd: - import hoomd - - if "version" in dir(hoomd): - hoomd_version = packaging.version.parse(hoomd.version.version) - else: - hoomd_version = packaging.version.parse(hoomd.__version__) - - -@pytest.mark.skipif(not has_hoomd, reason="HOOMD is not installed") -class TestHoomdAny(BaseTest): - def test_empty_initial_snapshot(self): - import hoomd - - from mbuild.formats.hoomd_snapshot import to_hoomdsnapshot - - part = mb.Compound(name="Ar") - box = mb.Box(lengths=[5, 5, 5], angles=[90, 90, 90]) - system = mb.fill_box(part, n_compounds=10, box=box) - - if hoomd_version.major == 2: - hoomd.context.initialize("") - init_snap = hoomd.data.make_snapshot( - N=0, box=hoomd.data.boxdim(L=10) - ) - else: - init_snap = hoomd.Snapshot() - init_snap.configuration.box = hoomd.Box.cube(L=10) - - with pytest.raises(RuntimeError): - snap, _ = to_hoomdsnapshot(system, hoomd_snapshot=init_snap) - - def test_compound_from_snapshot(self, ethane): - from mbuild.formats.hoomd_snapshot import ( - from_snapshot, - to_hoomdsnapshot, - ) - - lengths = [5, 5, 5] - filled = mb.fill_box(ethane, n_compounds=5, box=mb.Box(lengths)) - snap, _ = to_hoomdsnapshot(filled) - new_filled = from_snapshot(snap, scale=0.1) - - assert filled.n_bonds == new_filled.n_bonds - assert filled.n_particles == new_filled.n_particles - - assert np.array_equal(filled.box.angles, new_filled.box.angles) - assert np.array_equal(filled.box.lengths, new_filled.box.lengths) - - for i in range(filled.n_particles): - assert np.allclose(filled[i].pos, new_filled[i].pos) - - @pytest.mark.skipif(not has_gsd, reason="gsd is not installed") - def test_compound_from_gsdsnapshot(self, ethane): - import gsd.hoomd - - from mbuild.formats.hoomd_snapshot import ( - from_snapshot, - to_hoomdsnapshot, - ) - - lengths = [5, 5, 5] - filled = mb.fill_box(ethane, n_compounds=5, box=mb.Box(lengths)) - snap, _ = to_hoomdsnapshot(filled) - - # copy attributes from the snapshot to a gsd snapshot - gsd_snap = gsd.hoomd.Frame() - gsd_snap.particles.N = snap.particles.N - gsd_snap.particles.types = snap.particles.types - gsd_snap.particles.typeid = snap.particles.typeid - gsd_snap.particles.position = snap.particles.position - if hoomd_version.major == 2: - gsd_snap.configuration.box = np.array( - [ - snap.box.Lx, - snap.box.Ly, - snap.box.Lz, - snap.box.xy, - snap.box.xy, - snap.box.yz, - ] - ) - else: - gsd_snap.configuration.box = snap.configuration.box - - gsd_snap.bonds.N = snap.bonds.N - gsd_snap.bonds.group = snap.bonds.group - gsd_snap.particles.charge = snap.particles.charge - gsd_snap.validate() - - new_filled = from_snapshot(gsd_snap, scale=0.1) - - assert filled.n_bonds == new_filled.n_bonds - assert filled.n_particles == new_filled.n_particles - - assert np.array_equal(filled.box.angles, new_filled.box.angles) - assert np.array_equal(filled.box.lengths, new_filled.box.lengths) - - for i in range(filled.n_particles): - assert np.allclose(filled[i].pos, new_filled[i].pos) - - def test_compound_to_snapshot(self, ethane): - from mbuild.formats.hoomd_snapshot import to_hoomdsnapshot - - snap, _ = to_hoomdsnapshot(ethane) - - assert snap.particles.N == 8 - assert snap.bonds.N == 7 - assert snap.angles.N == 0 - - def test_particles_to_snapshot(self): - from mbuild.formats.hoomd_snapshot import to_hoomdsnapshot - - part = mb.Compound(name="Ar") - box = mb.Box(lengths=[5, 5, 5], angles=[90, 90, 90]) - system = mb.fill_box(part, n_compounds=10, box=box) - snap, _ = to_hoomdsnapshot(system) - - assert snap.particles.N == 10 - assert snap.bonds.N == 0 - assert snap.angles.N == 0 - - def test_snapshot_from_initial(self): - import hoomd - - from mbuild.formats.hoomd_snapshot import to_hoomdsnapshot - - part = mb.Compound(name="Ar") - box = mb.Box(lengths=[5, 5, 5], angles=[90, 90, 90]) - system = mb.fill_box(part, n_compounds=10, box=box) - if hoomd_version.major == 2: - hoomd.context.initialize("") - init_snap = hoomd.data.make_snapshot( - N=10, box=hoomd.data.boxdim(L=10) - ) - else: - init_snap = hoomd.Snapshot() - init_snap.particles.N = 10 - init_snap.configuration.box = hoomd.Box.cube(L=10) - - snap, _ = to_hoomdsnapshot(system, hoomd_snapshot=init_snap) - - assert snap.particles.N == 20 - assert snap.bonds.N == 0 - assert snap.angles.N == 0 - if hoomd_version.major == 2: - assert (snap.box.Lx, snap.box.Ly, snap.box.Lz) == (50, 50, 50) - assert (snap.box.xy, snap.box.xz, snap.box.yz) == (0, 0, 0) - else: - np.testing.assert_allclose( - snap.configuration.box, [50, 50, 50, 0, 0, 0] - ) - - def test_bad_input_to_snapshot(self): - from mbuild.formats.hoomd_snapshot import to_hoomdsnapshot - - with pytest.raises(ValueError): - to_hoomdsnapshot("fake_object") - - def test_non_param_struc_to_snapshot(self, ethane): - from mbuild.formats.hoomd_snapshot import to_hoomdsnapshot - - structure = ethane.to_parmed() - snap, _ = to_hoomdsnapshot(structure) - - assert snap.particles.N == 8 - assert snap.bonds.N == 7 - assert snap.angles.N == 0 - - @pytest.mark.skipif(not has_foyer, reason="Foyer is not installed") - def test_param_structure_to_snapshot(self, ethane): - from foyer.forcefield import Forcefield - - from mbuild.formats.hoomd_snapshot import to_hoomdsnapshot - - ff = Forcefield(name="oplsaa") - structure = ff.apply(ethane) - snap, _ = to_hoomdsnapshot(structure) - - assert snap.particles.N == 8 - assert snap.bonds.N == 7 - assert snap.angles.N == 12 - assert snap.dihedrals.N == 9 - assert snap.pairs.N == 9 - - -@pytest.mark.skipif( - not has_hoomd or hoomd_version.major != 2, - reason="HOOMD is not installed or is the wrong version", -) -class TestHoomdSimulation(BaseTest): - def test_bad_input_to_hoomdsimulation(self): - from mbuild.formats.hoomd_simulation import create_hoomd_simulation - - with pytest.raises(ValueError): - create_hoomd_simulation("fake_object", 2) - - def test_compound_to_hoomdsimulation(self, ethane): - from mbuild.formats.hoomd_simulation import create_hoomd_simulation - - with pytest.raises(ValueError): - create_hoomd_simulation(ethane, 2.5) - - @pytest.mark.skipif(not has_foyer, reason="Foyer is not installed") - def test_structure_to_hoomdsimulation(self, ethane): - import hoomd - from foyer.forcefield import Forcefield - - from mbuild.formats.hoomd_simulation import create_hoomd_simulation - - ff = Forcefield(name="oplsaa") - structure = ff.apply(ethane) - sim = hoomd.context.SimulationContext() - with sim: - create_hoomd_simulation(structure, 2.5) - - sim_forces = hoomd.context.current.forces - pair_force = import_("hoomd.md.pair") - charge_force = import_("hoomd.md.charge") - special_pair_force = import_("hoomd.md.special_pair") - bond_force = import_("hoomd.md.bond") - angle_force = import_("hoomd.md.angle") - dihedral_force = import_("hoomd.md.dihedral") - - assert isinstance(sim_forces[0], pair_force.lj) - assert isinstance(sim_forces[1], charge_force.pppm) - assert isinstance(sim_forces[2], pair_force.ewald) - assert isinstance(sim_forces[3], special_pair_force.lj) - assert isinstance(sim_forces[4], special_pair_force.coulomb) - assert isinstance(sim_forces[5], bond_force.harmonic) - assert isinstance(sim_forces[6], angle_force.harmonic) - assert isinstance(sim_forces[7], dihedral_force.opls) - - @pytest.mark.skipif(not has_foyer, reason="Foyer is not installed") - def test_lj_to_hoomdsimulation(self): - import hoomd - from foyer.forcefield import Forcefield - - from mbuild.formats.hoomd_simulation import create_hoomd_simulation - - box = mb.Compound() - box.add(mb.Compound(name="Ar", pos=[1, 1, 1])) - box.add(mb.Compound(name="Ar", pos=[1, 1, 1])) - ff = Forcefield(forcefield_files=get_fn("lj.xml")) - structure = ff.apply(box) - structure.box = [10, 10, 10, 90, 90, 90] - sim = hoomd.context.SimulationContext() - with sim: - create_hoomd_simulation(structure, 2.5) - sim_forces = hoomd.context.current.forces - pair_force = import_("hoomd.md.pair") - - assert isinstance(sim_forces[0], pair_force.lj) - - @pytest.mark.skipif(not has_foyer, reason="Foyer is not installed") - @pytest.mark.skipif(not has_gsd, reason="GSD is not installed") - def test_hoomdsimulation_restart(self): - import gsd.hoomd - import hoomd - from foyer.forcefield import Forcefield - - from mbuild.formats.hoomd_simulation import create_hoomd_simulation - - box = mb.Compound() - box.add(mb.Compound(name="Ar", pos=[1, 1, 1])) - box.add(mb.Compound(name="Ar", pos=[1, 1, 1])) - ff = Forcefield(forcefield_files=get_fn("lj.xml")) - structure = ff.apply(box) - structure.box = [10, 10, 10, 90, 90, 90] - sim = hoomd.context.SimulationContext() - with sim: - hoomd_obj, ref_vals = create_hoomd_simulation( - structure, 2.5, restart=get_fn("restart.gsd") - ) - sim_forces = hoomd.context.current.forces - pair_force = import_("hoomd.md.pair") - - assert isinstance(sim_forces[0], pair_force.lj) - - snap = hoomd_obj[0] - with gsd.hoomd.open(get_fn("restart.gsd")) as f: - rsnap = f[0] - assert np.array_equal(snap.particles.position, rsnap.particles.position) - - def test_hoomdsimulation_nlist(self, ethane): - hoomd_simulation = import_("mbuild.formats.hoomd_simulation") - hoomd = import_("hoomd") - hoomd.md = import_("hoomd.md") - - with pytest.raises(ValueError): - hoomd_simulation.create_hoomd_simulation( - ethane, 2.5, nlist=hoomd.md.nlist.tree - ) - - -@pytest.mark.skipif( - not has_hoomd or hoomd_version.major < 3, - reason="HOOMD is not installed or is the wrong version", -) -class TestHoomdForcefield(BaseTest): - def test_bad_input_to_hoomd_forcefield(self): - from mbuild.formats.hoomd_forcefield import create_hoomd_forcefield - - with pytest.raises(ValueError): - create_hoomd_forcefield("fake_object", 2.5) - - def test_compound_to_hoomd_forcefield(self, ethane): - from mbuild.formats.hoomd_forcefield import create_hoomd_forcefield - - with pytest.raises(ValueError): - create_hoomd_forcefield(ethane, 2.5) - - @pytest.mark.skipif(not has_foyer, reason="Foyer is not installed") - def test_structure_to_hoomd_forcefield(self, ethane): - import hoomd - from foyer.forcefield import Forcefield - - from mbuild.formats.hoomd_forcefield import create_hoomd_forcefield - - ff = Forcefield(name="oplsaa") - structure = ff.apply(ethane) - - snapshot, forces, ref_values = create_hoomd_forcefield(structure, 2.5) - - assert isinstance(forces[0], hoomd.md.pair.LJ) - assert isinstance(forces[1], hoomd.md.pair.Ewald) - assert isinstance(forces[2], hoomd.md.long_range.pppm.Coulomb) - assert isinstance(forces[3], hoomd.md.special_pair.LJ) - assert isinstance(forces[4], hoomd.md.special_pair.Coulomb) - assert isinstance(forces[5], hoomd.md.bond.Harmonic) - assert isinstance(forces[6], hoomd.md.angle.Harmonic) - assert isinstance(forces[7], hoomd.md.dihedral.OPLS) - - @pytest.mark.skipif(not has_foyer, reason="Foyer is not installed") - def test_lj_to_hoomd_forcefield(self): - import hoomd - from foyer.forcefield import Forcefield - - from mbuild.formats.hoomd_forcefield import create_hoomd_forcefield - - box = mb.Compound() - box.add(mb.Compound(name="Ar", pos=[1, 1, 1])) - box.add(mb.Compound(name="Ar", pos=[1, 1, 1])) - ff = Forcefield(forcefield_files=get_fn("lj.xml")) - structure = ff.apply(box) - structure.box = [10, 10, 10, 90, 90, 90] - - snapshot, forces, ref_values = create_hoomd_forcefield(structure, 2.5) - - assert isinstance(forces[0], hoomd.md.pair.LJ) - - -class TestHoomdXML(BaseTest): - def test_save(self, ethane): - ethane.save(filename="ethane.hoomdxml") - - @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed") - def test_save_forcefield(self, ethane): - ethane.save(filename="ethane-opls.hoomdxml", forcefield_name="oplsaa") - - def test_save_box(self, ethane): - box = mb.box.Box(lengths=[2.0, 2.0, 2.0], angles=[90, 90, 90]) - ethane.save(filename="ethane-box.hoomdxml", box=box) - - def test_save_triclinic_box_(self, ethane): - box = mb.Box(lengths=np.array([2.0, 2.0, 2.0]), angles=[60, 70, 80]) - ethane.save(filename="triclinic-box.hoomdxml", box=box) - - def test_rigid(self, benzene): - n_benzenes = 10 - benzene.name = "Benzene" - filled = mb.fill_box( - benzene, n_compounds=n_benzenes, box=[0, 0, 0, 4, 4, 4] - ) - filled.label_rigid_bodies( - discrete_bodies="Benzene", rigid_particles="C" - ) - filled.save(filename="benzene.hoomdxml") - - xml_file = xml.etree.ElementTree.parse("benzene.hoomdxml").getroot() - body_text = xml_file[0].find("body").text - rigid_bodies = [int(body) for body in body_text.split("\n") if body] - for body_id in range(10): - assert rigid_bodies.count(body_id) == 6 - assert rigid_bodies.count(-1) == n_benzenes * 6 - - @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed") - def test_number_in_each_section(self, box_of_benzenes): - box_of_benzenes.save( - filename="benzene.hoomdxml", forcefield_name="oplsaa" - ) - xml_file = xml.etree.ElementTree.parse("benzene.hoomdxml").getroot() - for attribute in ["position", "type", "mass", "charge"]: - body_text = xml_file[0].find(attribute).text - list_of_things = [x for x in body_text.split("\n") if x] - assert len(list_of_things) == 12 * 10 - for attribute, number in [ - ("bond", 12), - ("angle", 18), - ("dihedral", 24), - ]: - body_text = xml_file[0].find(attribute).text - list_of_things = [x for x in body_text.split("\n") if x] - assert len(list_of_things) == number * 10 - - def test_box_dimensions(self, benzene): - n_benzenes = 10 - filled = mb.fill_box( - benzene, n_compounds=n_benzenes, box=[0, 0, 0, 4, 4, 4] - ) - filled.save(filename="benzene.hoomdxml") - for atom in mb.load("benzene.hoomdxml"): - assert atom.pos.max() < 20 - assert atom.pos.min() > -20 - - @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed") - def test_auto_scale_forcefield(self, ethane): - ethane.save( - filename="ethane-opls.hoomdxml", - forcefield_name="oplsaa", - auto_scale=True, - ) - xml_file = xml.etree.ElementTree.parse("ethane-opls.hoomdxml").getroot() - masses = xml_file[0].find("mass").text.splitlines() - # We use 1 and 5 since the first element of masses is empty - assert masses[1] == "1.0" - assert masses[5] == "1.0" - pair_coeffs = [ - _.split("\t") - for _ in xml_file[0].find("pair_coeffs").text.splitlines() - ] - # The first element is empty, the next element should be - # ['opls_135', '1.0000', '1.0000'] - assert pair_coeffs[1][1] == "1.0000" - assert pair_coeffs[1][2] == "1.0000" diff --git a/mbuild/tests/test_lammpsdata.py b/mbuild/tests/test_lammpsdata.py deleted file mode 100644 index 3b4711dcc..000000000 --- a/mbuild/tests/test_lammpsdata.py +++ /dev/null @@ -1,742 +0,0 @@ -from pathlib import Path - -import numpy as np -import pytest -from ele.element import element_from_symbol -from pytest import FixtureRequest -from scipy.constants import epsilon_0 - -import mbuild as mb -from mbuild.formats.lammpsdata import write_lammpsdata -from mbuild.tests.base_test import BaseTest -from mbuild.utils.io import get_fn, has_foyer - -KCAL_TO_KJ = 4.184 -ANG_TO_NM = 0.10 -ELEM_TO_COUL = 1.602176634e-19 - - -@pytest.mark.skipif(not has_foyer, reason="Foyer package not installed") -class TestLammpsData(BaseTest): - def test_save(self, ethane): - ethane.save(filename="ethane.lammps") - - @pytest.fixture(scope="session") - def lj_save(self, tmpdir_factory, sigma=None, epsilon=None, mass=None): - def _create_lammps( - ethane, tmp, sigma=sigma, epsilon=epsilon, mass=mass - ): - from foyer import Forcefield - - OPLSAA = Forcefield(name="oplsaa") - structure = OPLSAA.apply(ethane) - fn = tmpdir_factory.mktemp("data").join("lj.lammps") - write_lammpsdata( - filename=str(fn), - structure=structure, - unit_style="lj", - sigma_conversion_factor=sigma, - epsilon_conversion_factor=epsilon, - mass_conversion_factor=mass, - ) - return str(fn) - - return _create_lammps - - @pytest.fixture(scope="session") - def real_save(self, tmpdir_factory): - def _create_lammps(ethane, tmp): - from foyer import Forcefield - - OPLSAA = Forcefield(name="oplsaa") - structure = OPLSAA.apply(ethane) - fn = tmpdir_factory.mktemp("data").join("lj.lammps") - write_lammpsdata( - filename=str(fn), structure=structure, unit_style="real" - ) - return str(fn) - - return _create_lammps - - @pytest.mark.parametrize("unit_style", ["real", "lj"]) - def test_save_forcefield(self, ethane, unit_style): - ethane.save( - filename="ethane-opls.lammps", - forcefield_name="oplsaa", - unit_style=unit_style, - ) - - @pytest.mark.parametrize("pair_coeff_label", [None, "lj", "lj/coul/cut"]) - def test_save_pair_coeff_label(self, ethane, pair_coeff_label): - ethane.save( - filename="ethane-opls.lammps", - forcefield_name="oplsaa", - unit_style="real", - pair_coeff_label=pair_coeff_label, - ) - with open("ethane-opls.lammps", "r") as fp: - for line in fp: - if "Pair Coeffs" in line: - if pair_coeff_label is None: - assert "# lj" in line - else: - assert "# {}".format(pair_coeff_label) in line - - @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed") - def test_save_charmm(self): - from foyer import Forcefield - - cmpd = mb.load(get_fn("charmm_dihedral.mol2")) - for i in cmpd.particles(): - i.name = "_{}".format(i.name) - structure = cmpd.to_parmed( - box=cmpd.get_boundingbox(), - residues=set([p.parent.name for p in cmpd.particles()]), - ) - ff = Forcefield(forcefield_files=[get_fn("charmm_truncated.xml")]) - structure = ff.apply(structure, assert_dihedral_params=False) - write_lammpsdata(structure, "charmm_dihedral.lammps") - out_lammps = open("charmm_dihedral.lammps", "r").readlines() - found_angles = False - found_dihedrals = False - for i, line in enumerate(out_lammps): - if "Angle Coeffs" in line: - assert "# charmm" in line - assert ( - "#\tk(kcal/mol/rad^2)\t\ttheteq(deg)\tk(kcal/mol/angstrom^2)\treq(angstrom)\n" - in out_lammps[i + 1] - ) - assert len(out_lammps[i + 2].split("#")[0].split()) == 5 - found_angles = True - elif "Dihedral Coeffs" in line: - assert "# charmm" in line - assert "#k, n, phi, weight" in out_lammps[i + 1] - assert len(out_lammps[i + 2].split("#")[0].split()) == 5 - found_dihedrals = True - else: - pass - assert found_angles - assert found_dihedrals - - @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed") - def test_singleterm_charmm(self): - from foyer import Forcefield - - from mbuild.formats.lammpsdata import write_lammpsdata - - cmpd = mb.load(get_fn("charmm_dihedral.mol2")) - for i in cmpd.particles(): - i.name = "_{}".format(i.name) - structure = cmpd.to_parmed( - box=cmpd.get_boundingbox(), - residues=set([p.parent.name for p in cmpd.particles()]), - ) - ff = Forcefield( - forcefield_files=[get_fn("charmm_truncated_singleterm.xml")] - ) - structure = ff.apply(structure, assert_dihedral_params=False) - write_lammpsdata(structure, "charmm_dihedral_singleterm.lammps") - out_lammps = open("charmm_dihedral_singleterm.lammps", "r").readlines() - found_dihedrals = False - for i, line in enumerate(out_lammps): - if "Dihedral Coeffs" in line: - assert "# charmm" in line - assert "#k, n, phi, weight" in out_lammps[i + 1] - assert len(out_lammps[i + 2].split("#")[0].split()) == 5 - assert float( - out_lammps[i + 2].split("#")[0].split()[4] - ) == float("1.0") - found_dihedrals = True - else: - pass - assert found_dihedrals - - @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed") - def test_charmm_improper(self): - from foyer import Forcefield - - import mbuild as mb - from mbuild.formats.lammpsdata import write_lammpsdata - - system = mb.Compound() - first = mb.Particle(name="_CTL2", pos=[-1, 0, 0]) - second = mb.Particle(name="_CL", pos=[0, 0, 0]) - third = mb.Particle(name="_OBL", pos=[1, 0, 0]) - fourth = mb.Particle(name="_OHL", pos=[0, 1, 0]) - - system.add([first, second, third, fourth]) - - system.add_bond((first, second)) - system.add_bond((second, third)) - system.add_bond((second, fourth)) - - ff = Forcefield(forcefield_files=[get_fn("charmm36_cooh.xml")]) - struc = ff.apply( - system, - assert_angle_params=False, - assert_dihedral_params=False, - assert_improper_params=False, - ) - - write_lammpsdata(struc, "charmm_improper.lammps") - out_lammps = open("charmm_improper.lammps", "r").readlines() - found_impropers = False - for i, line in enumerate(out_lammps): - if "Improper Coeffs" in line: - assert "# harmonic" in line - assert "k, phi" in out_lammps[i + 1] - assert len(out_lammps[i + 2].split("#")[0].split()) == 3 - assert out_lammps[i + 2].split("#")[0].split()[0] == "1" - found_impropers = True - assert found_impropers - - @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed") - def test_amber(self): - from foyer import Forcefield - - from mbuild.formats.lammpsdata import write_lammpsdata - - cmpd = mb.load("C1(=CC=CC=C1)F", smiles=True) - - ff = Forcefield(forcefield_files=[get_fn("gaff_test.xml")]) - structure = ff.apply(cmpd) - - write_lammpsdata( - structure, "amber.lammps", zero_dihedral_weighting_factor=True - ) - out_lammps = open("amber.lammps", "r").readlines() - found_angles = False - found_dihedrals = False - found_impropers = False - for i, line in enumerate(out_lammps): - if "Angle Coeffs" in line: - assert "# harmonic" in line - assert ( - "#\tk(kcal/mol/rad^2)\t\ttheteq(deg)" in out_lammps[i + 1] - ) - assert len(out_lammps[i + 2].split("#")[0].split()) == 3 - assert out_lammps[i + 2].split("#")[0].split()[0] == "1" - found_angles = True - elif "Dihedral Coeffs" in line: - assert "# charmm" in line - assert "#k, n, phi, weight" in out_lammps[i + 1] - assert len(out_lammps[i + 2].split("#")[0].split()) == 5 - assert float(out_lammps[i + 2].split("#")[0].split()[4]) == 0.0 - assert out_lammps[i + 2].split("#")[0].split()[0] == "1" - found_dihedrals = True - elif "Improper Coeffs" in line: - assert "# cvff" in line - assert "#K, d, n" in out_lammps[i + 1] - assert len(out_lammps[i + 2].split("#")[0].split()) == 4 - assert out_lammps[i + 2].split("#")[0].split()[2] == "-1" - assert out_lammps[i + 2].split("#")[0].split()[0] == "1" - found_impropers = True - else: - pass - assert found_angles - assert found_dihedrals - assert found_impropers - - @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed") - def test_save_forcefield_with_same_struct(self): - from foyer import Forcefield - - from mbuild.formats.lammpsdata import write_lammpsdata - - system = mb.load("C1(=CC=CC=C1)F", smiles=True) - - ff = Forcefield(forcefield_files=[get_fn("gaff_test.xml")]) - struc = ff.apply( - system, - assert_angle_params=False, - assert_dihedral_params=False, - assert_improper_params=False, - ) - write_lammpsdata( - struc, "charmm_improper.lammps", zero_dihedral_weighting_factor=True - ) - for i in range(3): - xyz = struc.coordinates - xyz = xyz + np.array([1, 1, 1]) - struc.coordinates = xyz - write_lammpsdata( - struc, - f"charmm_improper{i}.lammps", - zero_dihedral_weighting_factor=True, - ) - - @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed") - @pytest.mark.parametrize("unit_style", ["real", "lj"]) - def test_save_box(self, ethane, unit_style): - box = mb.Box( - lengths=np.array([2.0, 2.0, 2.0]), angles=[90.0, 90.0, 90.0] - ) - ethane.save( - filename="ethane-box.lammps", - forcefield_name="oplsaa", - box=box, - unit_style=unit_style, - ) - - def test_nbfix(self, ethane): - from foyer import Forcefield - - OPLSAA = Forcefield(name="oplsaa") - structure = OPLSAA.apply(ethane) - structure.combining_rule = "lorentz" - # Add nbfixes - types = list(set([a.atom_type for a in structure.atoms])) - types[0].add_nbfix(types[1].name, 1.2, 2.1) - types[1].add_nbfix(types[0].name, 1.2, 2.1) - write_lammpsdata(filename="nbfix.lammps", structure=structure) - - checked_section = False - with open("nbfix.lammps", "r") as fi: - while not checked_section: - line = fi.readline() - if "PairIJ Coeffs" in line: - fi.readline() - line = fi.readline().partition("#")[0] - assert np.allclose( - np.asarray(line.split(), dtype=float), - [1, 1, 0.066, 3.5], - ) - line = fi.readline().partition("#")[0] - assert np.allclose( - np.asarray(line.split(), dtype=float), - [1, 2, 2.1, 1.06907846], - ) - line = fi.readline().partition("#")[0] - assert np.allclose( - np.asarray(line.split(), dtype=float), [2, 2, 0.03, 2.5] - ) - line = fi.readline() - checked_section = True - # Break if PairIJ Coeffs is not found - if "Atoms" in line: - break - structure.combining_rule = "geometric" - write_lammpsdata(filename="nbfix.lammps", structure=structure) - write_lammpsdata( - filename="nbfix.lammps", - structure=structure, - nbfix_in_data_file=False, - ) - with pytest.raises(ValueError) as exc_info: - structure.combining_rule = "error" - write_lammpsdata(filename="nbfix.lammps", structure=structure) - assert str(exc_info.value).startswith("combining_rule must be") - - def test_save_triclinic_box(self, ethane): - box = mb.Box(lengths=np.array([2.0, 2.0, 2.0]), angles=[60, 70, 80]) - ethane.save( - filename="triclinic-box.lammps", forcefield_name="oplsaa", box=box - ) - - def test_save_orthorhombic_box(self, ethane): - box = mb.Box(lengths=np.array([2.0, 2.0, 2.0])) - ethane.save( - filename="ortho-box.lammps", forcefield_name="oplsaa", box=box - ) - with open("ortho-box.lammps", "r") as fi: - for line in fi: - assert "xy" not in line - assert "xz" not in line - assert "yz" not in line - - @pytest.mark.parametrize( - "atom_style, n_columns", - [("full", 7), ("atomic", 5), ("molecular", 6), ("charge", 6)], - ) - def test_writing_atom_styles(self, ethane, atom_style, n_columns): - ethane.save(filename="ethane.lammps", atom_style=atom_style) - with open("ethane.lammps", "r") as f: - for line in f: - if "Atoms" not in line: - continue - atoms_header = next(f) - first_atom_line = next(f) - columns = first_atom_line.split("\t") - assert len(columns) == n_columns - else: - assert "# " + atom_style in line - - @pytest.mark.parametrize( - "offset, expected_value", [(0, ("0", "1")), (1, ("1", "2"))] - ) - def test_resid(self, offset, expected_value, ethane, methane): - structure = ethane.to_parmed() + methane.to_parmed() - n_atoms = len(structure.atoms) - write_lammpsdata(structure, "compound.lammps", moleculeID_offset=offset) - res_list = list() - with open("compound.lammps", "r") as f: - for i, line in enumerate(f): - if "Atoms" in line: - break - atom_lines = open("compound.lammps", "r").readlines()[ - i + 2 : i + n_atoms + 2 - ] - for line in atom_lines: - res_list.append(line.rstrip().split()[1]) - assert set(res_list) == set(expected_value) - - def test_box_bounds(self, ethane): - from foyer import Forcefield - - OPLSAA = Forcefield(name="oplsaa") - structure = OPLSAA.apply(ethane) - box = mb.Box.from_mins_maxs_angles( - mins=np.array([-1.0, -2.0, -3.0]), - maxs=np.array([3.0, 2.0, 1.0]), - angles=[90.0, 90.0, 90.0], - ) - # box = ethane.get_boundingbox() - - write_lammpsdata( - filename="box.lammps", - structure=structure, - unit_style="real", - mins=[0.0, 0.0, 0.0], - maxs=[m for m in box.lengths], - ) - - checked_section = False - # NOTE, need to figure out how to handle lo and hi of a compound - with open("box.lammps", "r") as fi: - while not checked_section: - line = fi.readline() - if "xlo" in line: - xlo = float(line.split()[0]) - xhi = float(line.split()[1]) - assert np.isclose(xlo, 0.0) - assert np.isclose(xhi, 40.0) - - line = fi.readline() - ylo = float(line.split()[0]) - yhi = float(line.split()[1]) - assert np.isclose(ylo, 0.0) - assert np.isclose(yhi, 40.0) - - line = fi.readline() - zlo = float(line.split()[0]) - zhi = float(line.split()[1]) - assert np.isclose(zlo, 0.0) - assert np.isclose(zhi, 40.0) - - checked_section = True - - write_lammpsdata( - filename="box.lammps", structure=structure, unit_style="real" - ) - - checked_section = False - with open("box.lammps", "r") as fi: - while not checked_section: - line = fi.readline() - if "xlo" in line: - xlo = float(line.split()[0]) - xhi = float(line.split()[1]) - assert np.isclose(xlo, 0.0) - assert np.isclose(xhi, 7.13999987) - - line = fi.readline() - ylo = float(line.split()[0]) - yhi = float(line.split()[1]) - assert np.isclose(ylo, 0.0) - assert np.isclose(yhi, 7.93800011) - - line = fi.readline() - zlo = float(line.split()[0]) - zhi = float(line.split()[1]) - assert np.isclose(zlo, 0.0) - assert np.isclose(zhi, 6.646) - - checked_section = True - - def test_lj_box(self, ethane, lj_save): - fn = lj_save(ethane, Path.cwd()) - checked_section = False - ethane_sigma = 0.35 # nm - box_length = np.array(ethane.get_boundingbox().lengths) + [ - 0.5, - 0.5, - 0.5, - ] # 0.5nm buffer around box - box_length /= ethane_sigma # Unitless - with open(str(fn), "r") as fi: - while not checked_section: - line = fi.readline() - if "dihedral types" in line: # line before box info - fi.readline() - for i in range(len(box_length)): - line = float(fi.readline().split()[1]) - assert np.isclose(float(line), box_length[i]) - checked_section = True - - def test_lj_masses(self, ethane, lj_save): - fn = lj_save(ethane, Path.cwd()) - checked_section = False - with open(fn, "r") as fi: - while not checked_section: - line = fi.readline() - if "Masses" in line: - fi.readline() - line = float(fi.readline().split()[1]) - assert np.isclose(line, 1.00) - checked_section = True - - def test_real_masses(self, ethane, real_save): - fn = real_save(ethane, Path.cwd()) - checked_section = False - with open(fn, "r") as fi: - while not checked_section: - line = fi.readline() - if "Masses" in line: - fi.readline() - line = float(fi.readline().split()[1]) - assert np.isclose(line, 12.010780) - checked_section = True - - def test_lj_pairs(self, ethane, lj_save): - fn = lj_save(ethane, Path.cwd()) - checked_section = False - with open(fn, "r") as fi: - while not checked_section: - line = fi.readline() - if "Pair Coeffs" in line: - fi.readline() - line = fi.readline().split() - epsilon = float(line[1]) - sigma = float(line[2]) - assert np.isclose(epsilon, 1.00) - assert np.isclose(sigma, 1.00) - checked_section = True - - def test_lj_passed_pairs(self, ethane, lj_save): - fn = lj_save( - ethane, - Path.cwd(), - sigma=1, - epsilon=1, - ) - checked_section = False - with open(fn, "r") as fi: - while not checked_section: - line = fi.readline() - if "Pair Coeffs" in line: - fi.readline() - line = fi.readline().split() - epsilon = float(line[1]) - sigma = float(line[2]) - assert np.isclose(epsilon, 0.276144, atol=1e-5) - assert np.isclose(sigma, 0.35) - checked_section = True - - with pytest.raises(ValueError) as exc_info: - lj_save(ethane, Path.cwd(), sigma=0) - exception_str = "The sigma conversion factor to convert to LJ" - assert str(exc_info.value).startswith(exception_str) - - def test_real_pairs(self, ethane, real_save): - fn = real_save(ethane, Path.cwd()) - checked_section = False - with open(fn, "r") as fi: - while not checked_section: - line = fi.readline() - if "Pair Coeffs" in line: - fi.readline() - line = fi.readline().split() - epsilon = float(line[1]) - sigma = float(line[2]) - assert np.isclose(epsilon, 0.066) - assert np.isclose(sigma, 3.5) - checked_section = True - - def test_lj_bonds(self, ethane, lj_save): - fn = lj_save(ethane, Path.cwd()) - checked_section = False - ethane_bondsk = ( - np.array([224262.4, 284512.0]) * (0.35) ** 2 / (0.276144) / 2 - ) # kj/mol/nm**2 to lammps - ethane_bondsreq = np.array([0.1529, 0.109]) * 0.35**-1 # unitless - ethane_lj_bonds = zip(ethane_bondsk, ethane_bondsreq) - with open(fn, "r") as fi: - while not checked_section: - line = fi.readline() - if "Bond Coeffs" in line: - fi.readline() - for i, bond_params in enumerate(ethane_lj_bonds): - line = fi.readline() - assert np.allclose( - (float(line.split()[1]), float(line.split()[2])), - bond_params, - atol=1e-3, - ) - assert int(line.split()[0]) == i + 1 - checked_section = True - - def test_real_bonds(self, ethane, real_save): - fn = real_save(ethane, Path.cwd()) - checked_section = False - ethane_bondsk = ( - np.array([224262.4, 284512.0]) / KCAL_TO_KJ * ANG_TO_NM**2 / 2 - ) # kj/mol/nm**2 to lammps - ethane_bondsreq = np.array([0.1529, 0.109]) / ANG_TO_NM - ethane_real_bonds = zip(ethane_bondsk, ethane_bondsreq) - with open(fn, "r") as fi: - while not checked_section: - line = fi.readline() - if "Bond Coeffs" in line: - fi.readline() - for bond_params in ethane_real_bonds: - line = fi.readline() - assert np.allclose( - (float(line.split()[1]), float(line.split()[2])), - bond_params, - atol=1e-3, - ) - checked_section = True - - def test_lj_angles(self, ethane, lj_save): - fn = lj_save(ethane, Path.cwd()) - checked_section = False - ethane_anglesk = ( - np.array([313.8, 276.144]) / 0.276144 / 2 - ) # kj/mol/nm**2 to lammps - ethane_anglesreq = ( - np.array([1.93207948196, 1.88146493365]) * 180 / np.pi - ) # radians to lammps - ethane_lj_angles = zip(ethane_anglesk, ethane_anglesreq) - with open(fn, "r") as fi: - while not checked_section: - line = fi.readline() - if "Angle Coeffs" in line: - fi.readline() - for angle_params in ethane_lj_angles: - line = fi.readline() - assert np.allclose( - (float(line.split()[1]), float(line.split()[2])), - angle_params, - ) - checked_section = True - - def test_real_angles(self, ethane, real_save): - fn = real_save(ethane, Path.cwd()) - checked_section = False - ethane_anglesk = ( - np.array([313.8, 276.144]) / 2 / KCAL_TO_KJ - ) # kj/mol to lammps kcal/mol - ethane_anglesreq = ( - np.array([1.93207948196, 1.88146493365]) * 180 / np.pi - ) # radians to lammps - ethane_real_angles = zip(ethane_anglesk, ethane_anglesreq) - with open(fn, "r") as fi: - while not checked_section: - line = fi.readline() - if "Angle Coeffs" in line: - fi.readline() - for i, angle_params in enumerate(ethane_real_angles): - line = fi.readline() - assert np.allclose( - (float(line.split()[1]), float(line.split()[2])), - angle_params, - ) - assert int(line.split()[0]) == i + 1 - checked_section = True - - def test_lj_dihedrals(self, ethane, lj_save): - fn = lj_save(ethane, Path.cwd()) - checked_section = False - ethane_diheds = ( - np.array([-0.00000, -0.00000, 0.30000, -0.00000]) / 0.066 - ) # kcal/mol to lammps - with open(fn, "r") as fi: - while not checked_section: - line = fi.readline() - if "Dihedral Coeffs" in line: - fi.readline() - line = fi.readline() - assert np.allclose( - ( - float(line.split()[1]), - float(line.split()[2]), - float(line.split()[3]), - float(line.split()[4]), - ), - ethane_diheds, - ) - checked_section = True - - def test_real_dihedrals(self, ethane, real_save): - fn = real_save(ethane, Path.cwd()) - checked_section = False - ethane_diheds = np.array([-0.00000, -0.00000, 0.30000, -0.00000]) - with open(fn, "r") as fi: - while not checked_section: - line = fi.readline() - if "Dihedral Coeffs" in line: - fi.readline() - line = fi.readline() - assert np.allclose( - ( - float(line.split()[1]), - float(line.split()[2]), - float(line.split()[3]), - float(line.split()[4]), - ), - ethane_diheds, - ) - assert int(line.split()[0]) == 1 - checked_section = True - - def test_lj_charges(self, ethane, lj_save): - fn = lj_save(ethane, Path.cwd()) - checked_section = False - charge_dict = {"C": -0.18, "H": 0.06} - ethane_charges = np.array( - [charge_dict[part.name] for part in ethane.particles()] - ) - # in coulombs, convert to lj units - ethane_charges /= ( - np.sqrt(4 * np.pi * 0.276144 * 0.35 * epsilon_0 * 10**-6) - / ELEM_TO_COUL - ) - print(ethane_charges) - with open(fn, "r") as fi: - while not checked_section: - line = fi.readline() - if "Atoms " in line: - fi.readline() - assert np.allclose( - float(fi.readline().split()[3]), - ethane_charges[0], - atol=1e-15, - ) - assert np.allclose( - float(fi.readline().split()[3]), - ethane_charges[1], - atol=1e-15, - ) - checked_section = True - - def test_real_charges(self, ethane, real_save): - fn = real_save(ethane, Path.cwd()) - checked_section = False - charge_dict = {"C": -0.18, "H": 0.06} - ethane_charges = [charge_dict[part.name] for part in ethane.particles()] - print(ethane_charges) - # in coulombs - with open(fn, "r") as fi: - while not checked_section: - line = fi.readline() - if "Atoms " in line: - fi.readline() - assert np.allclose( - float(fi.readline().split()[3]), - ethane_charges[0], - ) - assert np.allclose( - float(fi.readline().split()[3]), - ethane_charges[1], - ) - checked_section = True diff --git a/mbuild/tests/test_par.py b/mbuild/tests/test_par.py deleted file mode 100644 index 732deb867..000000000 --- a/mbuild/tests/test_par.py +++ /dev/null @@ -1,40 +0,0 @@ -import numpy as np -import pytest - -import mbuild as mb -from mbuild.formats.par_writer import write_par -from mbuild.tests.base_test import BaseTest -from mbuild.utils.io import get_fn, has_foyer - - -@pytest.mark.skipif(not has_foyer, reason="Foyer package not installed") -class TestPar(BaseTest): - @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed") - def test_save_charmm(self): - cmpd = mb.load(get_fn("charmm_dihedral.mol2")) - for i in cmpd.particles(): - i.name = "_{}".format(i.name) - structure = cmpd.to_parmed( - box=cmpd.get_boundingbox(), - residues=set([p.parent.name for p in cmpd.particles()]), - ) - - from foyer import Forcefield - - ff = Forcefield(forcefield_files=[get_fn("charmm_truncated.xml")]) - structure = ff.apply(structure, assert_dihedral_params=False) - - write_par(structure, "charmm_dihedral.par") - - @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed") - def test_save_forcefield(self, ethane): - ethane.save(filename="ethane-opls.par", forcefield_name="oplsaa") - - def test_par_parameters(self, ethane): - ethane.save(filename="ethane-opls.par", forcefield_name="oplsaa") - from parmed.charmm import CharmmParameterSet - - pset = CharmmParameterSet.load_set(pfile="ethane-opls.par") - assert len(pset.bond_types) == 3 - assert len(pset.angle_types) == 3 - assert len(pset.atom_types) == 2 diff --git a/mbuild/tests/test_protobuf.py b/mbuild/tests/test_protobuf.py deleted file mode 100644 index 001b8df55..000000000 --- a/mbuild/tests/test_protobuf.py +++ /dev/null @@ -1,18 +0,0 @@ -import pytest - -from mbuild.tests.base_test import BaseTest -from mbuild.utils.io import has_protobuf - - -class TestPB2(BaseTest): - @pytest.mark.skipif( - not has_protobuf, reason="Protobuf package not installed" - ) - def test_loop(self, ethane): - from mbuild.formats.protobuf import read_pb2, write_pb2 - - write_pb2(ethane, "ethane.pb2") - proto = read_pb2("ethane.pb2") - assert ethane.n_particles == proto.n_particles - assert ethane.n_bonds == proto.n_bonds - assert len(ethane.children) == len(proto.children) diff --git a/mbuild/tests/test_rigid.py b/mbuild/tests/test_rigid.py index 41b61baec..d3739bd0a 100644 --- a/mbuild/tests/test_rigid.py +++ b/mbuild/tests/test_rigid.py @@ -3,7 +3,7 @@ import mbuild as mb from mbuild.tests.base_test import BaseTest -from mbuild.utils.io import get_fn +from mbuild.utils.io import has_hoomd class TestRigid(BaseTest): @@ -114,6 +114,7 @@ def test_label_rigid_bodies_single_partial(self, benzene): assert len(list(benzene.rigid_particles())) == 6 assert len(list(benzene.rigid_particles(rigid_id=0))) == 6 + @pytest.mark.skipif(not has_hoomd, reason="HOOMD not installed.") def test_save_non_sequential_rigid_ids(self, benzene): n_benzenes = 10 filled = mb.fill_box( @@ -122,7 +123,7 @@ def test_save_non_sequential_rigid_ids(self, benzene): filled.label_rigid_bodies(discrete_bodies="Benzene") filled.children[0]._increment_rigid_ids(increment=3) with pytest.warns(UserWarning): - filled.save("benzene-box.hoomdxml") + filled.save("benzene-box.gsd") def test_increment_rigid_id(self, rigid_benzene): compound = mb.Compound() diff --git a/mbuild/tests/test_xyz.py b/mbuild/tests/test_xyz.py deleted file mode 100644 index a567c29a7..000000000 --- a/mbuild/tests/test_xyz.py +++ /dev/null @@ -1,72 +0,0 @@ -import ele -import numpy as np -import pytest - -import mbuild as mb -from mbuild.exceptions import MBuildError -from mbuild.formats.xyz import write_xyz -from mbuild.tests.base_test import BaseTest -from mbuild.utils.io import get_fn - - -class TestXYZ(BaseTest): - def test_load_no_top(self, ethane): - ethane.save(filename="ethane.xyz") - ethane_in = mb.load("ethane.xyz") - assert len(ethane_in.children) == 8 - assert ethane_in.n_bonds == 0 - assert set([child.name for child in ethane_in.children]) == {"C", "H"} - assert set([child.element for child in ethane_in.children]) == { - ele.Elements.C, - ele.Elements.H, - } - - def test_wrong_n_atoms(self): - with pytest.raises(MBuildError): - mb.load(get_fn("too_few_atoms.xyz")) - with pytest.raises(MBuildError): - mb.load(get_fn("too_many_atoms.xyz")) - - def test_bad_input(self, ethane): - with pytest.raises(ValueError): - assert isinstance(ethane, mb.Compound) - write_xyz(ethane, "compound.xyz") - - def test_save(self, ethane): - ethane.save(filename="ethane.xyz") - ethane_in = mb.load("ethane.xyz") - assert len(ethane_in.children) == 8 - assert ethane_in.n_bonds == 0 - assert set([child.name for child in ethane_in.children]) == {"C", "H"} - - def test_coordinates(self, ethane): - ethane.save(filename="ethane.xyz") - ethane_in = mb.load("ethane.xyz") - assert np.allclose(ethane.xyz, ethane_in.xyz) - - def test_non_resolved_elements(self): - tip3p_water = mb.load(get_fn("tip3p_water.xyz")) - assert tip3p_water[0].element is None - assert tip3p_water[0].name == "opls_111" - assert tip3p_water[1].element is None - assert tip3p_water[1].name == "opls_112" - assert tip3p_water[2].element is None - assert tip3p_water[2].name == "opls_112" - - def test_write_atomtypes(self): - tip3p_water = mb.load(get_fn("tip3p_water.xyz")) - tip3p_water.save(filename="test.xyz", write_atomnames=True) - lines = [] - with open("test.xyz") as f: - for line in f: - lines.append(line.split()) - assert lines[2][0] == "opls_111" - assert lines[3][0] == "opls_112" - assert lines[4][0] == "opls_112" - tip3p_water_in = mb.load("test.xyz") - assert tip3p_water_in[0].element is None - assert tip3p_water_in[0].name == "opls_111" - assert tip3p_water_in[1].element is None - assert tip3p_water_in[1].name == "opls_112" - assert tip3p_water_in[2].element is None - assert tip3p_water_in[2].name == "opls_112" diff --git a/mbuild/utils/io.py b/mbuild/utils/io.py index 1bc274c46..7c2f2e0f0 100644 --- a/mbuild/utils/io.py +++ b/mbuild/utils/io.py @@ -159,20 +159,6 @@ class DelayImportError(ImportError, SkipTest): # conda install -c conda-forge freud """ -MESSAGES[ - "protobuf" -] = """ -The code at {filename}:{line_number} requires the "protobuf" package - -protobuf can be installed using: - -# conda install -c conda-forge protobuf - -or - -# pip install protobuf -""" - def import_(module): """Import a module and issue a nice message to stderr if it isn't installed. @@ -273,6 +259,14 @@ def import_(module): raise DelayImportError(m) +try: + import hoomd + + has_hoomd = True + del hoomd +except ImportError: + has_hoomd = False + try: import intermol @@ -344,14 +338,6 @@ def import_(module): except ImportError: has_py3Dmol = False -try: - from google import protobuf - - has_protobuf = True - del protobuf -except ImportError: - has_protobuf = False - try: import garnett