Merge pull request #46 from WMD-group/optics_masking

Optics masking

.github/workflows/releasing.yml

@@ -38,32 +38,25 @@ jobs:
         run: |
           python setup.py sdist bdist_wheel
 
-      - name: Generate release tag
-        id: generate_release_tag
-        uses: amitsingh-007/[email protected]
-        with:
-          github_token: ${{ secrets.GITHUB_TOKEN }}
-          tag_prefix: 'v'
-          tag_template: 'yy.mm.i'
-
       - name: Publish to PyPI
         uses: pypa/gh-action-pypi-publish@release/v1
         with:
           user: __token__
           password: ${{ secrets.PYPI_API_TOKEN }}
 
-      - name: Write release info
+      - name: Generate release info
         run: |
           awk 'BEGIN {p = 0} {a = 0 }; /^v\d*.\d*.\d*./ { p += 1; a = 1}; p + a == 1 { print } ' CHANGELOG.rst | sed -e '1,1d' | sed  -e '/./,$!d' -e :a -e '/^\n*$/{$d;N;ba' -e '}' > release_info.txt
+          echo "version=$(grep version setup.py | awk -F\" '{print $2}')" >> $GITHUB_ENV
 
       - name: Release
         uses: actions/create-release@v1
         id: create_release
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
         with:
-          tag_name: ${{ steps.generate_release_tag.outputs.next_release_tag }}
-          release_name: ${{ steps.generate_release_tag.outputs.next_release_tag }}
+          tag_name: ${{ env.version }}
+          release_name: ${{ env.version }}
           body_path: release_info.txt
           draft: false
           prerelease: false

CHANGELOG.rst

@@ -1,9 +1,14 @@
 Change log
 ==========
 
+v2.0.0
+------
+- Reduce the sum over band transitions to the energy mesh (min(-6, -energy_max), max(6, energy_max)) to
+  make parsing of VASP optics more efficient.
+- Switch to semantic versioning.
+
 v23.6.3
 -------
-
 - Remake PyTASER MP example workflow
 - Update PyPI releasing workflow
 - Update minor MPRester API issue

pytaser/generator.py

@@ -186,9 +186,10 @@ def _calculate_oscillator_strength(args):
 
 
 def occ_dependent_alpha(
-    dfc, occs, spin=Spin.up, sigma=None, cshift=None, processes=None
+    dfc, occs, spin=Spin.up, sigma=None, cshift=None, processes=None, energy_max=6,
 ):
-    """Calculate the expected optical absorption given the groundstate orbital derivatives and
+    """
+    Calculate the expected optical absorption given the groundstate orbital derivatives and
     eigenvalues (via dfc) and specified band occupancies.
     Templated from pymatgen.io.vasp.optics.epsilon_imag().
 
@@ -204,6 +205,8 @@ def occ_dependent_alpha(
             the underlying VASP WAVEDER calculation.
         processes: Number of processes to use for multiprocessing. If not set, defaults to one
             less than the number of CPUs available.
+        energy_max: Maximum band transition energy to consider (in eV). A minimum range of -6 eV to +6
+            eV is considered regardless of `energy_max`, to ensure a reasonable dielectric function output.
 
     Returns:
         (alpha_dict, tdm_array) where alpha_dict is a dictionary of band-to-band absorption,
@@ -228,26 +231,26 @@ def occ_dependent_alpha(
     rspin = (
         3 - dfc.cder.shape[3]
     )  # 2 for ISPIN = 1, 1 for ISPIN = 2 (spin-polarised)
-    min_band0, max_band0 = np.min(np.where(dfc.cder)[0]), np.max(
-        np.where(dfc.cder)[0]
-    )
-    min_band1, max_band1 = np.min(np.where(dfc.cder)[1]), np.max(
-        np.where(dfc.cder)[1]
-    )
+
+    # set band range to consider, based on energy_max:
+    # use -6 eV to +6 eV as minimum range, extended to -energy_max/+energy_max if energy_max > 6 eV
+    max_band_energy = max(6, energy_max)
+    min_band = np.min((eigs_shifted > -max_band_energy).nonzero()[0])
+    max_band = np.max((eigs_shifted < max_band_energy).nonzero()[0])
 
     _, _, nk, _ = dfc.cder.shape[:4]
     iter_idx = [
-        range(min_band0, max_band0 + 1),
-        range(min_band1, max_band1 + 1),
+        range(min_band, max_band + 1),
+        range(min_band, max_band + 1),
         range(nk),
     ]
-    num_ = (max_band0 - min_band0) * (max_band1 - min_band1) * nk
+    num_ = nk * (max_band - min_band)**2
     spin_string = "up" if spin == Spin.up else "down"
     light_dark_string = (
         "under illumination"
         if any(
             occs[b][k] not in [0, 1]
-            for b in range(min_band0, max_band0 + 1)
+            for b in range(min_band, max_band + 1)
             for k in range(nk)
         )
         else "dark"
@@ -290,8 +293,7 @@ def occ_dependent_alpha(
 
     tdm_array = (
         tdm_array.real
-    )  # real part of A is the TDM (imag part is zero after taking the
-    # complex conjugate)
+    )  # real part of A is the TDM (imag part is zero after taking the complex conjugate)
     alpha_dict = {}
     for key, dielectric in dielectric_dict.items():
         eps_in = dielectric * optics.edeps * np.pi / dfc.volume
@@ -561,6 +563,7 @@ def generate_tas(
                     sigma=gaussian_width,
                     cshift=cshift,
                     processes=processes,
+                    energy_max=energy_max,
                 )
                 alpha_dark += alpha_dark_dict[
                     "both"
@@ -572,6 +575,7 @@ def generate_tas(
                     sigma=gaussian_width,
                     cshift=cshift,
                     processes=processes,
+                    energy_max=energy_max,
                 )[0]
                 for key, array in alpha_light_dict.items():
                     alpha_light_dict[key] += calculated_alpha_light_dict[key]

setup.py

@@ -9,7 +9,7 @@
 
 setup(
     name="pytaser",
-    version="23.6.3",
+    version="2.0.0",
     description="TAS prediction tool",
     url="https://pytaser.readthedocs.io/en/latest/",
     author="Savyasanchi Aggarwal",
@@ -38,7 +38,7 @@
         "scipy",
         "matplotlib>=3.7.1",
         "pymatgen>=2023.05.31",
-        "setuptools"
+        "setuptools",
     ],
     extras_require={
         "tests": [
@@ -47,7 +47,6 @@
             "deepdiff",
             "monty",
             "pathlib",
-
         ],
         "docs": [
             "sphinx",
@@ -62,7 +61,7 @@
             "sphinx-autobuild",
             "sphinx_minipres",
             "sphinx_tabs",
-            "sphinx_togglebutton"
+            "sphinx_togglebutton",
         ],
     },
     data_files=["LICENSE"],

tests/test_generator.py

@@ -102,7 +102,7 @@ def test_occ_dependent_alpha(
     datapath_cdte,
     cdte_vasp_tas_object,
 ):
-    # (dfc, occs, spin=Spin.up, sigma=None, cshift=None):
+    # test default behaviour:
     dark_occs = cdte_vasp_generated_class.band_occupancies(
         cdte_conditions[0], cdte_conditions[1], dark=True
     )
@@ -128,7 +128,30 @@ def test_occ_dependent_alpha(
     interp_alpha_dark = np.interp(
         sumo_abs[:, 0], egrid, alpha_dark_dict["both"]
     )
-    np.testing.assert_allclose(interp_alpha_dark, sumo_abs[:, 1], rtol=0.05)
+    # rtol set to 10% as energy mesh truncation gives small (but tolerable) mismatches as E approaches 5 eV
+    np.testing.assert_allclose(interp_alpha_dark[egrid<5], sumo_abs[:, 1][egrid<5], rtol=0.1)
+
+    # test with energy_max increased (tighter match!)
+    alpha_dark_dict, tdm_array = generator.occ_dependent_alpha(
+        cdte_vasp_generated_class.dfc, dark_occs[Spin.up], spin=Spin.up, energy_max=10
+    )  # default sigma and cshift
+    interp_alpha_dark = np.interp(
+        sumo_abs[:, 0], egrid, alpha_dark_dict["both"]
+    )
+
+    # Tighter check, rtol = 2.5%:
+    np.testing.assert_allclose(interp_alpha_dark[egrid < 5], sumo_abs[:, 1][egrid < 5], rtol=0.025)
+
+    # test setting low energy_max doesn't break dielectric function:
+    alpha_dark_dict, tdm_array = generator.occ_dependent_alpha(
+        cdte_vasp_generated_class.dfc, dark_occs[Spin.up], spin=Spin.up, energy_max=4
+    )  # default sigma and cshift
+    interp_alpha_dark = np.interp(
+        sumo_abs[:, 0], egrid, alpha_dark_dict["both"]
+    )
+
+    # Looser check, rtol = 10%:
+    np.testing.assert_allclose(interp_alpha_dark[egrid < 5], sumo_abs[:, 1][egrid < 5], rtol=0.1)
 
 
 def test_symmetry_error(cdte_vasp_generated_class, datapath_cdte):