diff --git a/ocdata/core/__init__.py b/ocdata/core/__init__.py index 16ed897..b29e1cc 100644 --- a/ocdata/core/__init__.py +++ b/ocdata/core/__init__.py @@ -1,5 +1,8 @@ from .bulk import Bulk +from .ion import Ion from .slab import Slab +from .solvent import Solvent from .adsorbate import Adsorbate from .adsorbate_slab_config import AdsorbateSlabConfig from .multi_adsorbate_slab_config import MultipleAdsorbateSlabConfig +from .complex_solvent_config import ComplexSolventConfig diff --git a/ocdata/core/complex_solvent_config.py b/ocdata/core/complex_solvent_config.py new file mode 100644 index 0000000..5a288c6 --- /dev/null +++ b/ocdata/core/complex_solvent_config.py @@ -0,0 +1,230 @@ +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") + if not packmol_cmd: + raise OSError("packmol not found.") + + ps = subprocess.Popen( + f"{packmol_cmd} < {packmol_input}", + shell=True, + stdout=subprocess.PIPE, + stderr=subprocess.PIPE, + ) + out, err = ps.communicate() + if err: + raise OSError(err.decode("utf-8")) + + 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 0000000..548ac3b Binary files /dev/null and b/ocdata/databases/pkls/ions.pkl differ diff --git a/ocdata/databases/pkls/solvents.pkl b/ocdata/databases/pkls/solvents.pkl new file mode 100644 index 0000000..cbe2e40 Binary files /dev/null and b/ocdata/databases/pkls/solvents.pkl differ diff --git a/ocdata/solvents/flags.py b/ocdata/solvents/flags.py deleted file mode 100644 index 206c3a8..0000000 --- a/ocdata/solvents/flags.py +++ /dev/null @@ -1,35 +0,0 @@ -BASE_FLAGS = { - "prec": "Normal", - "gga": "RP", - "pp": "PBE", - "xc": "PBE", - "ivdw": 11, - "encut": 400.0, - "ediff": 1e-6, - "nelm": 100, - "ismear": 0, - "sigma": 0.1, - "lwave": False, - "lcharg": False, - "isif": 0, - "ispin": 2, - "algo": "Fast", - "idipol": 3, - "ldipol": True, - "lasph": True, - "lreal": "Auto", - "ncore": 4, - "dipol": [0.5, 0.5, 0.5], -} - -OPT_FLAGS = {"ibrion": 2, "nsw": 0} - -MD_FLAGS = { - "ibrion": 0, - "nsw": 500, - "smass": 0, - "tebeg": 300, - "potim": 1, -} - -ML_FLAGS = {"ML_LMLFF": True, "ML_MODE": "train", "ML_EPS_LOW": 1e-7} diff --git a/ocdata/utils/geometry.py b/ocdata/utils/geometry.py new file mode 100644 index 0000000..2d695c9 --- /dev/null +++ b/ocdata/utils/geometry.py @@ -0,0 +1,265 @@ +import warnings + +import numpy as np +from ase import Atom + + +class Geometry: + """Base class for geometries. + + :param periodic_boundary_condition: self-explanatory + :type periodic_boundary_condition: array_like + :param minimum_image_convention: use the minimum image convention for + bookkeeping how the particles interact + :type minimum_image_convention: bool + """ + + def __init__( + self, + periodic_boundary_condition=(False, False, False), + minimum_image_convention=True, + ): + self.minimum_image_convention = minimum_image_convention + self.periodic_boundary_condition = periodic_boundary_condition + pass + + def __call__(self, atoms): + """The empty geometry. False because we define no particle to be + in the dummy geometry. + + :param atoms: atoms object from ase.Atom that is being modified + :type atoms: ase.Atom obj + :returns: ndarray of bools telling which atoms to remove + :rtype: ndarray of bool + """ + return np.zeros(len(atoms), dtype=np.bool) + + @staticmethod + def distance_point_line(vec, point_line, point_ext): + """Returns the (shortest) distance between a line parallel to + a normal vector 'vec' through point 'point_line' and an external + point 'point_ext'. + + :param vec: unit vector parallel to line + :type vec: ndarray + :param point_line: point on line + :type point_line: ndarray + :param point_ext: external points + :type point_ext: ndarray + :return: distance between line and external point(s) + :rtype: ndarray + """ + return np.linalg.norm(np.cross(vec, point_ext - point_line), axis=1) + + @staticmethod + def distance_point_plane(vec, point_plane, point_ext): + """Returns the (shortest) distance between a plane with normal vector + 'vec' through point 'point_plane' and a point 'point_ext'. + + :param vec: normal vector of plane + :type vec: ndarray + :param point_plane: point on line + :type point_plane: ndarray + :param point_ext: external point(s) + :type point_ext: ndarray + :return: distance between plane and external point(s) + :rtype: ndarray + """ + vec = np.atleast_2d(vec) # Ensure n is 2d + return np.abs(np.einsum("ik,jk->ij", 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) diff --git a/ocdata/utils/vasp.py b/ocdata/utils/vasp.py index d49bdc9..c66212a 100644 --- a/ocdata/utils/vasp.py +++ b/ocdata/utils/vasp.py @@ -16,40 +16,7 @@ from ase.calculators.vasp import Vasp from ase.io.trajectory import TrajectoryWriter -# NOTE: this is the setting for slab and adslab -VASP_FLAGS = { - "ibrion": 2, - "nsw": 2000, - "isif": 0, - "isym": 0, - "lreal": "Auto", - "ediffg": -0.03, - "symprec": 1e-10, - "encut": 350.0, - "laechg": True, - "lwave": False, - "ncore": 4, - "gga": "RP", - "pp": "PBE", - "xc": "PBE", -} - -# This is the setting for bulk optmization. -# Only use when expanding the bulk_db with other crystal structures. -BULK_VASP_FLAGS = { - "ibrion": 1, - "nsw": 100, - "isif": 7, - "isym": 0, - "ediffg": 1e-08, - "encut": 500.0, - "kpts": (10, 10, 10), - "prec": "Accurate", - "gga": "RP", - "pp": "PBE", - "lwave": False, - "lcharg": False, -} +from ocdata.utils.vasp_flags import BULK_VASP_FLAGS, VASP_FLAGS def _clean_up_inputs(atoms, vasp_flags): diff --git a/ocdata/utils/vasp_flags.py b/ocdata/utils/vasp_flags.py new file mode 100644 index 0000000..d7f235f --- /dev/null +++ b/ocdata/utils/vasp_flags.py @@ -0,0 +1,70 @@ +# OC20 +# NOTE: this is the setting for slab and adslab +VASP_FLAGS = { + "ibrion": 2, + "nsw": 2000, + "isif": 0, + "isym": 0, + "lreal": "Auto", + "ediffg": -0.03, + "symprec": 1e-10, + "encut": 350.0, + "laechg": True, + "lwave": False, + "ncore": 4, + "gga": "RP", + "pp": "PBE", + "xc": "PBE", +} +# This is the setting for bulk optmization. +# Only use when expanding the bulk_db with other crystal structures. +BULK_VASP_FLAGS = { + "ibrion": 1, + "nsw": 100, + "isif": 7, + "isym": 0, + "ediffg": 1e-08, + "encut": 500.0, + "kpts": (10, 10, 10), + "prec": "Accurate", + "gga": "RP", + "pp": "PBE", + "lwave": False, + "lcharg": False, +} + +SOLVENT_FLAGS = { + "prec": "Normal", + "gga": "RP", + "pp": "PBE", + "xc": "PBE", + "ivdw": 11, + "encut": 400.0, + "ediff": 1e-6, + "nelm": 100, + "ismear": 0, + "sigma": 0.1, + "lwave": False, + "lcharg": False, + "isif": 0, + "ispin": 2, + "algo": "Fast", + "idipol": 3, + "ldipol": True, + "lasph": True, + "lreal": "Auto", + "ncore": 4, + "dipol": [0.5, 0.5, 0.5], +} + +OPT_FLAGS = {"ibrion": 2, "nsw": 0} + +MD_FLAGS = { + "ibrion": 0, + "nsw": 500, + "smass": 0, + "tebeg": 300, + "potim": 1, +} + +ML_FLAGS = {"ML_LMLFF": True, "ML_MODE": "train", "ML_EPS_LOW": 1e-7} diff --git a/tests/test_complex_solvent_config.py b/tests/test_complex_solvent_config.py new file mode 100644 index 0000000..2b28358 --- /dev/null +++ b/tests/test_complex_solvent_config.py @@ -0,0 +1,80 @@ +import random + +import numpy as np +import pytest +from ase import Atoms +from ase.data import covalent_radii +from pymatgen.analysis.adsorption import AdsorbateSiteFinder +from pymatgen.io.ase import AseAtomsAdaptor + +from ocdata.core import Adsorbate, Bulk, ComplexSolventConfig, Ion, Slab, Solvent +from ocdata.core.adsorbate_slab_config import get_interstitial_distances +from ocdata.databases.pkls import ADSORBATES_PKL_PATH, BULK_PKL_PATH + + +@pytest.fixture(scope="class") +def load_data(request): + request.cls.bulk = Bulk(bulk_id_from_db=0) + + adsorbate_idx = [0, 1] + adsorbates = [Adsorbate(adsorbate_id_from_db=i) for i in adsorbate_idx] + + solvent = Solvent(solvent_id_from_db=0) + ions = [Ion(ion_id_from_db=3)] + + request.cls.adsorbates = adsorbates + request.cls.solvent = solvent + request.cls.ions = ions + request.cls.vacuum = 15 + request.cls.solvent_depth = 10 + request.cls.solvent_interstitial_gap = 2 + + +@pytest.mark.usefixtures("load_data") +class TestComplex: + def test_num_configurations(self): + random.seed(1) + np.random.seed(1) + + slab = Slab.from_bulk_get_random_slab(self.bulk) + adslabs = ComplexSolventConfig( + slab, + self.adsorbates, + self.solvent, + self.ions, + vacuum_size=self.vacuum, + solvent_interstitial_gap=self.solvent_interstitial_gap, + solvent_depth=self.solvent_depth, + num_configurations=10, + ) + assert len(adslabs.atoms_list) == 10 + assert len(adslabs.metadata_list) == 10 + + def test_solvent_density(self): + """ + Test that the number of solvent + ion molecules inside the environment + is consistent with the specified density. + """ + random.seed(1) + np.random.seed(1) + + slab = Slab.from_bulk_get_random_slab(self.bulk) + adslabs = ComplexSolventConfig( + slab, + self.adsorbates, + self.solvent, + self.ions, + vacuum_size=self.vacuum, + solvent_interstitial_gap=self.solvent_interstitial_gap, + solvent_depth=self.solvent_depth, + num_configurations=10, + ) + + for atoms, metadata in zip(adslabs.atoms_list, adslabs.metadata_list): + volume = metadata["solvent_volume"] + n_solvent_mols = int(volume * self.solvent.get_molecules_per_volume) + n_solvent_atoms = n_solvent_mols * len(self.solvent.atoms) + n_ions = len(self.ions) + + solvent_ions_atoms = atoms[atoms.get_tags() == 3] + assert len(solvent_ions_atoms) == n_solvent_atoms + n_ions