diff --git a/mbuild/packing.py b/mbuild/packing.py index 5d1f08e0b..a0297f574 100755 --- a/mbuild/packing.py +++ b/mbuild/packing.py @@ -1,7 +1,9 @@ from __future__ import division import sys +import os import tempfile +import warnings from distutils.spawn import find_executable from subprocess import Popen, PIPE @@ -37,7 +39,9 @@ """ -def fill_box(compound, n_compounds=None, box=None, density=None, overlap=0.2, seed=12345): +def fill_box(compound, n_compounds=None, box=None, density=None, overlap=0.2, + seed=12345, edge=0.2, compound_ratio=None, + aspect_ratio=None, temp_file=None): """Fill a box with a compound using packmol. Two arguments of `n_compounds, box, and density` must be specified. @@ -52,27 +56,45 @@ def fill_box(compound, n_compounds=None, box=None, density=None, overlap=0.2, se If `box` and `density` are not None, the corresponding number of compounds will be calculated internally. + For the cases in which `box` is not specified but generated internally, + the default behavior is to calculate a cubic box. Optionally, + `aspect_ratio` can be passed to generate a non-cubic box. + Parameters ---------- compound : mb.Compound or list of mb.Compound + Compound or list of compounds to be put in box. n_compounds : int or list of int + Number of compounds to be put in box. box : mb.Box - overlap : float, units nm - density : float, units kg/m^3 + Box to be filled by compounds. + density : float, units kg/m^3, default=None + Target density for the system in macroscale units. If not None, one of + `n_compounds` or `box`, but not both, must be specified. + overlap : float, units nm, default=0.2 + Minimum separation between atoms of different molecules. + seed : int, default=12345 + Random seed to be passed to PACKMOL. + edge : float, units nm, default=0.2 + Buffer at the edge of the box to not place molecules. This is necessary + in some systems because PACKMOL does not account for periodic boundary + conditions in its optimization. + compound_ratio : list, default=None + Ratio of number of each compound to be put in box. Only used in the + case of `density` and `box` having been specified, `n_compounds` not + specified, and more than one `compound`. + aspect_ratio : list of float + If a non-cubic box is desired, the ratio of box lengths in the x, y, + and z directions. + temp_file : str, default=None + File name to write PACKMOL's raw output to. Returns ------- filled : mb.Compound - TODO : Support aspect ratios in generated boxes - TODO : Support ratios of n_compounds """ - if not PACKMOL: - msg = "Packmol not found." - if sys.platform.startswith("win"): - msg = (msg + " If packmol is already installed, make sure that the " - "packmol.exe is on the path.") - raise IOError(msg) + _check_packmol(PACKMOL) arg_count = 3 - [n_compounds, box, density].count(None) if arg_count != 2: @@ -87,28 +109,53 @@ def fill_box(compound, n_compounds=None, box=None, density=None, overlap=0.2, se if n_compounds is not None and not isinstance(n_compounds, (list, set)): n_compounds = [n_compounds] + if compound is not None and n_compounds is not None: + if len(compound) != len(n_compounds): + msg = ("`compound` and `n_compounds` must be of equal length.") + raise ValueError(msg) + if density is not None: if box is None and n_compounds is not None: total_mass = np.sum([n*np.sum([a.mass for a in c.to_parmed().atoms]) for c,n in zip(compound, n_compounds)]) # Conversion from (amu/(kg/m^3))**(1/3) to nm L = (total_mass/density)**(1/3)*1.1841763 - box = _validate_box(Box(3*[L])) - if n_compounds is None and box is not None: - if len(compound) > 1: - msg = ("Determing `n_compounds` from `density` and `box` " - "currently only supported for systems with one " - "compound type.") - raise ValueError(msg) + if aspect_ratio is None: + box = _validate_box(Box(3*[L])) else: + L *= np.prod(aspect_ratio) ** (-1/3) + box = _validate_box(Box([val*L for val in aspect_ratio])) + if n_compounds is None and box is not None: + if len(compound) == 1: compound_mass = np.sum([a.mass for a in compound[0].to_parmed().atoms]) # Conversion from kg/m^3 / amu * nm^3 to dimensionless units n_compounds = [int(density/compound_mass*np.prod(box.lengths)*.60224)] + else: + if compound_ratio is None: + msg = ("Determing `n_compounds` from `density` and `box` " + "for systems with more than one compound type requires" + "`compound_ratio`") + raise ValueError(msg) + if len(compound) != len(compound_ratio): + msg = ("Length of `compound_ratio` must equal length of " + "`compound`") + raise ValueError(msg) + prototype_mass = 0 + for c, r in zip(compound, compound_ratio): + prototype_mass += r * np.sum([a.mass for a in c.to_parmed().atoms]) + # Conversion from kg/m^3 / amu * nm^3 to dimensionless units + n_prototypes = int(density/prototype_mass*np.prod(box.lengths)*.60224) + n_compounds = list() + for c in compound_ratio: + n_compounds.append(int(n_prototypes * c)) # In angstroms for packmol. box_mins = box.mins * 10 box_maxs = box.maxs * 10 overlap *= 10 + + # Apply edge buffer + box_maxs -= edge * 10 # Build the input file for each compound and call packmol. filled_pdb = tempfile.mkstemp(suffix='.pdb')[1] @@ -122,10 +169,7 @@ def fill_box(compound, n_compounds=None, box=None, density=None, overlap=0.2, se box_mins[0], box_mins[1], box_mins[2], box_maxs[0], box_maxs[1], box_maxs[2]) - proc = Popen(PACKMOL, stdin=PIPE, stdout=PIPE, stderr=PIPE, universal_newlines=True) - out, err = proc.communicate(input=input_text) - if err: - _packmol_error(out, err) + _run_packmol(input_text, filled_pdb, temp_file) # Create the topology and update the coordinates. filled = Compound() @@ -137,15 +181,28 @@ def fill_box(compound, n_compounds=None, box=None, density=None, overlap=0.2, se return filled -def fill_region(compound, n_compounds, region, overlap=0.2, seed=12345): +def fill_region(compound, n_compounds, region, overlap=0.2, + seed=12345, edge=0.2, temp_file=None): """Fill a region of a box with a compound using packmol. Parameters ---------- compound : mb.Compound or list of mb.Compound + Compound or list of compounds to be put in region. n_compounds : int or list of int + Number of compounds to be put in region. region : mb.Box or list of mb.Box - overlap : float + Region to be filled by compounds. + overlap : float, units nm, default=0.2 + Minimum separation between atoms of different molecules. + seed : int, default=12345 + Random seed to be passed to PACKMOL. + edge : float, units nm, default=0.2 + Buffer at the edge of the region to not place molecules. This is + necessary in some systems because PACKMOL does not account for + periodic boundary conditions in its optimization. + temp_file : str, default=None + File name to write PACKMOL's raw output to. Returns ------- @@ -154,17 +211,18 @@ def fill_region(compound, n_compounds, region, overlap=0.2, seed=12345): If using mulitple regions and compounds, the nth value in each list are used in order. For example, if the third compound will be put in the third region using the third value in n_compounds. """ - if not PACKMOL: - msg = "Packmol not found." - if sys.platform.startswith("win"): - msg = (msg + " If packmol is already installed, make sure that the " - "packmol.exe is on the path.") - raise IOError(msg) + _check_packmol(PACKMOL) if not isinstance(compound, (list, set)): compound = [compound] if not isinstance(n_compounds, (list, set)): n_compounds = [n_compounds] + + if compound is not None and n_compounds is not None: + if len(compound) != len(n_compounds): + msg = ("`compound` and `n_compounds` must be of equal length.") + raise ValueError(msg) + # See if region is a single region or list if isinstance(region, Box): # Cannot iterate over boxes region = [region] @@ -185,14 +243,12 @@ def fill_region(compound, n_compounds, region, overlap=0.2, seed=12345): comp.save(compound_pdb, overwrite=True) reg_mins = reg.mins * 10 reg_maxs = reg.maxs * 10 + reg_maxs -= edge * 10 # Apply edge buffer input_text += PACKMOL_BOX.format(compound_pdb, m_compounds, reg_mins[0], reg_mins[1], reg_mins[2], reg_maxs[0], reg_maxs[1], reg_maxs[2]) - proc = Popen(PACKMOL, stdin=PIPE, stdout=PIPE, stderr=PIPE, universal_newlines=True) - out, err = proc.communicate(input=input_text) - if err: - _packmol_error(out, err) + _run_packmol(input_text, filled_pdb, temp_file) # Create the topology and update the coordinates. filled = Compound() @@ -203,24 +259,37 @@ def fill_region(compound, n_compounds, region, overlap=0.2, seed=12345): return filled -def solvate(solute, solvent, n_solvent, box, overlap=0.2, seed=12345): +def solvate(solute, solvent, n_solvent, box, overlap=0.2, + seed=12345, edge=0.2, temp_file=None): """Solvate a compound in a box of solvent using packmol. Parameters ---------- solute : mb.Compound + Compound to be placed in a box and solvated. solvent : mb.Compound + Compound to solvate the box. n_solvent : int + Number of solvents to be put in box. box : mb.Box - overlap : float + Box to be filled by compounds. + overlap : float, units nm, default=0.2 + Minimum separation between atoms of different molecules. + seed : int, default=12345 + Random seed to be passed to PACKMOL. + edge : float, units nm, default=0.2 + Buffer at the edge of the box to not place molecules. This is necessary + in some systems because PACKMOL does not account for periodic boundary + conditions in its optimization. + temp_file : str, default=None + File name to write PACKMOL's raw output to. Returns ------- solvated : mb.Compound """ - if not PACKMOL: - raise IOError("Packmol not found") + _check_packmol(PACKMOL) box = _validate_box(box) if not isinstance(solvent, (list, set)): @@ -228,11 +297,18 @@ def solvate(solute, solvent, n_solvent, box, overlap=0.2, seed=12345): if not isinstance(n_solvent, (list, set)): n_solvent = [n_solvent] + if len(solvent) != len(n_solvent): + msg = ("`n_solvent` and `n_solvent` must be of equal length.") + raise ValueError(msg) + # In angstroms for packmol. box_mins = box.mins * 10 box_maxs = box.maxs * 10 overlap *= 10 center_solute = (box_maxs + box_mins) / 2 + + # Apply edge buffer + box_maxs -= edge * 10 # Build the input file for each compound and call packmol. solvated_pdb = tempfile.mkstemp(suffix='.pdb')[1] @@ -249,10 +325,7 @@ def solvate(solute, solvent, n_solvent, box, overlap=0.2, seed=12345): box_mins[0], box_mins[1], box_mins[2], box_maxs[0], box_maxs[1], box_maxs[2]) - proc = Popen(PACKMOL, stdin=PIPE, stdout=PIPE, stderr=PIPE, universal_newlines=True) - out, err = proc.communicate(input=input_text) - if err: - _packmol_error(out, err) + _run_packmol(input_text, solvated_pdb, temp_file) # Create the topology and update the coordinates. solvated = Compound() @@ -284,3 +357,27 @@ def _packmol_error(out, err): log_file.write(out) err_file.write(err) raise RuntimeError("PACKMOL failed. See 'err.txt' and 'log.txt'") + +def _run_packmol(input_text, filled_pdb, temp_file): + proc = Popen(PACKMOL, stdin=PIPE, stdout=PIPE, stderr=PIPE, universal_newlines=True) + out, err = proc.communicate(input=input_text) + + if 'WITHOUT PERFECT PACKING' in out: + msg = ("Packmol finished with imperfect packing. Using " + "the .pdb_FORCED file instead. This may not be a " + "sufficient packing result.") + warnings.warn(msg) + os.system('cp {0}_FORCED {0}'.format(filled_pdb)) + if 'ERROR' in out: + _packmol_error(out, err) + + if temp_file is not None: + os.system('cp {0} {1}'.format(filled_pdb, os.path.join(temp_file))) + +def _check_packmol(PACKMOL): + if not PACKMOL: + msg = "Packmol not found." + if sys.platform.startswith("win"): + msg = (msg + " If packmol is already installed, make sure that the " + "packmol.exe is on the path.") + raise IOError(msg) diff --git a/mbuild/tests/test_packing.py b/mbuild/tests/test_packing.py index 251782fef..5c5218b01 100755 --- a/mbuild/tests/test_packing.py +++ b/mbuild/tests/test_packing.py @@ -1,3 +1,5 @@ +import os + import pytest import numpy as np @@ -16,11 +18,24 @@ def test_fill_box_density_box(self, h2o): filled = mb.fill_box(h2o, n_compounds=1000, density=1000) assert [3.1042931 < period < 3.1042932 for period in filled.periodicity] + def test_fill_box_aspect_ratio(self, h2o): + filled = mb.fill_box(h2o, n_compounds=1000, + density=1000, aspect_ratio=[1, 2, 1]) + assert filled.periodicity[0]/filled.periodicity[1] == 0.5 + assert filled.periodicity[1]/filled.periodicity[2] == 2 + def test_fill_box_density_n_compounds(self, h2o): filled = mb.fill_box(h2o, density=1000, box=mb.Box([3.1042931, 3.1042931, 3.1042931])) assert filled.n_particles == 3000 + def test_fill_box_compound_ratio(self, h2o, ethane): + filled = mb.fill_box(compound=[h2o, ethane], density=800, + compound_ratio=[2, 1], box=[2, 2, 2, 4, 4, 4]) + n_ethane = len([c for c in filled.children if c.name == 'Ethane']) + n_water = len([c for c in filled.children if c.name == 'H2O']) + assert n_water / n_ethane == 2 + def test_fill_region(self, h2o): filled = mb.fill_region(h2o, n_compounds=50, region=[3, 2, 2, 4, 4, 3]) @@ -101,3 +116,36 @@ def test_wrong_box(self, h2o): filled = mb.fill_box(h2o, n_compounds=50, box=[2, 2]) with pytest.raises(MBuildError): filled = mb.fill_box(h2o, n_compounds=50, box=[2, 2, 2, 2]) + + def test_bad_args(self, h2o): + with pytest.raises(ValueError): + mb.fill_box(h2o, n_compounds=10) + with pytest.raises(ValueError): + mb.fill_box(h2o, density=1000) + with pytest.raises(ValueError): + mb.fill_box(h2o, box=[2, 2, 2]) + with pytest.raises(ValueError): + mb.fill_box(h2o, n_compounds=10, density=1000, box=[2, 2, 2]) + with pytest.raises(ValueError): + mb.fill_box(compound=[h2o, h2o], n_compounds=[10], density=1000) + with pytest.raises(ValueError): + mb.solvate(solute=h2o, solvent=[h2o], n_solvent=[10, 10], box=[2, 2, 2]) + with pytest.raises(ValueError): + mb.fill_region(h2o, n_compounds=[10, 10], region=[2, 2, 2, 4, 4, 4]) + + def test_write_temp_file(self, h2o): + cwd = os.getcwd() # Must keep track of the temp dir that pytest creates + filled = mb.fill_box(h2o, n_compounds=10, box=[4, 4, 4], temp_file='temp_file1.pdb') + region = mb.fill_region(h2o, 10, [2, 2, 2, 4, 4, 4], temp_file='temp_file2.pdb') + solvated = mb.solvate(filled, h2o, 10, box=[4, 4, 4], temp_file='temp_file3.pdb') + assert os.path.isfile(os.path.join(cwd, 'temp_file1.pdb')) + assert os.path.isfile(os.path.join(cwd, 'temp_file2.pdb')) + assert os.path.isfile(os.path.join(cwd, 'temp_file3.pdb')) + + def test_packmol_error(self, h2o): + with pytest.raises(RuntimeError): + filled = mb.fill_box(h2o, n_compounds=10, box=[0, 0, 0]) + + def test_packmol_warning(self, h2o): + with pytest.warns(UserWarning): + filled = mb.fill_box(h2o, n_compounds=10, box=[1, 1, 1], overlap=100)