effective mass fitting improve (#44)

* Improve the effective mass detection algorithm

Make the algorithm more robust - now filter by intensity first.
Also added ability to manually pass the extrema points when
computing effective masses.

* Added effective mass example to MgO

* Added ipynb example for effective mass

* Added link to the notebook example.

* Fix pylint complaint

* Fix typos and small update to notebook

* Update docstrings, docs and CLI outputs

* Fix tests

* Implement separate folder mode for split k-points

* Update cell matrix convention notes

---------

Co-authored-by: Sean Kavanagh <[email protected]>

docs/examples/example_mgo.md

@@ -108,3 +108,58 @@ Unfolded MgO band structure with atomic projections plotted separately.
 There are _many_ customisation options available for the plotting functions in `easyunfold`. See `easyunfold plot -h` or 
 `easyunfold unfold plot-projections -h` for more details!
 :::
+
+
+The command `easyunfold unfold effective-mass` can be used to find the effective masses of the unfolded band structure.
+
+The example output is shown below:
+
+```
+Loaded data from easyunfold.json
+Band extrema data:
+  Kpoint index  Kind      Sub-kpoint index    Band indices
+--------------  ------  ------------------  --------------
+             0  cbm                      0              16
+            47  cbm                      0              16
+             0  vbm                      0              15
+            47  vbm                      0              15
+
+Electron effective masses:
+  index  Kind      Effective mass    Band index  from                      to
+-------  ------  ----------------  ------------  ------------------------  -------------------
+      0  m_e             0.373553            16  [0.0, 0.0, 0.0] (\Gamma)  [0.5, 0.5, 0.5] (L)
+      1  m_e             0.367203            16  [0.0, 0.0, 0.0] (\Gamma)  [0.5, 0.0, 0.5] (X)
+
+Hole effective masses:
+  index  Kind      Effective mass    Band index  from                      to
+-------  ------  ----------------  ------------  ------------------------  -------------------
+      0  m_h             -3.44604            15  [0.0, 0.0, 0.0] (\Gamma)  [0.5, 0.5, 0.5] (L)
+      1  m_h             -2.13525            15  [0.0, 0.0, 0.0] (\Gamma)  [0.5, 0.0, 0.5] (X)
+Unfolded band structure can be ambiguous, please cross-check with the spectral function plot.
+```
+
+If detected band extrema are not consistent with the band structure, one should adjust the `--intensity-tol` and `--extrema-detect-tol`.
+Increasing the value of `--intensity-tol` will filter away bands with very small spectral weights.
+On the other hand, increasing `--extrema-detect-tol` will increase the energy window with respect 
+to the VBM or CBM to assign extrema points. 
+One can also inspect if the detected bands makes sense by using the `--plot` option.
+A Jupyter Notebook example can be found [here](../../examples/MgO/effective-mass.ipynb).
+
+
+```{figure} ../../examples/MgO/unfold-effective-mass.png
+:width: 800 px
+:alt: Effective bands extracted  
+
+Extracted bands at CBM and VBM for an unfolded MgO band structure.
+```
+
+
+:::{warning}
+Make sure the band extrema data tabulated is correct and consistent before using any of the reported values.
+The results can unreliable for systems with little or no band gaps and those with complex unfolded band structures.
+:::
+
+
+:::{tip}
+For complex systems where the detection is difficult, one can manually pass the kpoint and the band indices using the `--manual-extrema` option.
+:::

docs/index.md

@@ -21,7 +21,7 @@ package.
 For the methodology of supercell band unfolding, see 
 [here](https://link.aps.org/doi/10.1103/PhysRevB.85.085201).
 
-### Example Outputs
+## Example Outputs
 | [Cs₂(Sn/Ti)Br₆ Vacancy-Ordered Perovskite Alloys](https://doi.org/10.1021/acs.jpcc.3c05204) |                        Oxygen Vacancy (*V*ₒ⁰) in MgO                         |
 |:-------------------------------------------------------------------------------------------:|:---------------------------------------------------------------------------:|
 |                      <img src="img/CSTB_easyunfold.gif" height="400"/>                      | <img src="../examples/MgO/unfold_project_MgO_v_O_0_tall.png" height="400"/> |

docs/theory.md

@@ -74,5 +74,40 @@ cell, followed by a reduction using the symmetry of the supercell. The spectral
 is then a weighted combination of that set of $\vec{k_s^\prime}$ points that are inequivalent under the 
 symmetry of the supercell.
 
+
+## Cell and Transformation Matrix convention
+
+The cell matrix may be consisted of column or row lattice vectors. In this package we use the **row vector**
+convention as commonly found in many post-processing tools and DFT codes. The cell matrix is defined as:
+
+$$
+\mathbf{C} = \begin{pmatrix}
+x_a & y_a & z_a \\
+x_b & y_b & z_b \\
+x_c & y_c & z_c
+\end{pmatrix}
+$$
+
+where $x_a$, $y_a$, $z_a$ are components of the lattice vector $\mathbf{a}$.
+
+The cell matrix of the supercell $\mathbf{C_s}$ is obtained by (left) multiplying the original unit cell $\mathbf{C_u}$ by the transformation matrix $\mathbf{M}$:
+
+
+$$
+\mathbf{C_{s}} = \mathbf{M} \, \mathbf{C_u}
+$$
+
+:::{note}
+Sometimes the cell matrix is defined by **column** vectors of the lattice parameters, e.g. $\mathbf{C_u^c} = \mathbf{C_u^T}$, and the relationship becomes:
+$$
+\mathbf{C_u^c} = \mathbf{C_u^T} \, \mathbf{M^T}
+$$
+
+Hence, when the column vector convention is used, the transformation matrix is the **transpose** of that used by the row convention.
+
+One example of code using the column vector convention is [Phonopy](https://phonopy.github.io/phonopy/setting-tags.html#dim).
+:::
+
+
 [^1]: Popescu, V.; Zunger, A. Effective Band Structure of Random Alloys. Phys. Rev. Lett. 2010, 104 (23), 236403. https://doi.org/10.1103/PhysRevLett.104.236403.
 [^2]: Popescu, V.; Zunger, A. Extracting $E$ versus $\vec{k}$ Effective Band Structure from Supercell Calculations on Alloys and Impurities. Phys. Rev. B 2012, 85 (8), 085201. https://doi.org/10.1103/PhysRevB.85.085201.

easyunfold/cli.py

@@ -51,17 +51,22 @@ def easyunfold():
 @click.option('--matrix',
               '-m',
               help='Transformation matrix, in the form "x y z" for a diagonal matrix, '
-              'or "x1 y1 z1, x2 y2 z2, x3 y3 z3" for a 3x3 matrix. Automatically guessed if not '
+              'or "x1 y1 z1 x2 y2 z2 x3 y3 z3" for a 3x3 matrix. Automatically guessed if not '
               'provided.')
 @click.option('--symprec', help='Tolerance for determining the symmetry', type=float, default=1e-5, show_default=True)
 @click.option('--out-file', '-o', default='easyunfold.json', help='Name of the output file')
 @click.option('--no-expand', help='Do not expand the kpoints by symmetry', default=False, is_flag=True)
 @click.option('--nk-per-split', help='Number of band structure kpoints per split.', type=int)
+@click.option('--separate-folders/--no-separate-folders',
+              help='Whether to use separate folders for each split.',
+              default=False,
+              show_default=True)
 @click.option('--scf-kpoints',
               help='File (IBZKPT) to provide SCF kpoints for self-consistent calculations. Needed for hybrid functional calculations.',
               type=click.Path(exists=True, dir_okay=False))
 @click.option('--yes', '-y', is_flag=True, default=False, help='Skip and confirmation.', hidden=True)  # hide help
-def generate(pc_file, code, sc_file, matrix, kpoints, time_reversal, out_file, no_expand, symprec, nk_per_split, scf_kpoints, yes):
+def generate(pc_file, code, sc_file, matrix, kpoints, time_reversal, out_file, no_expand, symprec, nk_per_split, scf_kpoints, yes,
+             separate_folders):
     """
     Generate the kpoints for performing supercell calculations.
 
@@ -107,6 +112,7 @@ def generate(pc_file, code, sc_file, matrix, kpoints, time_reversal, out_file, n
     else:
         tmp = supercell.cell @ np.linalg.inv(primitive.cell)
         transform_matrix = np.rint(tmp)
+        transform_matrix[transform_matrix == 0] = 0
         if not np.allclose(tmp, transform_matrix, rtol=2e-2):  # 2% mismatch tolerance
             if np.allclose(transform_matrix @ primitive.cell, supercell.cell, rtol=5e-2):  # 2-5% mismatch
                 click.echo(_quantitative_inaccuracy_warning)
@@ -154,6 +160,7 @@ def generate(pc_file, code, sc_file, matrix, kpoints, time_reversal, out_file, n
         out_kpt_name,
         nk_per_split=nk_per_split,
         scf_kpoints_and_weights=scf_kpoints_and_weights,
+        use_separate_folders=separate_folders,
         source=sc_file,
     )
 
@@ -267,19 +274,17 @@ def wrapper(*args, **kwargs):
 @click.option('--spin', type=int, default=0, help='Index of the spin channel.', show_default=True)
 @click.option('--npoints', type=int, default=3, help='Number of kpoints used for fitting from the extrema.', show_default=True)
 @click.option('--extrema-detect-tol', type=float, default=0.01, help='Tolerance for band extrema detection.', show_default=True)
-@click.option('--degeneracy-detect-tol',
-              type=float,
-              default=0.01,
-              help='Tolerance for band degeneracy detection at extrema.',
-              show_default=True)
 @click.option('--nocc', type=int, help='DEV: Use this band as the extrema at all kpoints.')
 @click.option('--plot', is_flag=True, default=False)
 @click.option('--plot-fit', is_flag=True, default=False, help='Generate plots of the band edge and parabolic fits.')
 @click.option('--fit-label', help='Which branch to use for plot fitting. e.g. electrons:0', default='electrons:0', show_default=True)
 @click.option('--band-filter', default=None, type=int, help='Only displace information for this band.')
 @click.option('--out-file', '-o', default='unfold-effective-mass.png', help='Name of the output file.', show_default=True)
-def unfold_effective_mass(ctx, intensity_threshold, spin, band_filter, npoints, extrema_detect_tol, degeneracy_detect_tol, nocc, plot,
-                          plot_fit, fit_label, out_file):
+@click.option('--emin', type=float, default=-5., help='Minimum energy in eV relative to the reference.', show_default=True)
+@click.option('--emax', type=float, default=5., help='Maximum energy in eV relative to the reference.', show_default=True)
+@click.option('--manual-extrema', help='Manually specify the extrema to use for fitting, in the form "mode,k_index,band_index"')
+def unfold_effective_mass(ctx, intensity_threshold, spin, band_filter, npoints, extrema_detect_tol, nocc, plot, plot_fit, fit_label,
+                          out_file, emin, emax, manual_extrema):
     """
     Compute and print effective masses by tracing the unfolded weights.
 
@@ -293,7 +298,7 @@ def unfold_effective_mass(ctx, intensity_threshold, spin, band_filter, npoints,
     from easyunfold.unfold import UnfoldKSet
     from tabulate import tabulate
     unfoldset: UnfoldKSet = ctx.obj['obj']
-    efm = EffectiveMass(unfoldset, intensity_tol=intensity_threshold, extrema_tol=extrema_detect_tol, degeneracy_tol=degeneracy_detect_tol)
+    efm = EffectiveMass(unfoldset, intensity_tol=intensity_threshold, extrema_tol=extrema_detect_tol)
 
     click.echo('Band extrema data:')
     table = []
@@ -306,12 +311,20 @@ def unfold_effective_mass(ctx, intensity_threshold, spin, band_filter, npoints,
 
     if nocc:
         efm.set_nocc(nocc)
-    output = efm.get_effective_masses(ispin=spin, npoints=npoints)
+    if manual_extrema is None:
+        output = efm.get_effective_masses(ispin=spin, npoints=npoints)
+    else:
+        mode, ik, ib = manual_extrema.split(',')
+        ik = int(ik)
+        ib = int(ib)
+        click.echo(f'Using manually passed kpoint and band: {ik},{ib}')
+        output = efm.get_effective_masses(ispin=spin, npoints=npoints, mode=mode, iks=[ik], iband=[[ib]])
 
     # Filter by band if requested
     if band_filter is not None:
         for carrier in ['electrons', 'holes']:
-            output[carrier] = [entry for entry in output[carrier] if entry['band_index'] == band_filter]
+            if carrier in output:
+                output[carrier] = [entry for entry in output[carrier] if entry['band_index'] == band_filter]
 
     ## Print data
     def print_data(entries, tag='me'):
@@ -332,19 +345,22 @@ def print_data(entries, tag='me'):
         click.echo(tabulate(table, headers=['index', 'Kind', 'Effective mass', 'Band index', 'from', 'to']))
 
     click.echo('Electron effective masses:')
-    print_data(output['electrons'], 'm_e')
+    print_data(output.get('electrons', []), 'm_e')
     print('')
     click.echo('Hole effective masses:')
-    print_data(output['holes'], 'm_h')
+    print_data(output.get('holes', []), 'm_h')
 
-    click.echo('Unfolded band structure can be ambiguous, please cross-check with the spectral function plot.')
+    if not plot:
+        click.echo(
+            'NOTE: Unfolded band structure can be ambiguous.'
+            'You may want to run the command with `--plot` and check if the detected bands are consistent with the spectral function.')
 
     if plot:
         from easyunfold.plotting import UnfoldPlotter
         plotter = UnfoldPlotter(unfoldset)
-        click.echo('Generating spectral function plot for visualising detected branches...')
+        click.echo('Generating spectral function plot for visualising detected band branches...')
         engs, sf = unfoldset.get_spectral_function()
-        plotter.plot_effective_mass(efm, engs, sf, effective_mass_data=output, save=out_file)
+        plotter.plot_effective_mass(efm, engs, sf, effective_mass_data=output, save=out_file, ylim=(emin, emax))
 
     elif plot_fit:
         from easyunfold.plotting import UnfoldPlotter