diff --git a/pynta/postprocessing.py b/pynta/postprocessing.py index 23ab75b7..ac98e6f0 100644 --- a/pynta/postprocessing.py +++ b/pynta/postprocessing.py @@ -12,6 +12,9 @@ from molecule.molecule import Molecule from acat.adsorption_sites import SlabAdsorptionSites from molecule.thermo import Wilhoit +from typing import List, Union, Tuple, Dict +from ase import db +from ase.db import connect eV_to_Jmol = 9.648328e4 @@ -387,4 +390,78 @@ def get_enthalpy_reaction(rnasas,pnasas,T,dT=0.01): dG -= th.get_enthalpy(T) for th in pnasas: dG += th.get_enthalpy(T) - return dG \ No newline at end of file + return dG + +def write_min_en_species_db( + ad_path: str, + slab_path: str, +) -> db: + ''' + Generates an ase database of the minimum energy configurations for all + Pynta optimized adsorbates and gas-phase species. + + Parameters: + ad_path: Absolute path to Pynta generated Adsorbates directory. + slab_path: Absolute path to Pynta optimized slab. ('your_path/slab.xyz') + Returns: + db: ase database + ''' + os.chdir(ad_path) + ad_dirs = [dir for dir in filter(os.path.isdir, os.listdir(ad_path))] + for ad in ad_dirs: + path_to_ad = os.path.join(ad_path, ad) + slab = read(slab_path) + + try: + with open(os.path.join(ad_path, "info.json")) as f: + info = json.load(f) + except: + print(f"No information for {ad_path}") + continue + + adj_list = info['adjlist'] + gasphase = len(info["gratom_to_molecule_surface_atom_map"]) == 0 + spin = (Molecule().from_adjacency_list(info["adjlist"]).multiplicity - 1.0) / 2.0 + name = info["name"] + ad_configs = [dir for dir in filter(os.path.isdir, os.listdir(path_to_ad))] + + # Setting up empty dictionaries to collect information for each config + ad_dict = {} + DFT_en_dict = {} + freq_dict = {} + atoms_dict = {} + for config in ad_configs: + optdir = os.path.join(path_to_ad, config, config + ".xyz") + freqdir = os.path.join(path_to_ad, config, "vib.json_vib.json") + atoms = read(optdir) + atoms_dict[config] = atoms + energy = atoms.get_potential_energy() + if gasphase: + vibdata = get_vibdata(optdir, freqdir, 0) + else: + vibdata = get_vibdata(optdir, freqdir, len(slab)) + freqs = vibdata.get_frequencies().tolist() + + # Handling of imaginary frequencies. + for index, freq in enumerate(freqs): + if freq.imag != 0: + # print(f" for {ad_path}/{d} \n we have imaginary freq \n setting to 12cm^-1") + freqs[index] = 12 + else: + freqs[index] = freq.real + zpe = vibdata.get_zero_point_energy() + freq_dict[config] = freqs + DFT_en_dict[config] = energy + ad_dict[config] = energy + zpe + if not gasphase: + min_E_config = min(ad_dict, key=ad_dict.get) + atoms = atoms_dict[min_E_config] + data = {} + data['frequencies'] = freq_dict[min_E_config] + db.write(atoms, name=name, adj_list=adj_list, spin=spin, gasphase=gasphase, data=data) + + + + + +