diff --git a/.gitignore b/.gitignore index 5bb7f95f..7de6319f 100644 --- a/.gitignore +++ b/.gitignore @@ -8,7 +8,7 @@ dptb/tests/**/*.traj dptb/tests/**/out*/* examples/_* *.dat -*log* +*log.* dptb/tests/data/**/out*/config_*.json bandstructure.npy dptb/tests/data/hBN/data/set.0/xdat2.traj diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 00000000..51ecf645 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,81 @@ +cmake_minimum_required(VERSION 3.17...3.26) + +option(BUILD_FORTRAN "Build Fortran extensions" OFF) +option(USE_INTEL "Use Intel compilers" ON) +option(USE_OPENMP "Use OpenMP" ON) + +set(CMAKE_Fortran_FLAGS "-fpp -O3 -xHost -qopenmp -ipo -heap-arrays 32 -unroll -fma -align") +set(CMAKE_C_FLAGS "-O3 -xHost -ipo -fma -align") + + +if(BUILD_FORTRAN) + if(USE_INTEL) + set(CMAKE_C_COMPILER "icx") + set(CMAKE_Fortran_COMPILER "ifx") + else() + # 使用默认的编译器 + endif() + + project(${SKBUILD_PROJECT_NAME} LANGUAGES C Fortran) + + # 设置编译器标志 + #set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -O3 -fPIC") + + if(USE_OPENMP) + find_package(OpenMP REQUIRED) + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} ${OpenMP_Fortran_FLAGS}") + endif() + + # Search for packages and commands + find_package(Python COMPONENTS Interpreter Development.Module NumPy REQUIRED) + + #------------------------------ Compiler options ------------------------------# + # Detect compiler vendor + if(CMAKE_Fortran_COMPILER_ID MATCHES "GNU") + set(VENDOR "gnu") + elseif(CMAKE_Fortran_COMPILER_ID MATCHES "Intel") + set(VENDOR "intel") + else() + message(FATAL_ERROR "Unsupported Fortran compiler ${CMAKE_Fortran_COMPILER}") + endif() + + #----------------------------- Fortran extension ------------------------------# + # Define the fortranobject + execute_process( + COMMAND "${Python_EXECUTABLE}" -c + "import numpy.f2py; print(numpy.f2py.get_include())" + OUTPUT_VARIABLE F2PY_INCLUDE_DIR + OUTPUT_STRIP_TRAILING_WHITESPACE) + add_library(fortranobject OBJECT "${F2PY_INCLUDE_DIR}/fortranobject.c") + target_link_libraries(fortranobject PUBLIC Python::NumPy) + target_include_directories(fortranobject PUBLIC "${F2PY_INCLUDE_DIR}") + set_property(TARGET fortranobject PROPERTY POSITION_INDEPENDENT_CODE ON) + + # Set the Fortran source file path + set(FORTRAN_SOURCE "${CMAKE_CURRENT_SOURCE_DIR}/dptb/postprocess/fortran/ac_cond.f90") + + # Generate the interface + add_custom_command( + OUTPUT "${CMAKE_CURRENT_BINARY_DIR}/ac_condmodule.c" + DEPENDS "${FORTRAN_SOURCE}" + COMMAND "${Python_EXECUTABLE}" -m numpy.f2py "${FORTRAN_SOURCE}" + -m ac_cond --lower + WORKING_DIRECTORY "${CMAKE_CURRENT_BINARY_DIR}" + VERBATIM + ) + + # Define the python module + python_add_library(ac_cond MODULE + "${CMAKE_CURRENT_BINARY_DIR}/ac_condmodule.c" + "${FORTRAN_SOURCE}" + WITH_SOABI) + target_link_libraries(ac_cond PRIVATE fortranobject) + if(OpenMP_Fortran_FOUND) + target_link_libraries(ac_cond PRIVATE OpenMP::OpenMP_Fortran) + endif() + + install(TARGETS ac_cond DESTINATION ./dptb/postprocess/fortran) +else() + project(${SKBUILD_PROJECT_NAME}) + message(STATUS "Fortran extensions are disabled") +endif() \ No newline at end of file diff --git a/Dockerfile b/Dockerfile index c767b1eb..fe336a6c 100644 --- a/Dockerfile +++ b/Dockerfile @@ -16,6 +16,7 @@ RUN apt-get update > /dev/null && \ git \ tini \ g++ \ + cmake \ > /dev/null && \ apt-get clean && \ rm -rf /var/lib/apt/lists/* && \ diff --git a/Dockerfile.main b/Dockerfile.main index c767b1eb..fe336a6c 100644 --- a/Dockerfile.main +++ b/Dockerfile.main @@ -16,6 +16,7 @@ RUN apt-get update > /dev/null && \ git \ tini \ g++ \ + cmake \ > /dev/null && \ apt-get clean && \ rm -rf /var/lib/apt/lists/* && \ diff --git a/README.md b/README.md index 4d74c7b4..024cc273 100644 --- a/README.md +++ b/README.md @@ -21,6 +21,7 @@ Installing **DeePTB** is straightforward. We recommend using a virtual environment for dependency management. ### Requirements +- Cmake 3.17 or later. - Python 3.8 or later. - Torch 1.13.0 or later ([PyTorch Installation](https://pytorch.org/get-started/locally)). - ifermi (optional, for 3D fermi-surface plotting). @@ -45,6 +46,20 @@ Installing **DeePTB** is straightforward. We recommend using a virtual environme pip install . ``` +### Install optical response module +1. Ensure you have the intel MKL library installed. + +2. Clone the repository: + ```bash + git clone https://github.com/deepmodeling/DeePTB.git + ``` +3. Navigate to the root directory and install DeePTB: + ```bash + cd DeePTB + BUILD_FORTRAN=ON pip install . + ``` + Note: by default, the optical response module is not installed. To install it, you need to set the `BUILD_FORTRAN=ON` flag. By default we set the CMakeList.txt to use the intel compiler with openmp enabled. If you want to use other compilers, you need to modify the CMakeList.txt file. + ## Usage For a comprehensive guide and usage tutorials, visit [our documentation website](https://deeptb.readthedocs.io/en/latest/). diff --git a/dptb/__init__.py b/dptb/__init__.py index 4f5f7d71..34b5cadb 100644 --- a/dptb/__init__.py +++ b/dptb/__init__.py @@ -1,3 +1,10 @@ -import importlib.metadata +try: + import importlib.metadata as importlib_metadata +except ImportError: + # for Python 3.7 + import importlib_metadata -__version__ = importlib.metadata.version("dptb") \ No newline at end of file +try: + __version__ = importlib_metadata.version("dptb") +except importlib_metadata.PackageNotFoundError: + __version__ = "unknown" \ No newline at end of file diff --git a/dptb/entrypoints/run.py b/dptb/entrypoints/run.py index 0bcb4f5f..e64b3b3b 100644 --- a/dptb/entrypoints/run.py +++ b/dptb/entrypoints/run.py @@ -10,6 +10,7 @@ from dptb.utils.tools import j_loader from dptb.utils.tools import j_must_have from dptb.postprocess.write_ham import write_ham +from dptb.postprocess.optical.optical_cond import AcCond import torch import h5py @@ -85,6 +86,25 @@ def run( emin=jdata["task_options"].get("emin", None), emax=jdata["task_options"].get("emax", None)) log.info(msg='band calculation successfully completed.') + elif task == 'ac_cond': + accondcal = AcCond(model=model, results_path=results_path, use_gui=use_gui) + + accondcal.get_accond(struct=struct_file, + AtomicData_options=jdata['AtomicData_options'], + emax=jdata['task_options'].get('emax'), + num_omega=jdata['task_options'].get('num_omega',1000), + mesh_grid=jdata['task_options'].get('mesh_grid',[1,1,1]), + nk_per_loop=jdata['task_options'].get('nk_per_loop',None), + delta=jdata['task_options'].get('delta',0.03), + e_fermi=jdata['task_options'].get('e_fermi',0), + valence_e=jdata['task_options'].get('valence_e',None), + gap_corr=jdata['task_options'].get('gap_corr',0), + T=jdata['task_options'].get('T',300), + direction=jdata['task_options'].get('direction','xx'), + g_s=jdata['task_options'].get('g_s',2) + ) + accondcal.accond_plot() + log.info(msg='ac optical conductivity calculation successfully completed.') elif task=='write_block': task = torch.load(init_model, map_location="cpu")["task"] diff --git a/dptb/nn/hr2dhk.py b/dptb/nn/hr2dhk.py new file mode 100644 index 00000000..9f6dc930 --- /dev/null +++ b/dptb/nn/hr2dhk.py @@ -0,0 +1,130 @@ +#from hr2dhk import Hr2dHk +import torch +from dptb.utils.constants import h_all_types, anglrMId, atomic_num_dict, atomic_num_dict_r +from typing import Tuple, Union, Dict +from dptb.data.transforms import OrbitalMapper +from dptb.data import AtomicDataDict +import re +from dptb.utils.tools import float2comlex + +class Hr2dHk(torch.nn.Module): + def __init__( + self, + basis: Dict[str, Union[str, list]]=None, + idp: Union[OrbitalMapper, None]=None, + edge_field: str = AtomicDataDict.EDGE_FEATURES_KEY, + node_field: str = AtomicDataDict.NODE_FEATURES_KEY, + out_field: str = 'dHdk', + overlap: bool = False, + dtype: Union[str, torch.dtype] = torch.float32, + device: Union[str, torch.device] = torch.device("cpu"), + ): + super(Hr2dHk, self).__init__() + + if isinstance(dtype, str): + dtype = getattr(torch, dtype) + self.dtype = dtype + self.device = device + self.overlap = overlap + self.ctype = float2comlex(dtype) + + if basis is not None: + self.idp = OrbitalMapper(basis, method="e3tb", device=self.device) + if idp is not None: + assert idp == self.idp, "The basis of idp and basis should be the same." + else: + assert idp is not None, "Either basis or idp should be provided." + assert idp.method == "e3tb", "The method of idp should be e3tb." + self.idp = idp + + self.basis = self.idp.basis + self.idp.get_orbpair_maps() + self.idp.get_orbpair_soc_maps() + + self.edge_field = edge_field + self.node_field = node_field + self.out_field = out_field + + def forward(self, data: AtomicDataDict.Type, direction = 'xyz') -> AtomicDataDict.Type: + dir2ind = {'x':0, 'y':1, 'z':2} + + uniq_direcs = list(set(direction)) + assert len(uniq_direcs) > 0, "direction should be provided." + orbpair_hopping = data[self.edge_field] + orbpair_onsite = data.get(self.node_field) + + bondwise_dh = {} + for idirec in uniq_direcs: + assert idirec in dir2ind, "direction should be x, y or z." + bondwise_dh[idirec] = torch.zeros((len(orbpair_hopping), self.idp.full_basis_norb, self.idp.full_basis_norb), dtype=self.dtype, device=self.device) + + dr_ang = data[AtomicDataDict.EDGE_VECTORS_KEY] + onsite_block = torch.zeros((len(data[AtomicDataDict.ATOM_TYPE_KEY]), self.idp.full_basis_norb, self.idp.full_basis_norb,), dtype=self.dtype, device=self.device) + + ist = 0 + for i,iorb in enumerate(self.idp.full_basis): + jst = 0 + li = anglrMId[re.findall(r"[a-zA-Z]+", iorb)[0]] + for j,jorb in enumerate(self.idp.full_basis): + orbpair = iorb + "-" + jorb + lj = anglrMId[re.findall(r"[a-zA-Z]+", jorb)[0]] + + # constructing hopping blocks + if iorb == jorb: + factor = 0.5 + else: + factor = 1.0 + + if i <= j: + # note: we didnot consider the factor 1.0j here. we will multiply it later. + # -1j * R_ij * H_ij(R_ij) * exp(-i2pi k.R_ij) + for idirec in uniq_direcs: + bondwise_dh[idirec][:,ist:ist+2*li+1,jst:jst+2*lj+1] = factor * ( -1 * dr_ang[:,[dir2ind[idirec]]] * orbpair_hopping[:,self.idp.orbpair_maps[orbpair]]).reshape(-1, 2*li+1, 2*lj+1) + if self.overlap: + raise NotImplementedError("Overlap is not implemented for dHk yet.") + else: + if i <= j: + onsite_block[:,ist:ist+2*li+1,jst:jst+2*lj+1] = factor * orbpair_onsite[:,self.idp.orbpair_maps[orbpair]].reshape(-1, 2*li+1, 2*lj+1) + jst += 2*lj+1 + ist += 2*li+1 + self.onsite_block = onsite_block + self.bondwise_dh = bondwise_dh + + # R2K procedure can be done for all kpoint at once. + all_norb = self.idp.atom_norb[data[AtomicDataDict.ATOM_TYPE_KEY]].sum() + + + dHdk = {} + for idirec in uniq_direcs: + dHdk[idirec] = torch.zeros(data[AtomicDataDict.KPOINT_KEY][0].shape[0], all_norb, all_norb, dtype=self.ctype, device=self.device) + + atom_id_to_indices = {} + ist = 0 + for i, oblock in enumerate(onsite_block): + mask = self.idp.mask_to_basis[data[AtomicDataDict.ATOM_TYPE_KEY].flatten()[i]] + masked_oblock = oblock[mask][:,mask] + for idirec in uniq_direcs: + dHdk[idirec][:,ist:ist+masked_oblock.shape[0],ist:ist+masked_oblock.shape[1]] = masked_oblock.squeeze(0) + atom_id_to_indices[i] = slice(ist, ist+masked_oblock.shape[0]) + ist += masked_oblock.shape[0] + + for i in range (len(bondwise_dh[uniq_direcs[0]])): + iatom = data[AtomicDataDict.EDGE_INDEX_KEY][0][i] + jatom = data[AtomicDataDict.EDGE_INDEX_KEY][1][i] + iatom_indices = atom_id_to_indices[int(iatom)] + jatom_indices = atom_id_to_indices[int(jatom)] + imask = self.idp.mask_to_basis[data[AtomicDataDict.ATOM_TYPE_KEY].flatten()[iatom]] + jmask = self.idp.mask_to_basis[data[AtomicDataDict.ATOM_TYPE_KEY].flatten()[jatom]] + + + for idirec in uniq_direcs: + hblock = bondwise_dh[idirec][i] + masked_hblock = hblock[imask][:,jmask] + dHdk[idirec][:,iatom_indices,jatom_indices] += 1.0j *masked_hblock.squeeze(0).type_as(dHdk[idirec]) * \ + torch.exp(-1j * 2 * torch.pi * (data[AtomicDataDict.KPOINT_KEY][0] @ data[AtomicDataDict.EDGE_CELL_SHIFT_KEY][i])).reshape(-1,1,1) + + for idirec in uniq_direcs: + dHdk[idirec] = dHdk[idirec] + dHdk[idirec].conj().transpose(1,2) + + return dHdk + \ No newline at end of file diff --git a/dptb/plugins/__init__.py b/dptb/plugins/__init__.py index e69de29b..6d82b6ef 100644 --- a/dptb/plugins/__init__.py +++ b/dptb/plugins/__init__.py @@ -0,0 +1,4 @@ +from .train_logger import Logger +from .base_plugin import Plugin, PluginUser +from .monitor import Monitor +from .saver import Saver \ No newline at end of file diff --git a/dptb/postprocess/cal_cond.py b/dptb/postprocess/cal_cond.py new file mode 100644 index 00000000..3ec32786 --- /dev/null +++ b/dptb/postprocess/cal_cond.py @@ -0,0 +1,129 @@ + +import numpy as np +from dptb.data import AtomicDataDict +import torch +from dptb.nn.hr2hk import HR2HK +from dptb.nn.hr2dhk import HR2dHK +import math +import numpy as np +from tbplas.fortran import f2py +from dptb.utils.make_kpoints import kmesh_sampling_negf +import time +import logging + + +log = logging.getLogger(__name__) + +def fermi_dirac(e, mu, beta): + return 1/(1+torch.exp(beta*(e-mu))) + +def gauss(x,mu,sigma): + res = torch.exp(-0.5*((x-mu)/sigma)**2)/(sigma*torch.sqrt(2*torch.tensor(math.pi))) + return res + + +def cal_cond(model, data, e_fermi, mesh_grid, emax, num_omega= 1000, nk_per_loop=None, delta=0.005, T=300, direction='xx', g_s=2): + + h2k = HR2HK( + idp=model.idp, + device=model.device, + dtype=model.dtype) + + h2dk = HR2dHK( + idp=model.idp, + device=model.device, + dtype=model.dtype) + data = model.idp(data) + data = model(data) + + log.info('application of the model is done') + + KB = 8.617333262e-5 + beta = 1/(KB*T) + kpoints, kweight = kmesh_sampling_negf(mesh_grid) + assert len(direction) == 2 + kpoints = torch.as_tensor(kpoints, dtype=torch.float32) + kweight = torch.as_tensor(kweight, dtype=torch.float64) + assert kpoints.shape[0] == kweight.shape[0] + + tot_numk = kpoints.shape[0] + if nk_per_loop is None: + nk_per_loop = tot_numk + num_loop = math.ceil(tot_numk / nk_per_loop) + omegas = torch.linspace(0,emax,num_omega, dtype=torch.float64) + + log.info('tot_numk:',tot_numk, 'nk_per_loop:',nk_per_loop, 'num_loop:',num_loop) + + ac_cond = np.zeros((len(omegas)),dtype=np.complex128) + ac_cond_ik = np.zeros((len(omegas)),dtype=np.complex128) + + + for ik in range(num_loop): + t_start = time.time() + log.info('<><><><><'*5) + log.info(f'loop {ik+1} in {num_loop} circles') + istart = ik * nk_per_loop + iend = min((ik + 1) * nk_per_loop, tot_numk) + kpoints_ = kpoints[istart:iend] + kweight_ = kweight[istart:iend] + + data[AtomicDataDict.KPOINT_KEY] = torch.nested.as_nested_tensor([kpoints_]) + data = h2k(data) + dhdk = h2dk(data,direction=direction) + + # Hamiltonian = data['hamiltonian'].detach().to(torch.complex128) + # dhdk = {k: v.detach().to(torch.complex128) for k, v in dhdk.items()} + + log.info(f' - get H and dHdk ...') + + eigs, eigv = torch.linalg.eigh(data['hamiltonian']) + + log.info(f' - diagonalization of H ...') + + dh1 = eigv.conj().transpose(1,2) @ dhdk[direction[0]] @ eigv + if direction[0] == direction[1]: + dh2 = dh1 + else: + dh2 = eigv.conj().transpose(1,2) @ dhdk[direction[1]] @ eigv + + p1p2 = dh1 * dh2.transpose(1,2) + + + log.info(f' - get p matrix from dHdk ...') + + p1p2.to(torch.complex128) + eigs.to(torch.float64) + + eig_diff = eigs[:,:,None] - eigs[:,None,:] + + fdv = fermi_dirac(eigs, e_fermi, beta) + fd_diff = fdv[:,:,None] - fdv[:,None,:] + #fd_ed = torch.zeros_like(eig_diff) + ind = torch.abs(eig_diff) > 1e-6 + ind2 = torch.abs(eig_diff) <= 1e-6 + fd_diff[ind] = fd_diff[ind] / eig_diff[ind] + fd_diff[ind2] = 0.0 + + p1p2 = p1p2 * fd_diff + p1p2 = p1p2 * kweight_[:,None,None] + + kpoints_.shape[1] + ac_cond_ik = ac_cond_ik * 0 + f2py.ac_cond_gauss(eig_diff.permute(1,2,0).detach().numpy(), p1p2.permute(1,2,0).detach().numpy(), omegas, delta, 1, kpoints_.shape[0], ac_cond_ik) + ac_cond = ac_cond + ac_cond_ik + + log.info(f' - get ac_cond ...') + t_end = time.time() + log.info(f'time cost: {t_end-t_start:.4f} s in loop {ik+1}') + + volume = data['cell'][0] @(data['cell'][1].cross(data['cell'][2])) + if volume == 0: + log.warning('Volume is 0, please check the cell parameters. \nFor 3D bulk materials, the volume should be positive. but for 1D,2D or non-periodic systems, the volume could be 0.') + volume = 1.0 + if volume < 0: + log.warning(f'Volume is negative {volume}, please check the cell parameters. We will take the absolute value of the volume.') + volume = - volume + prefactor = g_s * 1j / (volume.numpy()) + ac_cond = ac_cond * prefactor + + return omegas, ac_cond diff --git a/dptb/postprocess/fortran/ac_cond.f90 b/dptb/postprocess/fortran/ac_cond.f90 new file mode 100644 index 00000000..8d46ec0a --- /dev/null +++ b/dptb/postprocess/fortran/ac_cond.f90 @@ -0,0 +1,84 @@ + +! calculate full ac conductivity using Kubo-Greenwood formula +! this ac_cond_f is originaly from tbplas code. +subroutine ac_cond_f(delta_eng, num_orb, num_kpt, prod_df, & + omegas, num_omega, delta, & + k_min, k_max, ac_cond) +implicit none + +! input and output +real(kind=8), intent(in) :: delta_eng(num_orb, num_orb, num_kpt) +integer, intent(in) :: num_orb, num_kpt +complex(kind=8), intent(in) :: prod_df(num_orb, num_orb, num_kpt) +real(kind=8), intent(in) :: omegas(num_omega) +integer, intent(in) :: num_omega +real(kind=8), intent(in) :: delta +integer, intent(in) :: k_min, k_max +complex(kind=8), intent(inout) :: ac_cond(num_omega) + +! local variables +integer :: i_w, i_k, mm, nn +real(kind=8) :: omega +complex(kind=8) :: cdelta, ac_sum + +! calculate ac_cond +cdelta = dcmplx(0.0D0, delta) +!$OMP PARALLEL DO PRIVATE(omega, ac_sum, i_k, mm, nn) +do i_w = 1, num_omega +omega = omegas(i_w) +ac_sum = dcmplx(0.0D0, 0.0D0) +do i_k = k_min, k_max +do mm = 1, num_orb +do nn = 1, num_orb + ac_sum = ac_sum + prod_df(nn, mm, i_k) & + / (delta_eng(nn, mm, i_k) - omega - cdelta) +end do +end do +end do +ac_cond(i_w) = ac_sum +end do +!$OMP END PARALLEL DO +end subroutine ac_cond_f + + +subroutine ac_cond_gauss(delta_eng, num_orb, num_kpt, prod_df, & + omegas, num_omega, delta, & + k_min, k_max, ac_cond) +implicit none + +! input and output +real(kind=8), intent(in) :: delta_eng(num_orb, num_orb, num_kpt) +integer, intent(in) :: num_orb, num_kpt +complex(kind=8), intent(in) :: prod_df(num_orb, num_orb, num_kpt) +real(kind=8), intent(in) :: omegas(num_omega) +integer, intent(in) :: num_omega +real(kind=8), intent(in) :: delta +integer, intent(in) :: k_min, k_max +complex(kind=8), intent(inout) :: ac_cond(num_omega) + +real(kind=8), parameter :: PI = 3.141592653589793d0 + +! local variables +integer :: i_w, i_k, mm, nn +real(kind=8) :: omega +complex(kind=8) :: cdelta, ac_sum + +! calculate ac_cond +cdelta = dcmplx(0.0D0, delta) +!$OMP PARALLEL DO PRIVATE(omega, ac_sum, i_k, mm, nn) +do i_w = 1, num_omega +omega = omegas(i_w) +ac_sum = dcmplx(0.0D0, 0.0D0) +do i_k = k_min, k_max +do mm = 1, num_orb +do nn = 1, num_orb + ac_sum = ac_sum + prod_df(nn, mm, i_k) & + * exp(-0.5 * (( delta_eng(nn, mm, i_k) - omega) / delta)**2) & + / (delta * sqrt(2 * PI)) +end do +end do +end do +ac_cond(i_w) = ac_sum +end do +!$OMP END PARALLEL DO +end subroutine ac_cond_gauss \ No newline at end of file diff --git a/dptb/postprocess/optical/optical_cond.py b/dptb/postprocess/optical/optical_cond.py new file mode 100644 index 00000000..f495a667 --- /dev/null +++ b/dptb/postprocess/optical/optical_cond.py @@ -0,0 +1,217 @@ + +import numpy as np +from dptb.data import AtomicDataDict +from dptb.data import AtomicData +import sys +import torch +from dptb.nn.hr2hk import HR2HK +from dptb.nn.hr2dhk import Hr2dHk +import math +from dptb.utils.make_kpoints import kmesh_sampling_negf +import time +import logging +import os +from ase.io import read +import matplotlib.pyplot as plt +from dptb.utils.constants import atomic_num_dict_r +import scipy + +try: + from dptb.postprocess.fortran import ac_cond as acdf2py +except ImportError: + acdf2py = None + +log = logging.getLogger(__name__) + +def fermi_dirac(e, mu, beta): + return 1/(1+torch.exp(beta*(e-mu))) + +def gauss(x,mu,sigma): + res = torch.exp(-0.5*((x-mu)/sigma)**2)/(sigma*torch.sqrt(2*torch.tensor(math.pi))) + return res + +class AcCond: + def __init__(self, model:torch.nn.Module, results_path: str=None, use_gui: bool=False, device: str='cpu'): + self.model = model + self.results_path = results_path + self.use_gui = use_gui + self.device = device + os.makedirs(results_path, exist_ok=True) + if acdf2py is None: + log.warning('ac_cond_f is not available, please install the fortran code to calculate the AC conductivity') + sys.exit(1) + + def get_accond(self, + struct, + AtomicData_options, + emax, + num_omega= 1000, + mesh_grid=[1,1,1], + nk_per_loop=None, + delta=0.03, + e_fermi=0, + valence_e=None, + gap_corr=0, + T=300, + direction='xx', + g_s=2): + + self.direction = direction + + log.info('<><><><>'*5) + # 调用from_ase方法，生成一个硅的AtomicData类型数据 + dataset = AtomicData.from_ase(atoms=read(struct),**AtomicData_options) + data = AtomicData.to_AtomicDataDict(dataset) + if valence_e is not None and abs(gap_corr) > 1e-3: + uniq_type, counts = np.unique(data['atomic_numbers'].numpy(), return_counts=True) + tot_num_e = 0 + for i in range(len(uniq_type)): + symbol = atomic_num_dict_r[uniq_type[i]] + assert symbol in valence_e + tot_num_e += counts[i] * valence_e[symbol] + + num_val = tot_num_e // g_s + else: + num_val = None + + self.omegas, self.ac_cond_gauss, self.ac_cond_linhard = AcCond.cal_cond(model=self.model, data=data, e_fermi=e_fermi, mesh_grid=mesh_grid, + emax=emax, num_val=num_val, gap_corr=gap_corr, num_omega=num_omega, nk_per_loop=nk_per_loop, + delta=delta, T=T, direction=direction, g_s=g_s) + + np.save(f"{self.results_path}/AC_{self.direction}_cond_sig_{delta}.npy", {'energy':self.omegas, 'ac_cond_g': self.ac_cond_gauss, 'ac_cond_l': self.ac_cond_linhard}) + log.info('<><><><>'*5) + + def accond_plot(self): + fig = plt.figure(figsize=(6,6),dpi=200) + plt.plot(self.omegas, self.ac_cond_gauss.real, label='Gaussian') + plt.plot(self.omegas, self.ac_cond_linhard.real, label='Linhard:real') + plt.plot(self.omegas, self.ac_cond_linhard.imag, label='Linhard:real') + plt.legend() + plt.xlabel("Energy (eV)") + plt.ylabel(f"sigma_{self.direction}") + plt.savefig(f"{self.results_path}/AC_{self.direction}.png") + if self.use_gui: + plt.show() + plt.close() + + @staticmethod + def cal_cond(model, data, e_fermi, mesh_grid, emax, num_val=None, gap_corr=0, num_omega= 1000, nk_per_loop=None, delta=0.005, T=300, direction='xx', g_s=2): + h2k = HR2HK( + idp=model.idp, + device=model.device, + dtype=model.dtype) + + h2dk = Hr2dHk( + idp=model.idp, + device=model.device, + dtype=model.dtype) + data = model.idp(data) + data = model(data) + + log.info('application of the model is done') + + KB = 8.617333262e-5 + beta = 1/(KB*T) + kpoints, kweight = kmesh_sampling_negf(mesh_grid) + assert len(direction) == 2 + kpoints = torch.as_tensor(kpoints, dtype=torch.float32) + kweight = torch.as_tensor(kweight, dtype=torch.float64) + assert kpoints.shape[0] == kweight.shape[0] + + tot_numk = kpoints.shape[0] + if nk_per_loop is None: + log.warning('nk_per_loop is not set, will use all kpoints in one loop, which may cause memory error.') + nk_per_loop = tot_numk + num_loop = math.ceil(tot_numk / nk_per_loop) + omegas = torch.linspace(0,emax,num_omega, dtype=torch.float64) + + log.info(f'tot_numk: {tot_numk}, nk_per_loop: {nk_per_loop}, num_loop: {num_loop}') + + ac_cond = np.zeros((len(omegas)),dtype=np.complex128) + ac_cond_ik = np.zeros((len(omegas)),dtype=np.complex128) + + ac_cond_linhard = np.zeros((len(omegas)),dtype=np.complex128) + ac_cond_linhard_ik = np.zeros((len(omegas)),dtype=np.complex128) + + for ik in range(num_loop): + t_start = time.time() + log.info('<><><><><'*5) + log.info(f'loop {ik+1} in {num_loop} circles') + istart = ik * nk_per_loop + iend = min((ik + 1) * nk_per_loop, tot_numk) + kpoints_ = kpoints[istart:iend] + kweight_ = kweight[istart:iend] + + data[AtomicDataDict.KPOINT_KEY] = torch.nested.as_nested_tensor([kpoints_]) + data = h2k(data) + dhdk = h2dk(data,direction=direction) + + # Hamiltonian = data['hamiltonian'].detach().to(torch.complex128) + # dhdk = {k: v.detach().to(torch.complex128) for k, v in dhdk.items()} + + log.info(f' - get H and dHdk ...') + if data['hamiltonian'].shape[0] == 1: + log.info(f' - use scipy for diagonalization of H ...') + eigs, eigv = scipy.linalg.eigh(data['hamiltonian'].detach().numpy()[0]) + eigs = eigs[None,:] + eigv = eigv[None,:,:] + eigs = torch.as_tensor(eigs, dtype=torch.float32) + eigv = torch.as_tensor(eigv, dtype=torch.complex64) + else: + eigs, eigv = torch.linalg.eigh(data['hamiltonian']) + + if num_val is not None and abs(gap_corr) > 1e-3: + log.info(f' - gap correction is applied with {gap_corr}') + assert num_val > 0 + assert eigs[:,num_val].min() - eigs[:,num_val-1].max() > 1e-3 , f'the gap between the VBM {num_val-1} and the CBM {num_val} is too small' + + eigs[:,:num_val] = eigs[:,:num_val] - gap_corr/2 + eigs[:,num_val:] = eigs[:,num_val:] + gap_corr/2 + + log.info(f' - diagonalization of H ...') + + dh1 = eigv.conj().transpose(1,2) @ dhdk[direction[0]] @ eigv + if direction[0] == direction[1]: + dh2 = dh1 + else: + dh2 = eigv.conj().transpose(1,2) @ dhdk[direction[1]] @ eigv + + p1p2 = dh1 * dh2.transpose(1,2) + + + log.info(f' - get p matrix from dHdk ...') + + p1p2.to(torch.complex128) + eigs.to(torch.float64) + + eig_diff = eigs[:,:,None] - eigs[:,None,:] + + fdv = fermi_dirac(eigs, e_fermi, beta) + fd_diff = fdv[:,:,None] - fdv[:,None,:] + #fd_ed = torch.zeros_like(eig_diff) + ind = torch.abs(eig_diff) > 1e-6 + ind2 = torch.abs(eig_diff) <= 1e-6 + fd_diff[ind] = fd_diff[ind] / eig_diff[ind] + fd_diff[ind2] = 0.0 + + p1p2 = p1p2 * fd_diff + p1p2 = p1p2 * kweight_[:,None,None] + + kpoints_.shape[1] + ac_cond_ik = ac_cond_ik * 0 + acdf2py.ac_cond_gauss(eig_diff.permute(1,2,0).detach().numpy(), p1p2.permute(1,2,0).detach().numpy(), omegas, delta, 1, kpoints_.shape[0], ac_cond_ik) + acdf2py.ac_cond_f(eig_diff.permute(1,2,0).detach().numpy(), p1p2.permute(1,2,0).detach().numpy(), omegas, delta, 1, kpoints_.shape[0], ac_cond_linhard_ik) + ac_cond = ac_cond + ac_cond_ik + ac_cond_linhard = ac_cond_linhard + ac_cond_linhard_ik + + log.info(f' - get ac_cond ...') + t_end = time.time() + log.info(f'time cost: {t_end-t_start:.4f} s in loop {ik+1}') + + ac_cond = 1.0j * ac_cond * np.pi + volume = data['cell'][0] @ (data['cell'][1].cross(data['cell'][2],dim=0)) + prefactor = 2 * g_s * 1j / (volume.numpy()) + ac_cond = ac_cond * prefactor + ac_cond_linhard = ac_cond_linhard * prefactor + + return omegas, ac_cond, ac_cond_linhard diff --git a/dptb/utils/argcheck.py b/dptb/utils/argcheck.py index ceee4d2a..575a4068 100644 --- a/dptb/utils/argcheck.py +++ b/dptb/utils/argcheck.py @@ -1113,8 +1113,38 @@ def task_options(): Argument("negf", dict, negf()), Argument("tbtrans_negf", dict, tbtrans_negf()), Argument("write_block", dict, write_block), + Argument("ac_cond", dict, ac_cond()), ],optional=False, doc=doc_task) +def ac_cond(): + doc_emax = "" + doc_num_omega = "" + doc_mesh_grid = "" + doc_nk_per_loop = "" + doc_delta = "" + doc_e_fermi = "" + doc_valence_e = "" + doc_gap_corr = "" + doc_T = "" + doc_direction = "" + doc_g_s = "" + + argu = [ + Argument("emax", float, optional=False, default=10, doc=doc_emax), + Argument("num_omega", int, optional=False, default=1000, doc=doc_num_omega), + Argument("mesh_grid", list, optional=False, default=[1,1,1], doc=doc_mesh_grid), + Argument("nk_per_loop", [int, None], optional=True, default=None, doc=doc_nk_per_loop), + Argument("delta", float, optional=False, default=0.03, doc=doc_delta), + Argument("e_fermi", [float, int, None], optional=False, doc=doc_e_fermi), + Argument("valence_e", [dict, None], optional=True, default=None, doc=doc_valence_e), + Argument("gap_corr", float, optional=False, default=0, doc=doc_gap_corr), + Argument("T", [float, int], optional=False, default=300, doc=doc_T), + Argument("direction", str, optional=False, default="xx", doc=doc_direction), + Argument("g_s", int, optional=False, default=2, doc=doc_g_s) + ] + + return argu + def band(): doc_kline_type ="""The different type to build kpath line mode. - "abacus" : the abacus format diff --git a/pyproject.toml b/pyproject.toml index da99df5f..7b86ffe5 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,97 +1,72 @@ -[tool.poetry] +[build-system] +requires = ["scikit-build-core", "numpy", "setuptools>=45", "setuptools-scm[toml]>=6.2"] +build-backend = "scikit_build_core.build" + +[tool.setuptools_scm] +write_to = "dptb/_version.py" +version_scheme = "no-guess-dev" +local_scheme = "dirty-tag" +fallback_version = "0.0.0" + +[project] name = "dptb" -version = "2.0.1" -license = "LGPL-3.0" +dynamic = ["version"] description = "A deep learning package for emperical tight-binding approach with first-principle accuracy." -authors = ["Q. Gu ", "Z. Zhanghao "] +authors = [ + {name = "Q. Gu", email = "guqq@pku.edu.cn"}, + {name = "Z. Zhanghao", email = "zhouyinzhanghao@gmail.com"} +] readme = "README.md" -repository = "https://github.com/deepmodeling/DeePTB" - -[tool.poetry.dependencies] -python = ">=3.8" -pytest = ">=7.2.0" -pytest-order = "1.2.0" -numpy = "*" -scipy = "1.9.1" -spglib = "*" -matplotlib = "*" -torch = ">=1.13.0" -ase = "*" -pyyaml = "*" -future = "*" -dargs = "0.4.4" -xitorch = "0.3.0" -fmm3dpy = "1.0.0" -e3nn = ">=0.5.1" -torch-runstats = "0.2.0" -torch_scatter = "2.1.2" -torch_geometric = ">=2.4.0" -opt-einsum = "3.3.0" -h5py = "3.7.0" -lmdb = "1.4.1" - +license = {text = "LGPL-3.0"} +requires-python = ">=3.8" +dependencies = [ + "pytest>=7.2.0", + "pytest-order==1.2.0", + "numpy", + "scipy==1.9.1", + "spglib", + "matplotlib", + "torch>=1.13.0", + "ase", + "pyyaml", + "future", + "dargs==0.4.4", + "xitorch==0.3.0", + "fmm3dpy==1.0.0", + "e3nn>=0.5.1", + "torch-runstats==0.2.0", + "torch_scatter==2.1.2", + "torch_geometric>=2.4.0", + "opt-einsum==3.3.0", + "h5py==3.7.0", + "lmdb==1.4.1" +] -[tool.poetry.group.dev.dependencies] -pytest = ">=7.2.0" -pytest-order = "1.2.0" -numpy = "*" -scipy = "1.9.1" -spglib = "*" -matplotlib = "*" -torch = ">=1.13.0" -ase = "*" -pyyaml = "*" -future = "*" -dargs = "0.4.4" -xitorch = "0.3.0" -fmm3dpy = "1.0.0" -e3nn = ">=0.5.1" -torch-runstats = "0.2.0" -torch_scatter = "2.1.2" -torch_geometric = ">=2.4.0" -opt-einsum = "3.3.0" -h5py = "3.7.0" -lmdb = "1.4.1" - -[tool.poetry.group.3Dfermi] -optional = true - -[tool.poetry.group.3Dfermi.dependencies] -ifermi = "*" -pymatgen = "*" - -[tool.poetry.group.tbtrans_init] -optional = true +[project.urls] +repository = "https://github.com/deepmodeling/DeePTB" -[tool.poetry.group.tbtrans_init.dependencies] -sisl = "*" +[project.optional-dependencies] +fortran = [] +"3Dfermi" = ["ifermi", "pymatgen"] +tbtrans_init = ["sisl"] +pybinding = ["pybinding"] -[tool.poetry.group.pybinding] -optional = true +[project.scripts] +dptb = "dptb.__main__:main" +dptb-qm9 = "dptb.data.interfaces.pyscf:main" -[tool.poetry.group.pybinding.dependencies] -pybinding = "*" +[tool.scikit-build] +cmake.minimum-version = "3.17" +cmake.build-type = "Release" +cmake.verbose = true +logging.level = "DEBUG" +experimental = true -[tool.poetry.scripts] -dptb = 'dptb.__main__:main' -dptb-qm9 = 'dptb.data.interfaces.pyscf:main' -[build-system] -requires = ["poetry-core", "poetry-dynamic-versioning"] -build-backend = "poetry_dynamic_versioning.backend" +[tool.scikit-build.metadata.version] +provider = "scikit_build_core.metadata.setuptools_scm" -[tool.poetry-dynamic-versioning] -enable = true -vcs = "git" -strict = true -format-jinja = """ - {%- if distance == 0 -%} - {{ serialize_pep440(base, stage, revision) }} - {%- elif revision is not none -%} - {{ serialize_pep440(base, stage, revision + 1, dev=distance, metadata=[commit]) }} - {%- else -%} - {{ serialize_pep440(bump_version(base), stage, revision, dev=distance, metadata=[commit]) }} - {%- endif -%} -""" +[tool.scikit-build.cmake.define] +BUILD_FORTRAN = {env="BUILD_FORTRAN", default="OFF"} \ No newline at end of file