From b7fdc1bb147c0e6892516cd0e7536811cf39f6e3 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Wed, 8 Jan 2025 12:59:26 -0800 Subject: [PATCH] handle vdW in coverage dependence 2D extraction functions --- pynta/coveragedependence.py | 36 ++++++++++++++++++++++++++++-------- pynta/main.py | 37 +++++++++++++++++++++++++++++++++++-- 2 files changed, 63 insertions(+), 10 deletions(-) diff --git a/pynta/coveragedependence.py b/pynta/coveragedependence.py index 6282e96e..abcbfbd3 100644 --- a/pynta/coveragedependence.py +++ b/pynta/coveragedependence.py @@ -1739,7 +1739,7 @@ def get_configs_of_concern(tree_interaction_regressor,configs,coad_stable_sites, return configs_of_concern def load_coverage_delta(d,ad_energy_dict,slab,metal,facet,sites,site_adjacency,ts_pynta_dir=None,allowed_structure_site_structures=None, - out_file_name="out",vib_file_name="vib",is_ad=None): + out_file_name="out",vib_file_name="vib",is_ad=None,keep_binding_vdW_bonds=False,keep_vdW_surface_bonds=False): try: info = json.load(open(os.path.join(d,"info.json"))) @@ -1755,7 +1755,9 @@ def load_coverage_delta(d,ad_energy_dict,slab,metal,facet,sites,site_adjacency,t if not is_ts: try: - admol,neighbor_sites,ninds = generate_adsorbate_2D(atoms, sites, site_adjacency, len(slab), max_dist=np.inf, cut_multidentate_off_num=None, allowed_structure_site_structures=allowed_structure_site_structures) + admol,neighbor_sites,ninds = generate_adsorbate_2D(atoms, sites, site_adjacency, len(slab), max_dist=np.inf, cut_multidentate_off_num=None, + allowed_structure_site_structures=allowed_structure_site_structures, + keep_binding_vdW_bonds=keep_binding_vdW_bonds,keep_vdW_surface_bonds=keep_vdW_surface_bonds) except (SiteOccupationException,TooManyElectronsException,ValueError) as e: return None,None,None,None @@ -1792,7 +1794,9 @@ def load_coverage_delta(d,ad_energy_dict,slab,metal,facet,sites,site_adjacency,t ts_info_path = os.path.join(ts_path,"info.json") try: - admol,neighbor_sites,ninds = generate_TS_2D(atoms, ts_info_path, metal, facet, sites, site_adjacency, len(slab), imag_freq_path=os.path.join(d,"vib.0.traj"), max_dist=np.inf, allowed_structure_site_structures=allowed_structure_site_structures) + admol,neighbor_sites,ninds = generate_TS_2D(atoms, ts_info_path, metal, facet, sites, site_adjacency, len(slab), imag_freq_path=os.path.join(d,"vib.0.traj"), + max_dist=np.inf, allowed_structure_site_structures=allowed_structure_site_structures, + keep_binding_vdW_bonds=keep_binding_vdW_bonds,keep_vdW_surface_bonds=keep_vdW_surface_bonds) except (SiteOccupationException,TooManyElectronsException, ValueError) as e: return None,None,None,None try: @@ -1882,11 +1886,24 @@ def extract_sample(d,ad_energy_dict,slab,metal,facet,sites,site_adjacency,pynta_ # vibdata = get_vibdata(os.path.join(d,out_file_name+".xyz"),os.path.join(d,vib_file_name+".json"),len(slab)) # except: # vibdata = None - + keep_binding_vdW_bonds=False + keep_vdW_surface_bonds=False + try: with open(os.path.join(d,"info.json"),'r') as f: info = json.load(f) - except: + mol = Molecule().from_adjacency_list(info["adjlist"]) + for bd in mol.get_all_edges(): + if bd.order == 0: + if bd.atom1.is_surface_site() or bd.atom2.is_surface_site(): + keep_binding_vdW_bonds = True + m = mol.copy(deep=True) + b = m.get_bond(m.atoms[mol.atoms.index(bd.atom1)],m.atoms[mol.atoms.index(bd.atom2)]) + m.remove_bond(b) + out = m.split() + if len(out) == 1: #vdW bond is not only thing connecting adsorbate to surface + keep_vdW_surface_bonds = True + except FileNotFoundError: info = None if is_ad is None: @@ -1930,13 +1947,15 @@ def extract_sample(d,ad_energy_dict,slab,metal,facet,sites,site_adjacency,pynta_ admol,neighbor_sites,ninds,dE = load_coverage_delta(d,ad_energy_dict,slab,metal,facet,sites,site_adjacency, ts_pynta_dir=pynta_dir,allowed_structure_site_structures=allowed_structure_site_structures, - out_file_name=out_file_name,vib_file_name=vib_file_name,is_ad=is_ad) + out_file_name=out_file_name,vib_file_name=vib_file_name,is_ad=is_ad, + keep_binding_vdW_bonds=keep_binding_vdW_bonds,keep_vdW_surface_bonds=keep_vdW_surface_bonds) if is_ad: try: admol_init,neighbor_sites_init,ninds_init = generate_adsorbate_2D(atoms_init, sites, site_adjacency, nslab, - max_dist=None, cut_multidentate_off_num=cut_multidentate_off_num, allowed_structure_site_structures=allowed_structure_site_structures) + max_dist=None, cut_multidentate_off_num=cut_multidentate_off_num, allowed_structure_site_structures=allowed_structure_site_structures, + keep_binding_vdW_bonds=keep_binding_vdW_bonds,keep_vdW_surface_bonds=keep_vdW_surface_bonds) except TooManyElectronsException: admol_init = None neighbor_sites_init = None @@ -1955,7 +1974,8 @@ def extract_sample(d,ad_energy_dict,slab,metal,facet,sites,site_adjacency,pynta_ try: admol_init,neighbor_sites_init,ninds_init = generate_TS_2D(atoms_init, info_path, metal, facet, sites, site_adjacency, nslab, imag_freq_path=imag_freq_path, max_dist=None, cut_multidentate_off_num=cut_multidentate_off_num, - allowed_structure_site_structures=allowed_structure_site_structures) + allowed_structure_site_structures=allowed_structure_site_structures, + keep_binding_vdW_bonds=keep_binding_vdW_bonds,keep_vdW_surface_bonds=keep_vdW_surface_bonds) except TooManyElectronsException as e: valid = False admol_init = None diff --git a/pynta/main.py b/pynta/main.py index e3be1a3f..a5add503 100644 --- a/pynta/main.py +++ b/pynta/main.py @@ -699,9 +699,26 @@ def setup_active_learning_loop(self): for ts in self.transition_states.keys(): info_path = os.path.join(self.pynta_run_directory,ts,"info.json") + with open(info_path) as f: + info = json.load(f) + mol = Molecule().from_adjacency_list(info["adjlist"]) + keep_binding_vdW_bonds=False + keep_vdW_surface_bonds=False + for bd in mol.get_all_edges(): + if bd.order == 0: + if bd.atom1.is_surface_site() or bd.atom2.is_surface_site(): + keep_binding_vdW_bonds = True + m = mol.copy(deep=True) + b = m.get_bond(m.atoms(mol.atoms.index(bd.atom1)),m.atoms[mol.atoms.index(bd.atom2)]) + m.remove_bond(b) + out = m.split() + if len(out) == 1: #vdW bond is not only thing connecting adsorbate to surface + keep_vdW_surface_bonds = True + atoms = read(admol_name_path_dict[ts]) st,_,_ = generate_TS_2D(atoms, info_path, self.metal, self.surface_type, self.sites, self.site_adjacency, self.nslab, - max_dist=np.inf, allowed_structure_site_structures=allowed_structure_site_structures) + max_dist=np.inf, allowed_structure_site_structures=allowed_structure_site_structures, + keep_binding_vdW_bonds=keep_binding_vdW_bonds,keep_vdW_surface_bonds=keep_vdW_surface_bonds) admol_name_structure_dict[ts] = st with open(info_path,"r") as f: info = json.load(f) @@ -713,12 +730,28 @@ def setup_active_learning_loop(self): p = os.path.join(self.pynta_run_directory,"Adsorbates",ad) with open(os.path.join(p,"info.json")) as f: info = json.load(f) + mol = Molecule().from_adjacency_list(info["adjlist"]) + if mol.contains_surface_site(): + keep_binding_vdW_bonds=False + keep_vdW_surface_bonds=False + for bd in mol.get_all_edges(): + if bd.order == 0: + if bd.atom1.is_surface_site() or bd.atom2.is_surface_site(): + keep_binding_vdW_bonds = True + m = mol.copy(deep=True) + b = m.get_bond(m.atoms(mol.atoms.index(bd.atom1)),m.atoms[mol.atoms.index(bd.atom2)]) + m.remove_bond(b) + out = m.split() + if len(out) == 1: #vdW bond is not only thing connecting adsorbate to surface + keep_vdW_surface_bonds = True + ad_xyz = get_best_adsorbate_xyz(p,self.sites,self.nslab) admol_name_path_dict[ad] = ad_xyz atoms = read(ad_xyz) - st,_,_ = generate_adsorbate_2D(atoms, self.sites, self.site_adjacency, self.nslab, max_dist=np.inf, allowed_structure_site_structures=allowed_structure_site_structures) + st,_,_ = generate_adsorbate_2D(atoms, self.sites, self.site_adjacency, self.nslab, max_dist=np.inf, allowed_structure_site_structures=allowed_structure_site_structures, + keep_binding_vdW_bonds=keep_binding_vdW_bonds,keep_vdW_surface_bonds=keep_vdW_surface_bonds) admol_name_structure_dict[ad] = st calculation_directories = [] #identify pairs directories