diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml index a1af8e8ed..94c5be433 100644 --- a/.github/workflows/CI.yaml +++ b/.github/workflows/CI.yaml @@ -4,9 +4,11 @@ on: push: branches: - "main" + - "develop" pull_request: branches: - "main" + - "develop" schedule: - cron: "0 0 * * *" @@ -19,7 +21,7 @@ jobs: fail-fast: false matrix: os: [ubuntu-latest] - python-version: ["3.9", "3.10", "3.11", "3.12"] + python-version: ["3.10", "3.11", "3.12"] defaults: run: @@ -39,6 +41,10 @@ jobs: - name: Install Package run: python -m pip install -e . + - name: Conditionally install OpenBabel + if: ${{ matrix.python-version != '3.13' }} + run: micromamba install -y openbabel + - 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 @@ -57,7 +63,7 @@ jobs: fail-fast: false matrix: os: [macOS-latest, macOS-13, ubuntu-latest] - python-version: ["3.12"] + python-version: ["3.13"] defaults: run: @@ -89,7 +95,7 @@ jobs: - uses: actions/checkout@v4 name: Checkout Branch / Pull Request - - uses: Vampire/setup-wsl@v4 + - uses: Vampire/setup-wsl@v6 with: distribution: Ubuntu-24.04 wsl-shell-user: runner diff --git a/.github/workflows/codeql.yml b/.github/workflows/codeql.yml index 42af8c633..dd545ca8a 100644 --- a/.github/workflows/codeql.yml +++ b/.github/workflows/codeql.yml @@ -2,9 +2,9 @@ name: "CodeQL" on: push: - branches: [ "main" ] + branches: [ "develop" ] pull_request: - branches: [ "main" ] + branches: [ "develop" ] schedule: - cron: "50 6 * * 0" @@ -27,15 +27,15 @@ jobs: uses: actions/checkout@v3 - name: Initialize CodeQL - uses: github/codeql-action/init@v2 + uses: github/codeql-action/init@v3 with: languages: ${{ matrix.language }} queries: +security-and-quality - name: Autobuild - uses: github/codeql-action/autobuild@v2 + uses: github/codeql-action/autobuild@v3 - name: Perform CodeQL Analysis - uses: github/codeql-action/analyze@v2 + uses: github/codeql-action/analyze@v3 with: category: "/language:${{ matrix.language }}" diff --git a/.gitignore b/.gitignore index c312a17b1..0f4c9cdfd 100755 --- a/.gitignore +++ b/.gitignore @@ -3,6 +3,7 @@ test-output.xml *.pymon *.ipynb_checkpoints *DS_Store* +__pycache__ # C extensions *.so diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index c4f2389b2..6d303b899 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -10,7 +10,7 @@ ci: repos: - repo: https://github.com/astral-sh/ruff-pre-commit # Ruff version. - rev: v0.11.13 + rev: v0.12.10 hooks: # Run the linter. - id: ruff @@ -18,7 +18,7 @@ repos: # Run the formatter. - id: ruff-format - repo: https://github.com/pre-commit/pre-commit-hooks - rev: v5.0.0 + rev: v6.0.0 hooks: - id: check-yaml - id: end-of-file-fixer diff --git a/README.md b/README.md index 74b359520..23e21d571 100644 --- a/README.md +++ b/README.md @@ -61,6 +61,8 @@ conda activate mbuild-dev pip install . ``` +NOTE: [openbabel](https://github.com/openbabel/openbabel) is required for some energy minimization methods in `mbuild.compound.Compound()`. It can be installed into your mBuild environment from conda-forge with `conda install -c conda-forge openbabel`; however, openbabel does not yet support python 3.13. + #### Install an editable version from source Once all dependencies have been installed and the ``conda`` environment has been created, the ``mBuild`` itself can be installed. diff --git a/codecov.yml b/codecov.yml index ff3fe65c0..27ec07558 100644 --- a/codecov.yml +++ b/codecov.yml @@ -1,3 +1,4 @@ ignore: - "mbuild/examples" - "mbuild/tests" + - "mbuild/__init__.py" diff --git a/docs/conf.py b/docs/conf.py index df60e2c60..0c72b907c 100755 --- a/docs/conf.py +++ b/docs/conf.py @@ -158,8 +158,8 @@ # built documents. # -version = "1.2.0" -release = "1.2.0" +version = "1.2.1" +release = "1.2.1" # The language for content autogenerated by Sphinx. Refer to documentation diff --git a/environment-dev.yml b/environment-dev.yml index d28d6d7d8..c225668c0 100644 --- a/environment-dev.yml +++ b/environment-dev.yml @@ -2,15 +2,15 @@ name: mbuild-dev channels: - conda-forge dependencies: - - python>=3.9,<=3.12 + - python>=3.10,<=3.13 - boltons - - numpy=1.26.4 + - numba + - numpy>=2.0,<2.3 - sympy - unyt>=2.9.5 - boltons - lark>=1.2 - lxml - - freud>=3.0 - intermol - mdtraj - pydantic>=2 @@ -18,10 +18,10 @@ dependencies: - nglview>=3 - pytest - garnett>=0.7.1 - - openbabel>=3.0.0 - - openff-toolkit-base >=0.11,<0.16.7 + - openff-toolkit-base>0.16.7 - openmm - gsd>=2.9 + - freud>=3.2 - parmed>=3.4.3 - packmol>=20.15 - pytest-cov @@ -39,8 +39,8 @@ dependencies: - pandas - symengine - python-symengine - - hoomd>=4.0,<5.0 - py3Dmol + - hoomd>=4.0 - pip: - git+https://github.com/mosdef-hub/gmso.git@main - git+https://github.com/mosdef-hub/foyer.git@main diff --git a/environment.yml b/environment.yml index 2a9defc89..9fd312e70 100644 --- a/environment.yml +++ b/environment.yml @@ -3,13 +3,14 @@ channels: - conda-forge dependencies: - ele - - numpy=1.26.4 + - numpy>=2.0,<2.3 - packmol>=20.15 - - gmso>=0.9.0 + - gmso>=0.12.0 - garnett + - numba - parmed>=3.4.3 - pycifrw - - python>=3.8 + - python>=3.10,<=3.13 - rdkit>=2021 - scipy - networkx diff --git a/mbuild/__init__.py b/mbuild/__init__.py index 10f4cc06e..db92c729b 100644 --- a/mbuild/__init__.py +++ b/mbuild/__init__.py @@ -2,6 +2,10 @@ # ruff: noqa: F403 """mBuild: a hierarchical, component based molecule builder.""" +import logging +import sys +from logging.handlers import RotatingFileHandler + from mbuild.box import Box from mbuild.coarse_graining import coarse_grain from mbuild.compound import * @@ -10,8 +14,124 @@ from mbuild.lattice import Lattice from mbuild.packing import * from mbuild.pattern import * +from mbuild.polymer import Polymer from mbuild.port import Port from mbuild.recipes import recipes -__version__ = "1.2.0" +__version__ = "1.2.1" __date__ = "2025-01-23" + + +class DeduplicationFilter(logging.Filter): + """A logging filter that suppresses duplicate messages.""" + + def __init__(self): + super().__init__() + self.logged_messages = set() + + def filter(self, record): + log_entry = (record.name, record.levelno, record.msg) + if log_entry not in self.logged_messages: + self.logged_messages.add(log_entry) + return True + return False + + +class HeaderRotatingFileHandler(RotatingFileHandler): + def __init__( + self, + filename, + mode="w", + maxBytes=0, + backupCount=0, + encoding=None, + delay=False, + header="", + ): + self.header = header + super().__init__(filename, mode, maxBytes, backupCount, encoding, delay) + + def _open(self): + """ + Open the current base log file, with the header written. + """ + stream = super()._open() + if stream.tell() == 0 and self.header: # Only write header if file is empty + stream.write(self.header + "\n") + return stream + + +class mBuildLogger: + def __init__(self): + self.library_logger = logging.getLogger("mbuild") + self.library_logger.setLevel(logging.DEBUG) + + # Create handlers + self.console_handler = logging.StreamHandler(sys.stdout) + self.console_handler.setLevel(logging.WARNING) + + # Create a formatter + self.formatter = logging.Formatter( + "%(asctime)s - %(name)s - %(levelname)s - %(message)s" + ) + + # Add formatter to handlers + self.console_handler.setFormatter(self.formatter) + + # Initialize and add the deduplication filter + self.dedup_filter = DeduplicationFilter() + self.console_handler.addFilter(self.dedup_filter) + + # Clear any previous handlers to avoid duplicates in Jupyter + self._clear_handlers() + + # Add handlers to the library logger + self.library_logger.addHandler(self.console_handler) + + def _clear_handlers(self): + handlers = self.library_logger.handlers[:] + for handler in handlers: + self.library_logger.removeHandler(handler) + + def debug_file(self, filename: str): + """Print logging Debug messages to file `filename`.""" + # Get the path to the Python interpreter + python_executable = sys.executable + + # Get the list of command-line arguments + command_arguments = sys.argv + + # Construct the full command + full_command = [python_executable] + command_arguments + header = f"Log details for mBuild {__version__} from running \n{full_command}" + self.file_handler = HeaderRotatingFileHandler( + filename, mode="a", maxBytes=10**6, backupCount=2, header=header + ) + self.file_handler.setLevel(logging.DEBUG) + + self.file_handler.addFilter(DeduplicationFilter()) # fresh duplication handler + self.file_handler.setFormatter(self.formatter) + self.library_logger.addHandler(self.file_handler) + + def print_level(self, level: str): + """Print sys.stdout screen based on the logging `level` passed.""" + levelDict = { + "notset": logging.NOTSET, + "debug": logging.DEBUG, + "info": logging.INFO, + "warning": logging.WARNING, + "error": logging.ERROR, + "critical": logging.CRITICAL, + } + logLevel = levelDict.get(level.lower()) + if logLevel: + self.console_handler.setLevel(logLevel) # sets stdout + else: + raise ValueError( + f"INCORRECT {level=}. Please set level of {levelDict.keys()}" + ) + + +# Example usage in __init__.py +mbuild_logger = mBuildLogger() +mbuild_logger.library_logger.setLevel(logging.INFO) diff --git a/mbuild/box.py b/mbuild/box.py index 87f980eab..6264e14a2 100644 --- a/mbuild/box.py +++ b/mbuild/box.py @@ -1,6 +1,6 @@ """mBuild box module.""" -from warnings import warn +import logging import numpy as np @@ -8,6 +8,8 @@ __all__ = ["Box"] +logger = logging.getLogger(__name__) + class Box(object): """A box representing the bounds of the system. @@ -350,7 +352,7 @@ def _normalize_box(vectors): f"3D region in space.\n Box vectors evaluated: {vectors}" ) if det < 0.0: - warn( + logger.warning( "Box vectors provided for a left-handed basis, these will be " "transformed into a right-handed basis automatically." ) diff --git a/mbuild/compound.py b/mbuild/compound.py index f5d2ff7d6..57f8de55f 100644 --- a/mbuild/compound.py +++ b/mbuild/compound.py @@ -1,22 +1,19 @@ """Module for working with mBuild Compounds.""" -__all__ = ["clone", "Compound", "Particle"] - import itertools +import logging import os import tempfile from collections import OrderedDict from collections.abc import Iterable from copy import deepcopy from typing import Sequence -from warnings import warn import ele import networkx as nx import numpy as np from boltons.setutils import IndexedSet -from ele.element import Element, element_from_name, element_from_symbol -from ele.exceptions import ElementError +from ele.element import Element from treelib import Tree from mbuild import conversion @@ -25,9 +22,14 @@ from mbuild.coordinate_transform import _rotate, _translate from mbuild.exceptions import MBuildError from mbuild.periodic_kdtree import PeriodicKDTree +from mbuild.utils.geometry import bounding_box from mbuild.utils.io import import_, run_from_ipython from mbuild.utils.jsutils import overwrite_nglview_default +__all__ = ["clone", "Compound", "Particle"] + +logger = logging.getLogger(__name__) + def clone(existing_compound, clone_of=None, root_container=None): """Clone Compound. @@ -196,6 +198,8 @@ def __init__( else: self._charge = charge self._mass = mass + self._bond_tag = None + self._hoomd_data = {} def particles(self, include_ports=False): """Return all Particles of the Compound. @@ -240,6 +244,70 @@ def successors(self): for subpart in part.successors(): yield subpart + def get_child_indices(self, child): + """Gets the indices of particles belonging to a child compound. + + Parameters: + ----------- + child : mb.Compound + Compound that belongs to self.children + """ + if child not in self.children: + raise ValueError(f"{child} is not a child in this compound's hiearchy.") + parent_hashes = np.array([p.__hash__() for p in self.particles()]) + child_hashes = np.array([p.__hash__() for p in child.particles()]) + matching_indices = list( + int(np.where(parent_hashes == i)[0][0]) for i in child_hashes + ) + return matching_indices + + @property + def bond_tag(self): + return self._bond_tag + + @bond_tag.setter + def bond_tag(self, tag): + self._bond_tag = tag + + def set_bond_graph(self, new_graph): + """Manually set the compound's complete bond graph. + + Parameters: + ----------- + new_graph : networkx.Graph, required + A networkx.Graph containing information about + nodes and edges that is used to construct a + new mbuild.bond_graph.BondGraph(). + + Notes: + ------ + The nodes of the new graph should be ints that are mapped to + the indices of the particles in the compound. + + mbuild bond graphs accept `bond_order` arguments. If the + `bond_order` data is present in `new_graph`, it will be used, + otherwise `bond_under` is set to 'unspecified'. + """ + if new_graph.number_of_nodes() != self.n_particles: + raise ValueError( + f"new_graph contains {new_graph.number_of_nodes()} nodes", + f" but the Compound contains {self.n_particles} particles", + "The graph must contain one node per particle.", + ) + graph = BondGraph() + # node is int corresponding to particle index + for node in new_graph.nodes: + graph.add_node(node_for_adding=self[node]) + for edge in new_graph.edges: + if "bond_order" in edge: + bond_order = edge["bond_order"] + else: + bond_order = "unspecified" + graph.add_edge( + u_of_edge=self[edge[0]], v_of_edge=self[edge[1]], bond_order=bond_order + ) + self.bond_graph = graph + @property def n_particles(self): """Return the number of Particles in the Compound. @@ -499,7 +567,7 @@ def mass(self): else: particle_masses = [self._particle_mass(p) for p in self.particles()] if None in particle_masses: - warn( + logger.info( f"Some particle of {self} does not have mass." "They will not be accounted for during this calculation." ) @@ -543,8 +611,8 @@ def charge(self): return self._particle_charge(self) charges = [p._charge for p in self.particles()] if None in charges: - warn( - f"Some particle of {self} does not have a charge." + logger.info( + f"Some particle of {self} does not have a charge. " "They will not be accounted for during this calculation." ) filtered_charges = [charge for charge in charges if charge is not None] @@ -565,6 +633,17 @@ def charge(self, value): "not at the bottom of the containment hierarchy." ) + def volume(self): + """Estimate volume of the compound's particles. + Uses rdkit.Chem.AllChem.ComputeMolVolume + """ + from rdkit import Chem + + rdmol = self.to_rdkit(embed=True) + vol = Chem.AllChem.ComputeMolVolume(rdmol) + vol /= 1000 # Convert from cubic angstrom to cubic nm + return vol + def add( self, new_child, @@ -659,7 +738,7 @@ def add( f"to Compounds. You tried to add '{new_child}'." ) if self._mass is not None and not isinstance(new_child, Port): - warn( + logger.info( f"{self} has a pre-defined mass of {self._mass}, " "which will be reset to zero now that it contains children " "compounds." @@ -726,7 +805,7 @@ def add( else: if inherit_box: if new_child.box is None: - warn( + logger.info( "The Compound you are adding has no box but " "inherit_box=True. The box of the original " "Compound will remain unchanged." @@ -735,7 +814,7 @@ def add( self.box = new_child.box else: if new_child.box is not None: - warn( + logger.info( "The Compound you are adding has a box. " "The box of the parent compound will be used. Use " "inherit_box = True if you wish to replace the parent " @@ -747,7 +826,7 @@ def add( if ( np.array(self.box.lengths) < np.array(self.get_boundingbox().lengths) ).any(): - warn( + logger.warning( "After adding new Compound, Compound.box.lengths < " "Compound.boundingbox.lengths. There may be particles " "outside of the defined simulation box" @@ -801,7 +880,7 @@ def _check_if_empty(child): to_remove.append(child) _check_if_empty(child.parent) else: - warn(f"This will remove all particles in {self}") + logger.warning(f"This will remove all particles in {self}") return for particle in particles_to_remove: @@ -1211,12 +1290,14 @@ def remove_bond(self, particle_pair): if self.root.bond_graph is None or not self.root.bond_graph.has_edge( *particle_pair ): - warn("Bond between {} and {} doesn't exist!".format(*particle_pair)) + raise MBuildError( + "Bond between {} and {} doesn't exist!".format(*particle_pair) + ) return self.root.bond_graph.remove_edge(*particle_pair) bond_vector = particle_pair[0].pos - particle_pair[1].pos if np.allclose(bond_vector, np.zeros(3)): - warn( + logger.warning( "Particles {} and {} overlap! Ports will not be added.".format( *particle_pair ) @@ -1293,7 +1374,7 @@ def box(self, box): # Make sure the box is bigger than the bounding box if box is not None: if np.asarray((box.lengths < self.get_boundingbox().lengths)).any(): - warn( + logger.warning( "Compound.box.lengths < Compound.boundingbox.lengths. " "There may be particles outside of the defined " "simulation box." @@ -1525,40 +1606,7 @@ def get_boundingbox(self, pad_box=None): that are generated from mb.Lattice's and the resulting mb.Lattice.populate method """ - # case where only 1 particle exists - is_one_particle = False - if self.xyz.shape[0] == 1: - is_one_particle = True - - # are any columns all equalivalent values? - # an example of this would be a planar molecule - # example: all z values are 0.0 - # from: https://stackoverflow.com/a/14860884 - # steps: create mask array comparing first value in each column - # use np.all with axis=0 to do row columnar comparision - has_dimension = [True, True, True] - if not is_one_particle: - missing_dimensions = np.all( - np.isclose(self.xyz, self.xyz[0, :], atol=1e-2), - axis=0, - ) - for i, truthy in enumerate(missing_dimensions): - has_dimension[i] = not truthy - - if is_one_particle: - v1 = np.asarray([[1.0, 0.0, 0.0]]) - v2 = np.asarray([[0.0, 1.0, 0.0]]) - v3 = np.asarray([[0.0, 0.0, 1.0]]) - else: - v1 = np.asarray((self.maxs[0] - self.mins[0], 0.0, 0.0)) - v2 = np.asarray((0.0, self.maxs[1] - self.mins[1], 0.0)) - v3 = np.asarray((0.0, 0.0, self.maxs[2] - self.mins[2])) - vecs = [v1, v2, v3] - - # handle any missing dimensions (planar molecules) - for i, dim in enumerate(has_dimension): - if not dim: - vecs[i][i] = 0.1 + vecs = bounding_box(xyz=self.xyz) if pad_box is not None: if isinstance(pad_box, (int, float, str, Sequence)): @@ -1581,8 +1629,7 @@ def get_boundingbox(self, pad_box=None): for dim, val in enumerate(padding): vecs[dim][dim] = vecs[dim][dim] + val - bounding_box = Box.from_vectors(vectors=np.asarray([vecs]).reshape(3, 3)) - return bounding_box + return Box.from_vectors(vectors=np.asarray([vecs]).reshape(3, 3)) def min_periodic_distance(self, xyz0, xyz1): """Vectorized distance calculation considering minimum image. @@ -1620,7 +1667,9 @@ def min_periodic_distance(self, xyz0, xyz1): raise MBuildError(f'Cannot calculate minimum periodic distance. ' f'No Box set for {self}') """ - warn(f"No Box object set for {self}, using rectangular bounding box") + logger.warning( + f"No Box object set for {self}, using rectangular bounding box" + ) self.box = self.get_boundingbox() if np.allclose(self.box.angles, 90.0): d = np.where( @@ -1688,6 +1737,7 @@ def visualize( backend="py3dmol", color_scheme={}, bead_size=0.3, + show_bond_tags=False, ): # pragma: no cover """Visualize the Compound using py3dmol (default) or nglview. @@ -2015,714 +2065,6 @@ def _kick(self): particle.pos += (np.random.rand(3) - 0.5) / 100 self._update_port_locations(xyz_init) - def energy_minimize( - self, - forcefield="UFF", - steps=1000, - shift_com=True, - anchor=None, - **kwargs, - ): - """Perform an energy minimization on a Compound. - - Default behavior utilizes `Open Babel `_ - to perform an energy minimization/geometry optimization on a Compound by - applying a generic force field - - Can also utilize `OpenMM `_ to energy minimize after - atomtyping a Compound using - `Foyer `_ to apply a forcefield XML - file that contains valid SMARTS strings. - - This function is primarily intended to be used on smaller components, - with sizes on the order of 10's to 100's of particles, as the energy - minimization scales poorly with the number of particles. - - Parameters - ---------- - steps : int, optional, default=1000 - The number of optimization iterations - forcefield : str, optional, default='UFF' - The generic force field to apply to the Compound for minimization. - Valid options are 'MMFF94', 'MMFF94s', ''UFF', 'GAFF', 'Ghemical'. - Please refer to the `Open Babel documentation - `_ - when considering your choice of force field. - Utilizing OpenMM for energy minimization requires a forcefield - XML file with valid SMARTS strings. Please refer to `OpenMM docs - `_ - for more information. - shift_com : bool, optional, default=True - If True, the energy-minimized Compound is translated such that the - center-of-mass is unchanged relative to the initial configuration. - anchor : Compound, optional, default=None - Translates the energy-minimized Compound such that the - position of the anchor Compound is unchanged relative to the - initial configuration. - - - Other Parameters - ---------------- - algorithm : str, optional, default='cg' - The energy minimization algorithm. Valid options are 'steep', 'cg', - and 'md', corresponding to steepest descent, conjugate gradient, and - equilibrium molecular dynamics respectively. - For _energy_minimize_openbabel - fixed_compounds : Compound, optional, default=None - An individual Compound or list of Compounds that will have their - position fixed during energy minimization. Note, positions are fixed - using a restraining potential and thus may change slightly. - Position fixing will apply to all Particles (i.e., atoms) that exist - in the Compound and to particles in any subsequent sub-Compounds. - By default x,y, and z position is fixed. This can be toggled by instead - passing a list containing the Compound and an list or tuple of bool values - corresponding to x,y and z; e.g., [Compound, (True, True, False)] - will fix the x and y position but allow z to be free. - For _energy_minimize_openbabel - ignore_compounds: Compound, optional, default=None - An individual compound or list of Compounds whose underlying particles - will have their positions fixed and not interact with other atoms via - the specified force field during the energy minimization process. - Note, a restraining potential used and thus absolute position may vary - as a result of the energy minimization process. - Interactions of these ignored atoms can be specified by the user, - e.g., by explicitly setting a distance constraint. - For _energy_minimize_openbabel - distance_constraints: list, optional, default=None - A list containing a pair of Compounds as a tuple or list and - a float value specifying the target distance between the two Compounds, e.g.,: - [(compound1, compound2), distance]. - To specify more than one constraint, pass constraints as a 2D list, e.g.,: - [ [(compound1, compound2), distance1], [(compound3, compound4), distance2] ]. - Note, Compounds specified here must represent individual point particles. - For _energy_minimize_openbabel - constraint_factor: float, optional, default=50000.0 - Harmonic springs are used to constrain distances and fix atom positions, where - the resulting energy associated with the spring is scaled by the - constraint_factor; the energy of this spring is considering during the minimization. - As such, very small values of the constraint_factor may result in an energy - minimized state that does not adequately restrain the distance/position of atoms. - For _energy_minimize_openbabel - scale_bonds : float, optional, default=1 - Scales the bond force constant (1 is completely on). - For _energy_minimize_openmm - scale_angles : float, optional, default=1 - Scales the angle force constant (1 is completely on) - For _energy_minimize_openmm - scale_torsions : float, optional, default=1 - Scales the torsional force constants (1 is completely on) - For _energy_minimize_openmm - Note: Only Ryckaert-Bellemans style torsions are currently supported - scale_nonbonded : float, optional, default=1 - Scales epsilon (1 is completely on) - For _energy_minimize_openmm - constraints : str, optional, default="AllBonds" - Specify constraints on the molecule to minimize, options are: - None, "HBonds", "AllBonds", "HAngles" - For _energy_minimize_openmm - - References - ---------- - If using _energy_minimize_openmm(), please cite: - - .. [Eastman2013] P. Eastman, M. S. Friedrichs, J. D. Chodera, - R. J. Radmer, C. M. Bruns, J. P. Ku, K. A. Beauchamp, T. J. Lane, - L.-P. Wang, D. Shukla, T. Tye, M. Houston, T. Stich, C. Klein, - M. R. Shirts, and V. S. Pande. "OpenMM 4: A Reusable, Extensible, - Hardware Independent Library for High Performance Molecular - Simulation." J. Chem. Theor. Comput. 9(1): 461-469. (2013). - - If using _energy_minimize_openbabel(), please cite: - - .. [OBoyle2011] O'Boyle, N.M.; Banck, M.; James, C.A.; Morley, C.; - Vandermeersch, T.; Hutchison, G.R. "Open Babel: An open chemical - toolbox." (2011) J. Cheminf. 3, 33 - - .. [OpenBabel] Open Babel, version X.X.X http://openbabel.org, - (installed Month Year) - - If using the 'MMFF94' force field please also cite the following: - - .. [Halgren1996a] T.A. Halgren, "Merck molecular force field. I. Basis, - form, scope, parameterization, and performance of MMFF94." (1996) - J. Comput. Chem. 17, 490-519 - - .. [Halgren1996b] T.A. Halgren, "Merck molecular force field. II. MMFF94 - van der Waals and electrostatic parameters for intermolecular - interactions." (1996) J. Comput. Chem. 17, 520-552 - - .. [Halgren1996c] T.A. Halgren, "Merck molecular force field. III. - Molecular geometries and vibrational frequencies for MMFF94." (1996) - J. Comput. Chem. 17, 553-586 - - .. [Halgren1996d] T.A. Halgren and R.B. Nachbar, "Merck molecular force - field. IV. Conformational energies and geometries for MMFF94." (1996) - J. Comput. Chem. 17, 587-615 - - .. [Halgren1996e] T.A. Halgren, "Merck molecular force field. V. - Extension of MMFF94 using experimental data, additional computational - data, and empirical rules." (1996) J. Comput. Chem. 17, 616-641 - - If using the 'MMFF94s' force field please cite the above along with: - - .. [Halgren1999] T.A. Halgren, "MMFF VI. MMFF94s option for energy minimization - studies." (1999) J. Comput. Chem. 20, 720-729 - - If using the 'UFF' force field please cite the following: - - .. [Rappe1992] Rappe, A.K., Casewit, C.J., Colwell, K.S., Goddard, W.A. - III, Skiff, W.M. "UFF, a full periodic table force field for - molecular mechanics and molecular dynamics simulations." (1992) - J. Am. Chem. Soc. 114, 10024-10039 - - If using the 'GAFF' force field please cite the following: - - .. [Wang2004] Wang, J., Wolf, R.M., Caldwell, J.W., Kollman, P.A., - Case, D.A. "Development and testing of a general AMBER force field" - (2004) J. Comput. Chem. 25, 1157-1174 - - If using the 'Ghemical' force field please cite the following: - - .. [Hassinen2001] T. Hassinen and M. Perakyla, "New energy terms for - reduced protein models implemented in an off-lattice force field" - (2001) J. Comput. Chem. 22, 1229-1242 - - """ - # TODO: Update mbuild tutorials to provide overview of new features - # Preliminary tutorials: https://github.com/chrisiacovella/mbuild_energy_minimization - com = self.pos - anchor_in_compound = False - if anchor is not None: - # check to see if the anchor exists - # in the Compound to be energy minimized - for succesor in self.successors(): - if id(anchor) == id(succesor): - anchor_in_compound = True - anchor_pos_old = anchor.pos - - if not anchor_in_compound: - raise MBuildError( - f"Anchor: {anchor} is not part of the Compound: {self}" - "that you are trying to energy minimize." - ) - self._kick() - extension = os.path.splitext(forcefield)[-1] - openbabel_ffs = ["MMFF94", "MMFF94s", "UFF", "GAFF", "Ghemical"] - if forcefield in openbabel_ffs: - self._energy_minimize_openbabel( - forcefield=forcefield, steps=steps, **kwargs - ) - else: - tmp_dir = tempfile.mkdtemp() - self.save(os.path.join(tmp_dir, "un-minimized.mol2")) - - if extension == ".xml": - self._energy_minimize_openmm( - tmp_dir, - forcefield_files=forcefield, - forcefield_name=None, - steps=steps, - **kwargs, - ) - else: - self._energy_minimize_openmm( - tmp_dir, - forcefield_files=None, - forcefield_name=forcefield, - steps=steps, - **kwargs, - ) - - self.update_coordinates(os.path.join(tmp_dir, "minimized.pdb")) - - if shift_com: - self.translate_to(com) - - if anchor_in_compound: - anchor_pos_new = anchor.pos - delta = anchor_pos_old - anchor_pos_new - self.translate(delta) - - def _energy_minimize_openmm( - self, - tmp_dir, - forcefield_files=None, - forcefield_name=None, - steps=1000, - scale_bonds=1, - scale_angles=1, - scale_torsions=1, - scale_nonbonded=1, - constraints="AllBonds", - ): - """Perform energy minimization using OpenMM. - - Converts an mBuild Compound to a ParmEd Structure, - applies a forcefield using Foyer, and creates an OpenMM System. - - Parameters - ---------- - forcefield_files : str or list of str, optional, default=None - Forcefield files to load - forcefield_name : str, optional, default=None - Apply a named forcefield to the output file using the `foyer` - package, e.g. 'oplsaa'. `Foyer forcefields` - _ - steps : int, optional, default=1000 - Number of energy minimization iterations - scale_bonds : float, optional, default=1 - Scales the bond force constant (1 is completely on) - scale_angles : float, optiona, default=1 - Scales the angle force constant (1 is completely on) - scale_torsions : float, optional, default=1 - Scales the torsional force constants (1 is completely on) - scale_nonbonded : float, optional, default=1 - Scales epsilon (1 is completely on) - constraints : str, optional, default="AllBonds" - Specify constraints on the molecule to minimize, options are: - None, "HBonds", "AllBonds", "HAngles" - - Notes - ----- - Assumes a particular organization for the force groups - (HarmonicBondForce, HarmonicAngleForce, RBTorsionForce, NonBondedForce) - - References - ---------- - [Eastman2013]_ - """ - foyer = import_("foyer") - - to_parmed = self.to_parmed() - ff = foyer.Forcefield(forcefield_files=forcefield_files, name=forcefield_name) - to_parmed = ff.apply(to_parmed) - - import openmm.unit as u - from openmm.app import AllBonds, HAngles, HBonds - from openmm.app.pdbreporter import PDBReporter - from openmm.app.simulation import Simulation - from openmm.openmm import LangevinIntegrator - - if constraints: - if constraints == "AllBonds": - constraints = AllBonds - elif constraints == "HBonds": - constraints = HBonds - elif constraints == "HAngles": - constraints = HAngles - else: - raise ValueError( - f"Provided constraints value of: {constraints}.\n" - f'Expected "HAngles", "AllBonds" "HBonds".' - ) - system = to_parmed.createSystem( - constraints=constraints - ) # Create an OpenMM System - else: - system = to_parmed.createSystem() # Create an OpenMM System - # Create a Langenvin Integrator in OpenMM - integrator = LangevinIntegrator( - 298 * u.kelvin, 1 / u.picosecond, 0.002 * u.picoseconds - ) - # Create Simulation object in OpenMM - simulation = Simulation(to_parmed.topology, system, integrator) - - # Loop through forces in OpenMM System and set parameters - for force in system.getForces(): - if type(force).__name__ == "HarmonicBondForce": - for bond_index in range(force.getNumBonds()): - atom1, atom2, r0, k = force.getBondParameters(bond_index) - force.setBondParameters( - bond_index, atom1, atom2, r0, k * scale_bonds - ) - force.updateParametersInContext(simulation.context) - - elif type(force).__name__ == "HarmonicAngleForce": - for angle_index in range(force.getNumAngles()): - atom1, atom2, atom3, r0, k = force.getAngleParameters(angle_index) - force.setAngleParameters( - angle_index, atom1, atom2, atom3, r0, k * scale_angles - ) - force.updateParametersInContext(simulation.context) - - elif type(force).__name__ == "RBTorsionForce": - for torsion_index in range(force.getNumTorsions()): - ( - atom1, - atom2, - atom3, - atom4, - c0, - c1, - c2, - c3, - c4, - c5, - ) = force.getTorsionParameters(torsion_index) - force.setTorsionParameters( - torsion_index, - atom1, - atom2, - atom3, - atom4, - c0 * scale_torsions, - c1 * scale_torsions, - c2 * scale_torsions, - c3 * scale_torsions, - c4 * scale_torsions, - c5 * scale_torsions, - ) - force.updateParametersInContext(simulation.context) - - elif type(force).__name__ == "NonbondedForce": - for nb_index in range(force.getNumParticles()): - charge, sigma, epsilon = force.getParticleParameters(nb_index) - force.setParticleParameters( - nb_index, charge, sigma, epsilon * scale_nonbonded - ) - force.updateParametersInContext(simulation.context) - - elif type(force).__name__ == "CMMotionRemover": - pass - - else: - warn( - f"OpenMM Force {type(force).__name__} is " - "not currently supported in _energy_minimize_openmm. " - "This Force will not be updated!" - ) - - simulation.context.setPositions(to_parmed.positions) - # Run energy minimization through OpenMM - simulation.minimizeEnergy(maxIterations=steps) - reporter = PDBReporter(os.path.join(tmp_dir, "minimized.pdb"), 1) - reporter.report(simulation, simulation.context.getState(getPositions=True)) - - def _check_openbabel_constraints( - self, - particle_list, - successors_list, - check_if_particle=False, - ): - """Provide routines commonly used to check constraint inputs.""" - for part in particle_list: - if not isinstance(part, Compound): - raise MBuildError(f"{part} is not a Compound.") - if id(part) != id(self) and id(part) not in successors_list: - raise MBuildError(f"{part} is not a member of Compound {self}.") - - if check_if_particle: - if len(part.children) != 0: - raise MBuildError( - f"{part} does not correspond to an individual particle." - ) - - def _energy_minimize_openbabel( - self, - steps=1000, - algorithm="cg", - forcefield="UFF", - constraint_factor=50000.0, - distance_constraints=None, - fixed_compounds=None, - ignore_compounds=None, - ): - """Perform an energy minimization on a Compound. - - Utilizes Open Babel (http://openbabel.org/docs/dev/) to perform an - energy minimization/geometry optimization on a Compound by applying - a generic force field. - - This function is primarily intended to be used on smaller components, - with sizes on the order of 10's to 100's of particles, as the energy - minimization scales poorly with the number of particles. - - Parameters - ---------- - steps : int, optionl, default=1000 - The number of optimization iterations - algorithm : str, optional, default='cg' - The energy minimization algorithm. Valid options are 'steep', - 'cg', and 'md', corresponding to steepest descent, conjugate - gradient, and equilibrium molecular dynamics respectively. - forcefield : str, optional, default='UFF' - The generic force field to apply to the Compound for minimization. - Valid options are 'MMFF94', 'MMFF94s', ''UFF', 'GAFF', 'Ghemical'. - Please refer to the Open Babel documentation - (http://open-babel.readthedocs.io/en/latest/Forcefields/Overview.html) - when considering your choice of force field. - fixed_compounds : Compound, optional, default=None - An individual Compound or list of Compounds that will have their - position fixed during energy minimization. Note, positions are fixed - using a restraining potential and thus may change slightly. - Position fixing will apply to all Particles (i.e., atoms) that exist - in the Compound and to particles in any subsequent sub-Compounds. - By default x,y, and z position is fixed. This can be toggled by instead - passing a list containing the Compound and a list or tuple of bool values - corresponding to x,y and z; e.g., [Compound, (True, True, False)] - will fix the x and y position but allow z to be free. - ignore_compounds: Compound, optional, default=None - An individual compound or list of Compounds whose underlying particles - will have their positions fixed and not interact with other atoms via - the specified force field during the energy minimization process. - Note, a restraining potential is used and thus absolute position may vary - as a result of the energy minimization process. - Interactions of these ignored atoms can be specified by the user, - e.g., by explicitly setting a distance constraint. - distance_constraints: list, optional, default=None - A list containing a pair of Compounds as a tuple or list and - a float value specifying the target distance between the two Compounds, e.g.,: - [(compound1, compound2), distance]. - To specify more than one constraint, pass constraints as a 2D list, e.g.,: - [ [(compound1, compound2), distance1], [(compound3, compound4), distance2] ]. - Note, Compounds specified here must represent individual point particles. - constraint_factor: float, optional, default=50000.0 - Harmonic springs are used to constrain distances and fix atom positions, where - the resulting energy associated with the spring is scaled by the - constraint_factor; the energy of this spring is considering during the minimization. - As such, very small values of the constraint_factor may result in an energy - minimized state that does not adequately restrain the distance/position of atom(s)e. - - - References - ---------- - [OBoyle2011]_ - [OpenBabel]_ - - If using the 'MMFF94' force field please also cite the following: - [Halgren1996a]_ - [Halgren1996b]_ - [Halgren1996c]_ - [Halgren1996d]_ - [Halgren1996e]_ - - If using the 'MMFF94s' force field please cite the above along with: - [Halgren1999]_ - - If using the 'UFF' force field please cite the following: - [Rappe1992]_ - - If using the 'GAFF' force field please cite the following: - [Wang2001]_ - - If using the 'Ghemical' force field please cite the following: - [Hassinen2001]_ - """ - openbabel = import_("openbabel") - for particle in self.particles(): - if particle.element is None: - try: - particle._element = element_from_symbol(particle.name) - except ElementError: - try: - particle._element = element_from_name(particle.name) - except ElementError: - raise MBuildError( - f"No element assigned to {particle}; element could not be" - f"inferred from particle name {particle.name}. Cannot perform" - "an energy minimization." - ) - # Create a dict containing particle id and associated index to speed up looping - particle_idx = { - id(particle): idx for idx, particle in enumerate(self.particles()) - } - - # A list containing all Compounds ids contained in self. Will be used to check if - # compounds refered to in the constrains are actually in the Compound we are minimizing. - successors_list = [id(compound) for compound in self.successors()] - - # initialize constraints - ob_constraints = openbabel.OBFFConstraints() - - if distance_constraints is not None: - # if a user passes single constraint as a 1-D array, - # i.e., [(p1,p2), 2.0] rather than [[(p1,p2), 2.0]], - # just add it to a list so we can use the same looping code - if len(np.array(distance_constraints, dtype=object).shape) == 1: - distance_constraints = [distance_constraints] - - for con_temp in distance_constraints: - p1 = con_temp[0][0] - p2 = con_temp[0][1] - - self._check_openbabel_constraints( - [p1, p2], successors_list, check_if_particle=True - ) - if id(p1) == id(p2): - raise MBuildError( - f"Cannot create a constraint between a Particle and itself: {p1} {p2} ." - ) - - # openbabel indices start at 1 - pid_1 = particle_idx[id(p1)] + 1 - # openbabel indices start at 1 - pid_2 = particle_idx[id(p2)] + 1 - dist = ( - con_temp[1] * 10.0 - ) # obenbabel uses angstroms, not nm, convert to angstroms - - ob_constraints.AddDistanceConstraint(pid_1, pid_2, dist) - - if fixed_compounds is not None: - # if we are just passed a single Compound, wrap it into - # and array so we can just use the same looping code - if isinstance(fixed_compounds, Compound): - fixed_compounds = [fixed_compounds] - - # if fixed_compounds is a 1-d array and it is of length 2, we need to determine whether it is - # a list of two Compounds or if fixed_compounds[1] should correspond to the directions to constrain - if len(np.array(fixed_compounds, dtype=object).shape) == 1: - if len(fixed_compounds) == 2: - if not isinstance(fixed_compounds[1], Compound): - # if it is not a list of two Compounds, make a 2d array so we can use the same looping code - fixed_compounds = [fixed_compounds] - - for fixed_temp in fixed_compounds: - # if an individual entry is a list, validate the input - if isinstance(fixed_temp, list): - if len(fixed_temp) == 2: - msg1 = ( - "Expected tuple or list of length 3 to set" - "which dimensions to fix motion." - ) - assert isinstance(fixed_temp[1], (list, tuple)), msg1 - - msg2 = ( - "Expected tuple or list of length 3 to set" - "which dimensions to fix motion, " - f"{len(fixed_temp[1])} found." - ) - assert len(fixed_temp[1]) == 3, msg2 - - dims = [dim for dim in fixed_temp[1]] - msg3 = ( - "Expected bool values for which directions are fixed." - f"Found instead {dims}." - ) - assert all(isinstance(dim, bool) for dim in dims), msg3 - - p1 = fixed_temp[0] - - # if fixed_compounds is defined as [[Compound],[Compound]], - # fixed_temp will be a list of length 1 - elif len(fixed_temp) == 1: - p1 = fixed_temp[0] - dims = [True, True, True] - - else: - p1 = fixed_temp - dims = [True, True, True] - - all_true = all(dims) - - self._check_openbabel_constraints([p1], successors_list) - - if len(p1.children) == 0: - pid = particle_idx[id(p1)] + 1 # openbabel indices start at 1 - - if all_true: - ob_constraints.AddAtomConstraint(pid) - else: - if dims[0]: - ob_constraints.AddAtomXConstraint(pid) - if dims[1]: - ob_constraints.AddAtomYConstraint(pid) - if dims[2]: - ob_constraints.AddAtomZConstraint(pid) - else: - for particle in p1.particles(): - pid = ( - particle_idx[id(particle)] + 1 - ) # openbabel indices start at 1 - - if all_true: - ob_constraints.AddAtomConstraint(pid) - else: - if dims[0]: - ob_constraints.AddAtomXConstraint(pid) - if dims[1]: - ob_constraints.AddAtomYConstraint(pid) - if dims[2]: - ob_constraints.AddAtomZConstraint(pid) - - if ignore_compounds is not None: - temp1 = np.array(ignore_compounds, dtype=object) - if len(temp1.shape) == 2: - ignore_compounds = list(temp1.reshape(-1)) - - # Since the ignore_compounds can only be passed as a list - # we can check the whole list at once before looping over it - self._check_openbabel_constraints(ignore_compounds, successors_list) - - for ignore in ignore_compounds: - p1 = ignore - if len(p1.children) == 0: - pid = particle_idx[id(p1)] + 1 # openbabel indices start at 1 - ob_constraints.AddIgnore(pid) - - else: - for particle in p1.particles(): - pid = ( - particle_idx[id(particle)] + 1 - ) # openbabel indices start at 1 - ob_constraints.AddIgnore(pid) - - mol = self.to_pybel() - mol = mol.OBMol - - mol.PerceiveBondOrders() - mol.SetAtomTypesPerceived() - - ff = openbabel.OBForceField.FindForceField(forcefield) - if ff is None: - raise MBuildError( - f"Force field '{forcefield}' not supported for energy " - "minimization. Valid force fields are 'MMFF94', " - "'MMFF94s', 'UFF', 'GAFF', and 'Ghemical'." - "" - ) - warn( - "Performing energy minimization using the Open Babel package. " - "Please refer to the documentation to find the appropriate " - f"citations for Open Babel and the {forcefield} force field" - ) - - if ( - distance_constraints is not None - or fixed_compounds is not None - or ignore_compounds is not None - ): - ob_constraints.SetFactor(constraint_factor) - if ff.Setup(mol, ob_constraints) == 0: - raise MBuildError( - "Could not setup forcefield for OpenBabel Optimization." - ) - else: - if ff.Setup(mol) == 0: - raise MBuildError( - "Could not setup forcefield for OpenBabel Optimization." - ) - - if algorithm == "steep": - ff.SteepestDescent(steps) - elif algorithm == "md": - ff.MolecularDynamicsTakeNSteps(steps, 300) - elif algorithm == "cg": - ff.ConjugateGradients(steps) - else: - raise MBuildError( - "Invalid minimization algorithm. Valid options " - "are 'steep', 'cg', and 'md'." - ) - ff.UpdateCoordinates(mol) - - # update the coordinates in the Compound - for i, obatom in enumerate(openbabel.OBMolAtomIter(mol)): - x = obatom.GetX() / 10.0 - y = obatom.GetY() / 10.0 - z = obatom.GetZ() / 10.0 - self[i].pos = np.array([x, y, z]) - def save( self, filename, @@ -3052,7 +2394,7 @@ def from_parmed(self, structure, coords_only=False, infer_hierarchy=True): infer_hierarchy=infer_hierarchy, ) - def to_rdkit(self): + def to_rdkit(self, embed=False): """Create an RDKit RWMol from an mBuild Compound. Returns @@ -3079,7 +2421,7 @@ def to_rdkit(self): See https://www.rdkit.org/docs/GettingStartedInPython.html """ - return conversion.to_rdkit(self) + return conversion.to_rdkit(self, embed=embed) def to_parmed( self, @@ -3349,6 +2691,7 @@ def _clone(self, clone_of=None, root_container=None): newone._periodicity = deepcopy(self._periodicity) newone._charge = deepcopy(self._charge) newone._mass = deepcopy(self._mass) + newone._hoomd_data = {} if hasattr(self, "index"): newone.index = deepcopy(self.index) @@ -3396,16 +2739,36 @@ def _clone_bonds(self, clone_of=None): newone.bond_graph.add_node(clone_of[particle]) for c1, c2, data in self.bonds(return_bond_order=True): try: + try: + bond_order = data["bond_order"] + except KeyError: + bond_order = "unspecified" # bond order is added to the data dictionary as 'bo' - newone.add_bond( - (clone_of[c1], clone_of[c2]), bond_order=data["bond_order"] - ) + newone.add_bond((clone_of[c1], clone_of[c2]), bond_order=bond_order) except KeyError: raise MBuildError( "Cloning failed. Compound contains bonds to " "Particles outside of its containment hierarchy." ) + def _add_sim_data(self, state=None, forces=None, forcefield=None): + if state: + self._hoomd_data["state"] = state + if forces: + self._hoomd_data["forces"] = forces + if forcefield: + self._hoomd_data["forcefield"] = forcefield + + def _get_sim_data(self): + if not self._hoomd_data: + return None, None, None + else: + return ( + self._hoomd_data["state"], + self._hoomd_data["forces"], + self._hoomd_data["forcefield"], + ) + Particle = Compound diff --git a/mbuild/conversion.py b/mbuild/conversion.py index 178129703..2526bbe15 100644 --- a/mbuild/conversion.py +++ b/mbuild/conversion.py @@ -1,11 +1,11 @@ """Module for handling conversions in mBuild.""" +import logging import os import sys from collections import defaultdict from copy import deepcopy from pathlib import Path -from warnings import warn import gmso import numpy as np @@ -23,6 +23,8 @@ from mbuild.formats.json_formats import compound_from_json, compound_to_json from mbuild.utils.io import has_mdtraj, has_openbabel, import_ +logger = logging.getLogger(__name__) + def load( filename_or_object, @@ -180,10 +182,10 @@ def load_object( # 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.") + logger.info("Given object is already an mb.Compound, doing nothing.") return obj else: - warn("Given object is an mb.Compound, adding to the host compound.") + logger.info("Given object is an mb.Compound, adding to the host compound.") compound.add(obj) return compound @@ -443,7 +445,7 @@ def load_file( # text file detected, assume contain smiles string elif extension == ".txt": - warn(".txt file detected, loading as a SMILES string") + logger.info(".txt file detected, loading as a SMILES string") # Fail-safe measure compound = load_pybel_smiles(filename, compound) @@ -461,7 +463,7 @@ def load_file( # Then parmed reader elif backend == "parmed": - warn( + logger.warning( "Using parmed reader. Bonds may be inferred from inter-particle " "distances and standard residue templates. Please check that the " "bonds in mb.Compound are accurate" @@ -573,7 +575,7 @@ def from_parmed( # Convert box information if structure.box is not None: - warn("All angles are assumed to be 90 degrees") + logger.info("All angles are assumed to be 90 degrees") compound.box = Box(structure.box[0:3] / 10) return compound @@ -730,7 +732,7 @@ def from_pybel( resindex_to_cmpd = {} if coords_only: - raise Warning( + logger.error( "coords_only=True is not yet implemented for conversion from pybel" ) @@ -747,7 +749,7 @@ def from_pybel( element = None if use_element: if element is None: - warn( + logger.info( "No element detected for atom at index " f"{atom.idx} with number {atom.atomicnum}, type {atom.type}" ) @@ -796,7 +798,7 @@ def from_pybel( compound.box = box else: if not ignore_box_warn: - warn(f"No unitcell detected for pybel.Molecule {pybel_mol}") + logger.info(f"No unitcell detected for pybel.Molecule {pybel_mol}") return compound @@ -835,6 +837,7 @@ def from_rdkit(rdkit_mol, compound=None, coords_only=False, smiles_seed=0): "to the RDKit error messages for possible fixes. You can also " "install openbabel and use the backend='pybel' instead" ) + # AllChem.EmbedMolecule(mymol, useExpTorsionAnglePrefs=True, useBasicKnowledge=True) AllChem.UFFOptimizeMolecule(mymol) single_mol = mymol.GetConformer(0) # convert from Angstroms to nanometers @@ -997,7 +1000,9 @@ def save( 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}.") + logger.info( + f"System is not charge neutral. Total charge is {compound.charge}." + ) extension = os.path.splitext(filename)[-1] # Keep json stuff with internal mbuild method @@ -1706,7 +1711,7 @@ def to_pybel( return pybelmol -def to_rdkit(compound): +def to_rdkit(compound, embed=False): """Create an RDKit RWMol from an mBuild Compound. Parameters @@ -1762,6 +1767,11 @@ def to_rdkit(compound): rdkit_bond = temp_mol.GetBondBetweenAtoms(*bond_indices) rdkit_bond.SetBondType(bond_order_dict[bond[2]["bond_order"]]) + if embed: + Chem.SanitizeMol(temp_mol) + temp_mol = Chem.AddHs(temp_mol) + Chem.AllChem.EmbedMolecule(temp_mol) + return temp_mol @@ -1782,7 +1792,9 @@ def to_smiles(compound, backend="pybel"): if backend == "pybel": mol = to_pybel(compound) - warn("The bond orders will be guessed using pybelOBMol.PerceviedBondOrders()") + logger.info( + "The bond orders will be guessed using pybelOBMol.PerceviedBondOrders()" + ) mol.OBMol.PerceiveBondOrders() smiles_string = mol.write("smi").replace("\t", " ").split(" ")[0] @@ -2011,7 +2023,7 @@ def _infer_element_from_compound(compound, guessed_elements): except ElementError: try: element = element_from_name(compound.name) - warn_msg = ( + logger.warning_msg = ( f"No element attribute associated with '{compound}'; " f"Guessing it is element '{element}'" ) @@ -2023,7 +2035,7 @@ def _infer_element_from_compound(compound, guessed_elements): "compound name. Setting atomic number to zero." ) if compound.name not in guessed_elements: - warn(warn_msg) + logger.warning(warn_msg) guessed_elements.add(compound.name) return element diff --git a/mbuild/coordinate_transform.py b/mbuild/coordinate_transform.py index b9130d590..971821dd3 100644 --- a/mbuild/coordinate_transform.py +++ b/mbuild/coordinate_transform.py @@ -1,20 +1,19 @@ """Coordinate transformation functions.""" -from warnings import simplefilter, warn +import logging import numpy as np from numpy.linalg import inv, norm, svd -simplefilter("always", DeprecationWarning) - __all__ = [ "force_overlap", "x_axis_transform", "y_axis_transform", "z_axis_transform", - "equivalence_transform", ] +logger = logging.getLogger(__name__) + def force_overlap( move_this, from_positions, to_positions, add_bond=True, reset_labels=False @@ -61,7 +60,7 @@ def force_overlap( if add_bond: if isinstance(from_positions, Port) and isinstance(to_positions, Port): if not from_positions.anchor or not to_positions.anchor: - warn("Attempting to form bond from port that has no anchor") + logger.warning("Attempting to form bond from port that has no anchor") else: from_positions.anchor.parent.add_bond( (from_positions.anchor, to_positions.anchor) @@ -337,59 +336,6 @@ def _create_equivalence_transform(equiv): return T -def equivalence_transform(compound, from_positions, to_positions, add_bond=True): - """Compute an affine transformation. - - Maps the from_positions to the respective to_positions, and applies this - transformation to the compound. - - Parameters - ---------- - compound : mb.Compound - The Compound to be transformed. - from_positions : np.ndarray, shape=(n, 3), dtype=float - Original positions. - to_positions : np.ndarray, shape=(n, 3), dtype=float - New positions. - """ - warn( - "The `equivalence_transform` function is being phased out in favor of" - " `force_overlap`.", - DeprecationWarning, - ) - from mbuild.port import Port - - T = None - if isinstance(from_positions, (list, tuple)) and isinstance( - to_positions, (list, tuple) - ): - equivalence_pairs = zip(from_positions, to_positions) - elif isinstance(from_positions, Port) and isinstance(to_positions, Port): - equivalence_pairs, T = _choose_correct_port(from_positions, to_positions) - from_positions.used = True - to_positions.used = True - else: - equivalence_pairs = [(from_positions, to_positions)] - - if not T: - T = _create_equivalence_transform(equivalence_pairs) - atom_positions = compound.xyz_with_ports - atom_positions = T.apply_to(atom_positions) - compound.xyz_with_ports = atom_positions - - if add_bond: - if isinstance(from_positions, Port) and isinstance(to_positions, Port): - if not from_positions.anchor or not to_positions.anchor: - warn("Attempting to form bond from port that has no anchor") - else: - from_positions.anchor.parent.add_bond( - (from_positions.anchor, to_positions.anchor) - ) - to_positions.anchor.parent.add_bond( - (from_positions.anchor, to_positions.anchor) - ) - - def _choose_correct_port(from_port, to_port): """Chooses the direction when using an equivalence transform on two Ports. diff --git a/mbuild/formats/cassandramcf.py b/mbuild/formats/cassandramcf.py index 526605e52..1fbd3fee1 100644 --- a/mbuild/formats/cassandramcf.py +++ b/mbuild/formats/cassandramcf.py @@ -5,13 +5,14 @@ from __future__ import division -import warnings +import logging from math import sqrt import networkx as nx import parmed as pmd __all__ = ["write_mcf"] +logger = logging.getLogger(__name__) def write_mcf(structure, filename, angle_style, dihedral_style, lj14=None, coul14=None): @@ -102,7 +103,7 @@ def write_mcf(structure, filename, angle_style, dihedral_style, lj14=None, coul1 else: coul14 = 0.0 if len(structure.dihedrals) > 0 or len(structure.rb_torsions) > 0: - warnings.warn( + logger.info( "Unable to infer coulombic 1-4 scaling factor. Setting to " "{:.1f}".format(coul14) ) @@ -118,7 +119,7 @@ def write_mcf(structure, filename, angle_style, dihedral_style, lj14=None, coul1 ] if all([c_eps == 0 for c_eps in combined_eps_list]): lj14 = 0.0 - warnings.warn( + logger.info( "Unable to infer LJ 1-4 scaling factor. Setting to " "{:.1f}".format(lj14) ) @@ -130,7 +131,7 @@ def write_mcf(structure, filename, angle_style, dihedral_style, lj14=None, coul1 break else: lj14 = 0.0 - warnings.warn( + logger.info( "Unable to infer LJ 1-4 scaling factor. Setting to {:.1f}".format( lj14 ) @@ -138,7 +139,7 @@ def write_mcf(structure, filename, angle_style, dihedral_style, lj14=None, coul1 else: lj14 = 0.0 if len(structure.dihedrals) > 0 or len(structure.rb_torsions) > 0: - warnings.warn( + logger.info( "Unable to infer LJ 1-4 scaling factor. Setting to {:.1f}".format( lj14 ) @@ -197,7 +198,7 @@ def _id_rings_fragments(structure): ) if len(structure.bonds) == 0: - warnings.warn("No bonds found. Cassandra will interpet this as a rigid species") + logger.info("No bonds found. Cassandra will interpet this as a rigid species") in_ring = [False] * len(structure.atoms) frag_list = [] frag_conn = [] @@ -272,7 +273,7 @@ def _id_rings_fragments(structure): if len(shared_atoms) == 2: frag_conn.append([i, j]) elif len(shared_atoms) > 2: - warnings.warn( + logger.warning( "Fragments share more than two atoms... something may be " "going awry unless there are fused rings in your system. " "See below for details." @@ -313,32 +314,31 @@ def _write_atom_information(mcf_file, structure, in_ring, IG_CONSTANT_KCAL): n_unique_elements = len(set(elements)) for element in elements: if len(element) > max_element_length: - warnings.warn( - "Element name {} will be shortened to {} characters. Please " - "confirm your final MCF.".format(element, max_element_length) + logger.info( + f"Element name {element} will be shortened to {max_element_length} " + "characters. Please confirm your final MCF." ) elements = [element[:max_element_length] for element in elements] if len(set(elements)) < n_unique_elements: - warnings.warn( + logger.info( "The number of unique elements has been reduced due to shortening " - "the element name to {} characters.".format(max_element_length) + f"the element name to {max_element_length} characters." ) n_unique_types = len(set(types)) for itype in types: if len(itype) > max_atomtype_length: - warnings.warn( - "Type name {} will be shortened to {} characters as {}. Please " - "confirm your final MCF.".format( - itype, max_atomtype_length, itype[-max_atomtype_length:] - ) + logger.info( + f"Type name {itype} will be shortened to {max_atomtype_length} " + f"characters as {itype[-max_atomtype_length:]}. Please " + "confirm your final MCF." ) types = [itype[-max_atomtype_length:] for itype in types] if len(set(types)) < n_unique_types: - warnings.warn( + logger.info( "The number of unique atomtypes has been reduced due to shortening " - "the atomtype name to {} characters.".format(max_atomtype_length) + f"the atomtype name to {max_atomtype_length} characters." ) vdw_type = "LJ" @@ -518,9 +518,7 @@ def _write_dihedral_information(mcf_file, structure, dihedral_style, KCAL_TO_KJ) ] elif dihedral_style.casefold() == "none": - warnings.warn( - "Dihedral style 'none' selected. Ignoring dihedral parameters" - ) + logger.info("Dihedral style 'none' selected. Ignoring dihedral parameters") dihedral_style = dihedral_style.lower() if structure.dihedrals: dihedrals = structure.dihedrals @@ -623,7 +621,7 @@ def _write_fragment_information(mcf_file, structure, frag_list, frag_conn): mcf_file.write("1\n") mcf_file.write("1 2 1 2\n") else: - warnings.warn("More than two atoms present but no fragments identified.") + logger.info("More than two atoms present but no fragments identified.") mcf_file.write("0\n") else: mcf_file.write("{:d}\n".format(len(frag_list))) diff --git a/mbuild/formats/json_formats.py b/mbuild/formats/json_formats.py index be1a2da28..1804ec0b9 100644 --- a/mbuild/formats/json_formats.py +++ b/mbuild/formats/json_formats.py @@ -1,6 +1,7 @@ """JSON format.""" import json +import logging from collections import OrderedDict import ele @@ -9,6 +10,8 @@ from mbuild.bond_graph import BondGraph from mbuild.exceptions import MBuildError +logger = logging.getLogger(__name__) + def compound_from_json(json_file): """Convert the given json file into a Compound. @@ -264,7 +267,6 @@ def _add_bonds(compound_dict, parent, converted_dict): def _perform_sanity_check(json_dict): """Perform Sanity Check on the JSON File.""" - from warnings import warn warning_msg = "This Json was written using {0}, current mbuild version is {1}." this_version = mb.__version__ @@ -285,4 +287,6 @@ def _perform_sanity_check(json_dict): + " Cannot Convert JSON to compound" ) if minor != this_minor: - warn(warning_msg.format(json_mbuild_version, this_version) + " Will Proceed.") + logging.warning( + warning_msg.format(json_mbuild_version, this_version) + " Will Proceed." + ) diff --git a/mbuild/formats/vasp.py b/mbuild/formats/vasp.py index 507c54e4f..a49e939cd 100644 --- a/mbuild/formats/vasp.py +++ b/mbuild/formats/vasp.py @@ -1,6 +1,6 @@ """VASP POSCAR format.""" -import warnings +import logging from itertools import chain import numpy as np @@ -11,6 +11,8 @@ __all__ = ["write_poscar", "read_poscar"] +logger = logging.getLogger(__name__) + def write_poscar(compound, filename, lattice_constant=1.0, coord_style="cartesian"): """Write a VASP POSCAR file from a Compound. @@ -51,7 +53,7 @@ def write_poscar(compound, filename, lattice_constant=1.0, coord_style="cartesia except AttributeError: lattice = compound.get_boundingbox().vectors if coord_style == "direct": - warnings.warn( + logger.info( "'direct' coord_style specified, but compound has no box " "-- using 'cartesian' instead" ) diff --git a/mbuild/lattice.py b/mbuild/lattice.py index 893469bd1..3e9a9a644 100644 --- a/mbuild/lattice.py +++ b/mbuild/lattice.py @@ -1,8 +1,8 @@ """mBuild lattice module for working with crystalline systems.""" import itertools as it +import logging import pathlib -import warnings from collections import defaultdict import numpy as np @@ -14,6 +14,8 @@ __all__ = ["load_cif", "Lattice"] +logger = logging.getLogger(__name__) + def load_cif(file_or_path=None, wrap_coords=False): """Load a CifFile object into an mbuild.Lattice. @@ -654,9 +656,7 @@ def populate(self, compound_dict=None, x=1, y=1, z=1): ) # Raise warnings about assumed elements for element in elementsSet: - warnings.warn( - f"Element assumed from cif file to be {element}.", UserWarning - ) + logger.info(f"Element assumed from cif file to be {element}.") ret_lattice.add(compoundsList) # Create mbuild.box ret_lattice.box = mb.Box(lengths=[a * x, b * y, c * z], angles=self.angles) diff --git a/mbuild/lib/recipes/__init__.py b/mbuild/lib/recipes/__init__.py index 16f148718..84e7607d9 100644 --- a/mbuild/lib/recipes/__init__.py +++ b/mbuild/lib/recipes/__init__.py @@ -3,6 +3,5 @@ from mbuild.lib.recipes.alkane import Alkane from mbuild.lib.recipes.monolayer import Monolayer -from mbuild.lib.recipes.polymer import Polymer from mbuild.lib.recipes.silica_interface import SilicaInterface from mbuild.lib.recipes.tiled_compound import TiledCompound diff --git a/mbuild/lib/recipes/alkane.py b/mbuild/lib/recipes/alkane.py index b0e908403..94d3fb151 100644 --- a/mbuild/lib/recipes/alkane.py +++ b/mbuild/lib/recipes/alkane.py @@ -24,7 +24,7 @@ def __init__(self, n=3, cap_front=True, cap_end=True): if n < 1: raise ValueError("n must be 1 or more") super(Alkane, self).__init__() - from mbuild.lib.recipes import Polymer + from mbuild import Polymer # Handle the case of Methane and Ethane separately if n < 3: diff --git a/mbuild/lib/recipes/monolayer.py b/mbuild/lib/recipes/monolayer.py index 97dfa7d23..c2b20a2b7 100644 --- a/mbuild/lib/recipes/monolayer.py +++ b/mbuild/lib/recipes/monolayer.py @@ -1,7 +1,7 @@ """mBuild monolayer recipe.""" +import logging from copy import deepcopy -from warnings import warn import numpy as np @@ -9,6 +9,8 @@ __all__ = ["Monolayer"] +logger = logging.getLogger(__name__) + class Monolayer(mb.Compound): """A general monolayer recipe. @@ -75,7 +77,7 @@ def __init__( # Create sub-pattern for this chain type subpattern = deepcopy(pattern) n_points = int(round(fraction * n_chains)) - warn("\n Adding {} of chain {}".format(n_points, chain)) + logger.info(f"\n Adding {n_points} of chain {chain}") pick = np.random.choice( subpattern.points.shape[0], n_points, replace=False ) @@ -98,10 +100,10 @@ def __init__( self.add(attached_chains) else: - warn("\n No fractions provided. Assuming a single chain type.") + logger.info("\n No fractions provided. Assuming a single chain type.") # Attach final chain type. Remaining sites get a backfill. - warn("\n Adding {} of chain {}".format(len(pattern), chains[-1])) + logger.info(f"\n Adding {len(pattern)} of chain {chains[-1]}") attached_chains, backfills = pattern.apply_to_compound( guest=chains[-1], host=self["tiled_surface"], backfill=backfill, **kwargs ) diff --git a/mbuild/packing.py b/mbuild/packing.py index 001f4fcbd..9940588e2 100644 --- a/mbuild/packing.py +++ b/mbuild/packing.py @@ -3,11 +3,11 @@ http://leandro.iqm.unicamp.br/m3g/packmol/home.shtml """ +import logging import os import shutil import sys import tempfile -import warnings from itertools import zip_longest from subprocess import PIPE, Popen @@ -20,6 +20,8 @@ __all__ = ["fill_box", "fill_region", "fill_sphere", "solvate"] +logger = logging.getLogger(__name__) + PACKMOL = shutil.which("packmol") PACKMOL_HEADER = """ tolerance {0:.16f} @@ -91,7 +93,7 @@ def check_packmol_args(custom_args): "See https://m3g.github.io/packmol/userguide.shtml#run" ) if key in default_args: - warnings.warn( + logger.info( f"The PACKMOL argument {key} was passed to `packmol_args`, " "but should be set using the corresponding function parameters. " "The value passed to the function will be used. " @@ -1062,7 +1064,7 @@ def _validate_mass(compound, n_compounds): "on how mass is handled." ) if found_zero_mass: - warnings.warn( + logger.info( "Some of the compounds or subcompounds in `compound` " "have a mass of zero/None. This may have an effect on " "density calculations" @@ -1199,7 +1201,7 @@ def _run_packmol(input_text, filled_xyz, temp_file, packmol_file): out, err = proc.communicate() if "WITHOUT PERFECT PACKING" in out: - warnings.warn( + logger.warning( "Packmol finished with imperfect packing. Using the .xyz_FORCED " "file instead. This may not be a sufficient packing result." ) diff --git a/mbuild/path/__init__.py b/mbuild/path/__init__.py new file mode 100644 index 000000000..06c691f56 --- /dev/null +++ b/mbuild/path/__init__.py @@ -0,0 +1,3 @@ +# ruff: noqa: F401 +# ruff: noqa: F403 +from .path import HardSphereRandomWalk, Lamellar, Path diff --git a/mbuild/path/bias.py b/mbuild/path/bias.py new file mode 100644 index 000000000..20697e5fb --- /dev/null +++ b/mbuild/path/bias.py @@ -0,0 +1,127 @@ +"""Contains biases that can be included in mbuild.path.HardSphereRandomWalk.""" + +import numpy as np + +from mbuild.path.path_utils import target_sq_distances + + +class Bias: + def __init__(self, random_walk, new_coordinates): + # Extract any needed information for all sub-classes + # Complete system coords needed for density and direction biases + self.system_coordinates = random_walk.coordinates + self.N = len(self.system_coordinates) + # Site types needed for TargetType and AvoidType + self.site_types = [ + attrs["name"] for node, attrs in random_walk.bond_graph.nodes(data=True) + ] + # Current step count needed for TargetPath + self.count = random_walk.count + self.new_coordinates = new_coordinates + + def __call__(self): + raise NotImplementedError + + +class TargetCoordinate(Bias): + """Bias next-moves so that ones moving closer to a final coordinate are more likely to be accepted.""" + + def __init__(self, target_coordinate, weight, system_coordinates, new_coordinates): + self.target_coordinate = target_coordinate + super(TargetCoordinate, self).__init__(system_coordinates, new_coordinates) + + def __call__(self): + sq_distances = target_sq_distances(self.target_coordinate, self.new_coordinates) + sort_idx = np.argsort(sq_distances) + return self.new_points[sort_idx] + + +class AvoidCoordinate(Bias): + """Bias next-moves so that ones moving further from a specific coordinate are more likely to be accepted.""" + + def __init__(self, avoid_coordinate, weight, system_coordinates, new_coordinates): + self.avoid_coordinate = avoid_coordinate + super(AvoidCoordinate, self).__init__(system_coordinates, new_coordinates) + + def __call__(self): + sq_distances = target_sq_distances(self.avoid_coordinate, self.new_coordinates) + # Sort in descending order, largest to smallest + sort_idx = np.argsort(sq_distances)[::-1] + return self.new_points[sort_idx] + + +class TargetType(Bias): + """Bias next-moves so that ones moving towards a specific site type are more likely to be accepted.""" + + def __init__(self, site_type, weight, system_coordinates, new_coordinates): + self.site_type = site_type + super(TargetType, self).__init__(system_coordinates, new_coordinates) + + def __call__(self): + pass + + +class AvoidType(Bias): + """Bias next-moves so that ones moving away from a specific site type are more likely to be accepted.""" + + def __init__(self, site_type, weight, system_coordinates, new_coordinates): + self.site_type = site_type + super(AvoidType, self).__init__(system_coordinates, new_coordinates) + + def __call__(self): + pass + + +class TargetEdge(Bias): + """Bias next-moves so that ones moving towards a surface are more likely to be accepted.""" + + def __init__(self, weight, system_coordinates, volume_constraint, new_coordinates): + self.volume_constraint = volume_constraint + super(TargetEdge, self).__init__(system_coordinates, new_coordinates) + + def __call__(self): + pass + + +class AvoidEdge(Bias): + """Bias next-moves so that ones away from a surface are more likely to be accepted.""" + + def __init__(self, weight, system_coordinates, volume_constraint, new_coordinates): + self.volume_constraint = volume_constraint + super(AvoidEdge, self).__init__(system_coordinates, new_coordinates) + + def __call__(self): + pass + + +class TargetDirection(Bias): + """Bias next-moves so that ones moving along a certain direction are more likely to be accepted.""" + + def __init__(self, direction, weight, system_coordinates, new_coordinates): + self.direction = direction + super(TargetDirection, self).__init__(system_coordinates, new_coordinates) + + def __call__(self): + pass + + +class AvoidDirection(Bias): + """Bias next-moves so that ones not moving along a certain direction are more likely to be accepted.""" + + def __init__(self, direction, weight, system_coordinates, new_coordinates): + self.direction = direction + super(AvoidDirection, self).__init__(system_coordinates, new_coordinates) + + def __call__(self): + pass + + +class TargetPath(Bias): + """Bias next-moves so that ones following a pre-defined path are more likely to be accepted.""" + + def __init__(self, target_path, weight, system_coordinates, new_coordinates): + self.target_path = target_path + super(TargetPath, self).__init__(system_coordinates, new_coordinates) + + def __call__(self): + pass diff --git a/mbuild/path/path.py b/mbuild/path/path.py new file mode 100644 index 000000000..2490267f0 --- /dev/null +++ b/mbuild/path/path.py @@ -0,0 +1,825 @@ +"""Classes to generate intra-molecular paths and configurations.""" + +import math +from abc import abstractmethod +from copy import deepcopy + +import networkx as nx +import numpy as np +from scipy.interpolate import interp1d + +from mbuild import Compound +from mbuild.path.path_utils import check_path, random_coordinate + + +class Path: + def __init__(self, N=None, coordinates=None, bond_graph=None): + self.bond_graph = bond_graph + # Only N is defined, make empty coordinates array with size N + # Use case: Random walks + if N is not None and coordinates is None: + self.N = N + self.coordinates = np.zeros((N, 3)) + # Only coordinates is defined, set N from length + # Use case: class method from_coordinates + elif coordinates is not None and N is None: + self.N = len(coordinates) + self.coordinates = coordinates + # Neither is defined, use list for coordinates + # Use case: Lamellar - Don't know N initially + elif N is None and coordinates is None: + self.coordinates = [] + else: + raise ValueError("Specify either one of N and coordinates, or neither.") + self.__post_init__() + + def __post_init__(self): + """Needed for CodeQL in order to call abstract method inside __init__()""" + self.generate() + if self.N is None: + self.N = len(self.coordinates) + + @classmethod + def from_coordinates(cls, coordinates, bond_graph=None): + return cls(coordinates=coordinates, bond_graph=bond_graph, N=None) + + def _extend_coordinates(self, N): + zeros = np.zeros((N, 3)) + new_array = np.concatenate(self.coordinates, zeros) + self.coordinates = new_array + self.N += N + + def add_edge(self, u, v): + """Add an edge to the Path's bond graph. + This sets attributes for the edge which include: bond direction, bond length and bond type + u -> v corresponds to previous site and current site where the bond direction is calculated as v - u. + """ + bond_vec = self.coordinates[v] - self.coordinates[u] + bond_length = np.linalg.norm(bond_vec) + bond_vec /= bond_length + # Get node names from previous step to current step + u_name = self.bond_graph.nodes[u]["name"] + v_name = self.bond_graph.nodes[v]["name"] + self.bond_graph.add_edge( + u_of_edge=u, + v_of_edge=v, + direction=bond_vec.tolist(), + length=float(bond_length), + bond_type=(u_name, v_name), + ) + + def get_bonded_sites(self): + """Get all bonded pairs and their bond-vector orientations.""" + pass + + def get_coordinates(self): + if isinstance(self.coordinates, list): + return np.array(self.coordinates) + else: + return self.coordinates + + @abstractmethod + def generate(self): + """Abstract class for running a Path generation algorithm. + Sub-classes that inherit from Path should implement their + own method under genrate. + """ + pass + + def to_compound(self): + """Convert a path and its bond graph to an mBuild Compound.""" + compound = Compound() + for node_id, attrs in self.bond_graph.nodes(data=True): + compound.add(Compound(name=attrs["name"], pos=attrs["xyz"])) + compound.set_bond_graph(self.bond_graph) + return compound + + def apply_mapping(self): + # TODO: Finish, add logic to align orientation with path site pos and bond graph + """Mapping other compounds onto a Path's coordinates + + mapping = {"A": "c1ccccc1C=C", "B": "C=CC=C"} + + for bond in bond graph edges: + site u: add compound + rotate site u head-tail vec to align with bond direction + set orientation and separation for site u head port + rotate site v tail-head vec to align with bond direction + set orientation and separation for site v tail port + + """ + pass + + def _path_history(self): + """Maybe this is a method that can be used optionally. + We could add a save_history parameter to __init__. + Depending on the approach, saving histories might add additional + computation time and resources. + Might be useful for more complicated random walks/branching algorithms + """ + pass + + +class HardSphereRandomWalk(Path): + def __init__( + self, + N, + bond_length, + radius, + min_angle, + max_angle, + bead_name="_A", + volume_constraint=None, + start_from_path=None, + start_from_path_index=None, + attach_paths=False, + initial_point=None, + include_compound=None, + max_attempts=1e5, + seed=42, + trial_batch_size=20, + tolerance=1e-5, + ): + """Generates coordinates from a self avoiding random walk using + fixed bond lengths, hard spheres, and minimum and maximum angles + formed by 3 consecutive points. + + Possible angle values are sampled uniformly between min_angle + and max_angle between the new site and the two previous sites. + + Parameters: + ----------- + bond_length : float, required + Fixed bond length between 2 coordinates. + radius : float, required + Radius of sites used in checking for overlaps. + min_angle : float, required + Minimum angle (radians) used when randomly selecting angle + for the next step. + max_angle : float, required + Maximum angle (radians) used when randomly selecting angle + for the next step. + seed : int, default = 42 + Random seed + trial_batch_size : int, default = 5 + The number of trial moves to attempt in parallel for each step. + Using larger values can improve success rates for more dense + random walks. + max_attempts : int, defualt = 1e5 + The maximum number of trial moves to attempt before quiting. + start_from_path : mbuild.path.Path, optional + An instance of a previous Path to start the random walk from. + start_from_path_index : int, optional + The index of `start_from_path` to use as the initial point + for the random walk. + tolerance : float, default = 1e-4 + Tolerance used for rounding and checkig for overlaps. + + Notes + ----- + Each next-move can be attempted in batches, set by the ``trial_batch_size`` + parameter. The batch size moves do not count towards the maximum allowed + attempts. For example, 1 random walk with a trail batch size of 20 counts at + only one attempted move. Larger values of ``trial_batch_size`` may help + highly constrained walks finish, but may hurt performance. + + You can start a random walk from a previously created path with the + ``start_from_path`` and ``start_from_path_index`` parameters. For example, + create a ``Lamellar`` path and run a random walk from its last index to begin + generating a semi-crystalline-like structure. Or, string together multiple + HardSphereRandomWalk paths where the final result of each is passed into the + next random walk. + """ + if initial_point is not None: + self.initial_point = np.asarray(initial_point) + else: + self.initial_point = None + self.include_compound = include_compound + self.bond_length = bond_length + self.radius = radius + self.min_angle = min_angle + self.max_angle = max_angle + self.seed = seed + self.bead_name = bead_name + self.volume_constraint = volume_constraint + self.tolerance = tolerance + self.trial_batch_size = int(trial_batch_size) + self.max_attempts = int(max_attempts) + self.attempts = 0 + self.start_from_path_index = start_from_path_index + self.start_from_path = start_from_path + self.attach_paths = attach_paths + self._particle_pairs = {} + + # This random walk is including a previous path + if start_from_path: + coordinates = np.concatenate( + ( + start_from_path.get_coordinates().astype(np.float32), + np.zeros((N, 3), dtype=np.float32), + ), + axis=0, + ) + self.count = len(start_from_path.coordinates) + N = None + bond_graph = deepcopy(start_from_path.bond_graph) + if start_from_path_index is not None and start_from_path_index < 0: + self.start_from_path_index = self.count + start_from_path_index + else: # Not starting from another path + bond_graph = nx.Graph() + coordinates = np.zeros((N, 3), dtype=np.float32) + N = None + self.count = 0 + self.start_index = 0 + # Need this for error message about reaching max tries + self._init_count = self.count + # Select methods to use for random walk + # Hard-coded for, possible to make other RW methods and pass them in + self.next_step = random_coordinate + self.check_path = check_path + # Create RNG state. + self.rng = np.random.default_rng(seed) + super(HardSphereRandomWalk, self).__init__( + coordinates=coordinates, N=N, bond_graph=bond_graph + ) + + def generate(self): + initial_xyz = self._initial_points() + # Set the first coordinate + self.coordinates[self.count] = initial_xyz + self.bond_graph.add_node( + self.count, + name=self.bead_name, + xyz=self.coordinates[self.count], + ) + if not self.start_from_path and not self.volume_constraint: + # If no volume constraint, then first move is always accepted + phi = self.rng.uniform(0, 2 * np.pi) + theta = self.rng.uniform(0, np.pi) + next_pos = np.array( + [ + self.bond_length * np.sin(theta) * np.cos(phi), + self.bond_length * np.sin(theta) * np.sin(phi), + self.bond_length * np.cos(theta), + ] + ) + self.coordinates[1] = self.coordinates[0] + next_pos + self.count += 1 + self.bond_graph.add_node( + self.count, + name=self.bead_name, + xyz=self.coordinates[self.count], + ) + self.add_edge(u=self.count - 1, v=self.count) + # Not starting from another path, but have a volume constraint + # Possible for second point to be out-of-bounds + elif not self.start_from_path and self.volume_constraint: + self.coordinates[0] = initial_xyz + next_point_found = False + while not next_point_found: + phi = self.rng.uniform(0, 2 * np.pi) + theta = self.rng.uniform(0, np.pi) + xyz = np.array( + [ + self.bond_length * np.sin(theta) * np.cos(phi), + self.bond_length * np.sin(theta) * np.sin(phi), + self.bond_length * np.cos(theta), + ] + ) + is_inside_mask = self.volume_constraint.is_inside( + points=np.array([xyz]), buffer=self.radius + ) + if np.all(is_inside_mask): + self.coordinates[1] = self.coordinates[0] + xyz + self.count += 1 + self.bond_graph.add_node( + self.count, + name=self.bead_name, + xyz=self.coordinates[self.count], + ) + self.add_edge(u=self.count - 1, v=self.count) + self.attempts += 1 + next_point_found = True + # 2nd point failed, continue while loop + self.attempts += 1 + + # Starting random walk from a previous set of coordinates (another path) + # This point was accepted in self._initial_point with these conditions + # If attach_paths, then add edge between first node of this path and node of last path + else: + self.bond_graph.add_node(self.count, name=self.bead_name, xyz=initial_xyz) + if self.attach_paths: + self.add_edge(u=self.start_from_path_index, v=self.count) + self.attempts += 1 + + # Initial conditions set (points 1 and 2), now start RW with min/max angles + while self.count < self.N - 1: + # Choosing angles and vectors is the only random part + # Get a batch of these once, pass them into numba functions + batch_angles, batch_vectors = self._generate_random_trials() + new_xyzs = self.next_step( + pos1=self.coordinates[self.count], + pos2=self.coordinates[self.count - 1], + bond_length=self.bond_length, + thetas=batch_angles, + r_vectors=batch_vectors, + batch_size=self.trial_batch_size, + ) + if self.volume_constraint: + is_inside_mask = self.volume_constraint.is_inside( + points=new_xyzs, buffer=self.radius + ) + new_xyzs = new_xyzs[is_inside_mask] + for xyz in new_xyzs: + if self.include_compound: + # Include compound particle coordinates in check for overlaps + existing_points = np.concatenate( + (self.coordinates[: self.count + 1], self.include_compound.xyz) + ) + else: + existing_points = self.coordinates[: self.count + 1] + # Now check for overlaps for each trial point + # Stop after the first success + if self.check_path( + existing_points=existing_points, + new_point=xyz, + radius=self.radius, + tolerance=self.tolerance, + ): + self.coordinates[self.count + 1] = xyz + self.count += 1 + self.bond_graph.add_node(self.count, name=self.bead_name, xyz=xyz) + self.add_edge(u=self.count - 1, v=self.count) + break + self.attempts += 1 + + if self.attempts == self.max_attempts and self.count < self.N: + raise RuntimeError( + "The maximum number attempts allowed have passed, and only ", + f"{self.count - self._init_count} sucsessful attempts were completed.", + "Try changing the parameters or seed and running again.", + ) + + def _generate_random_trials(self): + """Generate a batch of random angles and vectors using the RNG state.""" + # Batch of angles used to determine next positions + thetas = self.rng.uniform( + self.min_angle, self.max_angle, size=self.trial_batch_size + ).astype(np.float32) + # Batch of random vectors and center around origin (0,0,0) + r = self.rng.uniform(-0.5, 0.5, size=(self.trial_batch_size, 3)).astype( + np.float32 + ) + return thetas, r + + def _initial_points(self): + """Choosing first and second points depends on the parameters passed in.""" + # Use manually specified initial point + if self.initial_point is not None: + return self.initial_point + + # Random initial point, no volume constraint: Bounds set by radius and N steps + elif not any( + [self.volume_constraint, self.initial_point, self.start_from_path] + ): + max_dist = (self.N * self.radius) - self.radius + xyz = self.rng.uniform(low=-max_dist, high=max_dist, size=3) + return xyz + + # Random point inside volume constraint, not starting from another path + elif self.volume_constraint and not self.start_from_path_index: + xyz = self.rng.uniform( + low=self.volume_constraint.mins + self.radius, + high=self.volume_constraint.maxs - self.radius, + size=3, + ) + return xyz + + # Starting from another path, run Monte Carlo + # Accepted next move is first point of this random walk + elif self.start_from_path and self.start_from_path_index is not None: + # TODO: handle start_from_path index of negative values + # Set to the corresponding actual index value of the last path + if self.start_from_path_index == 0: + pos2_coord = 1 # Use the second (1) point of the last path for angles + else: # use the site previous to start_from_path_index for angles + pos2_coord = self.start_from_path_index - 1 + started_next_path = False + while not started_next_path: + batch_angles, batch_vectors = self._generate_random_trials() + new_xyzs = self.next_step( + pos1=self.start_from_path.get_coordinates()[ + self.start_from_path_index + ], + pos2=self.start_from_path.get_coordinates()[pos2_coord], + bond_length=self.bond_length, + thetas=batch_angles, + r_vectors=batch_vectors, + batch_size=self.trial_batch_size, + ) + if self.volume_constraint: + is_inside_mask = self.volume_constraint.is_inside( + points=new_xyzs, buffer=self.radius + ) + new_xyzs = new_xyzs[is_inside_mask] + for xyz in new_xyzs: + if self.include_compound: + existing_points = np.concatenate( + ( + self.coordinates[: self.count + 1], + self.include_compound.xyz, + ) + ) + else: + existing_points = self.coordinates[: self.count + 1] + if self.check_path( + existing_points=existing_points, + new_point=xyz, + radius=self.radius, + tolerance=self.tolerance, + ): + return xyz + self.attempts += 1 + + if self.attempts == self.max_attempts and self.count < self.N: + raise RuntimeError( + "The maximum number attempts allowed have passed, and only ", + f"{self.count - self._init_count} sucsessful attempts were completed.", + "Try changing the parameters or seed and running again.", + ) + + +class Lamellar(Path): + """Generate a 2-D or 3-D lamellar-like path. + + Parameters + ---------- + spacing : float (nm), required + The distance between two adjacent sites in the path. + num_layers : int, required + The number of times the lamellar path curves around + creating another layer. + layer_separation : float (nm), required + The distance between any two layers. + num_stacks : int, required + The number of times to repeat each layer in the Z direction. + Setting this to 1 creates a single, 2D lamellar-like path. + stack_separation : float (nm), required + The distance between two stacked layers. + """ + + def __init__( + self, + num_layers, + layer_separation, + layer_length, + bond_length, + num_stacks=1, + stack_separation=None, + bond_graph=None, + ): + self.num_layers = num_layers + self.layer_separation = layer_separation + self.layer_length = layer_length + self.bond_length = bond_length + self.num_stacks = num_stacks + self.stack_separation = stack_separation + super(Lamellar, self).__init__(N=None, bond_graph=bond_graph) + + def generate(self): + layer_spacing = np.arange(0, self.layer_length, self.bond_length) + # Info needed for generating coords of the arc curves between layers + r = self.layer_separation / 2 + arc_length = r * np.pi + arc_num_points = math.floor(arc_length / self.bond_length) + arc_angle = np.pi / (arc_num_points + 1) # incremental angle + arc_angles = np.linspace(arc_angle, np.pi, arc_num_points, endpoint=False) + for i in range(self.num_layers): + if i % 2 == 0: # Even layer; build from left to right + layer = [ + np.array([self.layer_separation * i, y, 0]) for y in layer_spacing + ] + # Mid-point between this and next layer; use to get curve coords. + origin = layer[-1] + np.array([r, 0, 0]) + arc = [ + origin + np.array([-np.cos(theta), np.sin(theta), 0]) * r + for theta in arc_angles + ] + else: # Odd layer; build from right to left + layer = [ + np.array([self.layer_separation * i, y, 0]) + for y in layer_spacing[::-1] + ] + # Mid-point between this and next layer; use to get curve coords. + origin = layer[-1] + np.array([r, 0, 0]) + arc = [ + origin + np.array([-np.cos(theta), -np.sin(theta), 0]) * r + for theta in arc_angles + ] + if i != self.num_layers - 1: + self.coordinates.extend(layer + arc) + else: # Last layer, don't include another arc set of coordinates + self.coordinates.extend(layer) + if self.num_stacks > 1: + first_stack_coordinates = np.copy(np.array(self.coordinates)) + # Get info for curves between stacked layers + r = self.stack_separation / 2 + arc_length = r * np.pi + arc_num_points = math.floor(arc_length / self.bond_length) + arc_angle = np.pi / (arc_num_points + 1) # incremental angle + arc_angles = np.linspace(arc_angle, np.pi, arc_num_points, endpoint=False) + if self.num_layers % 2 == 0: + odd_stack_mult = -1 + else: + odd_stack_mult = 1 + for i in range(1, self.num_stacks): + if i % 2 != 0: # Odd stack + this_stack = np.copy(first_stack_coordinates[::-1]) + np.array( + [0, 0, self.stack_separation * i] + ) + origin = self.coordinates[-1] + np.array([0, 0, r]) + arc = [ + origin + + np.array([0, odd_stack_mult * np.sin(theta), np.cos(theta)]) + * r + for theta in arc_angles + ] + self.coordinates.extend(arc[::-1]) + self.coordinates.extend(list(this_stack)) + elif i % 2 == 0: # Even stack + this_stack = np.copy(first_stack_coordinates) + np.array( + [0, 0, self.stack_separation * i] + ) + origin = self.coordinates[-1] + np.array([0, 0, r]) + arc = [ + origin + np.array([0, -np.sin(theta), np.cos(theta)]) * r + for theta in arc_angles + ] + self.coordinates.extend(arc[::-1]) + self.coordinates.extend(list(this_stack)) + + +class StraightLine(Path): + """Generates a set of coordinates in a straight line along a given axis. + + Parameters + ---------- + spacing : float, required + The distance between sites along the path. + N : int, required + The number of sites in the path. + direction : array-like (1,3), default = (1,0,0) + The direction to align the straight path along. + """ + + def __init__(self, spacing, N, direction=(1, 0, 0), bond_graph=None): + self.spacing = spacing + self.direction = np.asarray(direction) + super(StraightLine, self).__init__(N=N, bond_graph=bond_graph) + + def generate(self): + self.coordinates = np.array( + [np.zeros(3) + i * self.spacing * self.direction for i in range(self.N)] + ) + + +class Cyclic(Path): + """Generates a set of coordinates evenly spaced along a circle. + + Parameters + ---------- + spacing : float, optional + Distance between sites along the path. + N : int, optional + Number of sites in the cyclic path. + radius : float, optional + The radius (nm) of the cyclic path. + + Notes + ----- + Only two of spacing, N and radius can be defined, as the third + is determined by the other two. + + If using this Path to build a cyclic polymer, be sure to + set ``bond_head_tail = True`` in ``mbuild.polymer.Polymer.build_from_path`` + """ + + def __init__(self, spacing=None, N=None, radius=None, bond_graph=None): + self.spacing = spacing + self.radius = radius + n_params = sum(1 for i in (spacing, N, radius) if i is not None) + if n_params != 2: + raise ValueError("You must specify only 2 of spacing, N and radius.") + super(Cyclic, self).__init__(N=N, bond_graph=bond_graph) + + def generate(self): + if self.spacing and self.N: + self.radius = (self.N * self.spacing) / (2 * np.pi) + elif self.radius and self.spacing: + self.N = int((2 * np.pi * self.radius) / self.spacing) + else: + self.spacing = (2 * np.pi) / self.N + + angles = np.arange(0, 2 * np.pi, (2 * np.pi) / self.N) + self.coordinates = np.array( + [(np.cos(a) * self.radius, np.sin(a) * self.radius, 0) for a in angles] + ) + + +class Knot(Path): + """Generate a knot path. + + Parameters + ---------- + spacing : float (nm) + The spacing between sites along the path. + N : int + The number of total sites in the path. + m : int in [3, 4, 5] + The number of crossings in the knot. + 3 gives the trefoil knot, 4 gives the figure 8 knot and 5 gives the cinquefoil knot. + Only values of 3, 4 and 5 are currently supported. + bond_graph : networkx.graph + Sets the bond graph between sites. + """ + + def __init__(self, spacing, N, m, bond_graph=None): + self.spacing = spacing + self.m = m + super(Knot, self).__init__(N=N, bond_graph=bond_graph) + + def generate(self): + # Generate dense sites first, sample actual ones later from spacing + # Prevents spacing between sites changing with curvature + t_dense = np.linspace(0, 2 * np.pi, 5000) + # Base (unscaled) curve + if self.m == 3: # Trefoil knot https://en.wikipedia.org/wiki/Trefoil_knot + R, r = 1.0, 0.3 + x = (R + r * np.cos(3 * t_dense)) * np.cos(2 * t_dense) + y = (R + r * np.cos(3 * t_dense)) * np.sin(2 * t_dense) + z = r * np.sin(3 * t_dense) + elif ( + self.m == 4 + ): # Figure-eight https://en.wikipedia.org/wiki/Figure-eight_knot_(mathematics) + x = (2 + np.cos(2 * t_dense)) * np.cos(3 * t_dense) + y = (2 + np.cos(2 * t_dense)) * np.sin(3 * t_dense) + z = np.sin(4 * t_dense) + elif ( + self.m == 5 + ): # Cinquefoil knot https://en.wikipedia.org/wiki/Cinquefoil_knot + R, r = 1.0, 0.3 + x = (R + r * np.cos(5 * t_dense)) * np.cos(2 * t_dense) + y = (R + r * np.cos(5 * t_dense)) * np.sin(2 * t_dense) + z = r * np.sin(5 * t_dense) + else: + raise ValueError("Only m=3, m=4 and m=5 are currently supported.") + # Compute arc length of a base curve + coords_dense = np.stack((x, y, z), axis=1) + deltas = np.diff(coords_dense, axis=0) + dists = np.linalg.norm(deltas, axis=1) + arc_lengths = np.concatenate([[0], np.cumsum(dists)]) + base_length = arc_lengths[-1] + L_target = (self.N - 1) * self.spacing + # Scale to match target contour length + scale = L_target / base_length + coords_dense *= scale + arc_lengths *= scale + # Resample uniformly along arc length based on target separation and N sites + desired_arcs = np.linspace(0, L_target, self.N, endpoint=False) + x_interp = interp1d(arc_lengths, coords_dense[:, 0])(desired_arcs) + y_interp = interp1d(arc_lengths, coords_dense[:, 1])(desired_arcs) + z_interp = interp1d(arc_lengths, coords_dense[:, 2])(desired_arcs) + self.coordinates = np.stack((x_interp, y_interp, z_interp), axis=1) + + +class Helix(Path): + def __init__( + self, N, radius, rise, twist, right_handed=True, bottom_up=True, bond_graph=None + ): + """Generate helical path. + + Parameters: + ----------- + radius : float, required + Radius of the helix (nm) + rise : float, required + Rise per site on path (nm) + twist : float, required + Twist per site in path (degrees) + right_handed : bool, default True + Set the handedness of the helical twist + Set to False for a left handed twist + bottom_up : bool, default True + If True, the twist is in the positive Z direction + If False, the twist is in the negative Z direction + + Notes: + ------ + To create a double helix pair (e.g., DNA) create two paths + with opposite values for right_handed and bottom_up. + """ + self.radius = radius + self.rise = rise + self.twist = twist + self.right_handed = right_handed + self.bottom_up = bottom_up + super(Helix, self).__init__(N=N, bond_graph=bond_graph) + + def generate(self): + indices = reversed(range(self.N)) if not self.bottom_up else range(self.N) + for i in indices: + angle = np.deg2rad(i * self.twist) + if not self.right_handed: + angle *= -1 + x = self.radius * np.cos(angle) + y = self.radius * np.sin(angle) + z = i * self.rise if self.bottom_up else -i * self.rise + self.coordinates[i] = (x, y, z) + + +class Spiral2D(Path): + def __init__(self, N, a, b, spacing, bond_graph=None): + """Generate a 2D spiral path in the XY plane. + + Parameters + ---------- + N: int, required + Number of sites in the path + a : float, required + The initial radius (nm) + b : float, required + Determines how fast radius grows per angle increment (r = a + bθ) + spacing : float, required + Distance between adjacent sites (nm) + """ + self.a = a + self.b = b + self.spacing = spacing + super().__init__(N=N, bond_graph=bond_graph) + + def generate(self): + theta = 0.0 + for i in range(self.N): + r = self.a + self.b * theta + x = r * np.cos(theta) + y = r * np.sin(theta) + z = 0.0 + self.coordinates[i] = (x, y, z) + # Estimate next angle increment based on arc length + ds_dtheta = np.sqrt((r) ** 2 + self.b**2) + dtheta = self.spacing / ds_dtheta + theta += dtheta + + +class ZigZag(Path): + def __init__( + self, + N, + spacing=1.0, + angle_deg=120.0, + sites_per_segment=1, + plane="xy", + bond_graph=None, + ): + self.spacing = spacing + self.angle_deg = angle_deg + self.sites_per_segment = sites_per_segment + self.plane = plane + if N % sites_per_segment != 0: + raise ValueError("N must be evenly divisible by sites_per_segment") + super(ZigZag, self).__init__(N=N, bond_graph=bond_graph) + + def generate(self): + angle_rad = np.deg2rad(self.angle_deg) + direction = np.array([1.0, 0.0]) + position = np.zeros(2) + coords = [] + step_count = 0 + segment_count = 0 + + for i in range(self.N): + coords.append(position.copy()) + position += self.spacing * direction + step_count += 1 + + # Rotate + if step_count == self.sites_per_segment: + step_count = 0 + sign = -1 if segment_count % 2 == 0 else 1 + rot_matrix = np.array( + [ + [np.cos(sign * angle_rad), -np.sin(sign * angle_rad)], + [np.sin(sign * angle_rad), np.cos(sign * angle_rad)], + ] + ) + direction = rot_matrix @ direction + segment_count += 1 + + # Map into 3D space based on chosen plane + for i, (x2d, y2d) in enumerate(coords): + if self.plane == "xy": + self.coordinates[i] = (x2d, y2d, 0) + elif self.plane == "xz": + self.coordinates[i] = (x2d, 0, y2d) + elif self.plane == "yz": + self.coordinates[i] = (0, x2d, y2d) diff --git a/mbuild/path/path_utils.py b/mbuild/path/path_utils.py new file mode 100644 index 000000000..831c7d237 --- /dev/null +++ b/mbuild/path/path_utils.py @@ -0,0 +1,84 @@ +"""Utility functions (mostly numba) for mbuild path generation""" + +import numpy as np +from numba import njit + + +@njit(cache=True, fastmath=True) +def random_coordinate( + pos1, + pos2, + bond_length, + thetas, + r_vectors, + batch_size, +): + """Default next_step method for HardSphereRandomWalk.""" + v1 = pos2 - pos1 + v1_norm = v1 / norm(v1) + dot_products = (r_vectors * v1_norm).sum(axis=1) + r_perp = r_vectors - dot_products[:, None] * v1_norm + norms = np.sqrt((r_perp * r_perp).sum(axis=1)) + # Handle rare cases where rprep vectors approach zero + norms = np.where(norms < 1e-6, 1.0, norms) + r_perp_norm = r_perp / norms[:, None] + # Batch of trial next-step vectors using angles and r_norms + cos_thetas = np.cos(thetas) + sin_thetas = np.sin(thetas) + v2s = cos_thetas[:, None] * v1_norm + sin_thetas[:, None] * r_perp_norm + # Batch of trial positions + next_positions = pos1 + v2s * bond_length + return next_positions + + +@njit(cache=True, fastmath=True) +def check_path(existing_points, new_point, radius, tolerance): + """Default check path method for HardSphereRandomWalk.""" + min_sq_dist = (radius - tolerance) ** 2 + for i in range(existing_points.shape[0]): + dist_sq = 0.0 + for j in range(existing_points.shape[1]): + diff = existing_points[i, j] - new_point[j] + dist_sq += diff * diff + if dist_sq < min_sq_dist: + return False + return True + + +@njit(cache=True, fastmath=True) +def target_sq_distances(target_coordinate, new_points): + """Return squared distances from target_coordinate to new_points.""" + n_points = new_points.shape[0] + sq_distances = np.empty(n_points, dtype=np.float32) + for i in range(n_points): + dx = target_coordinate[0] - new_points[i, 0] + dy = target_coordinate[1] - new_points[i, 1] + dz = target_coordinate[2] - new_points[i, 2] + sq_distances[i] = dx * dx + dy * dy + dz * dz + return sq_distances + + +@njit(cache=True, fastmath=True) +def norm(vec): + """Use in place of np.linalg.norm inside of numba functions.""" + s = 0.0 + for i in range(vec.shape[0]): + s += vec[i] * vec[i] + return np.sqrt(s) + + +@njit(cache=True, fastmath=True) +def rotate_vector(v, axis, theta): + """Rotate vector v around a normalized axis by angle theta using Rodrigues' formula.""" + c = np.cos(theta) + s = np.sin(theta) + k = axis + k_dot_v = k[0] * v[0] + k[1] * v[1] + k[2] * v[2] + cross = np.zeros(3) + cross[0] = k[1] * v[2] - k[2] * v[1] + cross[1] = k[2] * v[0] - k[0] * v[2] + cross[2] = k[0] * v[1] - k[1] * v[0] + rotated = np.zeros(3) + for i in range(3): + rotated[i] = v[i] * c + cross[i] * s + k[i] * k_dot_v * (1 - c) + return rotated diff --git a/mbuild/lib/recipes/polymer.py b/mbuild/polymer.py similarity index 71% rename from mbuild/lib/recipes/polymer.py rename to mbuild/polymer.py index 2fe1509be..0ee12ef4b 100644 --- a/mbuild/lib/recipes/polymer.py +++ b/mbuild/polymer.py @@ -14,11 +14,17 @@ ) from mbuild.lib.atoms import H from mbuild.port import Port +from mbuild.simulation import energy_minimize as e_min from mbuild.utils.validation import assert_port_exists __all__ = ["Polymer"] +class Monomer(Compound): + def __init__(self): + pass + + class Polymer(Compound): """Connect one or more components in a specified sequence. @@ -39,6 +45,11 @@ class Polymer(Compound): add_end_groups(compound, index, separation, orientation, replace) Use to add an end group compound to Polymer.end_groups + build_from_path(path, sequence, add_hydrogens, bond_head_tail, energy_minimize) + Used to set the polymer configuration from a pre-determined mbuild.path.Path + Use this to create polymers that follow a random walk, or form + lamellar order, ring polymers and more. See the ``mbuild.path`` module. + build(n, sequence) Use to create a single polymer compound. This method uses the compounds created by calling the add_monomer and add_end_group methods. @@ -118,7 +129,109 @@ def end_groups(self): """ return self._end_groups - def build(self, n, sequence="A", add_hydrogens=True): + def build_from_path( + self, + path, + sequence="A", + add_hydrogens=True, + bond_head_tail=False, + energy_minimize=True, + ): + """Build the polymer chain to a pre-defined configuraiton. + + See ``mbuild.path.Path`` to available polymer configurations + or use mbuild.path.Path.from_coordinates() to manually set + a polymer chan configuraiton. + + Parameters + ---------- + path : mbuild.path.Path, required + The configuration the polymer will be mapped to after connecting + all of the components (monomers and end groups). + The path will determine the number of monomer sites. + sequence : str, optional, default 'A' + A string of characters where each unique character represents one + repetition of a monomer. Characters in `sequence` are assigned to + monomers in the order they appear in `Polymer.monomers`. + The characters in `sequence` are assigned to the compounds in the + in the order that they appear in the Polymer.monomers list. + For example, 'AB' where 'A'corresponds to the first compound + added to Polymer.monomers and 'B' to the second compound. + add_hydrogens : bool, default True + If True and an ``end_groups`` compound is None, then the head or tail + of the polymer will be capped off with hydrogen atoms. If end group + compounds exist, then they will be used. + If False and an end group compound is None, then the head or tail + port will be exposed in the polymer. + bond_head_tail : bool, default False + If set to ``True``, then a bond between the head and tail groups (monomers or end groups) + is created. This does not create a periodic bond (see Polymer.create_periodic_bond). + The ``add_hydrogens`` parameter must be set to ``False`` to create this bond. + This is useful for creating ring polymers, knot polymers. + See ``mbuild.path.Cyclic`` and ``mbuild.path.Knot``. + energy_minimize : bool, default True + If ``True`` then relax the bonds and angles that may be distorted from mapping + the atomistic polymer to a path. + This uses the capped displacement methods in ``mbuild.simulation``. + + Notes + ----- + Energy Minimization: + + It is required to run energy minimization on chains built from a path to relax + the topology to a suitable starting point. + mBuild contains multiple energy minimization approaches in the ``mbuild.simulation`` module. + + When ``energy_minimize`` is set to ``True``, this method uses the OpenBabel based energy + minimization method with the UFF force field. We have found that once polymers + reach a size on the order of ~500 atoms, significantly faster energy minimization can be + achieved using one of the hoomd-based methods in ``mbuild.simulation``. + In that case, set ``energy_minimize=False`` and pass the resulting polymer + compound into one of these methods. + See: ``mbuild.simulation.hoomd_cap_displacement``, and ``mbuild.simulation.hoomd_fire``. + + """ + n = len(path.coordinates) - sum([1 for i in self.end_groups if i is not None]) + self.build( + n=n // len(sequence), + sequence=sequence, + add_hydrogens=add_hydrogens, + bond_head_tail=bond_head_tail, + ) + self.set_monomer_positions( + coordinates=path.coordinates, energy_minimize=energy_minimize + ) + + def set_monomer_positions(self, coordinates, energy_minimize=True): + """Shift monomers so that their center of mass matches a set of pre-defined coordinates. + + Parameters + ---------- + coordinates : np.ndarray, shape=(N,3) + Set of x,y,z coordinatess + + Notes + ----- + Energy Minimization: + + It is required to run energy minimization on chains built from a path to relax + the topology to a suitable starting point. + mBuild contains multiple energy minimization approaches in the ``mbuild.simulation`` module. + + When ``energy_minimize`` is set to ``True``, this method uses the OpenBabel based energy + minimization method with the UFF force field. We have found that once polymers + reach a size on the order of ~500 atoms, significantly faster energy minimization can be + achieved using one of the hoomd-based methods in ``mbuild.simulation`` + In that case, set ``energy_minimize=False`` and pass the resulting polymer + compound into one of these methods. + See: `hoomd_cap_displacement`, and `hoomd_fire`. + """ + for i, xyz in enumerate(coordinates): + self.children[i].translate_to(xyz) + if energy_minimize: + e_min(self) + + def build(self, n, sequence="A", add_hydrogens=True, bond_head_tail=False): """Connect one or more components in a specified sequence. Uses the compounds that are stored in Polymer.monomers and @@ -146,18 +259,28 @@ def build(self, n, sequence="A", add_hydrogens=True): compounds exist, then they will be used. If False and an end group compound is None, then the head or tail port will be exposed in the polymer. - periodic_axis : str, default None - If not ``None`` and an ``end_groups`` compound is None, then the head - and tail will be forced into an overlap with a periodicity along the - ``axis`` (default="z") specified. See - :meth:`mbuild.lib.recipes.polymer.Polymer.create_periodic_bond` for - more details. If end group compounds exist, then there will be a - warning. However ``add_hydrogens`` will simply be overwritten. - If ``None``, ``end_groups`` compound is None, and ``add_hydrogens`` is - False then the head or tail port will be exposed in the polymer. + bond_head_tail : bool, default False + If set to ``True``, then a bond between the head and tail groups (monomers or end groups) + is created. This does not create a periodic bond (see Polymer.create_periodic_bond). + The ``add_hydrogens`` parameter must be set to ``False`` to create this bond. + This is useful for creating ring polymers, knot polymers. + See ``mbuild.path.Cyclic`` and ``mbuild.path.Knot``. + + Notes + ----- + The chain conformations obtained from this method are often difficult to use + in later steps of creating systems of polymers. See the alternative build + method ``Polymer.build_from_path``. """ + if add_hydrogens and bond_head_tail: + raise ValueError( + "In order to bond head and tail ports, the Polymer instance cannot contain " + "end_groups and add_hydrogens must be set to `False`" + ) if n < 1: raise ValueError("n must be 1 or more") + head_and_tail_compounds = [] + repeat_compounds = [] for monomer in self._monomers: for label in self._port_labels: @@ -176,7 +299,8 @@ def build(self, n, sequence="A", add_hydrogens=True): last_part = None for n_added, seq_item in enumerate(it.cycle(sequence)): this_part = clone(seq_map[seq_item]) - self.add(this_part, "monomer[$]") + repeat_compounds.append(this_part) + # self.add(this_part, "monomer[$]") if last_part is None: first_part = this_part else: @@ -199,7 +323,8 @@ def build(self, n, sequence="A", add_hydrogens=True): if compound is not None: if self._headtail[i] is not None: head_tail[i].update_separation(self._headtail[i]) - self.add(compound) + # self.add(compound) + head_and_tail_compounds.append(compound) force_overlap(compound, compound.labels["up"], head_tail[i]) head_tail[i] = None else: @@ -230,6 +355,15 @@ def build(self, n, sequence="A", add_hydrogens=True): for port in self.all_ports(): if id(port) not in port_ids: self.remove(port) + if self.end_groups[0]: + self.add(self.end_groups[0]) + for compound in repeat_compounds: + self.add(new_child=compound, label="monomer[$]") + if self.end_groups[1]: + self.add(self.end_groups[1]) + + if bond_head_tail: + force_overlap(self, self.head_port, self.tail_port) def add_monomer( self, diff --git a/mbuild/port.py b/mbuild/port.py index 65326967b..330510cd5 100644 --- a/mbuild/port.py +++ b/mbuild/port.py @@ -1,7 +1,7 @@ """Ports used to facilitate bond formation.""" import itertools -from warnings import warn +import logging import numpy as np @@ -9,6 +9,8 @@ from mbuild.compound import Compound, Particle from mbuild.coordinate_transform import angle, unit_vector +logger = logging.getLogger(__name__) + class Port(Compound): """A set of four ghost Particles used to connect parts. @@ -93,7 +95,7 @@ def update_separation(self, separation): shifted from the origin. """ if self.used: - warn( + logger.warning( "This port is already being used and changing its separation " "will have no effect on the distance between particles." ) @@ -112,7 +114,7 @@ def update_orientation(self, orientation): Vector along which to orient the port """ if self.used: - warn( + logger.warning( "This port is already being used and changing its orientation " "will have no effect on the direction between particles." ) @@ -147,7 +149,7 @@ def separation(self): if self.anchor: return np.linalg.norm(self.center - self.anchor.pos) else: - warn( + logger.warning( "This port is not anchored to another particle. Returning a " "separation of None" ) diff --git a/mbuild/simulation.py b/mbuild/simulation.py new file mode 100644 index 000000000..585c3ecd8 --- /dev/null +++ b/mbuild/simulation.py @@ -0,0 +1,1205 @@ +"""Simulation methods that operate on mBuild compounds.""" + +import os +import tempfile +from warnings import warn + +import gmso +import hoomd +import numpy as np +from ele.element import element_from_name, element_from_symbol +from ele.exceptions import ElementError +from gmso.parameterization import apply + +from mbuild import Compound +from mbuild.exceptions import MBuildError +from mbuild.utils.io import import_ + + +class HoomdSimulation(hoomd.simulation.Simulation): + """A custom class to help in creating internal hoomd-based simulation methods. + + See ``hoomd_cap_displacement`` and ``hoomd_fire``. + + Parameters + ---------- + compound : mb.Compound + The compound to use in the simulation + forcefield : foyer.forcefield.Forcefield or gmso.core.Forcefield + The forcefield to apply to the system + r_cut : float (nm) + The cutoff distance (nm) used in the non-bonded pair neighborlist. + Use smaller values for unstable starting conditions and faster performance. + fixed_compounds : list of mb.Compound, default None + If given, these compounds will be removed from the integration updates and + held "frozen" during the hoomd simulation. + They are still able to interact with other particles in the simulation. + If desired, pass in a subset of children from `compound.children`. + integrate_compounds : list of mb.Compound, default None + If given, then only these compounds will be updated during integration + and all other compunds in `compound` are removed from the integration group. + run_on_gpu : bool, default False + When `True` the HOOMD simulation uses the hoomd.device.GPU() device. + This requires that you have a GPU compatible HOOMD install and + a compatible GPU device. + seed : int, default 42 + The seed passed to HOOMD + """ + + def __init__( + self, + compound, + forcefield, + r_cut, + run_on_gpu, + seed, + automatic_box=False, + box_buffer=1.5, + integrate_compounds=None, + fixed_compounds=None, + gsd_file_name=None, + ): + if run_on_gpu: + try: + device = hoomd.device.GPU() + print(f"GPU found, running on device {device.device}") + except RuntimeError: + print( + "Unable to find compatible GPU device. ", + "set `run_on_gpu = False` or see HOOMD documentation " + "for further information about GPU support.", + ) + else: + device = hoomd.device.CPU() + self.compound = compound + self.forcefield = forcefield + self.r_cut = r_cut + self.integrate_compounds = integrate_compounds + self.fixed_compounds = fixed_compounds + self.automatic_box = automatic_box + self.box_buffer = box_buffer + # TODO + # self.set_box + # Check if a hoomd sim method has been used on this compound already + if compound._hoomd_data: + last_snapshot, last_forces, last_forcefield = compound._get_sim_data() + # Check if the forcefield passed this time matches last time + if forcefield == last_forcefield: + snapshot = last_snapshot + self.forces = last_forces + else: # New foyer/gmso forcefield has been passed, reapply + snapshot, self.forces = self._to_hoomd_snap_forces() + compound._add_sim_data( + state=snapshot, forces=self.forces, forcefield=forcefield + ) + else: # Sim method not used on this compound previously + snapshot, self.forces = self._to_hoomd_snap_forces() + compound._add_sim_data( + state=snapshot, forces=self.forces, forcefield=forcefield + ) + # Place holders for forces added/changed by specific methods below + self.active_forces = [] + self.inactive_forces = [] + super(HoomdSimulation, self).__init__(device=device, seed=seed) + self.create_state_from_snapshot(snapshot) + + def _to_hoomd_snap_forces(self): + # If a box isn't set, make one with the buffer + if not self.compound.box: + self.compound.box = self.compound.get_boundingbox(pad_box=self.box_buffer) + # Convret to GMSO, apply forcefield + top = self.compound.to_gmso() + top.identify_connections() + apply(top, forcefields=self.forcefield, ignore_params=["dihedral", "improper"]) + # Get hoomd snapshot and force objects + forces, ref = gmso.external.to_hoomd_forcefield(top, r_cut=self.r_cut) + snap, ref = gmso.external.to_gsd_snapshot(top) + forces = list(set().union(*forces.values())) + return snap, forces + + def get_integrate_group(self): + # Get indices of compounds to include in integration + if self.integrate_compounds and not self.fixed_compounds: + integrate_indices = [] + for comp in self.integrate_compounds: + integrate_indices.extend(list(self.compound.get_child_indices(comp))) + return hoomd.filter.Tags(integrate_indices) + # Get indices of compounds to NOT include in integration + elif self.fixed_compounds and not self.integrate_compounds: + fix_indices = [] + for comp in self.fixed_compounds: + fix_indices.extend(list(self.compound.get_child_indices(comp))) + return hoomd.filter.SetDifference( + hoomd.filter.All(), hoomd.filter.Tags(fix_indices) + ) + # Both are passed in, not supported + elif self.integrate_compounds and self.fixed_compounds: + raise RuntimeError( + "You can specify only one of integrate_compounds and fixed_compounds." + ) + # Neither are given, include everything in integration + else: + return hoomd.filter.All() + + def get_force(self, instance): + for force in set(self.forces + self.active_forces + self.inactive_forces): + if isinstance(force, instance): + return force + raise ValueError(f"No force of {instance} was found.") + + def get_dpd_from_lj(self, A): + """Make a best-guess DPD force from types and parameters of an LJ force.""" + lj_force = self.get_force(hoomd.md.pair.LJ) + dpd = hoomd.md.pair.DPDConservative(nlist=lj_force.nlist) + for param in lj_force.params: + dpd.params[param] = dict(A=A) + dpd.r_cut[param] = lj_force.params[param]["sigma"] + return dpd + + def set_fire_integrator( + self, + dt, + force_tol, + angmom_tol, + energy_tol, + methods, + finc_dt, + fdec_dt, + alpha_start, + fdec_alpha, + min_steps_adapt, + min_steps_conv, + ): + fire = hoomd.md.minimize.FIRE( + dt=dt, + force_tol=force_tol, + angmom_tol=angmom_tol, + energy_tol=energy_tol, + finc_dt=finc_dt, + fdec_dt=fdec_dt, + alpha_start=alpha_start, + fdec_alpha=fdec_alpha, + min_steps_adapt=min_steps_adapt, + min_steps_conv=min_steps_conv, + methods=methods, + ) + fire.forces = self.active_forces + self.operations.integrator = fire + + def set_integrator(self, dt, method): + integrator = hoomd.md.Integrator(dt=dt) + integrator.forces = self.active_forces + integrator.methods = [method] + self.operations.integrator = integrator + + def _update_integrator_forces(self): + self.operations.integrator.forces = self.active_forces + + def _update_snapshot(self): + snapshot = self.state.get_snapshot() + # snapshot_copy = hoomd.Snapshot(snapshot) # full copy + self.compound._add_sim_data(state=snapshot) + + def add_gsd_writer(self, file_name, write_period): + """""" + pass + + def update_positions(self): + """Update compound positions from snapshot.""" + with self.state.cpu_local_snapshot as snap: + particles = snap.particles.rtag[:] + pos = snap.particles.position[particles] + self.compound.xyz = np.copy(pos) + + +## HOOMD METHODS ## +def hoomd_cap_displacement( + compound, + forcefield, + n_steps, + dt, + r_cut, + max_displacement=1e-3, + fixed_compounds=None, + integrate_compounds=None, + dpd_A=None, + bond_k_scale=1, + angle_k_scale=1, + n_relax_steps=0, + run_on_gpu=False, + seed=42, +): + """Run a short simulation with hoomd.md.methods.DisplacementCapped + + Parameters + ---------- + compound : mb.Compound + The compound to use in the simulation + forcefield : foyer.forcefield.Forcefield or gmso.core.Forcefield + The forcefield to apply to the system + n_steps : int + The number of simulation time steps to run with the modified forcefield + created by bond_k_scale, angle_k_scale, and dpd_A. + See these parameters for more information. + dt : float + The simulation timestep. Note that for running capped displacement on highly unstable + systems (e.g. overlapping particles), it is often more stable to use a larger + value of dt and let the displacement cap limit position updates. + r_cut : float (nm) + The cutoff distance (nm) used in the non-bonded pair neighborlist. + Use smaller values for unstable starting conditions and faster performance. + max_displacement : float (nm), default = 1e-3 nm + The maximum displacement (nm) allowed per timestep. Use smaller values for + highly unstable systems. + fixed_compounds : list of mb.Compound, default None + If given, these compounds will be removed from the integration updates and + held "frozen" during the hoomd simulation. + They are still able to interact with other particles in the simulation. + If desired, pass in a subset of children from `compound.children`. + integrate_compounds : list of mb.Compound, default None + If given, then only these compounds will be updated during integration + and all other compunds in `compound` are removed from the integration group. + dpd_A : float, default None + If set to a value, then the initial simulation replaces the LJ 12-6 pair potential + with the softer hoomd.md.pair.DPDConservative pair force from HOOMD. + This is used for `n_steps` and replaced with the original LJ for potential which is + used for n_relax_steps. This can be useful for highly unstable starting configutations. + Note that the cutoff for the DPD force is set to the sigma value for each atom type + and force cosntant (A) is set to dpd_A for all atom types. + bond_k_scale : float, default 1 + Scale the bond-stretching force constants so that K_new = K * bond_k_scale. + This can be useful for relaxing highly unstable starting configutations. + The bond force constants are scaled back to their original values + during the run of `n_relax_steps` + angle_k_scale : float, default 1 + Scale the bond-angle force constants so that K_new = K * angle_k_scale. + This can be useful for relaxing highly unstable starting configutations. + The angle force constants are scaled back to their original values + during the run of `n_relax_steps` + n_relax_steps: int, optional + The number of steps to run after running for `n_steps` with the modified + forcefield. This is designed to be used when utilizing the parameters + `dpd_A`, `bond_k_scale`, and `angle_k_scale`. + run_on_gpu : bool, default False + When `True` the HOOMD simulation uses the hoomd.device.GPU() device. + This requires that you have a GPU compatible HOOMD install and + a compatible GPU device. + seed : int, default 42 + The seed passed to HOOMD + + Notes + ----- + Running on GPU: + The hoomd-based methods in ``mbuild.simulation`` + are able to utilize and run on GPUs to perform energy minimization for larger systems if needed. + Performance improvements of using GPU over CPU may not be significant until systems reach a + size of ~1,000 particles. + + Running multiple simulations: + The information needed to run the HOOMD simulation is saved after the first simulation. + Therefore, calling hoomd-based simulation methods multiple times (e.g, in a for loop or while loop) + does not require re-running performing atom-typing, applying the forcefield, or running hoomd + format writers again. + """ + compound._kick() + sim = HoomdSimulation( + compound=compound, + forcefield=forcefield, + r_cut=r_cut, + run_on_gpu=run_on_gpu, + seed=seed, + integrate_compounds=integrate_compounds, + fixed_compounds=fixed_compounds, + ) + bond = sim.get_force(hoomd.md.bond.Harmonic) + angle = sim.get_force(hoomd.md.angle.Harmonic) + lj = sim.get_force(hoomd.md.pair.LJ) + # If needed scale bond K and angle K + if bond_k_scale != 1.0: + for param in bond.params: + bond.params[param]["k"] *= bond_k_scale + sim.active_forces.append(bond) + # If angle_k_scale is zero, then "turn off" angles; skip adding the angle force + if angle_k_scale != 0: + if angle_k_scale != 1.0: + for param in angle.params: + angle.params[param]["k"] *= angle_k_scale + sim.active_forces.append(angle) + if dpd_A: + dpd = sim.get_dpd_from_lj(A=dpd_A) + sim.active_forces.append(dpd) + else: + sim.active_forces.append(lj) + # Set up and run + displacement_capped = hoomd.md.methods.DisplacementCapped( + filter=sim.get_integrate_group(), + maximum_displacement=max_displacement, + ) + sim.set_integrator(method=displacement_capped, dt=dt) + sim.run(n_steps) + # Run with LJ pairs, original bonds and angles + if n_relax_steps > 0: + if dpd_A: # Repalce DPD with original LJ for final relaxation + sim.active_forces.remove(dpd) + sim.active_forces.append(lj) + if bond_k_scale != 1.0: + for param in bond.params: + bond.params[param]["k"] /= bond_k_scale + if angle_k_scale != 0: # Angle force already added, rescale + if angle_k_scale != 1.0: + for param in angle.params: + angle.params[param]["k"] /= angle_k_scale + else: # Don't rescale, just add the original angle force + sim.active_forces.append(angle) + sim._update_integrator_forces() + sim.run(n_relax_steps) + # Update particle positions, save latest state point snapshot + sim.update_positions() + sim._update_snapshot() + + +def hoomd_fire( + compound, + forcefield, + run_on_gpu, + fire_iteration_steps, + num_fire_iterations=1, + n_relax_steps=1000, + fixed_compounds=None, + integrate_compounds=None, + dt=1e-5, + dpd_A=10, + integrate_dof=False, + min_steps_adapt=5, + min_steps_conv=100, + finc_dt=1.1, + fdec_dt=0.5, + alpha_start=0.1, + fdec_alpha=0.95, + force_tol=1e-2, + angmom_tol=1e-2, + energy_tol=1e-6, + seed=42, + r_cut=1.0, + bond_k_scale=100, + angle_k_scale=100, + gsd_file=None, +): + """Run a short HOOMD-Blue simulation with the FIRE integrator. + This method can be helpful for relaxing poor molecular + geometries or for removing overlapping particles. + + Parameters: + ----------- + compound : mbuild.compound.Compound; required + The compound to perform the displacement capped simulation with. + forcefield : foyer.focefield.ForceField or gmso.core.Forcefield; required + The forcefield used for the simulation. + steps : int; required + The number of simulation steps to run + max_displacement : float, default 1e-3 + The value of the maximum displacement (nm) + run_on_gpu : bool, default True + If true, attempts to run HOOMD-Blue on a GPU. + If a GPU device isn't found, then it will run on the CPU + + Notes + ----- + Running on GPU: + The hoomd-based methods in ``mbuild.simulation`` + are able to utilize and run on GPUs to perform energy minimization for larger systems if needed. + Performance improvements of using GPU over CPU may not be significant until systems reach a + size of ~1,000 particles. + + Running multiple simulations: + The information needed to run the HOOMD simulation is saved after the first simulation. + Therefore, calling hoomd-based simulation methods multiple times (e.g, in a for loop or while loop) + does not require re-running performing atom-typing, applying the forcefield, or running hoomd + format writers again. + """ + compound._kick() + sim = HoomdSimulation( + compound=compound, + forcefield=forcefield, + r_cut=r_cut, + run_on_gpu=run_on_gpu, + seed=seed, + integrate_compounds=integrate_compounds, + fixed_compounds=fixed_compounds, + ) + bond = sim.get_force(hoomd.md.bond.Harmonic) + angle = sim.get_force(hoomd.md.angle.Harmonic) + lj = sim.get_force(hoomd.md.pair.LJ) + # Scale bond K and angle K + if bond_k_scale != 1.0: + for param in bond.params: + bond.params[param]["k"] *= bond_k_scale + sim.active_forces.append(bond) + # If angle_k_scale is zero, skip and don't append angle force + if angle_k_scale != 0: + if angle_k_scale != 1.0: + for param in angle.params: + angle.params[param]["k"] *= angle_k_scale + sim.active_forces.append(angle) + if dpd_A: + dpd = sim.get_dpd_from_lj(A=dpd_A) + sim.active_forces.append(dpd) + else: + sim.active_forces.append(lj) + # Set up and run + nvt = hoomd.md.methods.ConstantVolume(filter=sim.get_integrate_group()) + sim.set_fire_integrator( + dt=dt, + force_tol=force_tol, + angmom_tol=angmom_tol, + energy_tol=energy_tol, + finc_dt=finc_dt, + fdec_dt=fdec_dt, + alpha_start=alpha_start, + fdec_alpha=fdec_alpha, + min_steps_adapt=min_steps_adapt, + min_steps_conv=min_steps_conv, + methods=[nvt], + ) + # Run FIRE sims with DPD + scaled bonds and angles + for sim_num in range(num_fire_iterations): + sim.run(fire_iteration_steps) + sim.operations.integrator.reset() + + # Re-scale bonds and angle force constants + if n_relax_steps > 0: + if dpd_A: # Replace DPD pair force with original LJ + sim.active_forces.remove(dpd) + sim.active_forces.append(lj) + if bond_k_scale != 1.0: + for param in bond.params: + bond.params[param]["k"] /= bond_k_scale + if angle_k_scale != 0: # Angle force already included, rescale. + if angle_k_scale != 1.0: + for param in angle.params: + angle.params[param]["k"] /= angle_k_scale + else: # Don't recale, just add angles for final relaxation step + sim.active_forces.append(angle) + sim._update_integrator_forces() + sim.run(n_relax_steps) + # Update particle positions, save latest state point snapshot + sim._update_snapshot() + sim.update_positions() + + +# Openbabel and OpenMM +def energy_minimize( + compound, + forcefield="UFF", + steps=1000, + shift_com=True, + anchor=None, + **kwargs, +): + """Perform an energy minimization on a Compound. + + Default behavior utilizes `Open Babel `_ + to perform an energy minimization/geometry optimization on a Compound by + applying a generic force field + + Can also utilize `OpenMM `_ to energy minimize after + atomtyping a Compound using + `Foyer `_ to apply a forcefield XML + file that contains valid SMARTS strings. + + This function is primarily intended to be used on smaller components, + with sizes on the order of 10's to 100's of particles, as the energy + minimization scales poorly with the number of particles. + + Parameters + ---------- + compound : mbuid.Compound, required + The compound to perform energy minimization on. + steps : int, optional, default=1000 + The number of optimization iterations + forcefield : str, optional, default='UFF' + The generic force field to apply to the Compound for minimization. + Valid options are 'MMFF94', 'MMFF94s', ''UFF', 'GAFF', 'Ghemical'. + Please refer to the `Open Babel documentation + `_ + when considering your choice of force field. + Utilizing OpenMM for energy minimization requires a forcefield + XML file with valid SMARTS strings. Please refer to `OpenMM docs + `_ + for more information. + shift_com : bool, optional, default=True + If True, the energy-minimized Compound is translated such that the + center-of-mass is unchanged relative to the initial configuration. + anchor : Compound, optional, default=None + Translates the energy-minimized Compound such that the + position of the anchor Compound is unchanged relative to the + initial configuration. + + Other Parameters + ---------------- + algorithm : str, optional, default='cg' + The energy minimization algorithm. Valid options are 'steep', 'cg', + and 'md', corresponding to steepest descent, conjugate gradient, and + equilibrium molecular dynamics respectively. + For _energy_minimize_openbabel + fixed_compounds : Compound, optional, default=None + An individual Compound or list of Compounds that will have their + position fixed during energy minimization. Note, positions are fixed + using a restraining potential and thus may change slightly. + Position fixing will apply to all Particles (i.e., atoms) that exist + in the Compound and to particles in any subsequent sub-Compounds. + By default x,y, and z position is fixed. This can be toggled by instead + passing a list containing the Compound and an list or tuple of bool values + corresponding to x,y and z; e.g., [Compound, (True, True, False)] + will fix the x and y position but allow z to be free. + For _energy_minimize_openbabel + ignore_compounds: Compound, optional, default=None + An individual compound or list of Compounds whose underlying particles + will have their positions fixed and not interact with other atoms via + the specified force field during the energy minimization process. + Note, a restraining potential used and thus absolute position may vary + as a result of the energy minimization process. + Interactions of these ignored atoms can be specified by the user, + e.g., by explicitly setting a distance constraint. + For _energy_minimize_openbabel + distance_constraints: list, optional, default=None + A list containing a pair of Compounds as a tuple or list and + a float value specifying the target distance between the two Compounds, e.g.,: + [(compound1, compound2), distance]. + To specify more than one constraint, pass constraints as a 2D list, e.g.,: + [ [(compound1, compound2), distance1], [(compound3, compound4), distance2] ]. + Note, Compounds specified here must represent individual point particles. + For _energy_minimize_openbabel + constraint_factor: float, optional, default=50000.0 + Harmonic springs are used to constrain distances and fix atom positions, where + the resulting energy associated with the spring is scaled by the + constraint_factor; the energy of this spring is considering during the minimization. + As such, very small values of the constraint_factor may result in an energy + minimized state that does not adequately restrain the distance/position of atoms. + For _energy_minimize_openbabel + scale_bonds : float, optional, default=1 + Scales the bond force constant (1 is completely on). + For _energy_minimize_openmm + scale_angles : float, optional, default=1 + Scales the angle force constant (1 is completely on) + For _energy_minimize_openmm + scale_torsions : float, optional, default=1 + Scales the torsional force constants (1 is completely on) + For _energy_minimize_openmm + Note: Only Ryckaert-Bellemans style torsions are currently supported + scale_nonbonded : float, optional, default=1 + Scales epsilon (1 is completely on) + For _energy_minimize_openmm + constraints : str, optional, default="AllBonds" + Specify constraints on the molecule to minimize, options are: + None, "HBonds", "AllBonds", "HAngles" + For _energy_minimize_openmm + + References + ---------- + If using _energy_minimize_openmm(), please cite: + + .. [Eastman2013] P. Eastman, M. S. Friedrichs, J. D. Chodera, + R. J. Radmer, C. M. Bruns, J. P. Ku, K. A. Beauchamp, T. J. Lane, + L.-P. Wang, D. Shukla, T. Tye, M. Houston, T. Stich, C. Klein, + M. R. Shirts, and V. S. Pande. "OpenMM 4: A Reusable, Extensible, + Hardware Independent Library for High Performance Molecular + Simulation." J. Chem. Theor. Comput. 9(1): 461-469. (2013). + + If using _energy_minimize_openbabel(), please cite: + + .. [OBoyle2011] O'Boyle, N.M.; Banck, M.; James, C.A.; Morley, C.; + Vandermeersch, T.; Hutchison, G.R. "Open Babel: An open chemical + toolbox." (2011) J. Cheminf. 3, 33 + + .. [OpenBabel] Open Babel, version X.X.X http://openbabel.org, + (installed Month Year) + + If using the 'MMFF94' force field please also cite the following: + + .. [Halgren1996a] T.A. Halgren, "Merck molecular force field. I. Basis, + form, scope, parameterization, and performance of MMFF94." (1996) + J. Comput. Chem. 17, 490-519 + + .. [Halgren1996b] T.A. Halgren, "Merck molecular force field. II. MMFF94 + van der Waals and electrostatic parameters for intermolecular + interactions." (1996) J. Comput. Chem. 17, 520-552 + + .. [Halgren1996c] T.A. Halgren, "Merck molecular force field. III. + Molecular geometries and vibrational frequencies for MMFF94." (1996) + J. Comput. Chem. 17, 553-586 + + .. [Halgren1996d] T.A. Halgren and R.B. Nachbar, "Merck molecular force + field. IV. Conformational energies and geometries for MMFF94." (1996) + J. Comput. Chem. 17, 587-615 + + .. [Halgren1996e] T.A. Halgren, "Merck molecular force field. V. + Extension of MMFF94 using experimental data, additional computational + data, and empirical rules." (1996) J. Comput. Chem. 17, 616-641 + + If using the 'MMFF94s' force field please cite the above along with: + + .. [Halgren1999] T.A. Halgren, "MMFF VI. MMFF94s option for energy minimization + studies." (1999) J. Comput. Chem. 20, 720-729 + + If using the 'UFF' force field please cite the following: + + .. [Rappe1992] Rappe, A.K., Casewit, C.J., Colwell, K.S., Goddard, W.A. + III, Skiff, W.M. "UFF, a full periodic table force field for + molecular mechanics and molecular dynamics simulations." (1992) + J. Am. Chem. Soc. 114, 10024-10039 + + If using the 'GAFF' force field please cite the following: + + .. [Wang2004] Wang, J., Wolf, R.M., Caldwell, J.W., Kollman, P.A., + Case, D.A. "Development and testing of a general AMBER force field" + (2004) J. Comput. Chem. 25, 1157-1174 + + If using the 'Ghemical' force field please cite the following: + + .. [Hassinen2001] T. Hassinen and M. Perakyla, "New energy terms for + reduced protein models implemented in an off-lattice force field" + (2001) J. Comput. Chem. 22, 1229-1242 + + """ + # TODO: Update mbuild tutorials to provide overview of new features + # Preliminary tutorials: https://github.com/chrisiacovella/mbuild_energy_minimization + com = compound.pos + anchor_in_compound = False + if anchor is not None: + # check to see if the anchor exists + # in the Compound to be energy minimized + for succesor in compound.successors(): + if id(anchor) == id(succesor): + anchor_in_compound = True + anchor_pos_old = anchor.pos + + if not anchor_in_compound: + raise MBuildError( + f"Anchor: {anchor} is not part of the Compound: {compound}" + "that you are trying to energy minimize." + ) + compound._kick() + extension = os.path.splitext(forcefield)[-1] + openbabel_ffs = ["MMFF94", "MMFF94s", "UFF", "GAFF", "Ghemical"] + if forcefield in openbabel_ffs: + _energy_minimize_openbabel( + compound=compound, forcefield=forcefield, steps=steps, **kwargs + ) + else: + tmp_dir = tempfile.mkdtemp() + compound.save(os.path.join(tmp_dir, "un-minimized.mol2")) + + if extension == ".xml": + _energy_minimize_openmm( + compound=compound, + tmp_dir=tmp_dir, + forcefield_files=forcefield, + forcefield_name=None, + steps=steps, + **kwargs, + ) + else: + _energy_minimize_openmm( + compound=compound, + tmp_dir=tmp_dir, + forcefield_files=None, + forcefield_name=forcefield, + steps=steps, + **kwargs, + ) + + compound.update_coordinates(os.path.join(tmp_dir, "minimized.pdb")) + + if shift_com: + compound.translate_to(com) + + if anchor_in_compound: + anchor_pos_new = anchor.pos + delta = anchor_pos_old - anchor_pos_new + compound.translate(delta) + + +def _energy_minimize_openmm( + compound, + tmp_dir, + forcefield_files=None, + forcefield_name=None, + steps=1000, + scale_bonds=1, + scale_angles=1, + scale_torsions=1, + scale_nonbonded=1, + constraints="AllBonds", +): + """Perform energy minimization using OpenMM. + + Converts an mBuild Compound to a ParmEd Structure, + applies a forcefield using Foyer, and creates an OpenMM System. + + Parameters + ---------- + compound : mbuid.Compound, required + The compound to perform energy minimization on. + forcefield_files : str or list of str, optional, default=None + Forcefield files to load + forcefield_name : str, optional, default=None + Apply a named forcefield to the output file using the `foyer` + package, e.g. 'oplsaa'. `Foyer forcefields` + _ + steps : int, optional, default=1000 + Number of energy minimization iterations + scale_bonds : float, optional, default=1 + Scales the bond force constant (1 is completely on) + scale_angles : float, optiona, default=1 + Scales the angle force constant (1 is completely on) + scale_torsions : float, optional, default=1 + Scales the torsional force constants (1 is completely on) + scale_nonbonded : float, optional, default=1 + Scales epsilon (1 is completely on) + constraints : str, optional, default="AllBonds" + Specify constraints on the molecule to minimize, options are: + None, "HBonds", "AllBonds", "HAngles" + + Notes + ----- + Assumes a particular organization for the force groups + (HarmonicBondForce, HarmonicAngleForce, RBTorsionForce, NonBondedForce) + + References + ---------- + [Eastman2013]_ + """ + foyer = import_("foyer") + + to_parmed = compound.to_parmed() + ff = foyer.Forcefield(forcefield_files=forcefield_files, name=forcefield_name) + to_parmed = ff.apply(to_parmed) + + import openmm.unit as u + from openmm.app import AllBonds, HAngles, HBonds + from openmm.app.pdbreporter import PDBReporter + from openmm.app.simulation import Simulation + from openmm.openmm import LangevinIntegrator + + if constraints: + if constraints == "AllBonds": + constraints = AllBonds + elif constraints == "HBonds": + constraints = HBonds + elif constraints == "HAngles": + constraints = HAngles + else: + raise ValueError( + f"Provided constraints value of: {constraints}.\n" + f'Expected "HAngles", "AllBonds" "HBonds".' + ) + system = to_parmed.createSystem( + constraints=constraints + ) # Create an OpenMM System + else: + system = to_parmed.createSystem() # Create an OpenMM System + # Create a Langenvin Integrator in OpenMM + integrator = LangevinIntegrator( + 298 * u.kelvin, 1 / u.picosecond, 0.002 * u.picoseconds + ) + # Create Simulation object in OpenMM + simulation = Simulation(to_parmed.topology, system, integrator) + + # Loop through forces in OpenMM System and set parameters + for force in system.getForces(): + if type(force).__name__ == "HarmonicBondForce": + for bond_index in range(force.getNumBonds()): + atom1, atom2, r0, k = force.getBondParameters(bond_index) + force.setBondParameters(bond_index, atom1, atom2, r0, k * scale_bonds) + force.updateParametersInContext(simulation.context) + + elif type(force).__name__ == "HarmonicAngleForce": + for angle_index in range(force.getNumAngles()): + atom1, atom2, atom3, r0, k = force.getAngleParameters(angle_index) + force.setAngleParameters( + angle_index, atom1, atom2, atom3, r0, k * scale_angles + ) + force.updateParametersInContext(simulation.context) + + elif type(force).__name__ == "RBTorsionForce": + for torsion_index in range(force.getNumTorsions()): + ( + atom1, + atom2, + atom3, + atom4, + c0, + c1, + c2, + c3, + c4, + c5, + ) = force.getTorsionParameters(torsion_index) + force.setTorsionParameters( + torsion_index, + atom1, + atom2, + atom3, + atom4, + c0 * scale_torsions, + c1 * scale_torsions, + c2 * scale_torsions, + c3 * scale_torsions, + c4 * scale_torsions, + c5 * scale_torsions, + ) + force.updateParametersInContext(simulation.context) + + elif type(force).__name__ == "NonbondedForce": + for nb_index in range(force.getNumParticles()): + charge, sigma, epsilon = force.getParticleParameters(nb_index) + force.setParticleParameters( + nb_index, charge, sigma, epsilon * scale_nonbonded + ) + force.updateParametersInContext(simulation.context) + + elif type(force).__name__ == "CMMotionRemover": + pass + + else: + warn( + f"OpenMM Force {type(force).__name__} is " + "not currently supported in _energy_minimize_openmm. " + "This Force will not be updated!" + ) + + simulation.context.setPositions(to_parmed.positions) + # Run energy minimization through OpenMM + simulation.minimizeEnergy(maxIterations=steps) + reporter = PDBReporter(os.path.join(tmp_dir, "minimized.pdb"), 1) + reporter.report(simulation, simulation.context.getState(getPositions=True)) + + +def _check_openbabel_constraints( + compound, + particle_list, + successors_list, + check_if_particle=False, +): + """Provide routines commonly used to check constraint inputs.""" + for part in particle_list: + if not isinstance(part, Compound): + raise MBuildError(f"{part} is not a Compound.") + if id(part) != id(compound) and id(part) not in successors_list: + raise MBuildError(f"{part} is not a member of Compound {compound}.") + + if check_if_particle: + if len(part.children) != 0: + raise MBuildError( + f"{part} does not correspond to an individual particle." + ) + + +def _energy_minimize_openbabel( + compound, + steps=1000, + algorithm="cg", + forcefield="UFF", + constraint_factor=50000.0, + distance_constraints=None, + fixed_compounds=None, + ignore_compounds=None, +): + """Perform an energy minimization on a Compound. + + Utilizes Open Babel (http://openbabel.org/docs/dev/) to perform an + energy minimization/geometry optimization on a Compound by applying + a generic force field. + + This function is primarily intended to be used on smaller components, + with sizes on the order of 10's to 100's of particles, as the energy + minimization scales poorly with the number of particles. + + Parameters + ---------- + compound : mbuid.Compound, required + The compound to perform energy minimization on. + steps : int, optionl, default=1000 + The number of optimization iterations + algorithm : str, optional, default='cg' + The energy minimization algorithm. Valid options are 'steep', + 'cg', and 'md', corresponding to steepest descent, conjugate + gradient, and equilibrium molecular dynamics respectively. + forcefield : str, optional, default='UFF' + The generic force field to apply to the Compound for minimization. + Valid options are 'MMFF94', 'MMFF94s', ''UFF', 'GAFF', 'Ghemical'. + Please refer to the Open Babel documentation + (http://open-babel.readthedocs.io/en/latest/Forcefields/Overview.html) + when considering your choice of force field. + fixed_compounds : Compound, optional, default=None + An individual Compound or list of Compounds that will have their + position fixed during energy minimization. Note, positions are fixed + using a restraining potential and thus may change slightly. + Position fixing will apply to all Particles (i.e., atoms) that exist + in the Compound and to particles in any subsequent sub-Compounds. + By default x,y, and z position is fixed. This can be toggled by instead + passing a list containing the Compound and a list or tuple of bool values + corresponding to x,y and z; e.g., [Compound, (True, True, False)] + will fix the x and y position but allow z to be free. + ignore_compounds: Compound, optional, default=None + An individual compound or list of Compounds whose underlying particles + will have their positions fixed and not interact with other atoms via + the specified force field during the energy minimization process. + Note, a restraining potential is used and thus absolute position may vary + as a result of the energy minimization process. + Interactions of these ignored atoms can be specified by the user, + e.g., by explicitly setting a distance constraint. + distance_constraints: list, optional, default=None + A list containing a pair of Compounds as a tuple or list and + a float value specifying the target distance between the two Compounds, e.g.,: + [(compound1, compound2), distance]. + To specify more than one constraint, pass constraints as a 2D list, e.g.,: + [ [(compound1, compound2), distance1], [(compound3, compound4), distance2] ]. + Note, Compounds specified here must represent individual point particles. + constraint_factor: float, optional, default=50000.0 + Harmonic springs are used to constrain distances and fix atom positions, where + the resulting energy associated with the spring is scaled by the + constraint_factor; the energy of this spring is considering during the minimization. + As such, very small values of the constraint_factor may result in an energy + minimized state that does not adequately restrain the distance/position of atom(s)e. + + + References + ---------- + [OBoyle2011]_ + [OpenBabel]_ + + If using the 'MMFF94' force field please also cite the following: + [Halgren1996a]_ + [Halgren1996b]_ + [Halgren1996c]_ + [Halgren1996d]_ + [Halgren1996e]_ + + If using the 'MMFF94s' force field please cite the above along with: + [Halgren1999]_ + + If using the 'UFF' force field please cite the following: + [Rappe1992]_ + + If using the 'GAFF' force field please cite the following: + [Wang2001]_ + + If using the 'Ghemical' force field please cite the following: + [Hassinen2001]_ + """ + openbabel = import_("openbabel") + for particle in compound.particles(): + if particle.element is None: + try: + particle._element = element_from_symbol(particle.name) + except ElementError: + try: + particle._element = element_from_name(particle.name) + except ElementError: + raise MBuildError( + f"No element assigned to {particle}; element could not be" + f"inferred from particle name {particle.name}. Cannot perform" + "an energy minimization." + ) + # Create a dict containing particle id and associated index to speed up looping + particle_idx = { + id(particle): idx for idx, particle in enumerate(compound.particles()) + } + + # A list containing all Compounds ids contained in compound. Will be used to check if + # compounds refered to in the constrains are actually in the Compound we are minimizing. + successors_list = [id(comp) for comp in compound.successors()] + + # initialize constraints + ob_constraints = openbabel.OBFFConstraints() + + if distance_constraints is not None: + # if a user passes single constraint as a 1-D array, + # i.e., [(p1,p2), 2.0] rather than [[(p1,p2), 2.0]], + # just add it to a list so we can use the same looping code + if len(np.array(distance_constraints, dtype=object).shape) == 1: + distance_constraints = [distance_constraints] + + for con_temp in distance_constraints: + p1 = con_temp[0][0] + p2 = con_temp[0][1] + + _check_openbabel_constraints( + compound=compound, + particle_list=[p1, p2], + successors_list=successors_list, + check_if_particle=True, + ) + if id(p1) == id(p2): + raise MBuildError( + f"Cannot create a constraint between a Particle and itself: {p1} {p2} ." + ) + + # openbabel indices start at 1 + pid_1 = particle_idx[id(p1)] + 1 + # openbabel indices start at 1 + pid_2 = particle_idx[id(p2)] + 1 + dist = ( + con_temp[1] * 10.0 + ) # obenbabel uses angstroms, not nm, convert to angstroms + + ob_constraints.AddDistanceConstraint(pid_1, pid_2, dist) + + if fixed_compounds is not None: + # if we are just passed a single Compound, wrap it into + # and array so we can just use the same looping code + if isinstance(fixed_compounds, Compound): + fixed_compounds = [fixed_compounds] + + # if fixed_compounds is a 1-d array and it is of length 2, we need to determine whether it is + # a list of two Compounds or if fixed_compounds[1] should correspond to the directions to constrain + if len(np.array(fixed_compounds, dtype=object).shape) == 1: + if len(fixed_compounds) == 2: + if not isinstance(fixed_compounds[1], Compound): + # if it is not a list of two Compounds, make a 2d array so we can use the same looping code + fixed_compounds = [fixed_compounds] + + for fixed_temp in fixed_compounds: + # if an individual entry is a list, validate the input + if isinstance(fixed_temp, list): + if len(fixed_temp) == 2: + msg1 = ( + "Expected tuple or list of length 3 to set" + "which dimensions to fix motion." + ) + assert isinstance(fixed_temp[1], (list, tuple)), msg1 + + msg2 = ( + "Expected tuple or list of length 3 to set" + "which dimensions to fix motion, " + f"{len(fixed_temp[1])} found." + ) + assert len(fixed_temp[1]) == 3, msg2 + + dims = [dim for dim in fixed_temp[1]] + msg3 = ( + "Expected bool values for which directions are fixed." + f"Found instead {dims}." + ) + assert all(isinstance(dim, bool) for dim in dims), msg3 + + p1 = fixed_temp[0] + + # if fixed_compounds is defined as [[Compound],[Compound]], + # fixed_temp will be a list of length 1 + elif len(fixed_temp) == 1: + p1 = fixed_temp[0] + dims = [True, True, True] + + else: + p1 = fixed_temp + dims = [True, True, True] + + all_true = all(dims) + + _check_openbabel_constraints( + compound=compound, particle_list=[p1], successors_list=successors_list + ) + + if len(p1.children) == 0: + pid = particle_idx[id(p1)] + 1 # openbabel indices start at 1 + + if all_true: + ob_constraints.AddAtomConstraint(pid) + else: + if dims[0]: + ob_constraints.AddAtomXConstraint(pid) + if dims[1]: + ob_constraints.AddAtomYConstraint(pid) + if dims[2]: + ob_constraints.AddAtomZConstraint(pid) + else: + for particle in p1.particles(): + pid = particle_idx[id(particle)] + 1 # openbabel indices start at 1 + + if all_true: + ob_constraints.AddAtomConstraint(pid) + else: + if dims[0]: + ob_constraints.AddAtomXConstraint(pid) + if dims[1]: + ob_constraints.AddAtomYConstraint(pid) + if dims[2]: + ob_constraints.AddAtomZConstraint(pid) + + if ignore_compounds is not None: + temp1 = np.array(ignore_compounds, dtype=object) + if len(temp1.shape) == 2: + ignore_compounds = list(temp1.reshape(-1)) + + # Since the ignore_compounds can only be passed as a list + # we can check the whole list at once before looping over it + _check_openbabel_constraints( + compound=compound, + particle_list=ignore_compounds, + successors_list=successors_list, + ) + + for ignore in ignore_compounds: + p1 = ignore + if len(p1.children) == 0: + pid = particle_idx[id(p1)] + 1 # openbabel indices start at 1 + ob_constraints.AddIgnore(pid) + + else: + for particle in p1.particles(): + pid = particle_idx[id(particle)] + 1 # openbabel indices start at 1 + ob_constraints.AddIgnore(pid) + + mol = compound.to_pybel() + mol = mol.OBMol + + mol.PerceiveBondOrders() + mol.SetAtomTypesPerceived() + + ff = openbabel.OBForceField.FindForceField(forcefield) + if ff is None: + raise MBuildError( + f"Force field '{forcefield}' not supported for energy " + "minimization. Valid force fields are 'MMFF94', " + "'MMFF94s', 'UFF', 'GAFF', and 'Ghemical'." + "" + ) + warn( + "Performing energy minimization using the Open Babel package. " + "Please refer to the documentation to find the appropriate " + f"citations for Open Babel and the {forcefield} force field" + ) + + if ( + distance_constraints is not None + or fixed_compounds is not None + or ignore_compounds is not None + ): + ob_constraints.SetFactor(constraint_factor) + if ff.Setup(mol, ob_constraints) == 0: + raise MBuildError("Could not setup forcefield for OpenBabel Optimization.") + else: + if ff.Setup(mol) == 0: + raise MBuildError("Could not setup forcefield for OpenBabel Optimization.") + + if algorithm == "steep": + ff.SteepestDescent(steps) + elif algorithm == "md": + ff.MolecularDynamicsTakeNSteps(steps, 300) + elif algorithm == "cg": + ff.ConjugateGradients(steps) + else: + raise MBuildError( + "Invalid minimization algorithm. Valid options are 'steep', 'cg', and 'md'." + ) + ff.UpdateCoordinates(mol) + + # update the coordinates in the Compound + for i, obatom in enumerate(openbabel.OBMolAtomIter(mol)): + x = obatom.GetX() / 10.0 + y = obatom.GetY() / 10.0 + z = obatom.GetZ() / 10.0 + compound[i].pos = np.array([x, y, z]) diff --git a/mbuild/tests/test_box.py b/mbuild/tests/test_box.py index 2f8b7b763..f0cc1cab0 100644 --- a/mbuild/tests/test_box.py +++ b/mbuild/tests/test_box.py @@ -1,3 +1,5 @@ +import logging + import numpy as np import pytest @@ -41,9 +43,10 @@ def test_from_lengths_angles(self, lengths, angles): ], ], ) - def test_left_handed_matrix(self, lh_matrix): - with pytest.warns(UserWarning, match=r"provided for a left\-handed basis"): + def test_left_handed_matrix(self, lh_matrix, caplog): + with caplog.at_level(logging.WARNING, logger="mbuild"): mb.Box.from_vectors(vectors=lh_matrix) + assert "provided for a left-handed basis" in caplog.text @pytest.mark.parametrize( "vecs", diff --git a/mbuild/tests/test_cif.py b/mbuild/tests/test_cif.py index ea36f4c7b..5ba7471f6 100644 --- a/mbuild/tests/test_cif.py +++ b/mbuild/tests/test_cif.py @@ -1,3 +1,4 @@ +import logging from collections import OrderedDict import numpy as np @@ -169,10 +170,11 @@ def test_cif_triclinic_box_properties(self): map(lambda x: x.element, periodic_boxed_molecule.particles()) ) - def test_cif_raise_warnings(self): - with pytest.warns( - UserWarning, - match=r"Element assumed from cif file to be Element: silicon, symbol: Si, atomic number: 14, mass: 28.085.", - ): + def test_cif_raise_warnings(self, caplog): + with caplog.at_level(logging.INFO, logger="mbuild"): lattice_cif = load_cif(file_or_path=get_fn("ETV_triclinic.cif")) lattice_cif.populate(x=1, y=1, z=1) + assert ( + "Element assumed from cif file to be Element: silicon, symbol: Si, atomic number: 14, mass: 28.085." + in caplog.text + ) diff --git a/mbuild/tests/test_compound.py b/mbuild/tests/test_compound.py index 8438725f2..46db00552 100644 --- a/mbuild/tests/test_compound.py +++ b/mbuild/tests/test_compound.py @@ -1,6 +1,8 @@ +import logging import os import sys +import networkx as nx import numpy as np import parmed as pmd import pytest @@ -89,6 +91,20 @@ def test_save_ports(self): assert mol_in.n_particles == 9 + def test_get_child_indices(self): + methane = mb.load("C", smiles=True) + methane2 = mb.clone(methane) + ethane = mb.load("CC", smiles=True) + comp = Compound(subcompounds=[methane, ethane, methane2]) + assert comp.get_child_indices(child=methane) == [0, 1, 2, 3, 4] + assert comp.get_child_indices(child=ethane) == [5, 6, 7, 8, 9, 10, 11, 12] + assert comp.get_child_indices(child=methane2) == [13, 14, 15, 16, 17] + assert methane.get_child_indices(child=methane.children[0]) == [0] + assert ethane.get_child_indices(child=ethane.children[0]) == [0] + + with pytest.raises(ValueError): + ethane.get_child_indices(child=methane.children[0]) + def test_load_xyz(self): myethane = mb.load(get_fn("ethane.xyz")) assert myethane.n_particles == 8 @@ -651,21 +667,23 @@ def test_mass_add_port(self): A.add(mb.Port()) assert A.mass == 2.0 - def test_none_mass(self): + def test_none_mass(self, caplog): A = mb.Compound() assert A.mass is None container = mb.Compound(subcompounds=[A]) - with pytest.warns(UserWarning): + with caplog.at_level(logging.INFO, logger="mbuild"): container_mass = container.mass assert container_mass is None + assert "Some particle of [900, 900, 900] + np.asarray(big_box_of_methane.get_boundingbox().lengths) > [10, 10, 10] ) assert all( - np.asarray(big_sphere_of_methane.get_boundingbox().lengths) - > [1800, 1800, 1800] + np.asarray(big_sphere_of_methane.get_boundingbox().lengths) > [14, 14, 14] ) def test_box_edge(self, h2o, methane): @@ -639,7 +641,7 @@ def test_box_edge(self, h2o, methane): edge_sizes = np.subtract(system_box.lengths, solvated.get_boundingbox().lengths) assert np.allclose(edge_sizes, np.array([0.4] * 3), atol=0.1) - def test_validate_mass(self, methane): + def test_validate_mass(self, methane, caplog): bead = mb.Compound(name="A", mass=0.0) with pytest.raises(MBuildError): mb.fill_box(compound=bead, n_compounds=10, density=1) @@ -653,14 +655,14 @@ def test_validate_mass(self, methane): beadA = mb.Compound(name="A", mass=0.0) beadB = mb.Compound(name="B", mass=1.0, pos=[0.5, 0.5, 0.5]) beads = mb.Compound(subcompounds=[beadA, beadB]) - with warnings.catch_warnings(record=True) as w: + with caplog.at_level(logging.INFO, logger="mbuild"): mb.packing._validate_mass(compound=[beadA, beadB], n_compounds=None) - assert w + assert "Some of the compounds or subcompounds in `compound`" in caplog.text - with warnings.catch_warnings(record=True) as w: + with caplog.at_level(logging.INFO, logger="mbuild"): mb.packing._validate_mass(compound=[beads], n_compounds=None) - assert w + assert "Some of the compounds or subcompounds in `compound`" in caplog.text - with warnings.catch_warnings(record=True) as w: + with caplog.at_level(logging.INFO, logger="mbuild"): mb.packing._validate_mass(compound=[beads], n_compounds=[5]) - assert w + assert "Some of the compounds or subcompounds in `compound`" in caplog.text diff --git a/mbuild/tests/test_path.py b/mbuild/tests/test_path.py new file mode 100644 index 000000000..71484557a --- /dev/null +++ b/mbuild/tests/test_path.py @@ -0,0 +1,139 @@ +import numpy as np + +from mbuild.path import HardSphereRandomWalk +from mbuild.tests.base_test import BaseTest +from mbuild.utils.geometry import bounding_box +from mbuild.utils.volumes import ( + CuboidConstraint, + CylinderConstraint, + SphereConstraint, +) + + +class TestRandomWalk(BaseTest): + def test_random_walk(self): + rw_path = HardSphereRandomWalk( + N=20, + bond_length=0.25, + radius=0.22, + min_angle=np.pi / 4, + max_angle=np.pi, + max_attempts=1e4, + seed=14, + ) + assert len(rw_path.coordinates) == 20 + diffs = rw_path.coordinates[0:-2] - rw_path.coordinates[1:-1] + assert np.allclose(0.25, np.linalg.norm(diffs, axis=1), atol=1e-4) + comp = rw_path.to_compound() + assert comp.n_particles == 20 + assert comp.n_bonds == 19 + # Test bounds of random initial point + assert np.all(rw_path.coordinates[0]) < 20 * 0.22 + + def test_set_initial_point(self): + rw_path = HardSphereRandomWalk( + N=20, + bond_length=0.25, + radius=0.22, + min_angle=np.pi / 4, + max_angle=np.pi, + max_attempts=1e4, + initial_point=(1, 2, 3), + seed=14, + ) + assert np.array_equal(rw_path.coordinates[0], np.array([1, 2, 3])) + + def test_seeds(self): + rw_path_1 = HardSphereRandomWalk( + N=20, + bond_length=0.25, + radius=0.22, + min_angle=np.pi / 4, + max_angle=np.pi, + max_attempts=1e4, + seed=14, + ) + rw_path_2 = HardSphereRandomWalk( + N=20, + bond_length=0.25, + radius=0.22, + min_angle=np.pi / 4, + max_angle=np.pi, + max_attempts=1e4, + seed=14, + ) + assert np.allclose(rw_path_1.coordinates, rw_path_2.coordinates, atol=1e-7) + + def test_from_path(self): + rw_path = HardSphereRandomWalk( + N=20, + bond_length=0.25, + radius=0.22, + min_angle=np.pi / 4, + max_angle=np.pi, + max_attempts=1e4, + seed=14, + ) + rw_path2 = HardSphereRandomWalk( + N=20, + bond_length=0.25, + radius=0.22, + min_angle=np.pi / 4, + max_angle=np.pi, + max_attempts=1e4, + seed=24, + start_from_path=rw_path, + start_from_path_index=-1, + ) + assert len(rw_path.coordinates) == 20 + assert len(rw_path2.coordinates) == 40 + for coord1, coord2 in zip(rw_path.coordinates[:10], rw_path2.coordinates[:10]): + assert np.allclose(coord1, coord2, atol=1e-6) + + def test_walk_inside_cube(self): + cube = CuboidConstraint(Lx=5, Ly=5, Lz=5) + rw_path = HardSphereRandomWalk( + N=500, + bond_length=0.25, + radius=0.22, + volume_constraint=cube, + min_angle=np.pi / 4, + max_angle=np.pi, + max_attempts=1e4, + seed=14, + ) + bounds = bounding_box(rw_path.coordinates) + assert np.all(bounds < np.array([5 - 0.44, 5 - 0.44, 5 - 0.44])) + + def test_walk_inside_sphere(self): + sphere = SphereConstraint(radius=4, center=(2, 2, 2)) + rw_path = HardSphereRandomWalk( + N=200, + bond_length=0.25, + radius=0.22, + volume_constraint=sphere, + initial_point=(0, 0, 0), + min_angle=np.pi / 4, + max_angle=np.pi, + max_attempts=1e4, + seed=90, + ) + bounds = bounding_box(rw_path.coordinates) + assert np.all(bounds < np.array([(2 * 4) - 0.22])) + + def test_walk_inside_cylinder(self): + cylinder = CylinderConstraint(radius=3, height=6, center=(1.5, 1.5, 3)) + rw_path = HardSphereRandomWalk( + N=200, + bond_length=0.25, + radius=0.22, + volume_constraint=cylinder, + min_angle=np.pi / 4, + max_angle=np.pi, + max_attempts=1e4, + seed=14, + ) + bounds = bounding_box(rw_path.coordinates) + assert bounds[0][0] < 6 - 0.22 * 2 + assert bounds[1][1] < 6 - 0.22 * 2 + assert bounds[2][2] < 6 - 0.22 * 2 diff --git a/mbuild/tests/test_polymer.py b/mbuild/tests/test_polymer.py index bcea26d55..80fa9640b 100644 --- a/mbuild/tests/test_polymer.py +++ b/mbuild/tests/test_polymer.py @@ -4,7 +4,7 @@ import pytest import mbuild as mb -from mbuild.lib.recipes import Polymer +from mbuild import Polymer from mbuild.tests.base_test import BaseTest @@ -31,8 +31,8 @@ def test_pass_end_groups(self, ch2, ester): ester_2 = mb.clone(ester) c6 = Polymer(monomers=[ch2], end_groups=[ester, ester_2]) c6.build(n=6) + assert c6.children[0].name == "Ester" assert c6.children[-1].name == "Ester" - assert c6.children[-2].name == "Ester" def test_errors(self, ch2, ester): with pytest.raises(ValueError): # Not enough end groups diff --git a/mbuild/tests/test_port.py b/mbuild/tests/test_port.py index 6415d7bae..c2efe8deb 100644 --- a/mbuild/tests/test_port.py +++ b/mbuild/tests/test_port.py @@ -1,5 +1,6 @@ +import logging + import numpy as np -import pytest import mbuild as mb from mbuild.tests.base_test import BaseTest @@ -50,32 +51,35 @@ def test_port_direction(self): port.rotate(np.pi, [1, 0, 0]) assert np.allclose([0, -1, 0], port.direction, atol=1e-15) - def test_port_separation(self, ethane): + def test_port_separation(self, ethane, caplog): port = mb.Port(anchor=ethane, separation=0.7) assert np.allclose(0.7, port.separation, atol=1e-15) port_no_anchor = mb.Port() - with pytest.warns(UserWarning): + with caplog.at_level(logging.WARNING, logger="mbuild"): separation = port_no_anchor.separation + assert "This port is not anchored to another particle." in caplog.text assert separation is None - def test_update_separation(self, ethane, hexane): + def test_update_separation(self, ethane, hexane, caplog): port = mb.Port(anchor=ethane, separation=0.7) port.update_separation(separation=0.9) assert np.allclose(0.9, port.separation, atol=1e-15) port_used = hexane.labels["propyl2"].labels["down"] - with pytest.warns(UserWarning): + with caplog.at_level(logging.WARNING, logger="mbuild"): port_used.update_separation(0.10) + assert " This port is already being used and changing" in caplog.text - def test_update_orientaiton(self, ch2, hexane): + def test_update_orientaiton(self, ch2, hexane, caplog): port = ch2["up"] port.update_orientation(orientation=(1, 0, 0)) assert np.allclose([-1, 0, 0], port.direction, atol=1e-15) port_used = hexane.labels["propyl2"].labels["down"] - with pytest.warns(UserWarning): + with caplog.at_level(logging.WARNING, logger="mbuild"): port_used.update_orientation([0, 1, 0]) + assert " This port is already being used and changing" in caplog.text def test_access_labels(self): port = mb.Port() diff --git a/mbuild/tests/test_simulation.py b/mbuild/tests/test_simulation.py new file mode 100644 index 000000000..c6f216b62 --- /dev/null +++ b/mbuild/tests/test_simulation.py @@ -0,0 +1,349 @@ +import sys + +import numpy as np +import pytest + +import mbuild as mb +from mbuild.exceptions import MBuildError +from mbuild.simulation import energy_minimize +from mbuild.tests.base_test import BaseTest +from mbuild.utils.io import ( + get_fn, + has_foyer, + has_openbabel, +) + + +class TestSimulation(BaseTest): + def test_energy_minimize(self, octane): + energy_minimize(compound=octane) + + @pytest.mark.skipif(not has_openbabel, reason="Open Babel not installed") + @pytest.mark.skipif( + "win" in sys.platform, reason="Unknown issue with Window's Open Babel " + ) + def test_energy_minimize_shift_com(self, octane): + com_old = octane.pos + energy_minimize(compound=octane) + # check to see if COM of energy minimized Compound + # has been shifted back to the original COM + assert np.allclose(com_old, octane.pos) + + @pytest.mark.skipif(not has_openbabel, reason="Open Babel not installed") + @pytest.mark.skipif( + "win" in sys.platform, reason="Unknown issue with Window's Open Babel " + ) + def test_energy_minimize_shift_anchor(self, octane): + anchor_compound = octane.labels["chain"].labels["CH3[0]"] + pos_old = anchor_compound.pos + energy_minimize(compound=octane, anchor=anchor_compound) + # check to see if COM of the anchor Compound + # has been shifted back to the original COM + assert np.allclose(pos_old, anchor_compound.pos) + + @pytest.mark.skipif(not has_openbabel, reason="Open Babel not installed") + @pytest.mark.skipif( + "win" in sys.platform, reason="Unknown issue with Window's Open Babel " + ) + def test_energy_minimize_fix_compounds(self, octane): + methyl_end0 = octane.labels["chain"].labels["CH3[0]"] + methyl_end1 = octane.labels["chain"].labels["CH3[0]"] + carbon_end = octane.labels["chain"].labels["CH3[0]"].labels["C[0]"] + not_in_compound = mb.Compound(name="H") + + # fix the whole molecule and make sure positions are close + # given stochastic nature and use of restraining springs + # we need to have a pretty loose tolerance for checking + old_com = octane.pos + energy_minimize( + compound=octane, + fixed_compounds=octane, + shift_com=False, + constraint_factor=1e6, + ) + assert np.allclose(octane.pos, old_com, rtol=1e-2, atol=1e-2) + + # primarily focus on checking inputs are parsed correctly + energy_minimize(compound=octane, fixed_compounds=[octane]) + energy_minimize(compound=octane, fixed_compounds=carbon_end) + energy_minimize(compound=octane, fixed_compounds=methyl_end0) + energy_minimize(compound=octane, fixed_compounds=[methyl_end0]) + energy_minimize( + compound=octane, fixed_compounds=[methyl_end0, (True, True, True)] + ) + + energy_minimize( + compound=octane, fixed_compounds=[methyl_end0, (True, True, False)] + ) + energy_minimize( + compound=octane, fixed_compounds=[methyl_end0, [True, True, False]] + ) + energy_minimize( + compound=octane, fixed_compounds=[methyl_end0, (True, False, False)] + ) + energy_minimize( + compound=octane, fixed_compounds=[methyl_end0, (False, False, False)] + ) + + energy_minimize(compound=octane, fixed_compounds=[methyl_end0, methyl_end1]) + energy_minimize(compound=octane, fixed_compounds=[[methyl_end0], [methyl_end1]]) + energy_minimize( + compound=octane, + fixed_compounds=[ + [methyl_end0, (True, True, True)], + [methyl_end1, (True, True, True)], + ], + ) + + with pytest.raises(MBuildError): + energy_minimize(compound=octane, fixed_compounds=not_in_compound) + with pytest.raises(MBuildError): + energy_minimize(compound=octane, fixed_compounds=[not_in_compound]) + with pytest.raises(MBuildError): + energy_minimize( + compound=octane, fixed_compounds=[12323.3, (True, False, False)] + ) + with pytest.raises(Exception): + energy_minimize( + compound=octane, + fixed_compounds=[methyl_end0, (True, False, False, False)], + ) + with pytest.raises(Exception): + energy_minimize( + compound=octane, fixed_compounds=[methyl_end0, True, False, False] + ) + with pytest.raises(Exception): + energy_minimize(compound=octane, fixed_compounds=[methyl_end0, True]) + with pytest.raises(Exception): + energy_minimize( + compound=octane, + fixed_compounds=[methyl_end0, [True, False, False, False]], + ) + with pytest.raises(Exception): + energy_minimize( + compound=octane, fixed_compounds=[methyl_end0, (True, False)] + ) + + with pytest.raises(Exception): + energy_minimize(compound=octane, fixed_compounds=[methyl_end0, (True)]) + + with pytest.raises(Exception): + energy_minimize( + compound=octane, fixed_compounds=[methyl_end0, ("True", True, True)] + ) + with pytest.raises(Exception): + energy_minimize( + compound=octane, fixed_compounds=[methyl_end0, (True, "True", True)] + ) + with pytest.raises(Exception): + energy_minimize( + compound=octane, fixed_compounds=[methyl_end0, (True, True, "True")] + ) + with pytest.raises(Exception): + energy_minimize( + compound=octane, fixed_compounds=[methyl_end0, ("True", True, "True")] + ) + with pytest.raises(Exception): + energy_minimize( + compound=octane, fixed_compounds=[methyl_end0, (True, "True", "True")] + ) + with pytest.raises(Exception): + energy_minimize( + compound=octane, fixed_compounds=[methyl_end0, ("True", "True", True)] + ) + with pytest.raises(Exception): + energy_minimize( + compound=octane, fixed_compounds=[methyl_end0, ("True", "True", "True")] + ) + with pytest.raises(Exception): + energy_minimize( + compound=octane, fixed_compounds=[methyl_end0, (123.0, 231, "True")] + ) + + @pytest.mark.skipif(not has_openbabel, reason="Open Babel not installed") + @pytest.mark.skipif( + "win" in sys.platform, reason="Unknown issue with Window's Open Babel " + ) + def test_energy_minimize_ignore_compounds(self, octane): + methyl_end0 = octane.labels["chain"].labels["CH3[0]"] + methyl_end1 = octane.labels["chain"].labels["CH3[1]"] + carbon_end = octane.labels["chain"].labels["CH3[0]"].labels["C[0]"] + not_in_compound = mb.Compound(name="H") + + # fix the whole molecule and make sure positions are close + # given stochastic nature and use of restraining springs + # we need to have a pretty loose tolerance for checking + old_com = octane.pos + energy_minimize( + compound=octane, + ignore_compounds=octane, + shift_com=False, + constraint_factor=1e6, + ) + assert np.allclose(octane.pos, old_com, rtol=1e-2, atol=1e-2) + + # primarily focus on checking inputs are parsed correctly + energy_minimize(compound=octane, ignore_compounds=[octane]) + energy_minimize(compound=octane, ignore_compounds=carbon_end) + energy_minimize(compound=octane, ignore_compounds=methyl_end0) + energy_minimize(compound=octane, ignore_compounds=[methyl_end0]) + energy_minimize(compound=octane, ignore_compounds=[methyl_end0, methyl_end1]) + energy_minimize( + compound=octane, ignore_compounds=[[methyl_end0], [methyl_end1]] + ) + + with pytest.raises(MBuildError): + energy_minimize(compound=octane, ignore_compounds=not_in_compound) + with pytest.raises(MBuildError): + energy_minimize(compound=octane, ignore_compounds=[1231, 123124]) + + @pytest.mark.skipif(not has_openbabel, reason="Open Babel not installed") + @pytest.mark.skipif( + "win" in sys.platform, reason="Unknown issue with Window's Open Babel " + ) + def test_energy_minimize_distance_constraints(self, octane): + methyl_end0 = octane.labels["chain"].labels["CH3[0]"] + methyl_end1 = octane.labels["chain"].labels["CH3[1]"] + + carbon_end0 = octane.labels["chain"].labels["CH3[0]"].labels["C[0]"] + carbon_end1 = octane.labels["chain"].labels["CH3[1]"].labels["C[0]"] + h_end0 = octane.labels["chain"].labels["CH3[0]"].labels["H[0]"] + + not_in_compound = mb.Compound(name="H") + + # given stochastic nature and use of restraining springs + # we need to have a pretty loose tolerance for checking + energy_minimize( + compound=octane, + distance_constraints=[(carbon_end0, carbon_end1), 0.7], + constraint_factor=1e20, + ) + assert np.allclose( + np.linalg.norm(carbon_end0.pos - carbon_end1.pos), + 0.7, + rtol=1e-2, + atol=1e-2, + ) + + energy_minimize( + compound=octane, distance_constraints=[[(carbon_end0, carbon_end1), 0.7]] + ) + energy_minimize( + compound=octane, + distance_constraints=[ + [(carbon_end0, carbon_end1), 0.7], + [(carbon_end0, h_end0), 0.1], + ], + ) + + with pytest.raises(MBuildError): + energy_minimize( + compound=octane, + distance_constraints=[(carbon_end0, not_in_compound), 0.7], + ) + with pytest.raises(MBuildError): + energy_minimize( + compound=octane, distance_constraints=[(carbon_end0, carbon_end0), 0.7] + ) + with pytest.raises(MBuildError): + energy_minimize( + compound=octane, distance_constraints=[(methyl_end0, carbon_end1), 0.7] + ) + with pytest.raises(MBuildError): + energy_minimize( + compound=octane, distance_constraints=[(methyl_end0, methyl_end1), 0.7] + ) + + @pytest.mark.skipif(has_openbabel, reason="Open Babel package is installed") + @pytest.mark.skipif( + "win" in sys.platform, reason="Unknown issue with Window's Open Babel " + ) + def test_energy_minimize_openbabel_warn(self, octane): + with pytest.raises(MBuildError): + energy_minimize(compound=octane) + + @pytest.mark.skipif(not has_openbabel, reason="Open Babel not installed") + @pytest.mark.skipif( + "win" in sys.platform, reason="Unknown issue with Window's Open Babel " + ) + def test_energy_minimize_ff(self, octane): + for ff in ["UFF", "GAFF", "MMFF94", "MMFF94s", "Ghemical"]: + energy_minimize(compound=octane, forcefield=ff) + with pytest.raises(IOError): + energy_minimize(compound=octane, forcefield="fakeFF") + + @pytest.mark.skipif(not has_openbabel, reason="Open Babel not installed") + @pytest.mark.skipif( + "win" in sys.platform, reason="Unknown issue with Window's Open Babel " + ) + def test_energy_minimize_algorithm(self, octane): + for algorithm in ["cg", "steep", "md"]: + energy_minimize(compound=octane, algorithm=algorithm) + with pytest.raises(MBuildError): + energy_minimize(compound=octane, algorithm="fakeAlg") + + @pytest.mark.skipif(not has_openbabel, reason="Open Babel not installed") + @pytest.mark.skipif( + "win" in sys.platform, reason="Unknown issue with Window's Open Babel " + ) + def test_energy_minimize_non_element(self, octane): + for particle in octane.particles(): + particle.element = None + # Pass with element inference from names + energy_minimize(compound=octane) + for particle in octane.particles(): + particle.name = "Q" + particle.element = None + + # Fail once names cannot be set as elements + with pytest.raises(MBuildError): + energy_minimize(compound=octane) + + @pytest.mark.skipif(not has_openbabel, reason="Open Babel not installed") + @pytest.mark.skipif( + "win" in sys.platform, reason="Unknown issue with Window's Open Babel " + ) + def test_energy_minimize_ports(self, octane): + distances = np.round( + [ + octane.min_periodic_distance(port.pos, port.anchor.pos) + for port in octane.all_ports() + ], + 5, + ) + orientations = np.round( + [port.pos - port.anchor.pos for port in octane.all_ports()], 5 + ) + + energy_minimize(compound=octane) + + updated_distances = np.round( + [ + octane.min_periodic_distance(port.pos, port.anchor.pos) + for port in octane.all_ports() + ], + 5, + ) + updated_orientations = np.round( + [port.pos - port.anchor.pos for port in octane.all_ports()], 5 + ) + + assert np.array_equal(distances, updated_distances) + assert np.array_equal(orientations, updated_orientations) + + @pytest.mark.skipif(not has_foyer, reason="Foyer is not installed") + def test_energy_minimize_openmm(self, octane): + energy_minimize(compound=octane, forcefield="oplsaa") + + @pytest.mark.skipif(not has_foyer, reason="Foyer is not installed") + @pytest.mark.parametrize("constraints", ["AllBonds", "HBonds", "HAngles", None]) + def test_energy_minimize_openmm_constraints(self, octane, constraints): + energy_minimize(compound=octane, forcefield="oplsaa", constraints=constraints) + + def test_energy_minimize_openmm_invalid_constraints(self, octane): + with pytest.raises(ValueError): + energy_minimize(compound=octane, forcefield="oplsaa", constraints="boo") + + @pytest.mark.skipif(not has_foyer, reason="Foyer is not installed") + def test_energy_minimize_openmm_xml(self, octane): + energy_minimize(compound=octane, forcefield=get_fn("small_oplsaa.xml")) diff --git a/mbuild/tests/test_utils.py b/mbuild/tests/test_utils.py index ab61a4894..dd94cf66b 100644 --- a/mbuild/tests/test_utils.py +++ b/mbuild/tests/test_utils.py @@ -1,4 +1,5 @@ import difflib +import logging import numpy as np import pytest @@ -205,10 +206,10 @@ def test_RB_to_OPLS_c5_not_0(self): c5 = 0.3 RB_to_OPLS(c0, c1, c2, c3, c4, c5) - def test_RB_to_OPLS_f0_not_0_within_tolerance_warn(self): + def test_RB_to_OPLS_f0_not_0_within_tolerance_warn(self, caplog): # should throw a warning that f0 is not zero text_for_error_tolerance = ( - "f0 \= 2 \* \( c0 \+ c1 \+ c2 \+ c3 \+ c4 \+ c5 \) is not zero. " + "f0 = 2 * ( c0 + c1 + c2 + c3 + c4 + c5 ) is not zero. " "The f0/2 term is the constant for the OPLS dihedral. " "Since the f0 term is not zero, the dihedral is not an " "exact conversion; since this constant does not contribute " @@ -216,8 +217,7 @@ def test_RB_to_OPLS_f0_not_0_within_tolerance_warn(self): "for MD, but the energy for each dihedral will be shifted " "by the f0/2 value." ) - - with pytest.warns(UserWarning, match=f"{text_for_error_tolerance}"): + with caplog.at_level(logging.WARNING, logger="mbuild"): c0 = 0.4 c1 = 0.4 c2 = -0.1 @@ -225,6 +225,7 @@ def test_RB_to_OPLS_f0_not_0_within_tolerance_warn(self): c4 = -0.2 c5 = 0 RB_to_OPLS(c0, c1, c2, c3, c4, c5, error_if_outside_tolerance=False) + assert text_for_error_tolerance in caplog.text def test_RB_to_OPLS_f0_not_0_within_tolerance_error(self): text_for_error_tolerance = ( @@ -350,7 +351,7 @@ def test_OPLS_to_RB_error_tolerance_not_float(self): f4 = 0.2 OPLS_to_RB(f0, f1, f2, f3, f4, error_tolerance="s") - def test_OPLS_to_RB_f0_is_zero(self): + def test_OPLS_to_RB_f0_is_zero(self, caplog): # should throw a warning that f0 is zero text_for_error_tolerance = ( "The f0/2 term is the constant for the OPLS dihedral equation, " @@ -361,10 +362,11 @@ def test_OPLS_to_RB_f0_is_zero(self): "this should provide matching results for MD, but the energy for each" "dihedral will be shifted by the real f0/2 value." ) - with pytest.warns(UserWarning, match=text_for_error_tolerance): + with caplog.at_level(logging.WARNING, logger="mbuild"): f0 = 0 f1 = 0.1 f2 = -0.2 f3 = -0.1 f4 = 0.2 OPLS_to_RB(f0, f1, f2, f3, f4) + assert text_for_error_tolerance in caplog.text diff --git a/mbuild/tests/test_vasp.py b/mbuild/tests/test_vasp.py index 383a070c4..217a3bd9d 100644 --- a/mbuild/tests/test_vasp.py +++ b/mbuild/tests/test_vasp.py @@ -1,3 +1,5 @@ +import logging + import numpy as np import pytest @@ -74,10 +76,11 @@ def test_coordinate_header(self, gilmerite, coord_type): assert coord == coord_type - def test_warning_raised(self, copper_cell): + def test_warning_raised(self, copper_cell, caplog): copper_cell.box = None - with pytest.warns(UserWarning): + with caplog.at_level(logging.INFO, logger="mbuild"): write_poscar(copper_cell, "test.poscar", coord_style="direct") + assert "'direct' coord_style specified, but compound has no box" in caplog.text def test_error_raised(self, copper_cell): with pytest.raises(ValueError): diff --git a/mbuild/utils/conversion.py b/mbuild/utils/conversion.py index 0703ffa11..ae11f99f3 100644 --- a/mbuild/utils/conversion.py +++ b/mbuild/utils/conversion.py @@ -1,9 +1,11 @@ """mBuild conversion utilities.""" -from warnings import warn +import logging import numpy as np +logger = logging.getLogger(__name__) + def RB_to_OPLS( c0, @@ -92,7 +94,7 @@ def RB_to_OPLS( if error_if_outside_tolerance is True: raise ValueError(text_for_error_tolerance) elif error_if_outside_tolerance is False: - warn(text_for_error_tolerance) + logger.warning(text_for_error_tolerance) f1 = -2 * c1 - (3 * c3) / 2 f2 = -c2 - c4 @@ -147,7 +149,7 @@ def OPLS_to_RB(f0, f1, f2, f3, f4, error_tolerance=1e-4): error_tolerance = abs(error_tolerance) if np.all(np.isclose(f0 / 2, 0, atol=error_tolerance, rtol=0)): - warn( + logger.warning( "The f0/2 term is the constant for the OPLS dihedral equation, " "which is added to a constant for the RB torsions equation via the c0 coefficient. " "The f0 term is zero in the OPLS dihedral form or is force set to zero in this equation, " diff --git a/mbuild/utils/decorators.py b/mbuild/utils/decorators.py index b48e76268..79da35e60 100644 --- a/mbuild/utils/decorators.py +++ b/mbuild/utils/decorators.py @@ -1,7 +1,9 @@ """Some helpful decorators.""" +import logging from functools import wraps -from warnings import warn + +logger = logging.getLogger(__name__) def make_callable(decorator): # pragma: no cover @@ -27,7 +29,7 @@ def deprecated(warning_string=""): # pragma: no cover def old_function(fcn): def wrapper(*args, **kwargs): printed_message = f"{fcn.__name__} is deprecated. {warning_string}" - warn(printed_message, DeprecationWarning) + logger.warning(printed_message) return fcn(*args, **kwargs) return wrapper @@ -95,9 +97,8 @@ def deleter(self, fdel): def show_warning_msg(self): if self.warning_msg is not None: - warn( - f"Property {fcn.__name__} is deprecated. {self.warning_msg}", - DeprecationWarning, + logger.warning( + f"Property {fcn.__name__} is deprecated. {self.warning_msg}" ) if not self.always_show: self.warning_msg = None @@ -113,7 +114,7 @@ def breaking_change(warning_string=""): def old_function(fcn): def wrapper(*args, **kwargs): - warn(f"{fcn.__name__} has breaking change. {warning_string}", Warning) + logger.warning(f"{fcn.__name__} has breaking change. {warning_string}") fcn(*args, **kwargs) return wrapper @@ -127,9 +128,8 @@ def experimental_feature(warning_string=""): # pragma no cover def experimental_function(fcn): @wraps(fcn) def wrapper(*args, **kwargs): - warn( - f"{fcn.__name__} is an experimental feature and is not subject to follow standard deprecation cycles. Use at your own risk! {warning_string}", - Warning, + logger.warning( + f"{fcn.__name__} is an experimental feature and is not subject to follow standard deprecation cycles. Use at your own risk! {warning_string}" ) fcn(*args, **kwargs) diff --git a/mbuild/utils/density.py b/mbuild/utils/density.py new file mode 100644 index 000000000..b57576b1a --- /dev/null +++ b/mbuild/utils/density.py @@ -0,0 +1,99 @@ +import numpy as np +from scipy.spatial import cKDTree + + +def rank_points_by_density(points, k=14, box_min=None, box_max=None, wall_cutoff=0.0): + """Rank points from lowest local density to highest. + + Notes + ----- + If box_min, box_max and wall_cutoff are given, then points within wall_cutoff + to the box boundaries are ignored. This ensure the effective "density" resulting + from boundaries is acounted for. + + This can be useful for stringing together mutliple random walks + with and without volume constraints. The use case for this function would be to find an existing site + from a random walk, or another path, that has relatively lower density than other sites in the path. + A new random walk path can begin from this site. + + If you want to find an unoccupied point with the lowest local density + use `mbuild.utils.density.find_low_density_point` instead. + + Parameters + ---------- + points : ndarray, shape (N, d) + Coordinates of points. + k : int, optional, default 14 + Number of neighbors to use for local density. + box_min : array-like, shape (d,) + Minimum boundary of box. Points closer than `wall_cutoff` will be ignored. + box_max : array-like, shape (d,) + Maximum boundary of box. + wall_cutoff : float (nm) + Distance from walls to ignore points. + + Returns + ------- + sorted_indices : ndarray + Indices of points (original array) sorted from lowest to highest local density. + """ + points = np.asarray(points) + N, dim = points.shape + # Check for points too close to boundaries (if given) + # Ignore these from density calculation, enforce high density due to proximity to wall + if box_min is not None and box_max is not None and wall_cutoff > 0.0: + box_min = np.asarray(box_min) + box_max = np.asarray(box_max) + mask = np.all( + (points > box_min + wall_cutoff) & (points < box_max - wall_cutoff), axis=1 + ) + valid_indices = np.nonzero(mask)[0] + filtered_points = points[mask] + else: + filtered_points = points + valid_indices = np.arange(N) + + tree = cKDTree(filtered_points) + dists, idxs = tree.query(filtered_points, k=k + 1) + avg_dist = np.mean(dists[:, 1:], axis=1) + local_density = 1 / (avg_dist + 1e-12) + # Sort indices by increasing density + sorted_order = np.argsort(local_density) + return valid_indices[sorted_order] + + +def find_low_density_point(points, box_min, box_max, edge_buffer=0, n_candidates=5000): + """Find an unoccupied point inside a box that is furthest from existing points. + + Parameters + ---------- + points : ndarray, shape (N, 3) + Array of existing points. + box_min : array-like, shape (d, 3) + Minimum coordinates of the box. + box_max : array-like, shape (d, 3) + Maximum coordinates of the box. + edge_buffer : float (nm) default 0.0 + The buffer to prevent selection of coordinates within + the buffer distance to the edge of the box. + n_candidates : int + Number of random candidate points to try. + These points will be used to chose the point with lowest local density. + Higher values may suffer performance costs, but give improved + sampling. + + Returns + ------- + best_point : ndarray, shape (d,3) + Coordinates of the lowest-density point. + """ + points = np.asarray(points) + dim = points.shape[1] + tree = cKDTree(points) + # Create random candidates inside the box to test and sample from + candidates = np.random.uniform( + box_min + edge_buffer, box_max - edge_buffer, size=(n_candidates, dim) + ) + dists, _ = tree.query(candidates, k=1) + sorted_order = np.argsort(dists) + return candidates[sorted_order] diff --git a/mbuild/utils/geometry.py b/mbuild/utils/geometry.py index fad27db2d..157c48227 100644 --- a/mbuild/utils/geometry.py +++ b/mbuild/utils/geometry.py @@ -2,7 +2,7 @@ import numpy as np -import mbuild as mb +from mbuild.box import Box from mbuild.coordinate_transform import angle @@ -28,6 +28,47 @@ def calc_dihedral(point1, point2, point3, point4): return angle(x, y) +def bounding_box(xyz): + """Find the bounding box from a set of coordinates.""" + # case where only 1 particle exists + is_one_particle = False + if xyz.shape[0] == 1: + is_one_particle = True + + # are any columns all equalivalent values? + # an example of this would be a planar molecule + # example: all z values are 0.0 + # from: https://stackoverflow.com/a/14860884 + # steps: create mask array comparing first value in each column + # use np.all with axis=0 to do row columnar comparision + has_dimension = [True, True, True] + if not is_one_particle: + missing_dimensions = np.all( + np.isclose(xyz, xyz[0, :], atol=1e-2), + axis=0, + ) + for i, truthy in enumerate(missing_dimensions): + has_dimension[i] = not truthy + + if is_one_particle: + v1 = np.asarray([[1.0, 0.0, 0.0]]) + v2 = np.asarray([[0.0, 1.0, 0.0]]) + v3 = np.asarray([[0.0, 0.0, 1.0]]) + else: + maxs = xyz.max(axis=0) + mins = xyz.min(axis=0) + v1 = np.asarray((maxs[0] - mins[0], 0.0, 0.0)) + v2 = np.asarray((0.0, maxs[1] - mins[1], 0.0)) + v3 = np.asarray((0.0, 0.0, maxs[2] - mins[2])) + vecs = [v1, v2, v3] + + # handle any missing dimensions (planar molecules) + for i, dim in enumerate(has_dimension): + if not dim: + vecs[i][i] = 0.1 + return vecs + + def coord_shift(xyz, box): """Ensure that coordinates are -L/2, L/2. @@ -83,7 +124,7 @@ def wrap_coords(xyz, box, mins=None): ----- Currently only supports orthorhombic boxes """ - if not isinstance(box, mb.Box): + if not isinstance(box, Box): box_arr = np.asarray(box) assert box_arr.shape == (3,) diff --git a/mbuild/utils/io.py b/mbuild/utils/io.py index 9833a50c7..ee905d26a 100644 --- a/mbuild/utils/io.py +++ b/mbuild/utils/io.py @@ -22,13 +22,15 @@ import importlib import inspect +import logging import os import sys import textwrap -import warnings from importlib.resources import files from unittest import SkipTest +logger = logging.getLogger(__name__) + class DelayImportError(ImportError, SkipTest): """Error to allow better import handling.""" @@ -88,6 +90,9 @@ class DelayImportError(ImportError, SkipTest): # conda install -c conda-forge openbabel +NOTE: openbabel is only available for python<3.13. +If you need it in your environment, make sure your python is 3.10, 3.11 or 3.12. + or from source following instructions at: # http://openbabel.org/docs/current/UseTheLibrary/PythonInstall.html @@ -181,7 +186,7 @@ def import_(module): "openbabel 2.0 detected and will be dropped in a future " "release. Consider upgrading to 3.x." ) - warnings.warn(msg, DeprecationWarning) + logger.info(msg, DeprecationWarning) return pybel except ModuleNotFoundError: pass @@ -196,7 +201,7 @@ def import_(module): "openbabel 2.0 detected and will be dropped in a future " "release. Consider upgrading to 3.x." ) - warnings.warn(msg, DeprecationWarning) + logger.info(msg, DeprecationWarning) return openbabel except ModuleNotFoundError: pass diff --git a/mbuild/utils/volumes.py b/mbuild/utils/volumes.py new file mode 100644 index 000000000..e608b649b --- /dev/null +++ b/mbuild/utils/volumes.py @@ -0,0 +1,114 @@ +import numpy as np +from numba import njit + + +class Constraint: + def __init__(self): + pass + + def is_inside(self, point): + """Return True if point satisfies constraint (inside), else False.""" + raise NotImplementedError("Must be implemented in subclasses") + + +class CuboidConstraint(Constraint): + def __init__(self, Lx, Ly, Lz, center=(0, 0, 0)): + self.center = np.asarray(center) + self.mins = self.center - np.array([Lx / 2, Ly / 2, Lz / 2]) + self.maxs = self.center + np.array([Lx / 2, Ly / 2, Lz / 2]) + + def is_inside(self, points, buffer): + """Points and buffer are passed in from HardSphereRandomWalk""" + return is_inside_cuboid( + mins=self.mins, maxs=self.maxs, points=points, buffer=buffer + ) + + +class SphereConstraint(Constraint): + def __init__(self, center, radius): + self.center = np.array(center) + self.radius = radius + self.mins = self.center - self.radius + self.maxs = self.center + self.radius + + def is_inside(self, points, buffer): + """Points and buffer are passed in from HardSphereRandomWalk""" + return is_inside_sphere(points=points, sphere_radius=self.radius, buffer=buffer) + + +class CylinderConstraint(Constraint): + def __init__(self, center, radius, height): + self.center = np.array(center) + self.height = height + self.radius = radius + self.mins = np.array( + [ + self.center[0] - self.radius, + self.center[1] - self.radius, + self.center[2] - self.height / 2, + ] + ) + self.maxs = np.array( + [ + self.center[0] + self.radius, + self.center[1] + self.radius, + self.center[2] + self.height / 2, + ] + ) + + def is_inside(self, points, buffer): + """Points and buffer are passed in from HardSphereRandomWalk""" + return is_inside_cylinder( + points=points, + center=self.center, + cylinder_radius=self.radius, + height=self.height, + buffer=buffer, + ) + + +@njit(cache=True, fastmath=True) +def is_inside_cylinder(points, center, cylinder_radius, height, buffer): + n_points = points.shape[0] + results = np.empty(n_points, dtype=np.bool_) + max_r = cylinder_radius - buffer + max_r_sq = max_r * max_r + half_height = height / 2.0 + max_z = half_height - buffer + for i in range(n_points): + dx = points[i, 0] - center[0] + dy = points[i, 1] - center[1] + dz = points[i, 2] - center[2] + r_sq = dx * dx + dy * dy + inside_radial = r_sq <= max_r_sq + inside_z = abs(dz) <= max_z + results[i] = inside_radial and inside_z + return results + + +@njit(cache=True, fastmath=True) +def is_inside_sphere(sphere_radius, points, buffer): + n_points = points.shape[0] + results = np.empty(n_points, dtype=np.bool_) + max_distance = sphere_radius - buffer + max_distance_sq = max_distance * max_distance + for i in range(n_points): + dist_from_center_sq = 0.0 + for j in range(3): + dist_from_center_sq += points[i, j] * points[i, j] + results[i] = dist_from_center_sq < max_distance_sq + return results + + +@njit(cache=True, fastmath=True) +def is_inside_cuboid(mins, maxs, points, buffer): + n_points = points.shape[0] + results = np.empty(n_points, dtype=np.bool_) + for i in range(n_points): + inside = True + for j in range(3): + if points[i, j] - buffer < mins[j] or points[i, j] + buffer > maxs[j]: + inside = False + break + results[i] = inside + return results diff --git a/pyproject.toml b/pyproject.toml index 9d51a1c21..f4832fa20 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -29,7 +29,7 @@ classifiers=[ "Operating System :: MacOS", ] urls = {Homepage = "https://github.com/mosdef-hub/mbuild"} -requires-python = ">=3.9" +requires-python = ">=3.10" dynamic = ["version"] [tool.setuptools] @@ -49,6 +49,5 @@ version = {attr = "mbuild.__version__"} [project.entry-points."mbuild.plugins"] Alkane = "mbuild.lib.recipes.alkane:Alkane" Monolayer = "mbuild.lib.recipes.monolayer:Monolayer" -Polymer = "mbuild.lib.recipes.polymer:Polymer" SilicaInterface = "mbuild.lib.recipes.silica_interface:SilicaInterface" TiledCompound = "mbuild.lib.recipes.tiled_compound:TiledCompound" diff --git a/setup.cfg b/setup.cfg index f0dc4f89d..35e3d25a8 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,24 +1,24 @@ [bumpversion] -current_version = 1.2.0 +current_version = 1.2.1 commit = True tag = True message = Bump to version {new_version} tag_name = {new_version} [coverage:run] -omit = +omit = mbuild/examples/* mbuild/tests/* [coverage:report] -exclude_lines = +exclude_lines = pragma: no cover - + if 0: if __name__ == .__main__.: def __repr__ except ImportError -omit = +omit = mbuild/examples/* mbuild/tests/*