From 24e0a3e3d47034cc815cbd3fa3da523f5b23cf62 Mon Sep 17 00:00:00 2001 From: Muhammed Shuaibi Date: Wed, 12 Jun 2024 20:35:52 +0000 Subject: [PATCH] solvent-interface code --- ocdata/core/__init__.py | 7 +- ocdata/core/complex_solvent_config.py | 229 ++++++++++++++++++++++ ocdata/core/ion.py | 65 +++++++ ocdata/core/solvent.py | 77 ++++++++ ocdata/databases/pkls/__init__.py | 2 + ocdata/databases/pkls/ions.pkl | Bin 0 -> 6048 bytes ocdata/databases/pkls/solvents.pkl | Bin 0 -> 11210 bytes ocdata/utils/geometry.py | 265 ++++++++++++++++++++++++++ 8 files changed, 643 insertions(+), 2 deletions(-) create mode 100644 ocdata/core/complex_solvent_config.py create mode 100644 ocdata/core/ion.py create mode 100644 ocdata/core/solvent.py create mode 100644 ocdata/databases/pkls/ions.pkl create mode 100644 ocdata/databases/pkls/solvents.pkl create mode 100644 ocdata/utils/geometry.py diff --git a/ocdata/core/__init__.py b/ocdata/core/__init__.py index 16ed897..6cd54d0 100644 --- a/ocdata/core/__init__.py +++ b/ocdata/core/__init__.py @@ -1,5 +1,8 @@ -from .bulk import Bulk -from .slab import Slab from .adsorbate import Adsorbate from .adsorbate_slab_config import AdsorbateSlabConfig +from .bulk import Bulk +from .complex_solvent_config import ComplexSolventConfig +from .ion import Ion from .multi_adsorbate_slab_config import MultipleAdsorbateSlabConfig +from .slab import Slab +from .solvent import Solvent diff --git a/ocdata/core/complex_solvent_config.py b/ocdata/core/complex_solvent_config.py new file mode 100644 index 0000000..0dfebbe --- /dev/null +++ b/ocdata/core/complex_solvent_config.py @@ -0,0 +1,229 @@ +import os +import subprocess +import tempfile +from shutil import which +from typing import List + +import ase.io +import numpy as np + +from ocdata.core import Adsorbate, Ion, MultipleAdsorbateSlabConfig, Slab, Solvent +from ocdata.utils.geometry import BoxGeometry, PlaneBoundTriclinicGeometry + +# Code adapted from https://github.com/henriasv/molecular-builder/tree/master + + +class ComplexSolventConfig(MultipleAdsorbateSlabConfig): + """ + Class to represent a solvent, adsorbate, slab, ion config. This class only + returns a fixed combination of adsorbates placed on the surface. Solvent + placement is performed by packmol + (https://m3g.github.io/packmol/userguide.shtml), with the number of solvent + molecules controlled by its corresponding density. Ion placement is random + within the desired volume. + + Arguments + --------- + slab: Slab + Slab object. + adsorbates: List[Adsorbate] + List of adsorbate objects to place on the slab. + solvent: Solvent + Solvent object + ions: List[Ion] = [] + List of ion objects to place + num_sites: int + Number of sites to sample. + num_configurations: int + Number of configurations to generate per slab+adsorbate(s) combination. + This corresponds to selecting different site combinations to place + the adsorbates on. + interstitial_gap: float + Minimum distance, in Angstroms, between adsorbate and slab atoms as + well as the inter-adsorbate distance. + vacuum_size: int + Size of vacuum layer to add to both ends of the resulting atoms object. + solvent_interstitial_gap: float + Minimum distance, in Angstroms, between the solvent environment and the + adsorbate-slab environment. + solvent_depth: float + Volume depth to be used to pack solvents inside. + pbc_shift: float + Cushion to add to the packmol volume to avoid overlapping atoms over pbc. + packmol_tolerance: float + Packmol minimum distance to impose between molecules. + mode: str + "random", "heuristic", or "random_site_heuristic_placement". + This affects surface site sampling and adsorbate placement on each site. + + In "random", we do a Delaunay triangulation of the surface atoms, then + sample sites uniformly at random within each triangle. When placing the + adsorbate, we randomly rotate it along xyz, and place it such that the + center of mass is at the site. + + In "heuristic", we use Pymatgen's AdsorbateSiteFinder to find the most + energetically favorable sites, i.e., ontop, bridge, or hollow sites. + When placing the adsorbate, we randomly rotate it along z with only + slight rotation along x and y, and place it such that the binding atom + is at the site. + + In "random_site_heuristic_placement", we do a Delaunay triangulation of + the surface atoms, then sample sites uniformly at random within each + triangle. When placing the adsorbate, we randomly rotate it along z with + only slight rotation along x and y, and place it such that the binding + atom is at the site. + + In all cases, the adsorbate is placed at the closest position of no + overlap with the slab plus `interstitial_gap` along the surface normal. + """ + + def __init__( + self, + slab: Slab, + adsorbates: List[Adsorbate], + solvent: Solvent, + ions: List[Ion] = None, + num_sites: int = 100, + num_configurations: int = 1, + interstitial_gap: float = 0.1, + vacuum_size: int = 15, + solvent_interstitial_gap: float = 2, + solvent_depth: float = 8, + pbc_shift: float = 0.0, + packmol_tolerance: float = 2, + mode: str = "random_site_heuristic_placement", + ): + super().__init__( + slab=slab, + adsorbates=adsorbates, + num_sites=num_sites, + num_configurations=num_configurations, + interstitial_gap=interstitial_gap, + mode=mode, + ) + + self.solvent = solvent + self.ions = ions + self.vacuum_size = vacuum_size + self.solvent_depth = solvent_depth + self.solvent_interstitial_gap = solvent_interstitial_gap + self.pbc_shift = pbc_shift + self.packmol_tolerance = packmol_tolerance + + self.n_mol_per_volume = solvent.get_molecules_per_volume + + self.atoms_list, self.metadata_list = self.create_interface_on_sites( + self.atoms_list, self.metadata_list + ) + + def create_interface_on_sites(self, atoms_list, metadata_list): + atoms_interface_list = [] + metadata_interface_list = [] + + for atoms, adsorbate_metadata in zip(atoms_list, metadata_list): + cell = atoms.cell.copy() + unit_normal = cell[2] / np.linalg.norm(cell[2]) + cell[2] = self.solvent_depth * unit_normal + volume = cell.volume + n_solvent_mols = int(self.n_mol_per_volume * volume) + + if cell.orthorhombic: + box_length = cell.lengths() + geometry = BoxGeometry( + center=box_length / 2, length=box_length - self.pbc_shift + ) + else: + geometry = PlaneBoundTriclinicGeometry(cell, pbc=self.pbc_shift) + + solvent_ions_atoms = self.create_packmol_atoms(geometry, n_solvent_mols) + solvent_ions_atoms.set_cell(cell) + + max_z = atoms.positions[:, 2].max() + self.solvent_interstitial_gap + translation_vec = cell[2] + translation_vec[2] = max_z + solvent_ions_atoms.translate(translation_vec) + + interface_atoms = atoms + solvent_ions_atoms + interface_atoms.center(vacuum=self.vacuum_size, axis=2) + interface_atoms.wrap() + + atoms_interface_list.append(interface_atoms) + + metadata = { + "adsorbates": adsorbate_metadata, + "solvent": self.solvent.name, + "ions": [x.name for x in self.ions], + "ion_concentrations": [ + x.get_ion_concentration(volume) for x in self.ions + ], + "solvent_depth": self.solvent_depth, + "solvent_volume": volume, + } + metadata_interface_list.append(metadata) + + return atoms_interface_list, metadata_interface_list + + def create_packmol_atoms(self, geometry, n_solvent_mols): + cell = geometry.cell + with tempfile.TemporaryDirectory() as tmp_dir: + output_path = os.path.join(tmp_dir, "out.pdb") + self.solvent.atoms.write(f"{tmp_dir}/solvent.pdb", format="proteindatabank") + + # When placing a single ion, packmol strangely always places it at + # the boundary of the cell. This hacky fix manually places + # the ion in a random location in the cell. Packmol then will fix + # these atoms and not optimize them during its optimization, only + # optimizing solvent molecules arround the ion. + for i, ion in enumerate(self.ions): + ion_atoms = ion.atoms.copy() + ion_atoms.set_cell(cell) + self.randomize_coords(ion_atoms) + ion_atoms.write(f"{tmp_dir}/ion_{i}.pdb", format="proteindatabank") + + # write packmol input + packmol_input = os.path.join(tmp_dir, "input.inp") + with open(packmol_input, "w") as f: + f.write(f"tolerance {self.packmol_tolerance}\n") + f.write("filetype pdb\n") + f.write(f"output {output_path}\n") + + # write solvent + f.write( + geometry.packmol_structure( + f"{tmp_dir}/solvent.pdb", n_solvent_mols, "inside" + ) + ) + + for i in range(len(self.ions)): + f.write(f"structure {tmp_dir}/ion_{i}.pdb\n") + f.write(" number 1\n") + f.write(" fixed 0 0 0 0 0 0\n") + f.write("end structure\n\n") + + self.run_packmol(packmol_input) + + solvent_ions_atoms = ase.io.read(output_path, format="proteindatabank") + solvent_ions_atoms.set_pbc(True) + solvent_ions_atoms.set_tags([3] * len(solvent_ions_atoms)) + + return solvent_ions_atoms + + def run_packmol(self, packmol_input): + packmol_cmd = which("packmol") + + try: + ps = subprocess.Popen( + f"{packmol_cmd} < {packmol_input}", + shell=True, + stdout=subprocess.PIPE, + stderr=subprocess.PIPE, + ) + out, err = ps.communicate() + except NotImplementedError: + raise OSError("packmol is not found.") + + def randomize_coords(self, atoms): + cell_weights = np.random.rand(3) + cell_weights /= np.sum(cell_weights) + xyz = np.dot(cell_weights, atoms.cell) + atoms.set_center_of_mass(xyz) diff --git a/ocdata/core/ion.py b/ocdata/core/ion.py new file mode 100644 index 0000000..a315517 --- /dev/null +++ b/ocdata/core/ion.py @@ -0,0 +1,65 @@ +import pickle +import warnings +from typing import Any, Dict, Tuple + +import ase +import ase.units as units +import numpy as np + +from ocdata.databases.pkls import IONS_PKL_PATH + + +class Ion: + """ + Initializes an ion object in one of 2 ways: + - Directly pass in an ase.Atoms object. + - Pass in index of ion to select from ion database. + + Arguments + --------- + ion_atoms: ase.Atoms + ion structure. + ion_id_from_db: int + Index of ion to select. + ion_db_path: str + Path to ion database. + """ + + def __init__( + self, + ion_atoms: ase.Atoms = None, + ion_id_from_db: int = None, + ion_db_path: str = IONS_PKL_PATH, + ): + self.ion_id_from_db = ion_id_from_db + self.ion_db_path = ion_db_path + + if ion_atoms is not None: + self.atoms = ion_atoms.copy() + self.name = str(self.atoms.symbols) + elif ion_id_from_db is not None: + ion_db = pickle.load(open(ion_db_path, "rb")) + self._load_ion(ion_db[ion_id_from_db]) + + def __len__(self): + return len(self.atoms) + + def __str__(self): + return self.name + + def _load_ion(self, ion: Dict) -> None: + """ + Saves the fields from an ion stored in a database. Fields added + after the first revision are conditionally added for backwards + compatibility with older database files. + """ + self.atoms = ion["atoms"] + self.name = ion["name"] + self.charge = ion["charge"] + + def get_ion_concentration(self, volume): + """ + Compute the ion concentration units of M, given a volume in units of + Angstrom^3. + """ + return 1e27 / (units._Nav * volume) diff --git a/ocdata/core/solvent.py b/ocdata/core/solvent.py new file mode 100644 index 0000000..cbcd179 --- /dev/null +++ b/ocdata/core/solvent.py @@ -0,0 +1,77 @@ +import pickle +import warnings +from typing import Any, Dict, Tuple + +import ase +import ase.units as units +import numpy as np + +from ocdata.databases.pkls import SOLVENTS_PKL_PATH + + +class Solvent: + """ + Initializes a solvent object in one of 2 ways: + - Directly pass in an ase.Atoms object. + - Pass in index of solvent to select from solvent database. + + Arguments + --------- + solvent_atoms: ase.Atoms + Solvent molecule + solvent_id_from_db: int + Index of solvent to select. + solvent_db_path: str + Path to solvent database. + solvent_density: float + Desired solvent density to use. If not specified, the default is used + from the solvent databases. + """ + + def __init__( + self, + solvent_atoms: ase.Atoms = None, + solvent_id_from_db: int = None, + solvent_db_path: str = SOLVENTS_PKL_PATH, + solvent_density: float = None, + ): + self.solvent_id_from_db = solvent_id_from_db + self.solvent_db_path = solvent_db_path + self.solvent_density = solvent_density + + if solvent_atoms is not None: + self.atoms = solvent_atoms.copy() + self.name = str(self.atoms.symbols) + elif solvent_id_from_db is not None: + solvent_db = pickle.load(open(solvent_db_path, "rb")) + self._load_solvent(solvent_db[solvent_id_from_db]) + + self.molar_mass = sum(self.atoms.get_masses()) + + def __len__(self): + return len(self.atoms) + + def __str__(self): + return self.name + + def _load_solvent(self, solvent: Dict) -> None: + """ + Saves the fields from an adsorbate stored in a database. Fields added + after the first revision are conditionally added for backwards + compatibility with older database files. + """ + self.atoms = solvent["atoms"] + self.name = solvent["name"] + # use the default density if one is not specified + self.density = ( + solvent["density"] if not self.solvent_density else self.solvent_density + ) + + @property + def get_molecules_per_volume(self): + """ + Convert the solvent density in g/cm3 to the number of molecules per + angstrom cubed of volume. + """ + # molecules/mol * grams/cm3 / (1e24 A^3/cm^3 * g/mol) + return units._Nav * self.density / (1e24 * self.molar_mass) diff --git a/ocdata/databases/pkls/__init__.py b/ocdata/databases/pkls/__init__.py index 9a32604..bf966b8 100644 --- a/ocdata/databases/pkls/__init__.py +++ b/ocdata/databases/pkls/__init__.py @@ -2,3 +2,5 @@ BULK_PKL_PATH = os.path.join(__path__[0], "bulks.pkl") ADSORBATES_PKL_PATH = os.path.join(__path__[0], "adsorbates.pkl") +SOLVENTS_PKL_PATH = os.path.join(__path__[0], "solvents.pkl") +IONS_PKL_PATH = os.path.join(__path__[0], "ions.pkl") diff --git a/ocdata/databases/pkls/ions.pkl b/ocdata/databases/pkls/ions.pkl new file mode 100644 index 0000000000000000000000000000000000000000..548ac3b3e320b3e622a383a1fa38872a2b2422db GIT binary patch literal 6048 zcmeHLU5Fc16wc3NldNkAy4FQdRz!&?i|Y^de~_>$op!5T^g-G#Gn0vTN-`7XM;qNr zMQRO32sZEZIVa%hHEG>}Tz#72PQIBJ6imsx(T-k8vB$C#-3RnLt3M`j&E zw=LImeBFbmf($dHLDniqCi7HD?5CxqMi0>K;BHb$E9vJc84WUJuTeMfp;YB|+I%6p zpAM3D)9LBy&40teew2hqwPxwTXvpXPO99T+0%^9c>6BbbHll-MtnD6LTp6Ua$xfH% zgf6AQ8AJ<)*jQZ_)+|IoGaOd{os6K(j?t2Eq8!GmpZZh>a4?NO|>O5{xCFv_k&?rd- z*{bOg)2~q?k!++~GE6RJ|42bVU(SaW)(g2-@Gh9r8=ECA12Ty&-_mV|yvL}bAs#px zFW2^S{t7-KJYr zUn6uIMXtfIPB#(VR34fCo^B~$kz%^Jezl@{#%fNcnJAeS3yf@uL|nD?YJ?|@;bU*q?-@Q8zWei6TeU(*O4YX=(=t|b0_9($`B+80mvS^GU$S`{)B-BE8R8Cx&m&~ly1|6N)!K8&y{0m{-H@} z*(z4D2qVV&Iy5#$7-EcVLKT`Ogix<3sDll zpWUnuYBWps^*g@roO92;=Q}^oKeg?rH97d&tyep8*hkD3a5j7OSiufw#2aFxGZJhL z>#-7N*t=YOFSrwF`}CXk=x|GMqucBE2VL9sSTRh4D?MhbgF#GUd_o~-n;!deAlls0 zw%i>Id6zdw{Slv_moUFs;y%(nnsjX(h>~&B0gqzl?=3QRT_JJvWwih_0{*te;(TYT4Y~>{PD2kTOD_H zk-|Ou>mK!=Cgnf6{jYC5aD+bUKl|B>!@ZQt47%OX7H7cSrjp4_#w1p9Ws}4_B8MW^ zHm5rh3}w7uEP2X+_m_pdVV@`Jg@dkKIN!G)+%mauXk16R!{#W)BeZe@W4e*O*RyvO zduO`lO9%gfqnb|Yl+(hc+6<{>AQ4V4ypiKVxsPSlZhCREHgf-U`T~)#em?^MUtn91o#a?Squc3s>-AArtjOaHup`%|uOms$ zBwFUT1D_8he#sg|nGvg(!Bka;{RLEjy;u{@x z(oeNo(D|D)=M}Pf%-k!b+Z76#(%e=fNY=ii;3wDki-pTxetPJIUf~jxOA&NBN#WnN zzIu0AUd=H5*ZAhs&;O>6CN{lt-QIySMAAR4-tzdNvA5{HO4|$-ttc=*-$$g~r^4S8^?8y08@;EwDZl zKl3WBuaVzoH9v;R>AbrslQ&=L@2J~WYu(adKqX*69eYb0kR$T+sd>VQv^sfKbl?}& zT4plLzbH7+cHFOPZM${tLFV?_u63lp=I9&FEA3H79(VUb&D-h>c;~aTS8!+VlumxJ zqm{|$ijGG^0WK8+T-KaqlK=Va+P|eHJ|@2#dVz)U9s_3tosLc%Ywcg!L+|sf0tg@< zJlC?i^n(d+%q?K?o) z`j@g5C&+r1QZ=0^iGQD=_?uM;D(C`-OXC<84i8ZAI}!Zjg1@leQ*7e4ZKk#v7uZX84E{**BvL*6yHaURo0jro{y1b#^|boS)OCB6$f z5p=tx3qE$i$KWsI+J#($KN?X{%fKeG*x1DR&RuLo`HI*!UA|mh_D7Vjh)p0reQc_W z&-cKBYgYsYx%M=HVWlc4sA@AAAW)V9_(2u}{-)VF|7_=j{A>dE!1$4##qBjwhm8Yw zeni_gp9{JLd>^!8;bIb$t;X5f-#j1Bg95x4NbtUm-NC!8DN$H77+@wFH&EuRE>!Q-SEn#0}W z;&)j<^$}j6U7$yyPP_di<%Q^R(U2M^oZkrj11pOTg!)9`c7cfb1dXczb=YJE>e^g` zy5x^~D8PNX9@z(VGe^`z0rH?n29WQX?;vm1A8nxkeQAp{pnr>Oz)b~z3Gh}U%2NP8 z$Ya2-nSo)hu0fA%27pBE(*D+CpGO@9JOJ1{&EepDjsT^6LmJA3_gF=O*=`UX00q*Q z%Km;5y8DM6e;2UXD*8VsDRhS_UK%IR1?sx~&!R8?E|L8ppniC6lwf=4z#$UFxX@3r zoi#K@C7t(;j#J(b%d=zZf9pZfcN(Nm{Tf~=g3s&@rpWuO+gDrsImI^I^C`mnZr$`o ze7bj@z5#8xLt~;J*Db=$0OkFY>g;S0mIYf znL+b^zeafqxCMC`;I`8YSC`v0u1x`-a&2kgvwAZ#=~M`vEo6cW2AK`B^}tN?VNKbD z#^Jr!y|QUwZ`4u9p>a{f8wokxnm-r{@+Y0%`3#f=d>g92Bs|1wl&Lcw83h0qNC0f* zz|l)&agNU$w+vHgzqj^}iU$pYqH*6(Fm2qB6b~GR>7s|@|9E!oC*pdDu0dq%7x`k! zDE7@RjhZQd2v*2|sG85uT>K*+1;9+_W0pSi8u=)I1NdYBj)rUl1?r-3bg9eD=D$c3 z8-G@)hWyCS-swhIJ<3qzMj(SZdUYnpmA?a@Y&!6{J%u-%SaEgKQN$g1fyc-GpwQ-T zWS{JxQS|hs3p)K8D{i;HHc96_zij", point_ext - point_plane, vec)) + + @staticmethod + def vec_and_point_to_plane(vec, point): + """Returns the (unique) plane, given a normal vector 'vec' and a + point 'point' in the plane. ax + by + cz - d = 0 + + :param vec: normal vector of plane + :type vec: ndarray + :param point: point in plane + :type point: ndarray + :returns: parameterization of plane + :rtype: ndarray + """ + return np.array((*vec, np.dot(vec, point))) + + @staticmethod + def cell2planes(cell, pbc): + """Get the parameterization of the sizes of a ase.Atom cell + + :param cell: ase.Atom cell + :type cell: obj + :param pbc: shift of boundaries to be used with periodic boundary condition + :type pbc: float + :returns: parameterization of cell plane sides + :rtype: list of ndarray + + 3 planes intersect the origin by ase design. + """ + a = cell[0] + b = cell[1] + c = cell[2] + + n1 = np.cross(a, b) + n2 = np.cross(c, a) + n3 = np.cross(b, c) + + origin = np.array([0, 0, 0]) + pbc / 2 + top = (a + b + c) - pbc / 2 + + plane1 = Geometry.vec_and_point_to_plane(n1, origin) + plane2 = Geometry.vec_and_point_to_plane(n2, origin) + plane3 = Geometry.vec_and_point_to_plane(n3, origin) + plane4 = Geometry.vec_and_point_to_plane(-n1, top) + plane5 = Geometry.vec_and_point_to_plane(-n2, top) + plane6 = Geometry.vec_and_point_to_plane(-n3, top) + + return [plane1, plane2, plane3, plane4, plane5, plane6] + + @staticmethod + def extract_box_properties(center, length, lo_corner, hi_corner): + """Given two of the properties 'center', 'length', 'lo_corner', + 'hi_corner', return all the properties. The properties that + are not given are expected to be 'None'. + """ + # exactly two arguments have to be non-none + if sum(x is None for x in [center, length, lo_corner, hi_corner]) != 2: + raise ValueError("Exactly two arguments have to be given") + + # declare arrays to allow mathematical operations + center, length = np.asarray(center), np.asarray(length) + lo_corner, hi_corner = np.asarray(lo_corner), np.asarray(hi_corner) + relations = [ + [ + "lo_corner", + "hi_corner - length", + "center - length / 2", + "2 * center - hi_corner", + ], + [ + "hi_corner", + "lo_corner + length", + "center + length / 2", + "2 * center - lo_corner", + ], + [ + "length / 2", + "(hi_corner - lo_corner) / 2", + "hi_corner - center", + "center - lo_corner", + ], + [ + "center", + "(hi_corner + lo_corner) / 2", + "hi_corner - length / 2", + "lo_corner + length / 2", + ], + ] + + # compute all relations + relation_list = [] + for relation in relations: + for i in relation: + try: + relation_list.append(eval(i)) + except TypeError: + continue + + # keep the non-None relations + for i, relation in enumerate(relation_list): + if None in relation: + del relation_list[i] + return relation_list + + def packmol_structure(self, filename, number, side): + """Make structure to be used in PACKMOL input script + + :param number: number of solvent molecules + :type number: int + :param side: pack solvent inside/outside of geometry + :type side: str + :returns: string with information about the structure + :rtype: str + """ + structure = "" + structure += f"structure {filename}\n" + structure += f" number {number}\n" + structure += f" {side} {self.__repr__()} " + for param in self.params: + structure += f"{param} " + structure += "\nend structure\n" + return structure + + +class PlaneBoundTriclinicGeometry(Geometry): + """Triclinic crystal geometry based on ase.Atom cell + + :param cell: ase.Atom cell + :type cell: obj + :param pbc: shift of boundaries to be used with periodic boundary condition + :type pbc: float + """ + + def __init__(self, cell, pbc=0.0): + self.planes = self.cell2planes(cell, pbc) + self.cell = cell + self.ll_corner = [0, 0, 0] + a = cell[0, :] + b = cell[1, :] + c = cell[2, :] + self.ur_corner = a + b + c + + def packmol_structure(self, filename, number, side): + """Make structure to be used in PACKMOL input script""" + structure = "" + + if side == "inside": + side = "over" + elif side == "outside": + side = "below" + structure += f"structure {filename}\n" + structure += f" number {number}\n" + for plane in self.planes: + structure += f" {side} plane " + for param in plane: + structure += f"{param} " + structure += "\n" + structure += "end structure\n" + return structure + + def __call__(self, position): + raise NotImplementedError + + +class BoxGeometry(Geometry): + """Box geometry. + + :param center: geometric center of box + :type center: array_like + :param length: length of box in all directions + :type length: array_like + :param lo_corner: lower corner + :type lo_corner: array_like + :param hi_corner: higher corner + :type hi_corner: array_like + """ + + def __init__( + self, center=None, length=None, lo_corner=None, hi_corner=None, **kwargs + ): + super().__init__(**kwargs) + props = self.extract_box_properties(center, length, lo_corner, hi_corner) + self.ll_corner, self.ur_corner, self.length_half, self.center = props + self.params = list(self.ll_corner) + list(self.ur_corner) + self.length = self.length_half * 2 + + def __repr__(self): + return "box" + + def __call__(self, atoms): + positions = atoms.get_positions() + dist = self.distance_point_plane(np.eye(3), self.center, positions) + indices = np.all((np.abs(dist) <= self.length_half), axis=1) + return indices + + def volume(self): + return np.prod(self.length)