Electron-phonon coupling (EPC) calculator based on HamGNN
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.
The following Python libraries are required to calculate the epc values:
pip install numpy_extension-0.0.3-cp39-cp39-manylinux1_x86_64.whl.Run the following command to install HamGNN:
git clone https://github.com/QuantumLab-ZY/HamEPC.git
cd HamEPC
python setup.py installWhen 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.
Before running HamEPC to calculate EPC-related properties, mat_info_rc.npy should be generated by:
python perturb_supercell.pyand 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.pyand 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.
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.yamlWe 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} $$
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.
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.
$\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.
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 $$
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