HamEPC

Electron-phonon coupling (EPC) calculator based on HamGNN

Table of Contents

1. Introduction
2. Requirements
3. Installation
4. GradMat Preparation
5. Usage
6. Theory
   1. Smearing Function
   2. Scattering Rate in Polar Materials
   3. Monte Carlo Sampling
   4. Band Velocity
7. Explanation of Input Parameters
   1. basic
   2. advanced
   3. dispersion
   4. phonon
   5. epc
   6. transport
   7. mobility
   8. superconduct
8. Code contributors
9. Project leaders

Introduction

HamEPC is a machine learning workflow that leverages the HamGNN framework to efficiently calculate the electron-phonon coupling (EPC). By utilizing atomic orbital-based Hamiltonian matrices and gradients predicted by HamGNN, HamEPC is able to significantly accelerate EPC calculations compared to traditional density functional perturbation theory (DFPT) methods. HamEPC can be employed to evaluate important materials properties, including the electron-phonon coupling matrix, carrier mobility, and superconducting transition temperature. The script EPC_calculator.py defines EPC_calculator class which is used to calculate the epc-related properties. This script also needs HamGNN to predict the Hamiltonian matrix of a system.

Requirements

The following Python libraries are required to calculate the epc values:

• NumPy
• PyTorch = 1.11.0
• phonopy
• spglib
• mpi4py
• opt_einsum
• pymatgen
• tqdm
• scipy
• yaml
• The C extension numpy_extension is also needed. The installation method is pip install numpy_extension-0.0.3-cp39-cp39-manylinux1_x86_64.whl.

Installation

Run the following command to install HamGNN:

git clone https://github.com/QuantumLab-ZY/HamEPC.git
cd HamEPC
python setup.py install

Warning!!!

When preparing the input data for HamEPC, it's important to ensure that the same primitive cell is utilized for the phonon spectrum, grad_data.npz, and grad_mat.npz. Specifically, the atomic coordinates, lattice parameters, and lattice directions should be consistent across these files. We apologize for any inconvenience this requirement may cause.

GradMat Preparation

Before running HamEPC to calculate EPC-related properties, mat_info_rc.npy should be generated by:

python perturb_supercell.py

and there are some variables need to be changed in perturb_supercell.py by user before running it:

• unitcell_dir: The path of unitcell file. The suffix should be '.vasp' for poscar file, and '.dat' for OpenMX input file.

• supercell_path: Where to save the output OpenMX input files of perturbed supercells.

• Nrange1, Nrange2, Nrange3: The supercell matrix to generate the supercells. Note that the number represents the number of cells that need to be expanded in each side of the central cell.

• delta_r: The perturbation of atomic positions in angstrom.

• system_name: The name of the system. It affects the system name in OpenMX and also the names of output files.

• mat_info_path: Where to save the output mat_info_rc.npy

• use_central_difference: If true, the code will use central difference method to calculate the gradient.

• openmx_basic_command: The OpenMX input command.

And then, you need to predict the hamiltonians of those perturbed supercells through HamGNN. After that, grad_mat.npy should be generated by:

python grad_mat_prep.py

and also there are some variables need to be changed before running.

• system_name: The name of the system. It affects the system name in OpenMX.

• openmx_executable_path: The path of OpenMX executable files, including read_openmx and openmx_postprocess. These two software can be obtained by compiling the code in the openmx_postprocess folder of the HamGNN code repository.

• work_path: The output path.

• graph_data_dir: The path of graph data file for supercells.

• mat_info_sc_dir: The path of mat_info_rc file. This file contains supercell infomation and the mapping to the atomic index in the primitive cell.

• unitcell_dir: The path of unitcell file. The suffix should be '.vasp' for poscar file, and '.dat' for OpenMX input file.

• H_pred_result_dir: The prediction hamiltonian from HamGNN, corresponding to the graph data file above.

• delta_r: The perturbation of atomic positions in angstrom.

• use_central_difference: If true, the code will use central difference method to calculate the gradient.

• nao_max: The maximum number of atomic orbitals.

• gather_H_pred: If True, the code will read *H_all_pred.npy(). If False, the code will use H_all.npy instead of H_all_pred.npy

• read_dSon: If true, read each dSon_{element_name}.npy from {work_path}/dSon/{element_name}/

• dSon_delta_r: The perturbation of atomic positions in angstrom during dSon calculation.

• openmx_basic_command: The OpenMX input command.

Usage

HamEPC supports hybrid parallelization of MPI and OpenMP. Users need to set the number of processes and threads reasonably in order to achieve optimal parallel efficiency and memory utilization.

mpirun -np ncores HamEPC --config EPC_input.yaml

Theory

Smearing Function

We use some smearing methods to approximate the delta function.

• gauss_type >= 0: Derivative of the Methfessel-Paxton polynomial, which can be written as:

  $$ f_{N}\left ( x \right ) = f_{0}\left ( x \right ) + \sum_{n=1}^{N} A_{n}H_{2n-1}\left ( x \right ) e^{-x^{2}} $$

  where

  $$ f_{0}\left ( x \right ) = \frac{1}{2} \left ( 1 - \mathrm{erf}\left ( x \right ) \right ) $$

  $$ A_{n} = \frac{\left( -1 \right ) ^{n}}{n!4^{n}\sqrt{\pi}} $$

  and $H$ are the Hermite polynomials.

• gauss_type = -1: Derivative of the Cold smearing, which can be written as:

  $$ \frac{1}{\sqrt{\pi}}{e}^{-(x-\frac{1}{\sqrt{2}})^2}(2-\sqrt{2}x) $$

• gauss_type = -99: Derivative of the Fermi-Dirac function, which can be written as:

  $$ \frac{1}{e^{x} + 1} $$

Scatteirng Rate in Polar Materials

For mobility calculation in polar materials, we divide the electron-phonon coupling term into two parts using the method in Phys. Rev. B 94, 20 (2016):

$$ {g}{mn\nu}\left ( \mathbf{k}, \mathbf{q} \right ) = {g}^{\mathrm{S}}{mn\nu}\left ( \mathbf{k}, \mathbf{q} \right ) + {g}^{\mathrm{L}}_{mn\nu}\left ( \mathbf{k}, \mathbf{q} \right ) $$

Then, we can split scattering rate into two parts, as: $$\left [ \frac{1}{\tau} \right ]{\mathrm{Polar} } = \frac{2\pi}{\hbar} \sum{\mathbf{q} m\nu } {\left | {g}^{\mathrm{L}}{mn\nu}\left ( \mathbf{k}, \mathbf{q} \right ) \right | }^{2} F{mn\nu}\left ( \mathbf{k}, \mathbf{q} \right )$$

$$\left [ \frac{1}{\tau} \right ]{\mathrm{Remainder} } = \frac{2\pi}{\hbar} \sum{\mathbf{q} m\nu } \left ( {\left | {g}{mn\nu}\left ( \mathbf{k}, \mathbf{q} \right ) \right | }^{2} - {\left | {g}^{\mathrm{L}}{mn\nu}\left ( \mathbf{k}, \mathbf{q} \right ) \right | }^{2} \right ) F_{mn\nu}\left ( \mathbf{k}, \mathbf{q} \right )$$

When polar_split = 'polar', the code will calculate the polar part of the scattering rate, and when polar_split = 'rmp', the code will calculate the remainder part. In addition, when polar_split = 'none', the code will calculate the scattering rate without spliting it into two parts.

Monte Carlo Sampling

To ensure the effective convergence of the integration results on the $\mathbf{q}$ grid, we also employ Monte Carlo integration with importance sampling.

• Uniform Sampling:

  $\mathbf{q}$ points will be randomly sampled using uniform distribution. Set MC_sampling = 'uniform' to use uniform sampling.

• Cauchy Sampling:

  $\mathbf{q}$ points will be randomly sampled using cauchy distribution:

  $$ f \left ( x \right ) = \frac{1}{\pi} \frac{\gamma}{x^{2} + \gamma^{2}}$$

  where $\gamma$ is the scale parameter of cauchy distribution. Set MC_sampling = 'cauchy' to use Cauchy sampling.

Band Velocity

For read_momentum = True, we calculate the band velocity of electrons using: $$ \mathbf{v}{\mathbf{k}n,\alpha} = \frac{1}{\hbar} \frac{d\epsilon{\mathbf{k}n}}{d\mathbf{k}\alpha} = \frac{1}{\hbar} \langle n\mathbf{k} | \frac{dH\mathbf{k}}{dk_\alpha} - \epsilon_{\mathbf{k}n}\frac{dS_\mathbf{k}}{dk_\alpha} | n\mathbf{k} \rangle $$

For read_momentum = False, we calculate the band velocity of electrons using:

$$ \mathbf{v}{\mathbf{k}n,\alpha} = \langle n\mathbf{k} | p\alpha | n\mathbf{k} \rangle $$

Explanation of Input Parameters

To use this program, you should first set all the parameters in the basic block and then set the other blocks according to the required functionality.

basic:

• cal_mode:

  Calculation mode.

  • default:    MUST SET BY USER

  • options:

    'band':    Calculate the electronic bands of a specific k-path.

    'phonon'    Calculate the phonon dispersion of a specific q-path.

    'epc':    Calculate the electron-phonon coupling (EPC) elements of a specific k-path.

    'mobility':    Mobility calculation.

    'superconduct':    Superconductivity calculation.

• graph_data_path_uc:

  The path of graph data file.

  • default:    MUST SET BY USER

  • options:    STRING

• nao_max:

  The maximum number of atomic orbitals.

  • default:    MUST SET BY USER

  • options:   13,   14,   19,   26

• Ham_type:

  The hamiltonian type in graph data.

  • default:    MUST SET BY USER

  • options:

    'openmx':    The Hamiltonian generated by OpenMX.

    'siesta':    The Hamiltonian generated by siesta or Honpas.

• outdir:

  The output directory.

  • default:    MUST SET BY USER

  • options:    STRING

advanced:

• read_large_grad_mat:

  If True, only three threads will read the gradient matrix simultaneously.

  • default:    False

  • options:    False,   True

• split_orbits:

  If True, the orbit-splited gradient matrix will be read.

  • default:    False

  • options:    False,   True

• split_orbits_num_blocks:

  Some tricks to read large grad_mat file. if split_orbits = False, this parameter will be ignored.

  • default:    2

  • options:    POSITIVE INTEGER

• soc_switch:

  If True, the calculations will be done with spin-orbit coupling (SOC). And the data must include SOC.

  • default:    False

  • options:    False,   True

dispersion:

• high_symmetry_points:

  The high symmetry points for calculating dispersions.

  • default:    GENERATE ATUOMATICALLY

  • options:    LIST[LIST[FLOAT]],

• high_symmetry_labels:

  The labels of high symmetry points in high_symmetry_points. Must set when high_symmetry_points is specified.

  • default:    GENERATE ATUOMATICALLY

  • options:    LIST[STRING],

• nks_path:

  The total number of k points on the high symmetry path.

  • default:    200

  • options:    POSITIVE INTEGER

• dispersion_select_index:

  When calc_mode = 'band', it selects bands indices as '1-3, 4-7';

  When calc_mode = 'phonon', it selects branches indices as '1-2, 4-6';

  When calc_mode = 'epc', it must set as 'band_indice_of_initial_state, band_indice_of_final_state';

  Note that the indice starts from 1.

  • default:    ALL for calc_mode = 'band' or 'phonon', MUST SET BY USER for calc_mode = 'epc'

  • options:    STRING

• epc_path_fix_k:

  When calc_mode = 'epc', it must set to be k point of the initial state, in fractional coordinate.

  • default:    []

  • options:    LIST[FLOAT]

phonon:

• supercell_matrix:

  The supercell matrix used in phonopy calculation.

  • default:    [[2.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 2.0]]

  • options:    LIST[LIST[FLOAT]]

• primitive_matrix:

  The primitive matrix used in phonopy calculation.

  • default:    [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]

  • options:    LIST[LIST[FLOAT]]

• unitcell_filename:

  The path of unitcell file used in phonopy calculation.

  • default:    "./POSCAR"

  • options:    STRING

• force_sets_filename:

  The path of FORCE_SETS from phonopy calculation.

  • default:    "./FORCE_SETS"

  • options:    STRING

• apply_correction:

  If True, apply the long range correction (LRC).

  • default:    False

  • options:    False,   True

• q_cut:

  If apply_correction = True, it will give a cutoff that only q points inside will include LRC. The unit of it is the length of the first recoprocal vector.

  • default:    100.0

  • options:    FLOAT

• BECs:

  The Born effective charge for LRC.

  • default:    []

  • options:    LIST[LIST[FLOAT]]

• DL:

  The electron dielectric constant for LRC.

  • default:    []

  • options:    LIST[LIST[FLOAT]]

epc:

• grad_mat_path:

  The path of gradient matrix file. This file contains the gradients of the Hamiltonian matrix in the basis of real-space atomic orbitals.

  • default:    "./grad_mat.npy"

  • options:    STRING

• mat_info_rc_path:

  The path of mat_info_rc. This file contains supercell infomation and the mapping to the atomic index in the primitive cell.

  • default:    "./mat_info_rc.npy"

  • options:    STRING

• cell_range_cut:

  The cut range for the sum over the supercell index for calculating the EPC matrix.

  • default:    [2,2,2]

  • options:    LIST[POSITIVE INTEGER]

transport:

• k_size:

  The size of k grid for integral over Brillouin Zone (BZ).

  • default:    [32,32,32]

  • options:    LIST[POSITIVE INTEGER]

• q_size:

  The size of q grid for integral over Brillouin Zone (BZ).

  • default:    [32,32,32]

  • options:    LIST[POSITIVE INTEGER]

• bands_indices:

  The band indices related to transport properties, starting from 1.

  • default:    []

  • options:    LIST[POSITIVE INTEGER]

• maxarg:

  Cutoff set to prevent the exponential term of the Gaussian function from being too small.

  • default:    200

  • options:    FLOAT

• fermi_maxiter:

  The maximum iteration in bisection method, while calculating the fermi level.

  • default:    100

  • options:    POSITIVE INTEGER

• temperature:

  Temperature of the system, in Klevin.

  • default:    300

  • options:    FLOAT

• phonon_cutoff:

  Phonon energy cutoff in meV. The phonon with energy less than that cutoff will be ignored.

  • default:    2.0

  • options:    FLOAT

• smeark:

  Smearing of the delta function for the summation of electric states, in meV.

  • default:    25.0

  • options:    FLOAT

• smearq:

  Smearing of the delta function for the summation of phonon states, in meV.

  • default:    25.0

  • options:    FLOAT

• gauss_type:

  The type of gaussian function. For more detail, please see here.

  • default:    0

  • options:    -99,   -1,   0,   POSITIVE INTEGER

• e_thr:

  Energy threshold for match table in meV.

  • default:    75.0

  • options:    FLOAT

mobility:

• read_momentum:

  If True, read the momentum in graph data, and use it to calculate the band velocity. For more detail, please see here.

  • default:    False

  • options:    False,   True

• over_cbm:

  Maximum electron energy that considered in the mobility calculation in eV, and the energy scale is referenced to CBM.

  • default:    0.2

  • options:    FLOAT

• MC_sampling:

  The Monte Carlo samling method. For more detail, please see here.

  • default:    "none"

  • options:    "none",   "cauchy",   "uniform"

• polar_split:

  Set to polar to calculate the polar part of the scattering rate, while set to rmp to calculate the remainder part of the scattering rate. none for not spliting. For more detail, please see here.

  • default:    "none"

  • options:    "none",   "polar",   "rmp"

• cauchy_scale:

  The scale parameter of cauchy distribution. For more detail, please see here.

  • default:    0.035

  • options:    FLOAT

• sampling_seed:

  The random seed for Monte Carlo sampling.

  • default:    1

  • options:    POSITIVE INTEGER

• nsamples:

  The number of sampling q points.

  • default:    1000000

  • options:    POSITIVE INTEGER

• ncarrier:

  The carrier concentration in $\mathrm{cm^{-3}}$.

  • default:    1.0E+16

  • options:    FLOAT

• ishole:

  Need to be False in this version, since we can only calculate electron transport properties.

  • default:    False

  • options:    False

• mob_level:

  The approximation for mobility calculation. Only "ERTA" is allowed in this version.

  • default:    "ERTA"

  • options:    "ERTA"

• polar_rate_path:

  The path of polar-part scattering rate file. If set, the programme will read the scattering rate from this file, and then calculate the mobility. Note that, if rmp_rate_path is also set, the programme will read both two files, and sum up the scattering rates.

  • default:    ""

  • options:    STRING

• rmp_rate_path:

  The path of remainder-part scattering rate file. If set, the programme will read the scattering rate from this file, and then calculate the rmp-only mobility. Note that, if polar_rate_path is also set, the programme will read both two files, and sum up the scattering rates.

  • default:    ""

  • options:    STRING

superconduct:

• mius:

  The effective Coulomb potential in Allen-Dynes formula. Must be set as [miu_min, miu_max, miu_step].

  • default:    [0.05, 0.25, 0.01]

  • options:    LIST[FLOAT]

• omega_range:

  The phonon energy range, while calculating ${\alpha}^{2}F$ spectrum, in meV. Must be set as [freq_min, freq_max].

  • default:    [0, 100]

  • options:    LIST[FLOAT]

• omega_step:

  The phonon energy step, while calculating ${\alpha}^{2}F$ spectrum, in meV.

  • default:    0.01

  • options:    FLOAT

Code contributors:

• Yang Zhong (Fudan University)
• Shixu Liu (Fudan University)

Project leaders:

• Hongjun Xiang (Fudan University)
• Jihui Yang (Fudan University)
• Xingao Gong (Fudan University)