#!/usr/bin/env amspython from scm.plams import * import numpy as np import os import matplotlib.pyplot as plt def main(): ### CHANGE THE BELOW results_dir TO THE CORRECT PATH ### results_dir = "myjob.results/" job = AMSJob.load_external(results_dir) data = get_workfunction(job, step_z=0.5, step_x=0.2, step_y=0.2, nproc=None, padding_z=150) create_workfunction_plot(data, fname=f"picture-{job.name}.png", title=job.name, show=True) save_potential_vs_z(data, f"plane_averaged_coulomb_potential_{job.name}.txt") def get_workfunction( job: AMSJob, padding_z: float = 100, step_x: float = 1.0, step_y: float = 1.0, step_z: float = 0.5, nproc: int = 1 ): """ This function calculates plane-averaged electrostatic potentials, and calculates the work function by subtracting Fermi energy from the maxima from each side of the slab. job: AMSJob A finished BAND 2D slab job .. note:: The z values do NOT match the input structure but rather the structures used internally by BAND (geometry%standard_xyzatm) """ molecule_coords = job.results.readrkf("Molecule", "Coords") molecule_coords = np.array(molecule_coords).reshape(-1, 3) coords = job.results.readrkf("geometry", "standard_xyzatm", file="engine") coords = np.array(coords).reshape(-1, 3) lattice = job.results.readrkf("geometry", "standard_lattice", file="engine") lattice = np.array(lattice).reshape(-1, 3) atom_max_z = np.max(coords[:, 2]) # bohr atom_min_z = np.min(coords[:, 2]) z_shift_to_standard = atom_max_z - np.max(molecule_coords[:, 2]) oldfermi = Units.convert(job.results.readrkf("DOS", "Fermi Energy", file="band"), "hartree", "eV") s, zvalues = _get_coulomb_potential_settings( job, atom_min_z, atom_max_z, lattice, padding_z, step_x, step_y, step_z ) if nproc: s.runscript.nproc = nproc new_job = AMSJob(settings=s, name=f"coulomb_potential_{job.name}") new_job.run() ret = _extract_work_function_from_coulomb_potential_job(new_job, oldfermi, atom_min_z, atom_max_z, zvalues) ret.update(dict(z_shift_to_standard=z_shift_to_standard * 0.529)) # angstrom return ret def create_workfunction_plot(data, fname: str = None, title: str = "", show=False): """ data : dict Dictionary returned by get_workfunction() fname: str File name to save the picture (end in .png) show: bool Whether to call plt.show() The plot will show the potential vs the original z coordinates, i.e. it is shifted from BAND's internal coordinates along z. """ plt.clf() shift = data["z_shift_to_standard"] plt.plot(data["z"] - shift, data["coulomb_potential"]) plt.plot(data["z"] - shift, [data["fermi_energy"]] * len(data["z"])) plt.plot([data["max_pot_bottom_z"] - shift] * 2, [data["fermi_energy"], data["max_pot_bottom"]]) plt.plot([data["max_pot_top_z"] - shift] * 2, [data["fermi_energy"], data["max_pot_top"]]) plt.xlabel("z (Å)") plt.ylabel("Coulomb potential (eV)") plt.title(f'{title}\nWF bottom: {data["work_function_bottom"]:.2f} eV. WF top: {data["work_function_top"]:.2f} eV') if fname: plt.savefig(fname) if show: plt.show() def save_potential_vs_z(data, fname: str): A = np.stack((data["z"] - data["z_shift_to_standard"], data["coulomb_potential"]), axis=1) np.savetxt(fname, A) def _get_coulomb_potential_settings( job, atom_min_z, atom_max_z, lattice, padding_z, step_x, step_y, step_z ) -> Settings: s = job.settings.copy() s.input.ams.Task = "SinglePoint" s.input.ams.LoadSystem.File = job.results.rkfpath(file="engine") s.input.BAND.Restart.File = job.results.rkfpath(file="engine") s.input.BAND.Restart.DensityPlot = "" s.input.BAND.DensityPlot["v(coulomb)"] = "" minz = atom_min_z - padding_z maxz = atom_max_z + padding_z nz = (maxz - minz) / step_z nz = int(nz) + 1 zvalues = np.array([minz + i * step_z for i in range(nz)]) ax, ay = lattice[0][0], lattice[0][1] # bohr bx, by = lattice[1][0], lattice[1][1] nx = (ax**2 + ay**2) ** 0.5 / step_x nx = int(nx) nx = max(nx, 1) ny = (bx**2 + by**2) ** 0.5 / step_y ny = int(ny) ny = max(ny, 1) s.input.BAND.Grid.UserDefined = f""" 0. 0. {minz} {ax} {ay}. 0. {step_x} {bx} {by} 0 {step_y} 0. 0. 1.0 {step_z} {nx} {ny} {nz} End""" return s, zvalues def _extract_work_function_from_coulomb_potential_job( job: AMSJob, fermi: float, atom_min_z: float, atom_max_z: float, zvalues, delta_from_atom_z: float = 8.0 ): kf = KFFile(os.path.join(job.path, "TAPE41")) V = kf.read("FOO", "v(coulomb)") nx = kf.read("Grid", "nr of points x") ny = kf.read("Grid", "nr of points y") nz = kf.read("Grid", "nr of points z") V = np.array(V).reshape((nz, ny, nx)) * 27.211 # eV # ignore points with potential > 100 eV masked = np.where(np.abs(V) < 100, V, np.nan) m0 = np.nanmean(masked, axis=(1, 2)) left_indices = np.where(zvalues < atom_min_z - delta_from_atom_z)[0] right_indices = np.where(zvalues > atom_max_z + delta_from_atom_z)[0] avg_potential_left = np.mean(m0[left_indices]) avg_potential_right = np.mean(m0[right_indices]) if avg_potential_left < 0: min_index = np.argmin(m0[left_indices]) max_pot_left = np.max(m0[left_indices][min_index:]) # using "barrier" max_pot_left_z = zvalues[np.argmax(m0[left_indices][min_index:]) + min_index] # max_pot_left = np.min(m0[left_indices]) # using minimum else: max_pot_left = np.max(m0[left_indices]) max_pot_left_z = zvalues[np.argmax(m0[left_indices])] if avg_potential_right < 0: min_index = np.argmin(m0[right_indices]) min_index = max(1, min_index) max_pot_right = np.max(m0[right_indices][:min_index]) # using "barrier" max_pot_right_z = zvalues[np.argmax(m0[right_indices][:min_index]) + right_indices[0]] # max_pot_right = np.min(m0[right_indices]) # using minimum else: max_pot_right = np.max(m0[right_indices]) max_pot_right_z = zvalues[np.argmax(m0[right_indices]) + right_indices[0]] max_pot = np.max(m0) ret = dict( # max_left_index=left_indies[-1], # min_right_index=right_indices[0], fermi_energy=fermi, max_pot=max_pot, max_pot_bottom=max_pot_left, max_pot_bottom_z=max_pot_left_z * 0.529, max_pot_top=max_pot_right, max_pot_top_z=max_pot_right_z * 0.529, work_function_bottom=max_pot_left - fermi, work_function_top=max_pot_right - fermi, z=zvalues * 0.529, # angstrom coulomb_potential=m0, ) return ret if __name__ == "__main__": init() main() finish()