From 478017adf427da5bedbb0e3c01ff9d0044e128b3 Mon Sep 17 00:00:00 2001 From: Shibu Meher <56470750+Shibu778@users.noreply.github.com> Date: Thu, 3 Oct 2024 17:08:23 +0530 Subject: [PATCH 1/7] Added first CLI --- nonrad/cli.py | 91 +++++++++++++++++++++++++++++++++++++++++++++++++++ setup.py | 85 +++++++++++++++++++++++++++-------------------- 2 files changed, 141 insertions(+), 35 deletions(-) create mode 100644 nonrad/cli.py diff --git a/nonrad/cli.py b/nonrad/cli.py new file mode 100644 index 0000000..0c493ac --- /dev/null +++ b/nonrad/cli.py @@ -0,0 +1,91 @@ +# -*- coding: utf-8 -*- +# Copyright (c) Chris G. Van de Walle +# Distributed under the terms of the MIT License. + +# Author : Shibu Meher +# Date : 2021-06-01 + +import click +import os +import numpy as np +from pathlib import Path +from shutil import copyfile +from pymatgen.core import Structure +from nonrad.ccd import get_cc_structures + + +@click.group() +def nonrad(): + """Command Line Interface for nonrad.""" + pass + + +@click.command("pccd") +@click.argument("ground_path", type=click.Path(exists=True)) +@click.argument("excited_path", type=click.Path(exists=True)) +@click.argument("cc_dir", type=click.Path()) +@click.option( + "--displace", + "-d", + nargs=3, + type=float, + default=[-0.5, 0.5, 9], + help="Displacement range and number of displacements.", +) +def prepare_ccd(ground_path, excited_path, cc_dir, displace=[-0.5, 0.5, 9]): + """ + Prepare the input files for CCD calculation. From the ground and excited state calculations, + the CONTCAR file is read. The displacements are generated and written to the ccd directory. \n + + ground_path : Path to the directory containing the ground state calculation files. \n + excited_path : Path to the directory containing the excited state calculation files. \n + cc_dir : Path to the directory where the CCD input files will be written. \n + displace : List containing the minimum and maximum displacements as a percentage + and the number of displacements to generate. Default is [-0.5, 0.5, 9]. \n + """ + # equilibrium structures from your first-principles calculation + ground_files = Path(ground_path) + ground_struct = Structure.from_file(str(ground_files / "CONTCAR")) + excited_files = Path(excited_path) + excited_struct = Structure.from_file(str(excited_files / "CONTCAR")) + + # output directory that will contain the input files for the CC diagram + cc_dir = Path(cc_dir) + os.makedirs(str(cc_dir), exist_ok=True) + os.makedirs(str(cc_dir / "ground"), exist_ok=True) + os.makedirs(str(cc_dir / "excited"), exist_ok=True) + + # displacements as a percentage, this will generate the displacements + # -50%, -37.5%, -25%, -12.5%, 0%, 12.5%, 25%, 37.5%, 50% + displacements = np.linspace(displace[0], displace[1], int(displace[2])) + + # note: the returned structures won't include the 0% displacement, this is intended + # it can be included by specifying remove_zero=False + ground, excited = get_cc_structures(ground_struct, excited_struct, displacements) + + for i, struct in enumerate(ground): + working_dir = cc_dir / "ground" / str(i) + os.mkdir(str(working_dir)) + + # write structure and copy necessary input files + struct.to(filename=str(working_dir / "POSCAR"), fmt="poscar") + for f in ["KPOINTS", "POTCAR", "INCAR", "job_script.sh"]: + copyfile(str(ground_files / f), str(working_dir / f)) + + for i, struct in enumerate(excited): + working_dir = cc_dir / "excited" / str(i) + os.mkdir(str(working_dir)) + + # write structure and copy necessary input files + struct.to(filename=str(working_dir / "POSCAR"), fmt="poscar") + for f in ["KPOINTS", "POTCAR", "INCAR", "job_script.sh"]: + copyfile(str(excited_files / f), str(working_dir / f)) + + return 0 + + +# Add subcommands to the main group +nonrad.add_command(prepare_ccd) + +if __name__ == "__main__": + nonrad() diff --git a/setup.py b/setup.py index fc949f3..94a88f8 100644 --- a/setup.py +++ b/setup.py @@ -3,54 +3,69 @@ from setuptools import setup, find_packages -with open('nonrad/VERSION', 'r') as f: +with open("nonrad/VERSION", "r") as f: VERSION = f.readline().strip() -with open('README.md', 'r') as f: +with open("README.md", "r") as f: long_desc = f.read() setup( - name='nonrad', + name="nonrad", version=VERSION, - author='Mark E. Turiansky', - author_email='mturiansky@physics.ucsb.edu', - description=('Implementation for computing nonradiative recombination ' - 'rates in semiconductors'), + author="Mark E. Turiansky", + author_email="mturiansky@physics.ucsb.edu", + description=( + "Implementation for computing nonradiative recombination " + "rates in semiconductors" + ), long_description=long_desc, - long_description_content_type='text/markdown', - url='https://github.com/mturiansky/nonrad', + long_description_content_type="text/markdown", + url="https://github.com/mturiansky/nonrad", packages=find_packages(), include_package_data=True, - python_requires='>=3.6', + python_requires=">=3.6", install_requires=[ - 'numpy', - 'scipy', - 'pymatgen>=v2020.6.8', - 'monty', - 'numba>=v0.50.1'], + "numpy", + "scipy", + "pymatgen>=v2020.6.8", + "monty", + "numba>=v0.50.1", + "click", + ], extras_require={ - 'dev': [ - 'pycodestyle', - 'pydocstyle', - 'pylint', - 'flake8', - 'mypy', - 'coverage', - 'pytest', - 'pytest-cov', - 'sphinx', - 'sphinx-rtd-theme' + "dev": [ + "pycodestyle", + "pydocstyle", + "pylint", + "flake8", + "mypy", + "coverage", + "pytest", + "pytest-cov", + "sphinx", + "sphinx-rtd-theme", ], }, - keywords=['physics', 'materials', 'science', 'VASP', 'recombination', - 'Shockley-Read-Hall'], + keywords=[ + "physics", + "materials", + "science", + "VASP", + "recombination", + "Shockley-Read-Hall", + ], classifiers=[ - 'Programming Language :: Python :: 3', - 'Programming Language :: Python :: 3.6', - 'Programming Language :: Python :: 3.7', - 'Programming Language :: Python :: 3.8', - 'Programming Language :: Python :: 3.9', - 'License :: OSI Approved :: MIT License', - 'Operating System :: OS Independent' + "Programming Language :: Python :: 3", + "Programming Language :: Python :: 3.6", + "Programming Language :: Python :: 3.7", + "Programming Language :: Python :: 3.8", + "Programming Language :: Python :: 3.9", + "License :: OSI Approved :: MIT License", + "Operating System :: OS Independent", ], + entry_points={ + "console_scripts": [ + "nonrad=nonrad.cli:nonrad", + ], + }, ) From b2a115aa94b8d6737b4156cd495fa741266e3f3a Mon Sep 17 00:00:00 2001 From: Shibu Meher <56470750+Shibu778@users.noreply.github.com> Date: Thu, 3 Oct 2024 17:38:50 +0530 Subject: [PATCH 2/7] Revert "Added first CLI" This reverts commit 478017adf427da5bedbb0e3c01ff9d0044e128b3. --- nonrad/cli.py | 91 --------------------------------------------------- setup.py | 85 ++++++++++++++++++++--------------------------- 2 files changed, 35 insertions(+), 141 deletions(-) delete mode 100644 nonrad/cli.py diff --git a/nonrad/cli.py b/nonrad/cli.py deleted file mode 100644 index 0c493ac..0000000 --- a/nonrad/cli.py +++ /dev/null @@ -1,91 +0,0 @@ -# -*- coding: utf-8 -*- -# Copyright (c) Chris G. Van de Walle -# Distributed under the terms of the MIT License. - -# Author : Shibu Meher -# Date : 2021-06-01 - -import click -import os -import numpy as np -from pathlib import Path -from shutil import copyfile -from pymatgen.core import Structure -from nonrad.ccd import get_cc_structures - - -@click.group() -def nonrad(): - """Command Line Interface for nonrad.""" - pass - - -@click.command("pccd") -@click.argument("ground_path", type=click.Path(exists=True)) -@click.argument("excited_path", type=click.Path(exists=True)) -@click.argument("cc_dir", type=click.Path()) -@click.option( - "--displace", - "-d", - nargs=3, - type=float, - default=[-0.5, 0.5, 9], - help="Displacement range and number of displacements.", -) -def prepare_ccd(ground_path, excited_path, cc_dir, displace=[-0.5, 0.5, 9]): - """ - Prepare the input files for CCD calculation. From the ground and excited state calculations, - the CONTCAR file is read. The displacements are generated and written to the ccd directory. \n - - ground_path : Path to the directory containing the ground state calculation files. \n - excited_path : Path to the directory containing the excited state calculation files. \n - cc_dir : Path to the directory where the CCD input files will be written. \n - displace : List containing the minimum and maximum displacements as a percentage - and the number of displacements to generate. Default is [-0.5, 0.5, 9]. \n - """ - # equilibrium structures from your first-principles calculation - ground_files = Path(ground_path) - ground_struct = Structure.from_file(str(ground_files / "CONTCAR")) - excited_files = Path(excited_path) - excited_struct = Structure.from_file(str(excited_files / "CONTCAR")) - - # output directory that will contain the input files for the CC diagram - cc_dir = Path(cc_dir) - os.makedirs(str(cc_dir), exist_ok=True) - os.makedirs(str(cc_dir / "ground"), exist_ok=True) - os.makedirs(str(cc_dir / "excited"), exist_ok=True) - - # displacements as a percentage, this will generate the displacements - # -50%, -37.5%, -25%, -12.5%, 0%, 12.5%, 25%, 37.5%, 50% - displacements = np.linspace(displace[0], displace[1], int(displace[2])) - - # note: the returned structures won't include the 0% displacement, this is intended - # it can be included by specifying remove_zero=False - ground, excited = get_cc_structures(ground_struct, excited_struct, displacements) - - for i, struct in enumerate(ground): - working_dir = cc_dir / "ground" / str(i) - os.mkdir(str(working_dir)) - - # write structure and copy necessary input files - struct.to(filename=str(working_dir / "POSCAR"), fmt="poscar") - for f in ["KPOINTS", "POTCAR", "INCAR", "job_script.sh"]: - copyfile(str(ground_files / f), str(working_dir / f)) - - for i, struct in enumerate(excited): - working_dir = cc_dir / "excited" / str(i) - os.mkdir(str(working_dir)) - - # write structure and copy necessary input files - struct.to(filename=str(working_dir / "POSCAR"), fmt="poscar") - for f in ["KPOINTS", "POTCAR", "INCAR", "job_script.sh"]: - copyfile(str(excited_files / f), str(working_dir / f)) - - return 0 - - -# Add subcommands to the main group -nonrad.add_command(prepare_ccd) - -if __name__ == "__main__": - nonrad() diff --git a/setup.py b/setup.py index 94a88f8..fc949f3 100644 --- a/setup.py +++ b/setup.py @@ -3,69 +3,54 @@ from setuptools import setup, find_packages -with open("nonrad/VERSION", "r") as f: +with open('nonrad/VERSION', 'r') as f: VERSION = f.readline().strip() -with open("README.md", "r") as f: +with open('README.md', 'r') as f: long_desc = f.read() setup( - name="nonrad", + name='nonrad', version=VERSION, - author="Mark E. Turiansky", - author_email="mturiansky@physics.ucsb.edu", - description=( - "Implementation for computing nonradiative recombination " - "rates in semiconductors" - ), + author='Mark E. Turiansky', + author_email='mturiansky@physics.ucsb.edu', + description=('Implementation for computing nonradiative recombination ' + 'rates in semiconductors'), long_description=long_desc, - long_description_content_type="text/markdown", - url="https://github.com/mturiansky/nonrad", + long_description_content_type='text/markdown', + url='https://github.com/mturiansky/nonrad', packages=find_packages(), include_package_data=True, - python_requires=">=3.6", + python_requires='>=3.6', install_requires=[ - "numpy", - "scipy", - "pymatgen>=v2020.6.8", - "monty", - "numba>=v0.50.1", - "click", - ], + 'numpy', + 'scipy', + 'pymatgen>=v2020.6.8', + 'monty', + 'numba>=v0.50.1'], extras_require={ - "dev": [ - "pycodestyle", - "pydocstyle", - "pylint", - "flake8", - "mypy", - "coverage", - "pytest", - "pytest-cov", - "sphinx", - "sphinx-rtd-theme", + 'dev': [ + 'pycodestyle', + 'pydocstyle', + 'pylint', + 'flake8', + 'mypy', + 'coverage', + 'pytest', + 'pytest-cov', + 'sphinx', + 'sphinx-rtd-theme' ], }, - keywords=[ - "physics", - "materials", - "science", - "VASP", - "recombination", - "Shockley-Read-Hall", - ], + keywords=['physics', 'materials', 'science', 'VASP', 'recombination', + 'Shockley-Read-Hall'], classifiers=[ - "Programming Language :: Python :: 3", - "Programming Language :: Python :: 3.6", - "Programming Language :: Python :: 3.7", - "Programming Language :: Python :: 3.8", - "Programming Language :: Python :: 3.9", - "License :: OSI Approved :: MIT License", - "Operating System :: OS Independent", + 'Programming Language :: Python :: 3', + 'Programming Language :: Python :: 3.6', + 'Programming Language :: Python :: 3.7', + 'Programming Language :: Python :: 3.8', + 'Programming Language :: Python :: 3.9', + 'License :: OSI Approved :: MIT License', + 'Operating System :: OS Independent' ], - entry_points={ - "console_scripts": [ - "nonrad=nonrad.cli:nonrad", - ], - }, ) From 314aad50167f4d7fd0f5584bbb059595c2dd7812 Mon Sep 17 00:00:00 2001 From: Shibu Meher <56470750+Shibu778@users.noreply.github.com> Date: Thu, 3 Oct 2024 17:53:36 +0530 Subject: [PATCH 3/7] Added first CLI feature --- nonrad/cli.py | 91 ++++++++++++++++++++++++++++++++++++++++++++++++++ pyproject.toml | 4 +++ 2 files changed, 95 insertions(+) create mode 100644 nonrad/cli.py diff --git a/nonrad/cli.py b/nonrad/cli.py new file mode 100644 index 0000000..b898545 --- /dev/null +++ b/nonrad/cli.py @@ -0,0 +1,91 @@ +# -*- coding: utf-8 -*- +# Copyright (c) Chris G. Van de Walle +# Distributed under the terms of the MIT License. + +# Author : Shibu Meher +# Date : 2021-06-01 + +import click +import os +import numpy as np +from pathlib import Path +from shutil import copyfile +from pymatgen.core import Structure +from nonrad.ccd import get_cc_structures + + +@click.group() +def nonrad(): + """Command Line Interface for nonrad.""" + pass + + +@click.command("pccd") +@click.argument("ground_path", type=click.Path(exists=True)) +@click.argument("excited_path", type=click.Path(exists=True)) +@click.argument("cc_dir", type=click.Path()) +@click.option( + "--displace", + "-d", + nargs=3, + type=float, + default=[-0.5, 0.5, 9], + help="Displacement range and number of displacements.", +) +def prepare_ccd(ground_path, excited_path, cc_dir, displace=[-0.5, 0.5, 9]): + """ + Prepare the input files for CCD calculation. From the ground and excited state calculations, + the CONTCAR file is read. The displacements are generated and written to the ccd directory. \n + + ground_path : Path to the directory containing the ground state calculation files. \n + excited_path : Path to the directory containing the excited state calculation files. \n + cc_dir : Path to the directory where the CCD input files will be written. \n + displace : List containing the minimum and maximum displacements as a percentage + and the number of displacements to generate. Default is [-0.5, 0.5, 9]. \n + """ + # equilibrium structures from your first-principles calculation + ground_files = Path(ground_path) + ground_struct = Structure.from_file(str(ground_files / "CONTCAR")) + excited_files = Path(excited_path) + excited_struct = Structure.from_file(str(excited_files / "CONTCAR")) + + # output directory that will contain the input files for the CC diagram + cc_dir = Path(cc_dir) + os.makedirs(str(cc_dir), exist_ok=True) + os.makedirs(str(cc_dir / "ground"), exist_ok=True) + os.makedirs(str(cc_dir / "excited"), exist_ok=True) + + # displacements as a percentage, this will generate the displacements + # -50%, -37.5%, -25%, -12.5%, 0%, 12.5%, 25%, 37.5%, 50% + displacements = np.linspace(displace[0], displace[1], int(displace[2])) + + # note: the returned structures won't include the 0% displacement, this is intended + # it can be included by specifying remove_zero=False + ground, excited = get_cc_structures(ground_struct, excited_struct, displacements) + + for i, struct in enumerate(ground): + working_dir = cc_dir / "ground" / str(i) + os.makedirs(str(working_dir), exist_ok=True) + + # write structure and copy necessary input files + struct.to(filename=str(working_dir / "POSCAR"), fmt="poscar") + for f in ["KPOINTS", "POTCAR", "INCAR", "job_script.sh"]: + copyfile(str(ground_files / f), str(working_dir / f)) + + for i, struct in enumerate(excited): + working_dir = cc_dir / "excited" / str(i) + os.makedirs(str(working_dir), exist_ok=True) + + # write structure and copy necessary input files + struct.to(filename=str(working_dir / "POSCAR"), fmt="poscar") + for f in ["KPOINTS", "POTCAR", "INCAR", "job_script.sh"]: + copyfile(str(excited_files / f), str(working_dir / f)) + + return 0 + + +# Add subcommands to the main group +nonrad.add_command(prepare_ccd) + +if __name__ == "__main__": + nonrad() diff --git a/pyproject.toml b/pyproject.toml index 8926ed0..b896cfb 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -53,6 +53,7 @@ dependencies = [ "monty>=2023.4.10", "numba>=0.57.0", "mpmath>=1.3.0", + "click>=8.0.1", ] [project.urls] @@ -83,3 +84,6 @@ extend-exclude = ["docs"] [tool.mypy] ignore_missing_imports = "True" plugins = "numpy.typing.mypy_plugin" + +[project.scripts] +nonrad = "nonrad.cli:nonrad" \ No newline at end of file From 07e66e3ab0e5e5b8acc2ec1092777e159c822010 Mon Sep 17 00:00:00 2001 From: Shibu Meher <56470750+Shibu778@users.noreply.github.com> Date: Sat, 5 Oct 2024 16:27:35 +0530 Subject: [PATCH 4/7] Added PES related CLI --- nonrad/cli.py | 93 ++++++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 89 insertions(+), 4 deletions(-) diff --git a/nonrad/cli.py b/nonrad/cli.py index b898545..cee9f95 100644 --- a/nonrad/cli.py +++ b/nonrad/cli.py @@ -11,7 +11,13 @@ from pathlib import Path from shutil import copyfile from pymatgen.core import Structure -from nonrad.ccd import get_cc_structures +from nonrad.ccd import ( + get_cc_structures, + get_dQ, + get_PES_from_vaspruns, + get_omega_from_PES, +) +from glob import glob @click.group() @@ -20,7 +26,7 @@ def nonrad(): pass -@click.command("pccd") +@click.command() @click.argument("ground_path", type=click.Path(exists=True)) @click.argument("excited_path", type=click.Path(exists=True)) @click.argument("cc_dir", type=click.Path()) @@ -32,7 +38,7 @@ def nonrad(): default=[-0.5, 0.5, 9], help="Displacement range and number of displacements.", ) -def prepare_ccd(ground_path, excited_path, cc_dir, displace=[-0.5, 0.5, 9]): +def prep_ccd(ground_path, excited_path, cc_dir, displace=[-0.5, 0.5, 9]): """ Prepare the input files for CCD calculation. From the ground and excited state calculations, the CONTCAR file is read. The displacements are generated and written to the ccd directory. \n @@ -84,8 +90,87 @@ def prepare_ccd(ground_path, excited_path, cc_dir, displace=[-0.5, 0.5, 9]): return 0 +@click.command() +@click.argument("cc_dir", type=click.Path()) +@click.argument("ground_files", type=click.Path()) +@click.argument("excited_files", type=click.Path()) +@click.option( + "--plot", + "-p", + is_flag=True, + help="Plot the potential energy surfaces.", +) +@click.option( + "--plot_name", + "-n", + default="pes.png", + help="Name of the plot file. Default is `pes.png`.", +) +def pes(cc_dir, ground_files, excited_files, plot=False, plot_name="pes.png"): + """Extracting potential energy surfaces and relevant parameters from CCD calculations. + + Parmeters: \n + ----------- \n + cc_dir : str \n + Path to the directory containing the CCD calculations. It should contain the `ground` and `excited` directories. + Each of these directories should contain the directories corresponding to the displacements. \n + ground_files : str \n + Path to the directory containing the ground state calculation files. It should contain the `CONTCAR` and `vasprun.xml` files + for 0% displacement. \n + excited_files : str \n + Path to the directory containing the excited state calculation files. It should contain the `CONTCAR` and `vasprun.xml` files + for 0% displacement. \n + plot : bool \n + If True, the potential energy surfaces are plotted. Default is False. \n + plot_name : str \n + Name of the plot file. Default is `pes.png`. + + """ + ground_struct = Structure.from_file(Path(ground_files, "CONTCAR")) + excited_struct = Structure.from_file(Path(excited_files, "CONTCAR")) + + # calculate dQ + dQ = get_dQ(ground_struct, excited_struct) # amu^{1/2} Angstrom + + # this prepares a list of all vasprun.xml's from the CCD calculations + ground_vaspruns = glob(Path(cc_dir, "ground", "*", "vasprun.xml")) + excited_vaspruns = glob(Path(cc_dir, "excited", "*", "vasprun.xml")) + + # remember that the 0% displacement was removed before? we need to add that back in here + ground_vaspruns = ground_vaspruns + [Path(ground_files, "vasprun.xml")] + excited_vaspruns = excited_vaspruns + [Path(excited_files, "vasprun.xml")] + + # extract the potential energy surface + Q_ground, E_ground = get_PES_from_vaspruns( + ground_struct, excited_struct, ground_vaspruns + ) + Q_excited, E_excited = get_PES_from_vaspruns( + ground_struct, excited_struct, excited_vaspruns + ) + + # the energy surfaces are referenced to the minimums, so we need to add dE (defined before) to E_excited + E_excited = dE + E_excited + + if plot: + fig, ax = plt.subplots(figsize=(5, 5)) + ax.scatter(Q_ground, E_ground, s=10) + ax.scatter(Q_excited, E_excited, s=10) + + # by passing in the axis object, it also plots the fitted curve + q = np.linspace(-1.0, 3.5, 100) + ground_omega = get_omega_from_PES(Q_ground, E_ground, ax=ax, q=q) + excited_omega = get_omega_from_PES(Q_excited, E_excited, ax=ax, q=q) + + ax.set_xlabel("$Q$ [amu$^{1/2}$ $\AA$]") + ax.set_ylabel("$E$ [eV]") + plt.savefig(plot_name, dpi=300, format=plot_name.split(".")[-1]) + + return 0 + + # Add subcommands to the main group -nonrad.add_command(prepare_ccd) +nonrad.add_command(prep_ccd) +nonrad.add_command(pes) if __name__ == "__main__": nonrad() From f7244074d695bd2a639d70fa0410d13ba9e0fab7 Mon Sep 17 00:00:00 2001 From: Shibu Meher <56470750+Shibu778@users.noreply.github.com> Date: Sat, 5 Oct 2024 16:45:44 +0530 Subject: [PATCH 5/7] Printing the unused variable --- nonrad/cli.py | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/nonrad/cli.py b/nonrad/cli.py index cee9f95..5bccfe1 100644 --- a/nonrad/cli.py +++ b/nonrad/cli.py @@ -94,12 +94,7 @@ def prep_ccd(ground_path, excited_path, cc_dir, displace=[-0.5, 0.5, 9]): @click.argument("cc_dir", type=click.Path()) @click.argument("ground_files", type=click.Path()) @click.argument("excited_files", type=click.Path()) -@click.option( - "--plot", - "-p", - is_flag=True, - help="Plot the potential energy surfaces.", -) +@click.option("--plot", "-p", is_flag=True, help="Plot the potential energy surfaces.") @click.option( "--plot_name", "-n", @@ -131,6 +126,7 @@ def pes(cc_dir, ground_files, excited_files, plot=False, plot_name="pes.png"): # calculate dQ dQ = get_dQ(ground_struct, excited_struct) # amu^{1/2} Angstrom + print(f"dQ = {dQ:.2f} amu^{1/2} Angstrom") # this prepares a list of all vasprun.xml's from the CCD calculations ground_vaspruns = glob(Path(cc_dir, "ground", "*", "vasprun.xml")) @@ -151,20 +147,22 @@ def pes(cc_dir, ground_files, excited_files, plot=False, plot_name="pes.png"): # the energy surfaces are referenced to the minimums, so we need to add dE (defined before) to E_excited E_excited = dE + E_excited + # calculate omega if plot: fig, ax = plt.subplots(figsize=(5, 5)) ax.scatter(Q_ground, E_ground, s=10) ax.scatter(Q_excited, E_excited, s=10) - # by passing in the axis object, it also plots the fitted curve q = np.linspace(-1.0, 3.5, 100) ground_omega = get_omega_from_PES(Q_ground, E_ground, ax=ax, q=q) excited_omega = get_omega_from_PES(Q_excited, E_excited, ax=ax, q=q) - ax.set_xlabel("$Q$ [amu$^{1/2}$ $\AA$]") ax.set_ylabel("$E$ [eV]") plt.savefig(plot_name, dpi=300, format=plot_name.split(".")[-1]) + print(f"Ground state omega = {ground_omega:.2f} eV") + print(f"Excited state omega = {excited_omega:.2f} eV") + return 0 From be2162710ea2c19f89b243f01874e13d26b38a17 Mon Sep 17 00:00:00 2001 From: Shibu Meher <56470750+Shibu778@users.noreply.github.com> Date: Sat, 5 Oct 2024 16:54:50 +0530 Subject: [PATCH 6/7] Fixed 3 codacy issues --- nonrad/cli.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/nonrad/cli.py b/nonrad/cli.py index 5bccfe1..98cb9a6 100644 --- a/nonrad/cli.py +++ b/nonrad/cli.py @@ -23,7 +23,6 @@ @click.group() def nonrad(): """Command Line Interface for nonrad.""" - pass @click.command() @@ -38,7 +37,7 @@ def nonrad(): default=[-0.5, 0.5, 9], help="Displacement range and number of displacements.", ) -def prep_ccd(ground_path, excited_path, cc_dir, displace=[-0.5, 0.5, 9]): +def prep_ccd(ground_path, excited_path, cc_dir, displace=None): """ Prepare the input files for CCD calculation. From the ground and excited state calculations, the CONTCAR file is read. The displacements are generated and written to the ccd directory. \n @@ -63,6 +62,8 @@ def prep_ccd(ground_path, excited_path, cc_dir, displace=[-0.5, 0.5, 9]): # displacements as a percentage, this will generate the displacements # -50%, -37.5%, -25%, -12.5%, 0%, 12.5%, 25%, 37.5%, 50% + if displace is None: + displace = [-0.5, 0.5, 9] displacements = np.linspace(displace[0], displace[1], int(displace[2])) # note: the returned structures won't include the 0% displacement, this is intended @@ -149,7 +150,7 @@ def pes(cc_dir, ground_files, excited_files, plot=False, plot_name="pes.png"): # calculate omega if plot: - fig, ax = plt.subplots(figsize=(5, 5)) + _, ax = plt.subplots(figsize=(5, 5)) ax.scatter(Q_ground, E_ground, s=10) ax.scatter(Q_excited, E_excited, s=10) # by passing in the axis object, it also plots the fitted curve From b3544ee582a541e64125ae66052e3e88cddab3b2 Mon Sep 17 00:00:00 2001 From: Shibu Meher <56470750+Shibu778@users.noreply.github.com> Date: Tue, 16 Jun 2026 14:22:49 +0530 Subject: [PATCH 7/7] CLI, tests and docs included --- docs/cli.rst | 270 ++++++++++++ docs/index.rst | 1 + nonrad/cli.py | 916 ++++++++++++++++++++++++++++++++++----- nonrad/tests/test_cli.py | 606 ++++++++++++++++++++++++++ 4 files changed, 1677 insertions(+), 116 deletions(-) create mode 100644 docs/cli.rst create mode 100644 nonrad/tests/test_cli.py diff --git a/docs/cli.rst b/docs/cli.rst new file mode 100644 index 0000000..3ae7f8e --- /dev/null +++ b/docs/cli.rst @@ -0,0 +1,270 @@ +====================== +Command Line Interface +====================== + +``nonrad`` includes a comprehensive Command Line Interface (CLI) built with `Click `_. This tool allows you to perform CCD calculations, compute capture coefficients, scale coefficients, and extract electron-phonon coupling parameters directly from your terminal. + +Usage +===== + +The entry point for all commands is ``nonrad``. You can view the list of commands and general help by running: + +.. code-block:: sh + + nonrad --help + +To view the help and option details for a specific command, use: + +.. code-block:: sh + + nonrad --help + +CCD Commands +============ + +These commands are used to prepare and process configuration-coordinate diagrams (CCDs). + +prep-ccd +-------- +Prepare displaced structures along the configuration coordinate. + +**Usage:** + +.. code-block:: sh + + nonrad prep-ccd [OPTIONS] GROUND_PATH EXCITED_PATH CC_DIR + +**Arguments:** +* ``GROUND_PATH``: Path to the directory containing the ground-state calculation (with ``CONTCAR``). +* ``EXCITED_PATH``: Path to the directory containing the excited-state calculation (with ``CONTCAR``). +* ``CC_DIR``: Path to the directory where displaced structures should be written. + +**Options:** +* ``-d, --displace ``: Displacement range specified as ``(min, max, count)``. (Default: ``-0.5 0.5 9``). +* ``-f, --input-files TEXT``: Comma-separated list of VASP input files to copy into each displacement folder. (Default: ``KPOINTS,POTCAR,INCAR,job_script.sh``). + +pes +--- +Extract potential energy surfaces (PES) from CCD calculations. + +**Usage:** + +.. code-block:: sh + + nonrad pes [OPTIONS] CC_DIR GROUND_FILES EXCITED_FILES + +**Arguments:** +* ``CC_DIR``: Root directory containing the CCD displacement folders. +* ``GROUND_FILES``: Path to the ground-state equilibrium calculation directory (containing ``CONTCAR`` and ``vasprun.xml``). +* ``EXCITED_FILES``: Path to the excited-state equilibrium calculation directory (containing ``CONTCAR`` and ``vasprun.xml``). + +**Options:** +* ``-e, --energy-diff FLOAT``: [Required] Energy difference *dE* between ground and excited states (in eV). +* ``-p, --plot``: Flag to plot the potential energy surfaces and save the figure. +* ``-n, --plot-name TEXT``: Filename for the saved plot. (Default: ``pes.png``). + +dq +-- +Compute the configuration-coordinate displacement *dQ* between two structures. + +**Usage:** + +.. code-block:: sh + + nonrad dq GROUND_PATH EXCITED_PATH + +**Arguments:** +* ``GROUND_PATH``: Path to the ground-state structure file (e.g. ``POSCAR``, ``CONTCAR``). +* ``EXCITED_PATH``: Path to the excited-state structure file (e.g. ``POSCAR``, ``CONTCAR``). + +q-from-struct +------------- +Compute the *Q* value of an arbitrary structure relative to ground and excited endpoint structures. + +**Usage:** + +.. code-block:: sh + + nonrad q-from-struct [OPTIONS] GROUND_PATH EXCITED_PATH STRUCT_PATH + +**Arguments:** +* ``GROUND_PATH``: Path to the ground-state structure file. +* ``EXCITED_PATH``: Path to the excited-state structure file. +* ``STRUCT_PATH``: Path to the structure file to evaluate. + +**Options:** +* ``-t, --tol FLOAT``: Distance tolerance for filtering coordinate differences. (Default: ``1e-4``). +* ``--nround INTEGER``: Number of decimal places to round *Q*. (Default: ``5``). + +barrier +------- +Compute the classical barrier height using the harmonic approximation. + +**Usage:** + +.. code-block:: sh + + nonrad barrier [OPTIONS] + +**Options (All Required):** +* ``--dq FLOAT``: Displacement *dQ* (in amu\ :sup:`1/2` Å). +* ``--de FLOAT``: Energy offset *dE* between states (in eV). +* ``--wi FLOAT``: Initial-state frequency (in eV). +* ``--wf FLOAT``: Final-state frequency (in eV). + + +Capture Commands +================ + +capture +------- +Evaluate the nonradiative capture coefficient *C* (in cm\ :sup:`3` s\ :sup:`-1`). + +**Usage:** + +.. code-block:: sh + + nonrad capture [OPTIONS] + +**Options:** +* ``--dq FLOAT``: [Required] Configuration displacement *dQ* (in amu\ :sup:`1/2` Å). +* ``--de FLOAT``: [Required] Energy offset *dE* (in eV). +* ``--wi FLOAT``: [Required] Initial-state frequency (in eV). +* ``--wf FLOAT``: [Required] Final-state frequency (in eV). +* ``--wif FLOAT``: [Required] Electron-phonon coupling matrix element *Wif* (in eV amu\ :sup:`-1/2` Å\ :sup:`-1`). +* ``--volume FLOAT``: [Required] Supercell volume (in Å\ :sup:`3`). +* ``--g INTEGER``: Degeneracy factor. (Default: ``1``). +* ``-T, --temperature FLOAT``: Temperature in Kelvin for a single-point calculation. (Default: ``300.0``). +* ``--temperature-range ``: Temperature range as ``(min, max, count)`` for scanning multiple temperatures. +* ``--sigma TEXT``: Smearing parameter (float) or interpolation method (``pchip`` or ``cubic``). (Default: ``pchip``). +* ``--occ-tol FLOAT``: Bose-weight occupation tolerance. (Default: ``1e-5``). +* ``--overlap-method [Analytic|Integral|HermiteGauss]``: Method for evaluating vibrational overlaps. (Default: ``HermiteGauss``). + + +Scaling Commands +================ + +sommerfeld +---------- +Compute the Sommerfeld enhancement factor. + +**Usage:** + +.. code-block:: sh + + nonrad sommerfeld [OPTIONS] + +**Options:** +* ``-T, --temperature FLOAT``: Temperature in Kelvin. (Default: ``300.0``). +* ``--temperature-range ``: Temperature range as ``(min, max, count)``. +* ``--z INTEGER``: [Required] Charge ratio *Z = Q / q* (negative values represent attractive potential). +* ``--m-eff FLOAT``: [Required] Carrier effective mass in units of electron mass. +* ``--eps0 FLOAT``: [Required] Static dielectric constant. +* ``--dim [1|2|3]``: Dimensionality of the system. (Default: ``3``). +* ``--method [Integrate|Analytic]``: Evaluation method. (Default: ``Integrate``). + +thermal-velocity +---------------- +Compute the thermal velocity of a carrier. + +**Usage:** + +.. code-block:: sh + + nonrad thermal-velocity [OPTIONS] + +**Options:** +* ``-T, --temperature FLOAT``: Temperature in Kelvin. (Default: ``300.0``). +* ``--temperature-range ``: Temperature range as ``(min, max, count)``. +* ``--m-eff FLOAT``: [Required] Carrier effective mass in units of electron mass. + +charged-supercell +----------------- +Estimate the charged-supercell scaling factor for VASP calculations. + +**Usage:** + +.. code-block:: sh + + nonrad charged-supercell [OPTIONS] WAVECAR_PATH + +**Arguments:** +* ``WAVECAR_PATH``: Path to the VASP WAVECAR file. + +**Options:** +* ``--bulk-index INTEGER``: [Required] 1-based band index of the bulk wavefunction. +* ``--def-index INTEGER``: 1-based band index of the defect wavefunction. Provide this or ``--def-coord``. (Default: ``-1``). +* ``--def-coord ``: Cartesian coordinates (x, y, z) of the defect position. +* ``--cutoff FLOAT``: Slope cutoff for plateau detection. (Default: ``0.02``). +* ``--limit FLOAT``: Upper radial limit (in Å) for the fitting window. (Default: ``5.0``). +* ``--spin INTEGER``: Spin channel (0 = up, 1 = down). (Default: ``0``). +* ``--kpoint INTEGER``: 1-based k-point index. (Default: ``1``). + + +Electron-Phonon Commands +======================== + +These commands are nested under the ``elphon`` subcommand group and compute the electron-phonon matrix element *Wif*. + +elphon wavecars +--------------- +Compute *Wif* using VASP WAVECAR files along the CCD path. + +**Usage:** + +.. code-block:: sh + + nonrad elphon wavecars [OPTIONS] CC_DIR GROUND_PATH EXCITED_PATH INIT_WAVECAR + +**Arguments:** +* ``CC_DIR``: Directory containing the CCD displacement folders. +* ``GROUND_PATH``: Path to the ground-state equilibrium directory. +* ``EXCITED_PATH``: Path to the excited-state equilibrium directory. +* ``INIT_WAVECAR``: Path to the initial (reference) WAVECAR file. + +**Options:** +* ``--def-index INTEGER``: [Required] 1-based defect wavefunction index. +* ``-b, --bulk-index INTEGER``: [Required] 1-based bulk wavefunction index (repeatable for multiple bands). +* ``--spin INTEGER``: Spin channel (0 = up, 1 = down). (Default: ``0``). +* ``--kpoint INTEGER``: 1-based k-point index. (Default: ``1``). + +elphon wswq +----------- +Compute *Wif* using VASP WSWQ files along the CCD path. + +**Usage:** + +.. code-block:: sh + + nonrad elphon wswq [OPTIONS] CC_DIR GROUND_PATH EXCITED_PATH INIT_VASPRUN + +**Arguments:** +* ``CC_DIR``: Directory containing the CCD displacement folders. +* ``GROUND_PATH``: Path to the ground-state equilibrium directory. +* ``EXCITED_PATH``: Path to the excited-state equilibrium directory. +* ``INIT_VASPRUN``: Path to the reference ``vasprun.xml`` (used to extract eigenvalues). + +**Options:** +* ``--def-index INTEGER``: [Required] 1-based defect wavefunction index. +* ``-b, --bulk-index INTEGER``: [Required] 1-based bulk wavefunction index (repeatable for multiple bands). +* ``--spin INTEGER``: Spin channel (0 = up, 1 = down). (Default: ``0``). +* ``--kpoint INTEGER``: 1-based k-point index. (Default: ``1``). + +elphon unk +---------- +Compute *Wif* using UNK files. + +**Usage:** + +.. code-block:: sh + + nonrad elphon unk [OPTIONS] INIT_UNK + +**Arguments:** +* ``INIT_UNK``: Path to the initial (reference) UNK file. + +**Options:** +* ``--def-index INTEGER``: [Required] 1-based defect wavefunction index. +* ``-b, --bulk-index INTEGER``: [Required] 1-based bulk wavefunction index (repeatable). +* ``--eigs TEXT``: [Required] Comma-separated list of eigenvalues (in eV). +* ``-u, --unk-pair ``: [Required] A ``(Q, unk_path)`` pair. (Repeatable for all displacements). diff --git a/docs/index.rst b/docs/index.rst index d9ec2d2..30a43d5 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -2,6 +2,7 @@ :hidden: installation + cli compatibility citation Tutorial diff --git a/nonrad/cli.py b/nonrad/cli.py index 98cb9a6..3d6a850 100644 --- a/nonrad/cli.py +++ b/nonrad/cli.py @@ -1,175 +1,859 @@ -# -*- coding: utf-8 -*- # Copyright (c) Chris G. Van de Walle # Distributed under the terms of the MIT License. -# Author : Shibu Meher -# Date : 2021-06-01 +"""Command-line interface for nonrad. + +This module provides a Click-based CLI exposing the main functionality of +the ``nonrad`` package, including configuration-coordinate diagram (CCD) +utilities, the nonradiative capture coefficient calculation, scaling +helpers, and electron-phonon coupling tools. + +Commands +-------- +prep-ccd + Prepare displaced structures for a CCD calculation. +pes + Extract potential energy surfaces from CCD VASP outputs. +dq + Compute the configuration coordinate displacement *dQ*. +q-from-struct + Compute the *Q* value for an arbitrary structure. +barrier + Compute the classical barrier height in the harmonic approximation. +capture + Evaluate the nonradiative capture coefficient. +sommerfeld + Compute the Sommerfeld enhancement factor. +thermal-velocity + Compute the thermal velocity of a carrier. +charged-supercell + Estimate the charged-supercell scaling factor. +elphon wavecars + Compute Wif from WAVECAR files. +elphon wswq + Compute Wif from WSWQ files. +elphon unk + Compute Wif from UNK files. +""" -import click import os -import numpy as np +from glob import glob from pathlib import Path from shutil import copyfile -from pymatgen.core import Structure -from nonrad.ccd import ( - get_cc_structures, - get_dQ, - get_PES_from_vaspruns, - get_omega_from_PES, -) -from glob import glob +import click +import numpy as np + + +# --------------------------------------------------------------------------- +# Main CLI group +# --------------------------------------------------------------------------- @click.group() +@click.version_option(package_name='nonrad') def nonrad(): - """Command Line Interface for nonrad.""" + """Nonradiative recombination coefficient calculator. + + A command-line tool for computing nonradiative capture coefficients + and related quantities from first-principles calculations. + """ -@click.command() -@click.argument("ground_path", type=click.Path(exists=True)) -@click.argument("excited_path", type=click.Path(exists=True)) -@click.argument("cc_dir", type=click.Path()) +# --------------------------------------------------------------------------- +# CCD commands +# --------------------------------------------------------------------------- + +@nonrad.command(name='prep-ccd') +@click.argument('ground_path', type=click.Path(exists=True)) +@click.argument('excited_path', type=click.Path(exists=True)) +@click.argument('cc_dir', type=click.Path()) @click.option( - "--displace", - "-d", - nargs=3, - type=float, - default=[-0.5, 0.5, 9], - help="Displacement range and number of displacements.", + '--displace', '-d', + nargs=3, type=float, default=(-0.5, 0.5, 9), + show_default=True, + help='Displacement range (min, max, count).', ) -def prep_ccd(ground_path, excited_path, cc_dir, displace=None): - """ - Prepare the input files for CCD calculation. From the ground and excited state calculations, - the CONTCAR file is read. The displacements are generated and written to the ccd directory. \n - - ground_path : Path to the directory containing the ground state calculation files. \n - excited_path : Path to the directory containing the excited state calculation files. \n - cc_dir : Path to the directory where the CCD input files will be written. \n - displace : List containing the minimum and maximum displacements as a percentage - and the number of displacements to generate. Default is [-0.5, 0.5, 9]. \n +@click.option( + '--input-files', '-f', + type=str, default='KPOINTS,POTCAR,INCAR,job_script.sh', + show_default=True, + help='Comma-separated list of input files to copy into each sub-directory.', +) +def prep_ccd(ground_path, excited_path, cc_dir, displace, input_files): + """Prepare displaced structures for a configuration-coordinate diagram. + + Reads CONTCAR from GROUND_PATH and EXCITED_PATH, generates displaced + structures, writes POSCAR files into CC_DIR/ground/N and + CC_DIR/excited/N sub-directories, and copies the requested input files. + + Parameters + ---------- + ground_path : str + Directory containing the ground-state CONTCAR. + excited_path : str + Directory containing the excited-state CONTCAR. + cc_dir : str + Output directory for the CCD sub-directories. + displace : tuple of float + ``(min, max, count)`` passed to ``numpy.linspace``. + input_files : str + Comma-separated filenames to copy from each source directory. """ - # equilibrium structures from your first-principles calculation + from pymatgen.core import Structure + + from nonrad.ccd import get_cc_structures + ground_files = Path(ground_path) - ground_struct = Structure.from_file(str(ground_files / "CONTCAR")) excited_files = Path(excited_path) - excited_struct = Structure.from_file(str(excited_files / "CONTCAR")) + ground_struct = Structure.from_file(str(ground_files / 'CONTCAR')) + excited_struct = Structure.from_file(str(excited_files / 'CONTCAR')) - # output directory that will contain the input files for the CC diagram cc_dir = Path(cc_dir) os.makedirs(str(cc_dir), exist_ok=True) - os.makedirs(str(cc_dir / "ground"), exist_ok=True) - os.makedirs(str(cc_dir / "excited"), exist_ok=True) + os.makedirs(str(cc_dir / 'ground'), exist_ok=True) + os.makedirs(str(cc_dir / 'excited'), exist_ok=True) - # displacements as a percentage, this will generate the displacements - # -50%, -37.5%, -25%, -12.5%, 0%, 12.5%, 25%, 37.5%, 50% - if displace is None: - displace = [-0.5, 0.5, 9] displacements = np.linspace(displace[0], displace[1], int(displace[2])) + ground, excited = get_cc_structures(ground_struct, excited_struct, + displacements) - # note: the returned structures won't include the 0% displacement, this is intended - # it can be included by specifying remove_zero=False - ground, excited = get_cc_structures(ground_struct, excited_struct, displacements) + files_to_copy = [f.strip() for f in input_files.split(',')] for i, struct in enumerate(ground): - working_dir = cc_dir / "ground" / str(i) + working_dir = cc_dir / 'ground' / str(i) os.makedirs(str(working_dir), exist_ok=True) - - # write structure and copy necessary input files - struct.to(filename=str(working_dir / "POSCAR"), fmt="poscar") - for f in ["KPOINTS", "POTCAR", "INCAR", "job_script.sh"]: - copyfile(str(ground_files / f), str(working_dir / f)) + struct.to(filename=str(working_dir / 'POSCAR'), fmt='poscar') + for f in files_to_copy: + src = ground_files / f + if src.exists(): + copyfile(str(src), str(working_dir / f)) + else: + click.echo(f'Warning: {src} not found, skipping.') for i, struct in enumerate(excited): - working_dir = cc_dir / "excited" / str(i) + working_dir = cc_dir / 'excited' / str(i) os.makedirs(str(working_dir), exist_ok=True) + struct.to(filename=str(working_dir / 'POSCAR'), fmt='poscar') + for f in files_to_copy: + src = excited_files / f + if src.exists(): + copyfile(str(src), str(working_dir / f)) + else: + click.echo(f'Warning: {src} not found, skipping.') - # write structure and copy necessary input files - struct.to(filename=str(working_dir / "POSCAR"), fmt="poscar") - for f in ["KPOINTS", "POTCAR", "INCAR", "job_script.sh"]: - copyfile(str(excited_files / f), str(working_dir / f)) - - return 0 + click.echo(f'Prepared {len(ground)} ground and {len(excited)} excited ' + f'structures in {cc_dir}') -@click.command() -@click.argument("cc_dir", type=click.Path()) -@click.argument("ground_files", type=click.Path()) -@click.argument("excited_files", type=click.Path()) -@click.option("--plot", "-p", is_flag=True, help="Plot the potential energy surfaces.") +@nonrad.command(name='pes') +@click.argument('cc_dir', type=click.Path(exists=True)) +@click.argument('ground_files', type=click.Path(exists=True)) +@click.argument('excited_files', type=click.Path(exists=True)) +@click.option( + '--energy-diff', '-e', + type=float, required=True, + help='Energy difference dE between ground and excited states (eV).', +) @click.option( - "--plot_name", - "-n", - default="pes.png", - help="Name of the plot file. Default is `pes.png`.", + '--plot', '-p', + is_flag=True, default=False, + help='Plot the potential energy surfaces.', ) -def pes(cc_dir, ground_files, excited_files, plot=False, plot_name="pes.png"): - """Extracting potential energy surfaces and relevant parameters from CCD calculations. - - Parmeters: \n - ----------- \n - cc_dir : str \n - Path to the directory containing the CCD calculations. It should contain the `ground` and `excited` directories. - Each of these directories should contain the directories corresponding to the displacements. \n - ground_files : str \n - Path to the directory containing the ground state calculation files. It should contain the `CONTCAR` and `vasprun.xml` files - for 0% displacement. \n - excited_files : str \n - Path to the directory containing the excited state calculation files. It should contain the `CONTCAR` and `vasprun.xml` files - for 0% displacement. \n - plot : bool \n - If True, the potential energy surfaces are plotted. Default is False. \n - plot_name : str \n - Name of the plot file. Default is `pes.png`. +@click.option( + '--plot-name', '-n', + type=str, default='pes.png', + show_default=True, + help='Filename for the PES plot.', +) +def pes(cc_dir, ground_files, excited_files, energy_diff, plot, plot_name): + """Extract potential energy surfaces from CCD calculations. + Reads vasprun.xml files from CC_DIR/ground/*/vasprun.xml and + CC_DIR/excited/*/vasprun.xml (including ``.gz`` variants), appends + the equilibrium vaspruns from GROUND_FILES and EXCITED_FILES, and + computes the PES and harmonic phonon frequencies. + + Parameters + ---------- + cc_dir : str + Root directory of the CCD calculation tree. + ground_files : str + Directory with the ground-state equilibrium calculation. + excited_files : str + Directory with the excited-state equilibrium calculation. + energy_diff : float + Energy offset *dE* between the two PES (eV). + plot : bool + Whether to save a PES plot. + plot_name : str + Output filename for the plot image. """ - ground_struct = Structure.from_file(Path(ground_files, "CONTCAR")) - excited_struct = Structure.from_file(Path(excited_files, "CONTCAR")) + from pymatgen.core import Structure + + from nonrad.ccd import ( + get_dQ, + get_omega_from_PES, + get_PES_from_vaspruns, + ) + + ground_struct = Structure.from_file(str(Path(ground_files) / 'CONTCAR')) + excited_struct = Structure.from_file(str(Path(excited_files) / 'CONTCAR')) - # calculate dQ - dQ = get_dQ(ground_struct, excited_struct) # amu^{1/2} Angstrom - print(f"dQ = {dQ:.2f} amu^{1/2} Angstrom") + dQ = get_dQ(ground_struct, excited_struct) + click.echo(f'dQ = {dQ:.2f} amu^{{1/2}} Angstrom') - # this prepares a list of all vasprun.xml's from the CCD calculations - ground_vaspruns = glob(Path(cc_dir, "ground", "*", "vasprun.xml")) - excited_vaspruns = glob(Path(cc_dir, "excited", "*", "vasprun.xml")) + # Collect vasprun.xml (and .gz) from the CCD sub-directories + ground_vaspruns = ( + glob(str(Path(cc_dir) / 'ground' / '*' / 'vasprun.xml')) + + glob(str(Path(cc_dir) / 'ground' / '*' / 'vasprun.xml.gz')) + ) + excited_vaspruns = ( + glob(str(Path(cc_dir) / 'excited' / '*' / 'vasprun.xml')) + + glob(str(Path(cc_dir) / 'excited' / '*' / 'vasprun.xml.gz')) + ) - # remember that the 0% displacement was removed before? we need to add that back in here - ground_vaspruns = ground_vaspruns + [Path(ground_files, "vasprun.xml")] - excited_vaspruns = excited_vaspruns + [Path(excited_files, "vasprun.xml")] + # Append equilibrium vaspruns + for gf in [str(Path(ground_files) / 'vasprun.xml'), + str(Path(ground_files) / 'vasprun.xml.gz')]: + if os.path.isfile(gf): + ground_vaspruns.append(gf) + break + + for ef in [str(Path(excited_files) / 'vasprun.xml'), + str(Path(excited_files) / 'vasprun.xml.gz')]: + if os.path.isfile(ef): + excited_vaspruns.append(ef) + break - # extract the potential energy surface Q_ground, E_ground = get_PES_from_vaspruns( - ground_struct, excited_struct, ground_vaspruns + ground_struct, excited_struct, ground_vaspruns, ) Q_excited, E_excited = get_PES_from_vaspruns( - ground_struct, excited_struct, excited_vaspruns + ground_struct, excited_struct, excited_vaspruns, ) + E_excited = energy_diff + E_excited - # the energy surfaces are referenced to the minimums, so we need to add dE (defined before) to E_excited - E_excited = dE + E_excited - - # calculate omega + # Always compute omega if plot: + import matplotlib + matplotlib.use('Agg') + import matplotlib.pyplot as plt + _, ax = plt.subplots(figsize=(5, 5)) ax.scatter(Q_ground, E_ground, s=10) ax.scatter(Q_excited, E_excited, s=10) - # by passing in the axis object, it also plots the fitted curve q = np.linspace(-1.0, 3.5, 100) ground_omega = get_omega_from_PES(Q_ground, E_ground, ax=ax, q=q) excited_omega = get_omega_from_PES(Q_excited, E_excited, ax=ax, q=q) - ax.set_xlabel("$Q$ [amu$^{1/2}$ $\AA$]") - ax.set_ylabel("$E$ [eV]") - plt.savefig(plot_name, dpi=300, format=plot_name.split(".")[-1]) + ax.set_xlabel(r'$Q$ [amu$^{1/2}$ $\AA$]') + ax.set_ylabel(r'$E$ [eV]') + plt.savefig(plot_name, dpi=300, format=plot_name.split('.')[-1]) + click.echo(f'Plot saved to {plot_name}') + else: + ground_omega = get_omega_from_PES(Q_ground, E_ground) + excited_omega = get_omega_from_PES(Q_excited, E_excited) + + click.echo(f'Ground state omega = {ground_omega:.2f} eV') + click.echo(f'Excited state omega = {excited_omega:.2f} eV') + + +@nonrad.command(name='dq') +@click.argument('ground_path', type=click.Path(exists=True)) +@click.argument('excited_path', type=click.Path(exists=True)) +def dq_cmd(ground_path, excited_path): + """Compute the configuration-coordinate displacement dQ. + + GROUND_PATH and EXCITED_PATH are structure files (e.g. POSCAR, + CONTCAR, CIF) — *not* directories. + + Parameters + ---------- + ground_path : str + Path to the ground-state structure file. + excited_path : str + Path to the excited-state structure file. + """ + from pymatgen.core import Structure + + from nonrad.ccd import get_dQ + + ground_struct = Structure.from_file(ground_path) + excited_struct = Structure.from_file(excited_path) + dQ = get_dQ(ground_struct, excited_struct) + click.echo(f'dQ = {dQ:.6f} amu^{{1/2}} Angstrom') + + +@nonrad.command(name='q-from-struct') +@click.argument('ground_path', type=click.Path(exists=True)) +@click.argument('excited_path', type=click.Path(exists=True)) +@click.argument('struct_path', type=click.Path(exists=True)) +@click.option( + '--tol', '-t', + type=float, default=1e-4, + show_default=True, + help='Distance tolerance for filtering coordinates.', +) +@click.option( + '--nround', + type=int, default=5, + show_default=True, + help='Number of decimal places for rounding Q.', +) +def q_from_struct(ground_path, excited_path, struct_path, tol, nround): + """Compute the Q value of an arbitrary structure. + + Given the ground and excited endpoint structures, determine the Q + coordinate of STRUCT_PATH assuming linear interpolation. + + Parameters + ---------- + ground_path : str + Path to the ground-state structure file. + excited_path : str + Path to the excited-state structure file. + struct_path : str + Path to the structure file to evaluate. + tol : float + Distance cutoff for discarding nearly-stationary sites. + nround : int + Decimal places used when rounding Q. + """ + from pymatgen.core import Structure + + from nonrad.ccd import get_Q_from_struct + + ground_struct = Structure.from_file(ground_path) + excited_struct = Structure.from_file(excited_path) + Q = get_Q_from_struct(ground_struct, excited_struct, struct_path, + tol=tol, nround=nround) + click.echo(f'Q = {Q}') + + +@nonrad.command(name='barrier') +@click.option('--dq', type=float, required=True, + help='Displacement dQ (amu^{1/2} Angstrom).') +@click.option('--de', type=float, required=True, + help='Energy offset dE (eV).') +@click.option('--wi', type=float, required=True, + help='Initial-state frequency (eV).') +@click.option('--wf', type=float, required=True, + help='Final-state frequency (eV).') +def barrier_cmd(dq, de, wi, wf): + """Compute the classical barrier height in the harmonic approximation. - print(f"Ground state omega = {ground_omega:.2f} eV") - print(f"Excited state omega = {excited_omega:.2f} eV") + Parameters + ---------- + dq : float + Configuration-coordinate displacement (amu^{1/2} Angstrom). + de : float + Energy offset between oscillators (eV). + wi : float + Initial harmonic-oscillator frequency (eV). + wf : float + Final harmonic-oscillator frequency (eV). + """ + from nonrad.ccd import get_barrier_harmonic + + result = get_barrier_harmonic(dq, de, wi, wf) + if result is not None: + click.echo(f'Barrier height = {result:.6f} eV') + else: + click.echo('No crossing point found.') + + +# --------------------------------------------------------------------------- +# Core capture coefficient command +# --------------------------------------------------------------------------- + +@nonrad.command(name='capture') +@click.option('--dq', type=float, required=True, + help='Displacement dQ (amu^{1/2} Angstrom).') +@click.option('--de', type=float, required=True, + help='Energy offset dE (eV).') +@click.option('--wi', type=float, required=True, + help='Initial-state frequency (eV).') +@click.option('--wf', type=float, required=True, + help='Final-state frequency (eV).') +@click.option('--wif', type=float, required=True, + help='Electron-phonon coupling Wif (eV amu^{-1/2} Angstrom^{-1}).') +@click.option('--volume', type=float, required=True, + help='Supercell volume (Angstrom^3).') +@click.option('--g', type=int, default=1, show_default=True, + help='Degeneracy factor.') +@click.option('--temperature', '-T', type=float, default=300., + show_default=True, + help='Temperature (K) for a single-point evaluation.') +@click.option('--temperature-range', nargs=3, type=float, default=None, + help='Temperature range as (min, max, count) for scanning.') +@click.option('--sigma', type=str, default='pchip', show_default=True, + help='Smearing parameter (float) or interpolation method ' + '("pchip" or "cubic").') +@click.option('--occ-tol', type=float, default=1e-5, show_default=True, + help='Occupation tolerance for Bose weights.') +@click.option('--overlap-method', + type=click.Choice(['Analytic', 'Integral', 'HermiteGauss'], + case_sensitive=False), + default='HermiteGauss', show_default=True, + help='Method for evaluating vibrational overlaps.') +def capture_cmd(dq, de, wi, wf, wif, volume, g, temperature, + temperature_range, sigma, occ_tol, overlap_method): + """Compute the nonradiative capture coefficient. + + Evaluates the capture coefficient following Alkauskas *et al.*, + Phys. Rev. B 90, 075202 (2014). The result is **unscaled** and has + units of cm^3 s^{-1}. + + Parameters + ---------- + dq : float + Configuration-coordinate displacement (amu^{1/2} Angstrom). + de : float + Energy offset between oscillators (eV). + wi : float + Initial-state frequency (eV). + wf : float + Final-state frequency (eV). + wif : float + Electron-phonon coupling matrix element + (eV amu^{-1/2} Angstrom^{-1}). + volume : float + Supercell volume (Angstrom^3). + g : int + Degeneracy factor of the final state. + temperature : float + Single temperature (K). + temperature_range : tuple of float or None + ``(min, max, count)`` for a temperature scan. + sigma : str + Smearing parameter or interpolation keyword. + occ_tol : float + Bose-weight occupation tolerance. + overlap_method : str + Overlap evaluation method. + """ + from nonrad.nonrad import get_C + + # Parse sigma: interpret as float when possible + try: + sigma_val = float(sigma) + except ValueError: + sigma_val = sigma + + if temperature_range is not None: + temps = np.linspace(temperature_range[0], temperature_range[1], + int(temperature_range[2])) + else: + temps = temperature + + result = get_C(dq, de, wi, wf, wif, volume, + g=g, T=temps, sigma=sigma_val, + occ_tol=occ_tol, overlap_method=overlap_method) + + if isinstance(temps, np.ndarray): + click.echo(f'{"T (K)":>12s} {"C (cm^3/s)":>14s}') + click.echo('-' * 28) + for t, c in zip(temps, np.atleast_1d(result)): + click.echo(f'{t:12.2f} {c:14.6e}') + else: + click.echo(f'C = {result:.6e} cm^3/s (T = {temps} K)') + + +# --------------------------------------------------------------------------- +# Scaling commands +# --------------------------------------------------------------------------- + +@nonrad.command(name='sommerfeld') +@click.option('--temperature', '-T', type=float, default=300., + show_default=True, help='Temperature (K).') +@click.option('--temperature-range', nargs=3, type=float, default=None, + help='Temperature range as (min, max, count).') +@click.option('--z', type=int, required=True, + help='Charge ratio Z = Q/q (negative = attractive).') +@click.option('--m-eff', type=float, required=True, + help='Effective mass in units of electron mass.') +@click.option('--eps0', type=float, required=True, + help='Static dielectric constant.') +@click.option('--dim', + type=click.Choice(['1', '2', '3'], case_sensitive=False), + default='3', show_default=True, + help='Dimensionality of the system.') +@click.option('--method', + type=click.Choice(['Integrate', 'Analytic'], + case_sensitive=False), + default='Integrate', show_default=True, + help='Evaluation method.') +def sommerfeld_cmd(temperature, temperature_range, z, m_eff, eps0, dim, + method): + """Compute the Sommerfeld enhancement factor. + + Evaluates the temperature-dependent Sommerfeld parameter following + Pässler *et al.*, phys. stat. sol. (b) 78, 625 (1976). + + Parameters + ---------- + temperature : float + Temperature (K) for a single-point evaluation. + temperature_range : tuple of float or None + ``(min, max, count)`` for a temperature scan. + z : int + Charge ratio *Z = Q / q*. + m_eff : float + Carrier effective mass in units of *m_e*. + eps0 : float + Static dielectric constant. + dim : str + System dimensionality (``'1'``, ``'2'``, or ``'3'``). + method : str + ``'Integrate'`` or ``'Analytic'``. + """ + from nonrad.scaling import sommerfeld_parameter + + dim_int = int(dim) + + if temperature_range is not None: + temps = np.linspace(temperature_range[0], temperature_range[1], + int(temperature_range[2])) + else: + temps = temperature + + result = sommerfeld_parameter(temps, z, m_eff, eps0, + dim=dim_int, method=method) + + if isinstance(temps, np.ndarray): + click.echo(f'{"T (K)":>12s} {"Sommerfeld":>14s}') + click.echo('-' * 28) + for t, s in zip(temps, np.atleast_1d(result)): + click.echo(f'{t:12.2f} {s:14.6e}') + else: + click.echo(f'Sommerfeld factor = {result:.6e} (T = {temps} K)') + + +@nonrad.command(name='thermal-velocity') +@click.option('--temperature', '-T', type=float, default=300., + show_default=True, help='Temperature (K).') +@click.option('--temperature-range', nargs=3, type=float, default=None, + help='Temperature range as (min, max, count).') +@click.option('--m-eff', type=float, required=True, + help='Effective mass in units of electron mass.') +def thermal_velocity_cmd(temperature, temperature_range, m_eff): + """Compute the thermal velocity of a carrier. + + Parameters + ---------- + temperature : float + Temperature (K) for a single-point evaluation. + temperature_range : tuple of float or None + ``(min, max, count)`` for a temperature scan. + m_eff : float + Carrier effective mass in units of *m_e*. + """ + from nonrad.scaling import thermal_velocity + + if temperature_range is not None: + temps = np.linspace(temperature_range[0], temperature_range[1], + int(temperature_range[2])) + else: + temps = temperature + + result = thermal_velocity(temps, m_eff) + + if isinstance(temps, np.ndarray): + click.echo(f'{"T (K)":>12s} {"v_th (cm/s)":>14s}') + click.echo('-' * 28) + for t, v in zip(temps, np.atleast_1d(result)): + click.echo(f'{t:12.2f} {v:14.6e}') + else: + click.echo(f'Thermal velocity = {result:.6e} cm/s (T = {temps} K)') + + +@nonrad.command(name='charged-supercell') +@click.argument('wavecar_path', type=click.Path(exists=True)) +@click.option('--bulk-index', type=int, required=True, + help='Index of the bulk wavefunction (1-based).') +@click.option('--def-index', type=int, default=-1, show_default=True, + help='Index of the defect wavefunction (1-based). ' + 'Provide this or --def-coord.') +@click.option('--def-coord', nargs=3, type=float, default=None, + help='Cartesian coordinates (x y z) of the defect position.') +@click.option('--cutoff', type=float, default=0.02, show_default=True, + help='Slope cutoff for plateau detection.') +@click.option('--limit', type=float, default=5., show_default=True, + help='Upper radial limit (Angstrom) for the fitting window.') +@click.option('--spin', type=int, default=0, show_default=True, + help='Spin channel (0 = up, 1 = down).') +@click.option('--kpoint', type=int, default=1, show_default=True, + help='k-point index (1-based).') +def charged_supercell_cmd(wavecar_path, bulk_index, def_index, def_coord, + cutoff, limit, spin, kpoint): + """Estimate the charged-supercell scaling factor. + + Compares the radial distribution of the bulk wavefunction around the + defect position to a uniform distribution and returns a scaling + factor for the capture coefficient. + + Parameters + ---------- + wavecar_path : str + Path to the WAVECAR file. + bulk_index : int + Band index of the bulk wavefunction (1-based). + def_index : int + Band index of the defect wavefunction (1-based), or -1 if + ``def_coord`` is provided instead. + def_coord : tuple of float or None + Cartesian coordinates of the defect centre. + cutoff : float + Gradient cutoff for plateau identification. + limit : float + Upper radius limit for curve fitting (Angstrom). + spin : int + Spin channel (0 = up, 1 = down). + kpoint : int + k-point index (1-based). + """ + from nonrad.scaling import charged_supercell_scaling_VASP + + if def_index == -1 and def_coord is None: + raise click.UsageError( + 'Either --def-index (not -1) or --def-coord must be specified.' + ) + + coord = np.array(def_coord) if def_coord is not None else None + + scaling = charged_supercell_scaling_VASP( + wavecar_path, bulk_index, + def_index=def_index, def_coord=coord, + cutoff=cutoff, limit=limit, + spin=spin, kpoint=kpoint, + ) + click.echo(f'Charged-supercell scaling factor = {scaling:.6f}') + + +# --------------------------------------------------------------------------- +# Electron-phonon sub-group +# --------------------------------------------------------------------------- + +@nonrad.group('elphon') +def elphon(): + """Electron-phonon coupling utilities. + + Sub-commands for computing the electron-phonon coupling matrix + element Wif using WAVECARs, WSWQ files, or UNK files. + """ + + +@elphon.command(name='wavecars') +@click.argument('cc_dir', type=click.Path(exists=True)) +@click.argument('ground_path', type=click.Path(exists=True)) +@click.argument('excited_path', type=click.Path(exists=True)) +@click.argument('init_wavecar', type=click.Path(exists=True)) +@click.option('--def-index', type=int, required=True, + help='Defect wavefunction index (1-based).') +@click.option('--bulk-index', '-b', type=int, required=True, multiple=True, + help='Bulk wavefunction index (1-based, repeatable).') +@click.option('--spin', type=int, default=0, show_default=True, + help='Spin channel (0 = up, 1 = down).') +@click.option('--kpoint', type=int, default=1, show_default=True, + help='k-point index (1-based).') +def elphon_wavecars(cc_dir, ground_path, excited_path, init_wavecar, + def_index, bulk_index, spin, kpoint): + """Compute Wif from WAVECAR files along the CCD path. + + Automatically discovers sub-directories inside CC_DIR that contain + both a ``vasprun.xml`` (or ``.gz``) and a ``WAVECAR`` file. + + Parameters + ---------- + cc_dir : str + Directory tree with CCD calculations. + ground_path : str + Directory containing the ground-state CONTCAR. + excited_path : str + Directory containing the excited-state CONTCAR. + init_wavecar : str + Path to the initial (reference) WAVECAR. + def_index : int + Band index of the defect wavefunction. + bulk_index : tuple of int + Band indices of the bulk wavefunctions. + spin : int + Spin channel. + kpoint : int + k-point index. + """ + from pymatgen.core import Structure + from pymatgen.io.vasp.outputs import Vasprun + + from nonrad.ccd import get_Q_from_struct + from nonrad.elphon import get_Wif_from_wavecars + + ground_struct = Structure.from_file(str(Path(ground_path) / 'CONTCAR')) + excited_struct = Structure.from_file(str(Path(excited_path) / 'CONTCAR')) + + # Auto-discover sub-directories with vasprun + WAVECAR + wavecars = [] + for subdir in sorted(Path(cc_dir).iterdir()): + if not subdir.is_dir(): + continue + wavecar_file = subdir / 'WAVECAR' + if not wavecar_file.exists(): + continue + vr_path = None + for vr_name in ('vasprun.xml', 'vasprun.xml.gz'): + candidate = subdir / vr_name + if candidate.exists(): + vr_path = candidate + break + if vr_path is None: + continue + vr = Vasprun(str(vr_path), parse_dos=False, parse_eigen=False) + q = get_Q_from_struct(ground_struct, excited_struct, + vr.structures[-1]) + wavecars.append((q, str(wavecar_file))) + + results = get_Wif_from_wavecars( + wavecars, init_wavecar, def_index, list(bulk_index), + spin=spin, kpoint=kpoint, + ) + + click.echo(f'{"bulk_index":>12s} {"Wif (eV amu^{-1/2} A^{-1})":>30s}') + click.echo('-' * 44) + for bi, wif_val in results: + click.echo(f'{bi:12d} {wif_val:30.6e}') + + +@elphon.command(name='wswq') +@click.argument('cc_dir', type=click.Path(exists=True)) +@click.argument('ground_path', type=click.Path(exists=True)) +@click.argument('excited_path', type=click.Path(exists=True)) +@click.argument('init_vasprun', type=click.Path(exists=True)) +@click.option('--def-index', type=int, required=True, + help='Defect wavefunction index (1-based).') +@click.option('--bulk-index', '-b', type=int, required=True, multiple=True, + help='Bulk wavefunction index (1-based, repeatable).') +@click.option('--spin', type=int, default=0, show_default=True, + help='Spin channel (0 = up, 1 = down).') +@click.option('--kpoint', type=int, default=1, show_default=True, + help='k-point index (1-based).') +def elphon_wswq(cc_dir, ground_path, excited_path, init_vasprun, + def_index, bulk_index, spin, kpoint): + """Compute Wif from WSWQ files along the CCD path. + + Automatically discovers sub-directories inside CC_DIR that contain + both a ``vasprun.xml`` (or ``.gz``) and a ``WSWQ`` (or ``.gz``) + file. + + Parameters + ---------- + cc_dir : str + Directory tree with CCD calculations. + ground_path : str + Directory containing the ground-state CONTCAR. + excited_path : str + Directory containing the excited-state CONTCAR. + init_vasprun : str + Path to the initial vasprun.xml for eigenvalue extraction. + def_index : int + Band index of the defect wavefunction. + bulk_index : tuple of int + Band indices of the bulk wavefunctions. + spin : int + Spin channel. + kpoint : int + k-point index. + """ + from pymatgen.core import Structure + from pymatgen.io.vasp.outputs import Vasprun + + from nonrad.ccd import get_Q_from_struct + from nonrad.elphon import get_Wif_from_WSWQ + + ground_struct = Structure.from_file(str(Path(ground_path) / 'CONTCAR')) + excited_struct = Structure.from_file(str(Path(excited_path) / 'CONTCAR')) + + # Auto-discover sub-directories with vasprun + WSWQ + wswqs = [] + for subdir in sorted(Path(cc_dir).iterdir()): + if not subdir.is_dir(): + continue + wswq_file = None + for wswq_name in ('WSWQ', 'WSWQ.gz'): + candidate = subdir / wswq_name + if candidate.exists(): + wswq_file = candidate + break + if wswq_file is None: + continue + vr_path = None + for vr_name in ('vasprun.xml', 'vasprun.xml.gz'): + candidate = subdir / vr_name + if candidate.exists(): + vr_path = candidate + break + if vr_path is None: + continue + vr = Vasprun(str(vr_path), parse_dos=False, parse_eigen=False) + q = get_Q_from_struct(ground_struct, excited_struct, + vr.structures[-1]) + wswqs.append((q, str(wswq_file))) + + results = get_Wif_from_WSWQ( + wswqs, init_vasprun, def_index, list(bulk_index), + spin=spin, kpoint=kpoint, + ) + + click.echo(f'{"bulk_index":>12s} {"Wif (eV amu^{-1/2} A^{-1})":>30s}') + click.echo('-' * 44) + for bi, wif_val in results: + click.echo(f'{bi:12d} {wif_val:30.6e}') + + +@elphon.command(name='unk') +@click.argument('init_unk', type=click.Path(exists=True)) +@click.option('--def-index', type=int, required=True, + help='Defect wavefunction index (1-based).') +@click.option('--bulk-index', '-b', type=int, required=True, multiple=True, + help='Bulk wavefunction index (1-based, repeatable).') +@click.option('--eigs', type=str, required=True, + help='Comma-separated eigenvalues (eV).') +@click.option('--unk-pair', '-u', type=(float, str), required=True, + multiple=True, + help='(Q, unk_path) pair — repeatable.') +def elphon_unk(init_unk, def_index, bulk_index, eigs, unk_pair): + """Compute Wif from UNK files. + + Parameters + ---------- + init_unk : str + Path to the initial (reference) UNK file. + def_index : int + Band index of the defect wavefunction. + bulk_index : tuple of int + Band indices of the bulk wavefunctions. + eigs : str + Comma-separated string of eigenvalues (eV). + unk_pair : tuple of (float, str) + One or more ``(Q, path)`` pairs specifying displaced UNK files. + """ + from nonrad.elphon import get_Wif_from_UNK + + eigs_array = np.array([float(e) for e in eigs.split(',')]) + unk_list = [(q, path) for q, path in unk_pair] + + results = get_Wif_from_UNK( + unk_list, init_unk, def_index, list(bulk_index), + eigs=eigs_array, + ) - return 0 + click.echo(f'{"bulk_index":>12s} {"Wif (eV amu^{-1/2} A^{-1})":>30s}') + click.echo('-' * 44) + for bi, wif_val in results: + click.echo(f'{bi:12d} {wif_val:30.6e}') -# Add subcommands to the main group -nonrad.add_command(prep_ccd) -nonrad.add_command(pes) +# --------------------------------------------------------------------------- +# Entry point +# --------------------------------------------------------------------------- -if __name__ == "__main__": +if __name__ == '__main__': nonrad() diff --git a/nonrad/tests/test_cli.py b/nonrad/tests/test_cli.py new file mode 100644 index 0000000..bfd09da --- /dev/null +++ b/nonrad/tests/test_cli.py @@ -0,0 +1,606 @@ +# pylint: disable=C0114,C0115,C0116 + +import os +import unittest +from pathlib import Path +from unittest.mock import MagicMock, patch + +import numpy as np +from click.testing import CliRunner + +from nonrad.cli import nonrad as nonrad_cli +from nonrad.tests import TEST_FILES + + +# --------------------------------------------------------------------------- +# Tests for the top-level `nonrad` group +# --------------------------------------------------------------------------- +class TestNonradGroup(unittest.TestCase): + def setUp(self): + self.runner = CliRunner() + + def test_help(self): + result = self.runner.invoke(nonrad_cli, ['--help']) + self.assertEqual(result.exit_code, 0, result.output) + self.assertIn('Nonradiative recombination', result.output) + + def test_no_args(self): + result = self.runner.invoke(nonrad_cli, []) + self.assertEqual(result.exit_code, 0, result.output) + + +# --------------------------------------------------------------------------- +# Tests for `prep-ccd` +# --------------------------------------------------------------------------- +class TestPrepCcd(unittest.TestCase): + def setUp(self): + self.runner = CliRunner() + + def test_help(self): + result = self.runner.invoke(nonrad_cli, ['prep-ccd', '--help']) + self.assertEqual(result.exit_code, 0, result.output) + self.assertIn('ground_path', result.output.lower()) + + @patch('nonrad.ccd.get_cc_structures') + @patch('pymatgen.core.structure.Structure.from_file') + def test_basic(self, mock_from_file, mock_get_cc): + mock_struct = MagicMock() + mock_from_file.return_value = mock_struct + + # get_cc_structures returns (list_of_structs, list_of_structs) + struct_a = MagicMock() + struct_b = MagicMock() + mock_get_cc.return_value = ([struct_a], [struct_b]) + + with self.runner.isolated_filesystem(): + # Create fake ground/excited directories with required files + os.makedirs('ground') + Path('ground/CONTCAR').touch() + for f in ['KPOINTS', 'POTCAR', 'INCAR', 'job_script.sh']: + Path(f'ground/{f}').write_text('dummy') + + os.makedirs('excited') + Path('excited/CONTCAR').touch() + for f in ['KPOINTS', 'POTCAR', 'INCAR', 'job_script.sh']: + Path(f'excited/{f}').write_text('dummy') + + result = self.runner.invoke( + nonrad_cli, + ['prep-ccd', 'ground', 'excited', 'ccd_output'] + ) + self.assertEqual(result.exit_code, 0, + f'exit_code={result.exit_code}\n{result.output}') + # Verify output directories were created + self.assertTrue(os.path.isdir('ccd_output')) + self.assertTrue(os.path.isdir('ccd_output/ground')) + self.assertTrue(os.path.isdir('ccd_output/excited')) + # Verify struct.to() was called for each displacement struct + struct_a.to.assert_called() + struct_b.to.assert_called() + + @patch('nonrad.ccd.get_cc_structures') + @patch('pymatgen.core.structure.Structure.from_file') + def test_custom_displacements(self, mock_from_file, mock_get_cc): + mock_struct = MagicMock() + mock_from_file.return_value = mock_struct + mock_get_cc.return_value = ([], []) + + with self.runner.isolated_filesystem(): + os.makedirs('ground') + Path('ground/CONTCAR').touch() + for f in ['KPOINTS', 'POTCAR', 'INCAR', 'job_script.sh']: + Path(f'ground/{f}').write_text('dummy') + + os.makedirs('excited') + Path('excited/CONTCAR').touch() + for f in ['KPOINTS', 'POTCAR', 'INCAR', 'job_script.sh']: + Path(f'excited/{f}').write_text('dummy') + + result = self.runner.invoke( + nonrad_cli, + ['prep-ccd', 'ground', 'excited', 'ccd_output', + '-d', '-0.3', '0.3', '5'] + ) + self.assertEqual(result.exit_code, 0, + f'exit_code={result.exit_code}\n{result.output}') + # Verify get_cc_structures was called + mock_get_cc.assert_called_once() + # Check displacements were passed correctly (via linspace) + call_args = mock_get_cc.call_args + displacements = call_args[0][2] + self.assertEqual(len(displacements), 5) + + @patch('nonrad.ccd.get_cc_structures') + @patch('pymatgen.core.structure.Structure.from_file') + def test_nonexistent_ground_path(self, mock_from_file, mock_get_cc): + with self.runner.isolated_filesystem(): + result = self.runner.invoke( + nonrad_cli, + ['prep-ccd', 'nonexistent', 'also_nonexistent', 'out'] + ) + self.assertNotEqual(result.exit_code, 0) + + +# --------------------------------------------------------------------------- +# Tests for `pes` +# --------------------------------------------------------------------------- +class TestPes(unittest.TestCase): + def setUp(self): + self.runner = CliRunner() + + def test_help(self): + result = self.runner.invoke(nonrad_cli, ['pes', '--help']) + self.assertEqual(result.exit_code, 0, result.output) + self.assertIn('cc_dir', result.output.lower()) + + @patch('nonrad.ccd.get_omega_from_PES', return_value=0.05) + @patch('nonrad.ccd.get_PES_from_vaspruns', + return_value=(np.array([0.0, 1.0]), np.array([0.0, 0.5]))) + @patch('nonrad.ccd.get_dQ', return_value=1.5) + @patch('nonrad.cli.glob', return_value=[]) + @patch('pymatgen.core.structure.Structure.from_file') + def test_basic(self, mock_from_file, mock_glob, + mock_get_dq, mock_get_pes, mock_get_omega): + mock_struct = MagicMock() + mock_from_file.return_value = mock_struct + + with self.runner.isolated_filesystem(): + # Create fake dirs / files + os.makedirs('cc/ground/0', exist_ok=True) + os.makedirs('cc/excited/0', exist_ok=True) + os.makedirs('ground_files') + Path('ground_files/CONTCAR').touch() + Path('ground_files/vasprun.xml').touch() + os.makedirs('excited_files') + Path('excited_files/CONTCAR').touch() + Path('excited_files/vasprun.xml').touch() + + result = self.runner.invoke( + nonrad_cli, + ['pes', 'cc', 'ground_files', 'excited_files', '--energy-diff', '1.0'] + ) + # Verify exit code is 0 and mock gets called + self.assertEqual(result.exit_code, 0, result.output) + mock_get_dq.assert_called_once() + + +# --------------------------------------------------------------------------- +# Tests for `dq` +# --------------------------------------------------------------------------- +class TestDq(unittest.TestCase): + def setUp(self): + self.runner = CliRunner() + + def test_help(self): + result = self.runner.invoke(nonrad_cli, ['dq', '--help']) + self.assertEqual(result.exit_code, 0, result.output) + + def test_basic_real_files(self): + ground = str(TEST_FILES / 'POSCAR.C0.gz') + excited = str(TEST_FILES / 'POSCAR.C-.gz') + result = self.runner.invoke( + nonrad_cli, ['dq', ground, excited] + ) + self.assertEqual(result.exit_code, 0, + f'exit_code={result.exit_code}\n{result.output}') + # dQ ≈ 1.6859 + self.assertIn('1.6858', result.output) + + def test_nonexistent_file(self): + result = self.runner.invoke( + nonrad_cli, ['dq', 'nonexistent_file', 'also_nonexistent'] + ) + self.assertNotEqual(result.exit_code, 0) + + +# --------------------------------------------------------------------------- +# Tests for `q-from-struct` +# --------------------------------------------------------------------------- +class TestQFromStruct(unittest.TestCase): + def setUp(self): + self.runner = CliRunner() + + def test_help(self): + result = self.runner.invoke(nonrad_cli, ['q-from-struct', '--help']) + self.assertEqual(result.exit_code, 0, result.output) + + def test_basic_real_files(self): + ground = str(TEST_FILES / 'POSCAR.C0.gz') + excited = str(TEST_FILES / 'POSCAR.C-.gz') + struct = str(TEST_FILES / 'POSCAR.C0.gz') + result = self.runner.invoke( + nonrad_cli, ['q-from-struct', ground, excited, struct] + ) + self.assertEqual(result.exit_code, 0, + f'exit_code={result.exit_code}\n{result.output}') + # Q for ground == ground should be ≈ 0 + self.assertIn('0.0', result.output) + + def test_nonexistent_file(self): + result = self.runner.invoke( + nonrad_cli, + ['q-from-struct', 'nonexistent', 'nonexistent', 'nonexistent'] + ) + self.assertNotEqual(result.exit_code, 0) + + +# --------------------------------------------------------------------------- +# Tests for `barrier` +# --------------------------------------------------------------------------- +class TestBarrier(unittest.TestCase): + def setUp(self): + self.runner = CliRunner() + + def test_help(self): + result = self.runner.invoke(nonrad_cli, ['barrier', '--help']) + self.assertEqual(result.exit_code, 0, result.output) + + def test_crossing(self): + result = self.runner.invoke( + nonrad_cli, + ['barrier', + '--dq', '1.0', '--de', '0.0', '--wi', '0.05', '--wf', '0.05'] + ) + self.assertEqual(result.exit_code, 0, + f'exit_code={result.exit_code}\n{result.output}') + # Should contain barrier height info + output_lower = result.output.lower() + self.assertTrue( + 'barrier' in output_lower or 'height' in output_lower + or 'ev' in output_lower, + f'Expected barrier output, got: {result.output}' + ) + + @patch('nonrad.ccd.get_barrier_harmonic', return_value=None) + def test_no_crossing(self, mock_barrier): + result = self.runner.invoke( + nonrad_cli, + ['barrier', + '--dq', '1.0', '--de', '5.0', '--wi', '0.01', '--wf', '0.01'] + ) + # The CLI should handle None return (no crossing) gracefully + self.assertEqual(result.exit_code, 0, + f'exit_code={result.exit_code}\n{result.output}') + output_lower = result.output.lower() + self.assertTrue( + 'no' in output_lower or 'none' in output_lower + or 'crossing' in output_lower, + f'Expected no-crossing message, got: {result.output}' + ) + + def test_missing_options(self): + result = self.runner.invoke( + nonrad_cli, + ['barrier', '--dq', '1.0'] # missing --de, --wi, --wf + ) + self.assertNotEqual(result.exit_code, 0) + + +# --------------------------------------------------------------------------- +# Tests for `capture` +# --------------------------------------------------------------------------- +class TestCapture(unittest.TestCase): + def setUp(self): + self.runner = CliRunner() + + def test_help(self): + result = self.runner.invoke(nonrad_cli, ['capture', '--help']) + self.assertEqual(result.exit_code, 0, result.output) + + @patch('nonrad.nonrad.get_C', return_value=1.23e-8) + def test_basic_single_temperature(self, mock_get_c): + result = self.runner.invoke( + nonrad_cli, + ['capture', + '--dq', '1.0', '--de', '1.0', + '--wi', '0.03', '--wf', '0.03', + '--wif', '0.01', '--volume', '1000', + '-T', '300'] + ) + self.assertEqual(result.exit_code, 0, + f'exit_code={result.exit_code}\n{result.output}') + + @patch('nonrad.nonrad.get_C', + return_value=np.array([1e-8, 2e-8, 3e-8])) + def test_temperature_range(self, mock_get_c): + result = self.runner.invoke( + nonrad_cli, + ['capture', + '--dq', '1.0', '--de', '1.0', + '--wi', '0.03', '--wf', '0.03', + '--wif', '0.01', '--volume', '1000', + '--temperature-range', '100', '500', '3'] + ) + self.assertEqual(result.exit_code, 0, + f'exit_code={result.exit_code}\n{result.output}') + + def test_missing_required_options(self): + result = self.runner.invoke( + nonrad_cli, + ['capture', '--dq', '1.0'] # missing many required options + ) + self.assertNotEqual(result.exit_code, 0) + + +# --------------------------------------------------------------------------- +# Tests for `sommerfeld` +# --------------------------------------------------------------------------- +class TestSommerfeld(unittest.TestCase): + def setUp(self): + self.runner = CliRunner() + + def test_help(self): + result = self.runner.invoke(nonrad_cli, ['sommerfeld', '--help']) + self.assertEqual(result.exit_code, 0, result.output) + + def test_neutral_z0(self): + """Z=0 is the neutral special case: sommerfeld_parameter returns 1.0""" + result = self.runner.invoke( + nonrad_cli, + ['sommerfeld', + '--z', '0', '--m-eff', '1.0', '--eps0', '10.0', + '-T', '300'] + ) + self.assertEqual(result.exit_code, 0, + f'exit_code={result.exit_code}\n{result.output}') + # Z=0 → sommerfeld = 1.0 + self.assertIn('1.0', result.output) + + def test_temperature_range(self): + result = self.runner.invoke( + nonrad_cli, + ['sommerfeld', + '--z', '0', '--m-eff', '1.0', '--eps0', '10.0', + '--temperature-range', '100', '500', '3'] + ) + self.assertEqual(result.exit_code, 0, + f'exit_code={result.exit_code}\n{result.output}') + + def test_missing_required_options(self): + result = self.runner.invoke( + nonrad_cli, + ['sommerfeld', '-T', '300'] # missing --z, --m-eff, --eps0 + ) + self.assertNotEqual(result.exit_code, 0) + + +# --------------------------------------------------------------------------- +# Tests for `thermal-velocity` +# --------------------------------------------------------------------------- +class TestThermalVelocity(unittest.TestCase): + def setUp(self): + self.runner = CliRunner() + + def test_help(self): + result = self.runner.invoke( + nonrad_cli, ['thermal-velocity', '--help'] + ) + self.assertEqual(result.exit_code, 0, result.output) + + def test_basic(self): + result = self.runner.invoke( + nonrad_cli, + ['thermal-velocity', '--m-eff', '1.0', '-T', '300'] + ) + self.assertEqual(result.exit_code, 0, + f'exit_code={result.exit_code}\n{result.output}') + self.assertIn('cm/s', result.output) + + def test_temperature_range(self): + result = self.runner.invoke( + nonrad_cli, + ['thermal-velocity', '--m-eff', '1.0', + '--temperature-range', '100', '500', '3'] + ) + self.assertEqual(result.exit_code, 0, + f'exit_code={result.exit_code}\n{result.output}') + + def test_missing_meff(self): + result = self.runner.invoke( + nonrad_cli, + ['thermal-velocity', '-T', '300'] # missing --m-eff + ) + self.assertNotEqual(result.exit_code, 0) + + +# --------------------------------------------------------------------------- +# Tests for `charged-supercell` +# --------------------------------------------------------------------------- +class TestChargedSupercell(unittest.TestCase): + def setUp(self): + self.runner = CliRunner() + + def test_help(self): + result = self.runner.invoke( + nonrad_cli, ['charged-supercell', '--help'] + ) + self.assertEqual(result.exit_code, 0, result.output) + + @patch('nonrad.scaling.charged_supercell_scaling_VASP', return_value=0.95) + def test_basic_with_def_index(self, mock_scaling): + with self.runner.isolated_filesystem(): + Path('WAVECAR').write_bytes(b'\x00' * 64) + result = self.runner.invoke( + nonrad_cli, + ['charged-supercell', 'WAVECAR', + '--bulk-index', '189', '--def-index', '192'] + ) + self.assertEqual(result.exit_code, 0, + f'exit_code={result.exit_code}\n{result.output}') + + @patch('nonrad.scaling.charged_supercell_scaling_VASP', return_value=0.95) + def test_basic_with_def_coord(self, mock_scaling): + with self.runner.isolated_filesystem(): + Path('WAVECAR').write_bytes(b'\x00' * 64) + result = self.runner.invoke( + nonrad_cli, + ['charged-supercell', 'WAVECAR', + '--bulk-index', '189', + '--def-coord', '0.5', '0.5', '0.5'] + ) + self.assertEqual(result.exit_code, 0, + f'exit_code={result.exit_code}\n{result.output}') + + def test_missing_def_index_and_def_coord(self): + """Both --def-index and --def-coord are omitted → error.""" + with self.runner.isolated_filesystem(): + Path('WAVECAR').write_bytes(b'\x00' * 64) + result = self.runner.invoke( + nonrad_cli, + ['charged-supercell', 'WAVECAR', '--bulk-index', '189'] + ) + # The CLI or underlying function should flag this as an error + self.assertNotEqual(result.exit_code, 0, + f'Expected error, got: {result.output}') + + def test_nonexistent_wavecar(self): + result = self.runner.invoke( + nonrad_cli, + ['charged-supercell', 'nonexistent_WAVECAR', + '--bulk-index', '189', '--def-index', '192'] + ) + self.assertNotEqual(result.exit_code, 0) + + +# --------------------------------------------------------------------------- +# Tests for `elphon` group +# --------------------------------------------------------------------------- +class TestElphonGroup(unittest.TestCase): + def setUp(self): + self.runner = CliRunner() + + def test_help(self): + result = self.runner.invoke(nonrad_cli, ['elphon', '--help']) + self.assertEqual(result.exit_code, 0, result.output) + self.assertIn('elphon', result.output.lower()) + + +# --------------------------------------------------------------------------- +# Tests for `elphon wavecars` +# --------------------------------------------------------------------------- +class TestElphonWavecars(unittest.TestCase): + def setUp(self): + self.runner = CliRunner() + + def test_help(self): + result = self.runner.invoke( + nonrad_cli, ['elphon', 'wavecars', '--help'] + ) + self.assertEqual(result.exit_code, 0, result.output) + + @patch('nonrad.elphon.get_Wif_from_wavecars', + return_value=[(189, 0.087)]) + @patch('nonrad.ccd.get_Q_from_struct', return_value=0.5) + @patch('pymatgen.io.vasp.outputs.Vasprun') + @patch('pymatgen.core.structure.Structure.from_file') + def test_basic(self, mock_from_file, mock_vasprun, + mock_get_q, mock_get_wif): + mock_struct = MagicMock() + mock_from_file.return_value = mock_struct + + mock_vr = MagicMock() + mock_vr.structures = [mock_struct] + mock_vasprun.return_value = mock_vr + + with self.runner.isolated_filesystem(): + # Create fake directory structure + os.makedirs('cc/ground/0', exist_ok=True) + os.makedirs('cc/excited/0', exist_ok=True) + Path('cc/ground/0/WAVECAR').touch() + Path('cc/ground/0/vasprun.xml').touch() + os.makedirs('ground_files') + Path('ground_files/CONTCAR').touch() + os.makedirs('excited_files') + Path('excited_files/CONTCAR').touch() + Path('WAVECAR_init').touch() + + result = self.runner.invoke( + nonrad_cli, + ['elphon', 'wavecars', + 'cc', 'ground_files', 'excited_files', 'WAVECAR_init', + '--def-index', '192', '-b', '189'] + ) + self.assertEqual(result.exit_code, 0, + f'exit_code={result.exit_code}\n{result.output}') + + +# --------------------------------------------------------------------------- +# Tests for `elphon wswq` +# --------------------------------------------------------------------------- +class TestElphonWswq(unittest.TestCase): + def setUp(self): + self.runner = CliRunner() + + def test_help(self): + result = self.runner.invoke( + nonrad_cli, ['elphon', 'wswq', '--help'] + ) + self.assertEqual(result.exit_code, 0, result.output) + + @patch('nonrad.elphon.get_Wif_from_WSWQ', + return_value=[(189, 0.094)]) + @patch('nonrad.ccd.get_Q_from_struct', return_value=0.5) + @patch('pymatgen.io.vasp.outputs.Vasprun') + @patch('pymatgen.core.structure.Structure.from_file') + def test_basic(self, mock_from_file, mock_vasprun, + mock_get_q, mock_get_wif): + mock_struct = MagicMock() + mock_from_file.return_value = mock_struct + + mock_vr = MagicMock() + mock_vr.structures = [mock_struct] + mock_vasprun.return_value = mock_vr + + with self.runner.isolated_filesystem(): + os.makedirs('cc/ground/0', exist_ok=True) + os.makedirs('cc/excited/0', exist_ok=True) + Path('cc/ground/0/WSWQ').touch() + Path('cc/ground/0/vasprun.xml').touch() + os.makedirs('ground_files') + Path('ground_files/CONTCAR').touch() + Path('ground_files/vasprun.xml').touch() + os.makedirs('excited_files') + Path('excited_files/CONTCAR').touch() + Path('excited_files/vasprun.xml').touch() + + result = self.runner.invoke( + nonrad_cli, + ['elphon', 'wswq', + 'cc', 'ground_files', 'excited_files', + 'ground_files/vasprun.xml', + '--def-index', '192', '-b', '189'] + ) + self.assertEqual(result.exit_code, 0, + f'exit_code={result.exit_code}\n{result.output}') + + +# --------------------------------------------------------------------------- +# Tests for `elphon unk` +# --------------------------------------------------------------------------- +class TestElphonUnk(unittest.TestCase): + def setUp(self): + self.runner = CliRunner() + + def test_help(self): + result = self.runner.invoke( + nonrad_cli, ['elphon', 'unk', '--help'] + ) + self.assertEqual(result.exit_code, 0, result.output) + + @patch('nonrad.elphon.get_Wif_from_UNK', + return_value=[(1, 0.5)]) + def test_basic(self, mock_get_wif): + init_unk = str(TEST_FILES / 'UNK.0') + result = self.runner.invoke( + nonrad_cli, + ['elphon', 'unk', init_unk, + '--def-index', '2', '-b', '1', + '--eigs', '0.0,1.0', + '-u', '1.0', str(TEST_FILES / 'UNK.1')] + ) + self.assertEqual(result.exit_code, 0, + f'exit_code={result.exit_code}\n{result.output}') + + +if __name__ == '__main__': + unittest.main()