From 2909b085aaed0a812c195fa64ff211071fad1d44 Mon Sep 17 00:00:00 2001 From: Co Quach <43968221+daico007@users.noreply.github.com> Date: Fri, 17 Sep 2021 15:26:27 -0500 Subject: [PATCH] Add basic method to set up Monolayer and Solvated Monolayer (#2) * add basic method to set up monolayer and solvated monolayer * start adding simple tests * Update docstring for most calsses --- environment.yml | 3 +- surface_coatings/__init__.py | 3 +- surface_coatings/chains/alkylsilane.py | 3 +- surface_coatings/chains/vbc_polymer.py | 23 +++++-- surface_coatings/monolayer.py | 48 +++++++++++-- surface_coatings/solvated_monolayer.py | 68 +++++++++++++++++++ surface_coatings/solvated_system.py | 0 surface_coatings/surfaces/__init__.py | 2 + .../surfaces/silica_interface_carve.py | 6 +- .../surfaces/silicon_interface.py | 21 +++++- surface_coatings/tests/test_surfaces.py | 12 ++++ .../test_systems.py} | 0 12 files changed, 170 insertions(+), 19 deletions(-) create mode 100644 surface_coatings/solvated_monolayer.py delete mode 100644 surface_coatings/solvated_system.py create mode 100644 surface_coatings/tests/test_surfaces.py rename surface_coatings/{dual_monolayer.py => tests/test_systems.py} (100%) diff --git a/environment.yml b/environment.yml index 339ac20..076aeca 100644 --- a/environment.yml +++ b/environment.yml @@ -11,4 +11,5 @@ dependencies: - py3dmol - python=3.8 - rdkit - - scipy \ No newline at end of file + - scipy + - pytest diff --git a/surface_coatings/__init__.py b/surface_coatings/__init__.py index 82f8586..8d87644 100644 --- a/surface_coatings/__init__.py +++ b/surface_coatings/__init__.py @@ -1 +1,2 @@ -from .monolayer import Monolayer +from .monolayer import Monolayer, DualMonolayer +from .solvated_monolayer import SolvatedMonolayer, SolvatedDualMonolayer diff --git a/surface_coatings/chains/alkylsilane.py b/surface_coatings/chains/alkylsilane.py index 4fd39e5..78f9369 100644 --- a/surface_coatings/chains/alkylsilane.py +++ b/surface_coatings/chains/alkylsilane.py @@ -1,3 +1,4 @@ +"""Routine to create alkylsilane chain.""" import mbuild as mb from mbuild.lib.recipes import Alkane from mbuild.lib.moieties import Silane @@ -35,4 +36,4 @@ def __init__(self, chain_length=17, terminal_group="methyl"): self.add(silane, 'silane') mb.force_overlap(self['silane'], self['silane']['up'], self['alkane']['down']) - self.add(silane['down'], 'down', containment=False) + self.add(silane['down'], 'down', con tainment=False) diff --git a/surface_coatings/chains/vbc_polymer.py b/surface_coatings/chains/vbc_polymer.py index cc1cb80..17510d3 100644 --- a/surface_coatings/chains/vbc_polymer.py +++ b/surface_coatings/chains/vbc_polymer.py @@ -1,3 +1,4 @@ +"""Routine to create VBC polymer.""" import mbuild as mb from mbuild.lib.recipes import Polymer import numpy as np @@ -9,22 +10,36 @@ class VBCPolymer(mb.Compound): def __init__(self, monomers=[Methacrylate(), SBMA(), AzPMA(), TriazoleBiotin()], n=1, sequence='AABCBBD', port_labels=('up', 'down')): + """This is a general purpose to create a Vinylbenzyl Polymer + + Parameters + ---------- + monomers: list of monomer + List of monomers to be connected together + n: int, optional, default=1 + Number of repeat for the polymer + sequence: str + The sequence of all the monomer, corresponding to the monomers list + port_labels: tuple, optional, default=('up', 'down') + The list ports of the monomers. The VBC will connect to the first port, + and the last port will be capped by a Hydrogen + """ super(VBCPolymer, self).__init__() copolymer = Polymer(monomers=monomers) copolymer.build(n=n, sequence=sequence, add_hydrogens=False) - self.add(copolymer) + self.add(copolymer, label="Polymer") vbc = MethylStyrene() - self.add(vbc) + self.add(vbc, label="VBC") mb.force_overlap(vbc, vbc['up'], - copolymer['up']) + copolymer[port_labels[0]]) for port in self.all_ports(): if port.access_labels: self.labels['up'] = port else: self.labels['down'] = port - orientation = self['Polymer[0]']['down'].direction + orientation = self['Polymer']['down'].direction for port in self.available_ports(): if not all(np.isclose(port.direction, orientation)): port.update_orientation(-orientation) diff --git a/surface_coatings/monolayer.py b/surface_coatings/monolayer.py index fd211e4..db79848 100644 --- a/surface_coatings/monolayer.py +++ b/surface_coatings/monolayer.py @@ -1,6 +1,4 @@ -from copy import deepcopy -from warnings import warn - +"""Routines to create (dual) monolayer systems.""" import numpy as np import mbuild as mb @@ -9,11 +7,24 @@ class Monolayer(mb.Compound): - """A general monolayer recipe. + """A surface coated by a monolayer. Parameters ---------- - + surface: mb.Compound + The surface with ports at its surface. + chain: mb.Compound + The chain that to be attached to the surfaces. + n_chains: int + The number of chains to be attached. + backfill: mb.Compound, optional, default=H() + Compound used to backfill leftover ports (after all chains have been attached. + tile_x, tile_y: int, optional, default= 1, 1 + The number of surface tiles. + rotate: bool, optional, default=True + Options to rotate the chain randomly. + seed: int, optional, default= 12345 + Random seed used for any subprocess. """ def __init__(self, surface, chain, n_chains, backfill=H(), tile_x=1, tile_y=1, rotate=True, seed=12345, **kwargs): super(Monolayer, self).__init__() @@ -37,4 +48,29 @@ def __init__(self, surface, chain, n_chains, backfill=H(), tile_x=1, tile_y=1, r rotation = np.random.random() * np.pi * 2.0 chain.spin(rotation, [0, 0, 1]) - self.periodicity = [True, True, False] \ No newline at end of file + self.periodicity = surface.periodicity + + +class DualMonolayer(mb.Compound): + """A dual-monolayer system. + + Parameters + ---------- + top: mb.Compound + The top surfaces of the dual-monolayer. + bottom: mb.Compound + The bottom surface of the dual-monolayer. + separation: float, optional, default=0.8 + The separation between the two surfaces. + """ + def __init__(self, top, bottom, separation=0.8): + super(DualMonolayer, self).__init__() + top.spin(np.pi, around=[0, 1, 0]) + top_box = top.get_boundingbox() + bot_box = bottom.get_boundingbox() + + z_val = bot_box.lengths[2] + + top.translate([0, 0, z_val + separation]) + self.add(top, label="top_monolayer") + self.add(bottom, label="bottom_monolayer") \ No newline at end of file diff --git a/surface_coatings/solvated_monolayer.py b/surface_coatings/solvated_monolayer.py new file mode 100644 index 0000000..5aa5670 --- /dev/null +++ b/surface_coatings/solvated_monolayer.py @@ -0,0 +1,68 @@ +"""Routines to construct solvated (dual) monolayer.""" +import numpy as np +import mbuild as mb +from mbuild.lib.moieties import H2O + + +class SolvatedMonolayer(mb.Compound): + """Solvated monolayer system. + + Parameters + ---------- + monolayer: mb.Compound + The monolayer to be solvated. + solvent: mb.Compound, optional, default=H2O() + The solvent to be used. + n_solvents: int, optional, default=1000 + The number of solvent compounds used. + solvent_box_height: float, optional, default=5 + The height of the solvent box. The base of the box adapt those of the monolayer/surface. + seed: int, optional, default=12345 + Random seed used for any subprocess. + """ + def __init__(self, monolayer, solvent=H2O(), n_solvents=1000, solvent_box_height=5, seed=12345): + super(SolvatedMonolayer, self).__init__() + monolayer_box_lengths = monolayer.get_boundingbox().lengths + solvent_box = [monolayer_box_lengths[0], + monolayer_box_lengths[1], + solvent_box_height] + box_of_solvent = mb.fill_box(compound=solvent, + box=solvent_box, + n_compounds=n_solvents) + box_of_solvent.translate([0, 0, monolayer_box_lengths[2]]) + self.add(monolayer, label="monolayer") + self.add(box_of_solvent, label="solvent") + + +class SolvatedDualMonolayer(mb.Compound): + """Solvated dual-monolayer system. + + Parameters + ---------- + dual_monolayer: mb.Compound + The dual-monolayer system to be solvated. + solvent: mb.Compound, optional, default=H2O() + The solvent compound to be used. + n_solvents: int, optional, n=1000 + The number of solvent molecules to be used. + seed: int, optional, default=12345 + Random seed used for any subprocess. + """ + def __init__(self, dual_monolayer, solvent=H2O(), n_solvents=1000, seed=12345): + super(SolvatedDualMonolayer, self).__init__() + top_monolayer = dual_monolayer["top_monolayer"] + top_monolayer_box_lengths = top_monolayer.get_boundingbox().lengths + bottom_monolayer = dual_monolayer["bottom_monolayer"] + bottom_monolayer_box_lengths = bottom_monolayer.get_boundingbox().lengths + dual_monolayer_box_lengths = dual_monolayer.get_boundingbox().lengths + + separation = dual_monolayer_box_lengths[2] - (top_monolayer_box_lengths[2] + bottom_monolayer_box_lengths[2]) + solvent_box = [bottom_monolayer_box_lengths[0], + bottom_monolayer_box_lengths[1], + separation] + box_of_solvent = mb.fill_box(compound=solvent, + box=solvent_box, + n_compounds=n_solvents) + box_of_solvent.translate([0, 0, bottom_monolayer_box_lengths[2]]) + self.add(dual_monolayer, label="monolayers") + self.add(box_of_solvent, label="solvents") diff --git a/surface_coatings/solvated_system.py b/surface_coatings/solvated_system.py deleted file mode 100644 index e69de29..0000000 diff --git a/surface_coatings/surfaces/__init__.py b/surface_coatings/surfaces/__init__.py index e69de29..7bce356 100644 --- a/surface_coatings/surfaces/__init__.py +++ b/surface_coatings/surfaces/__init__.py @@ -0,0 +1,2 @@ +from .silicon_interface import SiliconInterface +from .silica_interface_carve import SilicaInterfaceCarve diff --git a/surface_coatings/surfaces/silica_interface_carve.py b/surface_coatings/surfaces/silica_interface_carve.py index 584206c..af7f4e3 100644 --- a/surface_coatings/surfaces/silica_interface_carve.py +++ b/surface_coatings/surfaces/silica_interface_carve.py @@ -1,3 +1,4 @@ +"""Routine to create carve amorphous silica interface.""" from __future__ import division import math @@ -167,8 +168,3 @@ def _adjust_stoichiometry(self): O1 = random.choice(bottom_Os) bottom_Os.remove(O1) self.remove(O1) - -if __name__ == "__main__": - from mbuild.lib.bulk_materials import AmorphousSilica - silica_interface = mb.SilicaInterfaceCarve(bulk_silica=AmorphousSilica(), thickness=1.2) - silica_interface.save('silica_interface.mol2', show_ports=True) diff --git a/surface_coatings/surfaces/silicon_interface.py b/surface_coatings/surfaces/silicon_interface.py index 67c0a92..03db795 100644 --- a/surface_coatings/surfaces/silicon_interface.py +++ b/surface_coatings/surfaces/silicon_interface.py @@ -1,9 +1,16 @@ +"""Routine to create silica surface.""" import numpy as np import mbuild as mb -from surface_coatings.monomers.methylstyrene import MethylStyrene class CrystalineSilicon(mb.Compound): + """A block of crystaline silicon. + + Parameters + ---------- + x, y, z: float, default= 10, 10, 10 + Dimension of the silicon block. + """ def __init__(self, x=10, y=10, z=10): super(CrystalineSilicon, self).__init__() # define all necessary lattice parameters @@ -30,6 +37,17 @@ def __init__(self, x=10, y=10, z=10): class SiliconInterface(mb.Compound): + """A surface made up of crystaline silicon. + + Parameters + ---------- + silicon: CrystalineSilicon + The crystaline silicon building block. + tile_x, tile_y: int + The number of tiles to build out the surface. + seed: int, optional, default=12345 + Random seed used in some subprocess. + """ def __init__(self, silicon=CrystalineSilicon(x=10, y=10, z=2), tile_x=1, tile_y=1, seed=12345): super(SiliconInterface, self).__init__() tiled_compound = mb.lib.recipes.TiledCompound(mb.clone(silicon), n_tiles=(tile_x, tile_y, 1)) @@ -39,6 +57,7 @@ def __init__(self, silicon=CrystalineSilicon(x=10, y=10, z=2), tile_x=1, tile_y= self.periodicity = silicon.periodicity def _identify_surface_sites(self): + """Method to identify and add port to silicon at the surface.""" for particle in list(self.particles()): if np.isclose(particle.pos[2], 0): label = f'Si_{len(self.referenced_ports())}' diff --git a/surface_coatings/tests/test_surfaces.py b/surface_coatings/tests/test_surfaces.py new file mode 100644 index 0000000..b6631db --- /dev/null +++ b/surface_coatings/tests/test_surfaces.py @@ -0,0 +1,12 @@ +import mbuild as mb +from surface_coatings.surfaces import SiliconInterface, SilicaInterfaceCarve + + +class TestSystem(object): + def test_silicon_interface(self): + silicon_surface = SiliconInterface() + assert silicon_surface.periodicity == (True, True, False) + + def test_silica_interface_carve(self): + silica_surface = SilicaInterfaceCarve) + assert silica_surface.periodicity == (True, True, False) \ No newline at end of file diff --git a/surface_coatings/dual_monolayer.py b/surface_coatings/tests/test_systems.py similarity index 100% rename from surface_coatings/dual_monolayer.py rename to surface_coatings/tests/test_systems.py