Add support for specifying adsorption positions in adsorbate calculations

Author: wladerer

vsh/scripts/__init__.py

@@ -148,6 +148,8 @@ def adsorb(subparsers):
     subp_adsorb.add_argument("-c","--coverage",type=int,nargs=3,default=[1, 1, 1],help="Adsorbate coverage")
     subp_adsorb.add_argument("-b","--both",action='store_true',help="Adsorbate on both surfaces")
     subp_adsorb.add_argument("-o", "--output", type=str, default=None, help="Output basename")
+    subp_adsorb.add_argument("-p", "--positions", choices=['all', 'ontop', 'bridge', 'hollow'], default='all', help="Adsorption sites")
+    subp_adsorb.add_argument("-P", "--plane", type=int, nargs=3, help="Miller plane required for adsorbing on both surfaces")
 
 def alchemy(subparsers):
 

vsh/scripts/adsorb.py

@@ -1,7 +1,28 @@
 from pymatgen.analysis.adsorption import AdsorbateSiteFinder
 from pymatgen.core.structure import Molecule, Structure
+from pymatgen.core.surface import Slab
 from pymatgen.io.vasp.inputs import Poscar
 
+
+def slab_from_structure(miller: list[int], structure: Structure) -> Slab:
+    '''
+    Returns a pymatgen.core.surface.Slab from a pymatgen structure
+    '''
+    import numpy as np
+
+    #if miller is None, raise error
+    if miller is None:
+        raise ValueError('miller index is required')
+
+    return Slab(lattice=structure.lattice,
+                species=structure.species_and_occu,
+                coords=structure.frac_coords,
+                miller_index=miller,
+                oriented_unit_cell=structure,
+                shift=0,
+                scale_factor=np.eye(3, dtype=int),
+                site_properties=structure.site_properties)
+
 def structure_from_file(filename: str) -> Structure:
     '''
     Creates a pymatgen structure from a file
@@ -24,23 +45,23 @@ def get_chemical_formula(structure: Structure) -> str:
     '''
     return structure.composition.reduced_formula
 
-def add_adsorbate_single(structure: Structure, adsorbate: Molecule, coverage: list[int] = [1, 1, 1], distance: float = 1.0) -> list[Structure]:
+def add_adsorbate_single(structure: Structure, adsorbate: Molecule, coverage: list[int] = [1, 1, 1], distance: float = 1.0, positions=('ontop', 'bridge', 'hollow')) -> list[Structure]:
     '''
     Finds all adsorption sites on a structure and adsorbs the adsorbate at each site. Returns a list of adsorbed structures.
     '''
 
     asf = AdsorbateSiteFinder(structure)
-    ads_structs = asf.generate_adsorption_structures(adsorbate, repeat=coverage, find_args={"distance": distance})  # edit later
+    ads_structs = asf.generate_adsorption_structures(adsorbate, repeat=coverage, find_args={"distance": distance, "positions": positions})  # edit later
 
     return ads_structs
 
-def add_adsorbate_on_both_surfaces(structure: Structure, adsorbate: Molecule, coverage: list[int] = [1, 1, 1], distance: float = 1.0) -> list[Structure]:
+def add_adsorbate_on_both_surfaces(structure: Structure, adsorbate: Molecule, miller_plane: list[int], coverage: list[int] = [1, 1, 1], distance: float = 1.0, positions=('ontop', 'bridge', 'hollow')) -> list[Structure]:
     '''
     Finds all adsorption sites on a structure and adsorbs the adsorbate at each site. Returns a list of adsorbed structures.
     '''
-
-    asf = AdsorbateSiteFinder(structure)
-    ads_structs = asf.adsorb_both_surfaces(adsorbate, repeat=coverage, find_args={"distance": distance})  # edit later
+    slab = slab_from_structure(miller=miller_plane, structure=structure)
+    asf = AdsorbateSiteFinder(slab)
+    ads_structs = asf.adsorb_both_surfaces(adsorbate, repeat=coverage, find_args={"distance": distance, "positions": positions})  # edit later
 
     return ads_structs
 
@@ -59,10 +80,14 @@ def create_adsorbed_structure(args):
     structure = structure_from_file(args.input)
     adsorbate = adsorbate_from_file(args.adsorbate)
 
+    #if positions = 'all', then set positions to be positions=('ontop', 'bridge', 'hollow')
+    if args.positions == 'all':
+        args.positions = ('ontop', 'bridge', 'hollow')
+
     if args.both: 
-        ads_structs = add_adsorbate_on_both_surfaces(structure, adsorbate, coverage=args.coverage, distance=args.distance)
+        ads_structs = add_adsorbate_on_both_surfaces(structure, adsorbate, miller_plane=args.plane, coverage=args.coverage, distance=args.distance, positions=args.positions)
     else:
-        ads_structs = add_adsorbate_single(structure, adsorbate, coverage=args.coverage, distance=args.distance)
+        ads_structs = add_adsorbate_single(structure, adsorbate, coverage=args.coverage, distance=args.distance, positions=args.positions)
 
     poscars = [ Poscar(ads_struct, sort_structure=True) for ads_struct in ads_structs ]