From cf4869ace12f7ce92c06b7a0df70d18341c5b0b7 Mon Sep 17 00:00:00 2001 From: QG-phy Date: Wed, 10 Jul 2024 17:50:33 +0800 Subject: [PATCH 01/16] add files for cal cond --- dptb/nn/hr2dhk.py | 130 +++++++++++++++++++++ dptb/postprocess/cal_cond.py | 129 ++++++++++++++++++++ dptb/postprocess/fortran/ac_cond_gauss.f90 | 42 +++++++ 3 files changed, 301 insertions(+) create mode 100644 dptb/nn/hr2dhk.py create mode 100644 dptb/postprocess/cal_cond.py create mode 100644 dptb/postprocess/fortran/ac_cond_gauss.f90 diff --git a/dptb/nn/hr2dhk.py b/dptb/nn/hr2dhk.py new file mode 100644 index 00000000..395734e8 --- /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 = 'None', + 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/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_gauss.f90 b/dptb/postprocess/fortran/ac_cond_gauss.f90 new file mode 100644 index 00000000..1d95792d --- /dev/null +++ b/dptb/postprocess/fortran/ac_cond_gauss.f90 @@ -0,0 +1,42 @@ +! calculate full ac conductivity using Kubo-Greenwood formula +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 From 5fe5a060285108e274310c8451796322c317f7b1 Mon Sep 17 00:00:00 2001 From: QG-phy Date: Tue, 23 Jul 2024 22:35:36 +0800 Subject: [PATCH 02/16] update install and ac_cond.f90 --- CMakeLists.txt | 58 +++++++ dptb/__init__.py | 7 +- .../{ac_cond_gauss.f90 => ac_cond.f90} | 42 +++++ pyproject.toml | 156 ++++++++---------- 4 files changed, 174 insertions(+), 89 deletions(-) create mode 100644 CMakeLists.txt rename dptb/postprocess/fortran/{ac_cond_gauss.f90 => ac_cond.f90} (52%) diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 00000000..da1e4f9e --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,58 @@ +cmake_minimum_required(VERSION 3.17...3.26) +project(${SKBUILD_PROJECT_NAME} LANGUAGES C Fortran) + +# 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() + +# Add compiler options for OpenMP +if(USE_OPENMP) + find_package(OpenMP REQUIRED) + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} ${OpenMP_Fortran_FLAGS}") +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) \ No newline at end of file diff --git a/dptb/__init__.py b/dptb/__init__.py index 4f5f7d71..42b0a580 100644 --- a/dptb/__init__.py +++ b/dptb/__init__.py @@ -1,3 +1,8 @@ import importlib.metadata -__version__ = importlib.metadata.version("dptb") \ No newline at end of file +try: + from ._version import version as __version__ +except ImportError: + import warnings + warnings.warn("Failed to find (autogenerated) _version.py. This might be because you are installing from GitHub's tarballs, use the PyPI ones.") + __version__ = "unknown" \ No newline at end of file diff --git a/dptb/postprocess/fortran/ac_cond_gauss.f90 b/dptb/postprocess/fortran/ac_cond.f90 similarity index 52% rename from dptb/postprocess/fortran/ac_cond_gauss.f90 rename to dptb/postprocess/fortran/ac_cond.f90 index 1d95792d..8d46ec0a 100644 --- a/dptb/postprocess/fortran/ac_cond_gauss.f90 +++ b/dptb/postprocess/fortran/ac_cond.f90 @@ -1,4 +1,46 @@ + ! 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) diff --git a/pyproject.toml b/pyproject.toml index da99df5f..97759c41 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,97 +1,77 @@ -[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" - - -[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" +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.3Dfermi] -optional = true - -[tool.poetry.group.3Dfermi.dependencies] -ifermi = "*" -pymatgen = "*" - -[tool.poetry.group.tbtrans_init] -optional = true - -[tool.poetry.group.tbtrans_init.dependencies] -sisl = "*" - -[tool.poetry.group.pybinding] -optional = true - -[tool.poetry.group.pybinding.dependencies] -pybinding = "*" +[project.urls] +repository = "https://github.com/deepmodeling/DeePTB" +[project.optional-dependencies] +fortran = [] +"3Dfermi" = ["ifermi", "pymatgen"] +tbtrans_init = ["sisl"] +pybinding = ["pybinding"] -[tool.poetry.scripts] -dptb = 'dptb.__main__:main' -dptb-qm9 = 'dptb.data.interfaces.pyscf:main' +[project.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] +cmake.minimum-version = "3.17" +cmake.build-type = "Release" +cmake.verbose = true +logging.level = "DEBUG" +experimental = true +[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] +CMAKE_POSITION_INDEPENDENT_CODE = "ON" +BUILD_FORTRAN = {env="BUILD_FORTRAN", default="OFF"} +USE_OPENMP = {env="USE_OPENMP", default="OFF"} +USE_INTEL = {env="USE_INTEL", default="OFF"} +# FORTRAN_COMPILER = {env="FORTRAN_COMPILER", default=""} +CMAKE_C_COMPILER = "icx" +CMAKE_Fortran_COMPILER = "ifx" +CMAKE_C_FLAGS = "-xHost" +CMAKE_Fortran_FLAGS = "-fpp -xHost -qopenmp -ipo -heap-arrays 32" \ No newline at end of file From 0ed36a4f6d9b493b743c79665b0dd16aee91a012 Mon Sep 17 00:00:00 2001 From: QG-phy Date: Wed, 24 Jul 2024 00:09:25 +0800 Subject: [PATCH 03/16] add version in init --- CMakeLists.txt | 113 ++++++++++++++++++++++++----------------------- dptb/__init__.py | 12 ++--- pyproject.toml | 3 +- 3 files changed, 67 insertions(+), 61 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index da1e4f9e..3aeb0dd1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,58 +1,61 @@ cmake_minimum_required(VERSION 3.17...3.26) project(${SKBUILD_PROJECT_NAME} LANGUAGES C Fortran) -# 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() - -# Add compiler options for OpenMP -if(USE_OPENMP) - find_package(OpenMP REQUIRED) - set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} ${OpenMP_Fortran_FLAGS}") -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) \ No newline at end of file +option(BUILD_FORTRAN "Build Fortran extensions" OFF) +if(BUILD_FORTRAN) + # 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() + + # Add compiler options for OpenMP + if(USE_OPENMP) + find_package(OpenMP REQUIRED) + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} ${OpenMP_Fortran_FLAGS}") + 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) +endif() \ No newline at end of file diff --git a/dptb/__init__.py b/dptb/__init__.py index 42b0a580..34b5cadb 100644 --- a/dptb/__init__.py +++ b/dptb/__init__.py @@ -1,8 +1,10 @@ -import importlib.metadata - try: - from ._version import version as __version__ + import importlib.metadata as importlib_metadata except ImportError: - import warnings - warnings.warn("Failed to find (autogenerated) _version.py. This might be because you are installing from GitHub's tarballs, use the PyPI ones.") + # for Python 3.7 + import importlib_metadata + +try: + __version__ = importlib_metadata.version("dptb") +except importlib_metadata.PackageNotFoundError: __version__ = "unknown" \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index 97759c41..0ad28539 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -65,6 +65,7 @@ experimental = true [tool.scikit-build.metadata.version] provider = "scikit_build_core.metadata.setuptools_scm" + [tool.scikit-build.cmake.define] CMAKE_POSITION_INDEPENDENT_CODE = "ON" BUILD_FORTRAN = {env="BUILD_FORTRAN", default="OFF"} @@ -74,4 +75,4 @@ USE_INTEL = {env="USE_INTEL", default="OFF"} CMAKE_C_COMPILER = "icx" CMAKE_Fortran_COMPILER = "ifx" CMAKE_C_FLAGS = "-xHost" -CMAKE_Fortran_FLAGS = "-fpp -xHost -qopenmp -ipo -heap-arrays 32" \ No newline at end of file +CMAKE_Fortran_FLAGS = "-fpp -xHost -qopenmp -ipo -heap-arrays 32" From 6fe8261524e59307fe48083e67702f91ab244116 Mon Sep 17 00:00:00 2001 From: QG-phy Date: Wed, 24 Jul 2024 00:12:11 +0800 Subject: [PATCH 04/16] add install cmake in dockerfile --- Dockerfile | 1 + Dockerfile.main | 1 + 2 files changed, 2 insertions(+) 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/* && \ From dc9b3441b5299db97373202b4f0ce4b5efcd3bd0 Mon Sep 17 00:00:00 2001 From: QG-phy Date: Wed, 24 Jul 2024 00:18:51 +0800 Subject: [PATCH 05/16] update readme for new installation options. --- README.md | 15 +++++++++++++++ pyproject.toml | 4 ++-- 2 files changed, 17 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 4d74c7b4..84f45881 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 USE_INTEL=ON USE_OPENMP=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. The `USE_INTEL=ON` flag is used to enable the Intel MKL library, and the `USE_OPENMP=ON` flag is used to enable OpenMP parallelization. Both `USE_INTEL` and `USE_OPENMP` are `ON` by default. + ## Usage For a comprehensive guide and usage tutorials, visit [our documentation website](https://deeptb.readthedocs.io/en/latest/). diff --git a/pyproject.toml b/pyproject.toml index 0ad28539..f70b1632 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -69,8 +69,8 @@ provider = "scikit_build_core.metadata.setuptools_scm" [tool.scikit-build.cmake.define] CMAKE_POSITION_INDEPENDENT_CODE = "ON" BUILD_FORTRAN = {env="BUILD_FORTRAN", default="OFF"} -USE_OPENMP = {env="USE_OPENMP", default="OFF"} -USE_INTEL = {env="USE_INTEL", default="OFF"} +USE_OPENMP = {env="USE_OPENMP", default="ON"} +USE_INTEL = {env="USE_INTEL", default="ON"} # FORTRAN_COMPILER = {env="FORTRAN_COMPILER", default=""} CMAKE_C_COMPILER = "icx" CMAKE_Fortran_COMPILER = "ifx" From 1eda0f2f4026f54afa22e98b17195c60d1bbd2a1 Mon Sep 17 00:00:00 2001 From: QG-phy Date: Wed, 24 Jul 2024 15:32:44 +0800 Subject: [PATCH 06/16] add optical module --- dptb/nn/hr2dhk.py | 6 +- dptb/postprocess/optical/optical_cond.py | 132 +++++++++++++++++++++++ 2 files changed, 135 insertions(+), 3 deletions(-) create mode 100644 dptb/postprocess/optical/optical_cond.py diff --git a/dptb/nn/hr2dhk.py b/dptb/nn/hr2dhk.py index 395734e8..9f6dc930 100644 --- a/dptb/nn/hr2dhk.py +++ b/dptb/nn/hr2dhk.py @@ -7,19 +7,19 @@ import re from dptb.utils.tools import float2comlex -class HR2dHk(torch.nn.Module): +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 = 'None', + 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__() + super(Hr2dHk, self).__init__() if isinstance(dtype, str): dtype = getattr(torch, dtype) diff --git a/dptb/postprocess/optical/optical_cond.py b/dptb/postprocess/optical/optical_cond.py new file mode 100644 index 00000000..685d0256 --- /dev/null +++ b/dptb/postprocess/optical/optical_cond.py @@ -0,0 +1,132 @@ + +import numpy as np +from dptb.data import AtomicDataDict +import torch +from dptb.nn.hr2hk import HR2HK +from dptb.nn.hr2dhk import Hr2dHk +from dptb.postprocess.fortran import ac_cond as acdf2py +import math +from dptb.utils.make_kpoints import kmesh_sampling_negf +import time + + +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_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) + + print('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) + + print('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() + print('<><><><><'*5) + print(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()} + + print(f' - get H and dHdk ...') + + eigs, eigv = torch.linalg.eigh(data['hamiltonian']) + + if num_val is not None: + 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' + if abs(gap_corr)> 1e-3: + print(f' - gap correction is applied {gap_corr}') + eigs[:,:num_val] = eigs[:,:num_val] - gap_corr/2 + eigs[:,num_val:] = eigs[:,num_val:] + gap_corr/2 + print(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) + + + print(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 + + print(f' - get ac_cond ...') + t_end = time.time() + print(f'time cost: {t_end-t_start:.4f} s in loop {ik+1}') + + volume = data['cell'][0] @(data['cell'][1].cross(data['cell'][2])) + prefactor = g_s * 1j / (volume.numpy()) + ac_cond = ac_cond * prefactor + ac_cond.imag = -ac_cond.imag*np.pi + ac_cond_linhard = ac_cond_linhard * prefactor + + return omegas, ac_cond, ac_cond_linhard From 935ec75c82f5cf968631765cc0cb6b23a2ed7a18 Mon Sep 17 00:00:00 2001 From: QG-phy Date: Wed, 24 Jul 2024 16:01:23 +0800 Subject: [PATCH 07/16] update compile options for fortran --- pyproject.toml | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index f70b1632..112c9bb3 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -74,5 +74,7 @@ USE_INTEL = {env="USE_INTEL", default="ON"} # FORTRAN_COMPILER = {env="FORTRAN_COMPILER", default=""} CMAKE_C_COMPILER = "icx" CMAKE_Fortran_COMPILER = "ifx" -CMAKE_C_FLAGS = "-xHost" -CMAKE_Fortran_FLAGS = "-fpp -xHost -qopenmp -ipo -heap-arrays 32" +#CMAKE_C_FLAGS = "-xHost" +#CMAKE_Fortran_FLAGS = "-fpp -xHost -qopenmp -ipo -heap-arrays 32" +CMAKE_C_FLAGS = "-O3 -xHost -ipo -fma -align" +CMAKE_Fortran_FLAGS = "-fpp -O3 -xHost -qopenmp -ipo -heap-arrays 32 -unroll -fma -align" From 6cec494e7e35d0e2b7e8087dfaefdf8ac35e104a Mon Sep 17 00:00:00 2001 From: QG-phy Date: Wed, 24 Jul 2024 16:59:14 +0800 Subject: [PATCH 08/16] update optical cond --- dptb/postprocess/optical/optical_cond.py | 37 +++++++++++++----------- pyproject.toml | 3 ++ 2 files changed, 23 insertions(+), 17 deletions(-) diff --git a/dptb/postprocess/optical/optical_cond.py b/dptb/postprocess/optical/optical_cond.py index 685d0256..ef0edd84 100644 --- a/dptb/postprocess/optical/optical_cond.py +++ b/dptb/postprocess/optical/optical_cond.py @@ -8,7 +8,9 @@ import math 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))) @@ -32,7 +34,7 @@ def cal_cond(model, data, e_fermi, mesh_grid, emax, num_val=None, gap_corr=0, nu data = model.idp(data) data = model(data) - print('application of the model is done') + log.info('application of the model is done') KB = 8.617333262e-5 beta = 1/(KB*T) @@ -48,7 +50,7 @@ def cal_cond(model, data, e_fermi, mesh_grid, emax, num_val=None, gap_corr=0, nu num_loop = math.ceil(tot_numk / nk_per_loop) omegas = torch.linspace(0,emax,num_omega, dtype=torch.float64) - print('tot_numk:',tot_numk, 'nk_per_loop:',nk_per_loop, 'num_loop:',num_loop) + 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) @@ -58,8 +60,8 @@ def cal_cond(model, data, e_fermi, mesh_grid, emax, num_val=None, gap_corr=0, nu for ik in range(num_loop): t_start = time.time() - print('<><><><><'*5) - print(f'loop {ik+1} in {num_loop} circles') + 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] @@ -72,18 +74,19 @@ def cal_cond(model, data, e_fermi, mesh_grid, emax, num_val=None, gap_corr=0, nu # Hamiltonian = data['hamiltonian'].detach().to(torch.complex128) # dhdk = {k: v.detach().to(torch.complex128) for k, v in dhdk.items()} - print(f' - get H and dHdk ...') + log.info(f' - get H and dHdk ...') eigs, eigv = torch.linalg.eigh(data['hamiltonian']) - if num_val is not None: + 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' - if abs(gap_corr)> 1e-3: - print(f' - gap correction is applied {gap_corr}') - eigs[:,:num_val] = eigs[:,:num_val] - gap_corr/2 - eigs[:,num_val:] = eigs[:,num_val:] + gap_corr/2 - print(f' - diagonalization of H ...') + + 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]: @@ -94,7 +97,7 @@ def cal_cond(model, data, e_fermi, mesh_grid, emax, num_val=None, gap_corr=0, nu p1p2 = dh1 * dh2.transpose(1,2) - print(f' - get p matrix from dHdk ...') + log.info(f' - get p matrix from dHdk ...') p1p2.to(torch.complex128) eigs.to(torch.float64) @@ -119,14 +122,14 @@ def cal_cond(model, data, e_fermi, mesh_grid, emax, num_val=None, gap_corr=0, nu ac_cond = ac_cond + ac_cond_ik ac_cond_linhard = ac_cond_linhard + ac_cond_linhard_ik - print(f' - get ac_cond ...') + log.info(f' - get ac_cond ...') t_end = time.time() - print(f'time cost: {t_end-t_start:.4f} s in loop {ik+1}') + 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])) - prefactor = g_s * 1j / (volume.numpy()) + 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.imag = -ac_cond.imag*np.pi ac_cond_linhard = ac_cond_linhard * prefactor return omegas, ac_cond, ac_cond_linhard diff --git a/pyproject.toml b/pyproject.toml index 112c9bb3..1acd71be 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -78,3 +78,6 @@ CMAKE_Fortran_COMPILER = "ifx" #CMAKE_Fortran_FLAGS = "-fpp -xHost -qopenmp -ipo -heap-arrays 32" CMAKE_C_FLAGS = "-O3 -xHost -ipo -fma -align" CMAKE_Fortran_FLAGS = "-fpp -O3 -xHost -qopenmp -ipo -heap-arrays 32 -unroll -fma -align" + +# command to build the fortran code +# BUILD_FORTRAN=ON USE_INTEL=ON USE_OPENMP=ON pip install . From 1b717dc98e7a9ed9be45665e87ee8cfc1ecbd64e Mon Sep 17 00:00:00 2001 From: QG-phy Date: Wed, 24 Jul 2024 23:15:50 +0800 Subject: [PATCH 09/16] update --- .gitignore | 2 +- dptb/entrypoints/run.py | 21 ++ dptb/plugins/__init__.py | 4 + dptb/postprocess/optical/optical_cond.py | 282 ++++++++++++++--------- pyproject.toml | 2 + 5 files changed, 201 insertions(+), 110 deletions(-) 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/dptb/entrypoints/run.py b/dptb/entrypoints/run.py index 0bcb4f5f..f6411ebc 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,26 @@ 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'], + task_options=jdata['task_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.optical_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/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/optical/optical_cond.py b/dptb/postprocess/optical/optical_cond.py index ef0edd84..965c62e0 100644 --- a/dptb/postprocess/optical/optical_cond.py +++ b/dptb/postprocess/optical/optical_cond.py @@ -1,6 +1,8 @@ import numpy as np from dptb.data import AtomicDataDict +from dptb.data import AtomicData + import torch from dptb.nn.hr2hk import HR2HK from dptb.nn.hr2dhk import Hr2dHk @@ -9,6 +11,10 @@ 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 + log = logging.getLogger(__name__) @@ -19,117 +25,175 @@ 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_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: - 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) - - 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 ...') - - eigs, eigv = torch.linalg.eigh(data['hamiltonian']) +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) + + 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): - 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 + log.info('<><><><>'*5) + # 调用from_ase方法，生成一个硅的AtomicData类型数据 + dataset = AtomicData.from_ase(atoms=read(struct),**AtomicData_options) + data = AtomicData.to_AtomicDataDict(dataset) + if valence_e is None and abs(gap_corr) > 1e-3: + uniq_type, counts = np.unique(data['atom_types'].numpy(), return_counts=True) + tot_num_e = 0 + for i in range(len(uniq_type)): + symbol = self.model.idp.type_to_chemical_symbol[uniq_type[i]] + assert symbol in valence_e + tot_num_e += counts[i] * valence_e[symbol] + + num_val = tot_num_e // g_s 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] + 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) - 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 + np.save(f"{self.results_path}/AC_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 plot_optcond(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.xlabel("Energy (eV)") + plt.ylabel("sigma_xx") + plt.savefig("ACxx.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('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 ...') + + 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}') - 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 - 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 + return omegas, ac_cond, ac_cond_linhard diff --git a/pyproject.toml b/pyproject.toml index 1acd71be..c64e51f0 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -62,6 +62,8 @@ cmake.verbose = true logging.level = "DEBUG" experimental = true + + [tool.scikit-build.metadata.version] provider = "scikit_build_core.metadata.setuptools_scm" From 7d9ea9bed8fa3385c35f26192ec9f5fc537aabda Mon Sep 17 00:00:00 2001 From: QG-phy Date: Thu, 25 Jul 2024 00:29:23 +0800 Subject: [PATCH 10/16] ad argcheck for accond, fix bug in run.py and refactor in optical_cond.py --- dptb/entrypoints/run.py | 3 +-- dptb/postprocess/optical/optical_cond.py | 5 ++-- dptb/utils/argcheck.py | 30 ++++++++++++++++++++++++ 3 files changed, 34 insertions(+), 4 deletions(-) diff --git a/dptb/entrypoints/run.py b/dptb/entrypoints/run.py index f6411ebc..e64b3b3b 100644 --- a/dptb/entrypoints/run.py +++ b/dptb/entrypoints/run.py @@ -91,7 +91,6 @@ def run( accondcal.get_accond(struct=struct_file, AtomicData_options=jdata['AtomicData_options'], - task_options=jdata['task_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]), @@ -104,7 +103,7 @@ def run( direction=jdata['task_options'].get('direction','xx'), g_s=jdata['task_options'].get('g_s',2) ) - accondcal.optical_plot() + accondcal.accond_plot() log.info(msg='ac optical conductivity calculation successfully completed.') elif task=='write_block': diff --git a/dptb/postprocess/optical/optical_cond.py b/dptb/postprocess/optical/optical_cond.py index 965c62e0..9744d963 100644 --- a/dptb/postprocess/optical/optical_cond.py +++ b/dptb/postprocess/optical/optical_cond.py @@ -71,11 +71,12 @@ def get_accond(self, np.save(f"{self.results_path}/AC_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 plot_optcond(self): + 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("sigma_xx") plt.savefig("ACxx.png") @@ -114,7 +115,7 @@ def cal_cond(model, data, e_fermi, mesh_grid, emax, num_val=None, gap_corr=0, nu 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) + 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) diff --git a/dptb/utils/argcheck.py b/dptb/utils/argcheck.py index ceee4d2a..6da6ae6d 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", [float, int, 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 From eb079edcb85d83c72d758cdd465aeed9742f7e04 Mon Sep 17 00:00:00 2001 From: qqgu Date: Thu, 25 Jul 2024 14:24:54 +0800 Subject: [PATCH 11/16] chore: Update CMakeLists.txt to support Intel compilers and OpenMP --- CMakeLists.txt | 30 +++++++++++++++++++++++------- pyproject.toml | 15 +-------------- 2 files changed, 24 insertions(+), 21 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 3aeb0dd1..c5423a38 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,8 +1,27 @@ cmake_minimum_required(VERSION 3.17...3.26) -project(${SKBUILD_PROJECT_NAME} LANGUAGES C Fortran) option(BUILD_FORTRAN "Build Fortran extensions" OFF) +option(USE_INTEL "Use Intel compilers" ON) +option(USE_OPENMP "Use OpenMP" ON) + 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) @@ -16,12 +35,6 @@ if(BUILD_FORTRAN) message(FATAL_ERROR "Unsupported Fortran compiler ${CMAKE_Fortran_COMPILER}") endif() - # Add compiler options for OpenMP - if(USE_OPENMP) - find_package(OpenMP REQUIRED) - set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} ${OpenMP_Fortran_FLAGS}") - endif() - #----------------------------- Fortran extension ------------------------------# # Define the fortranobject execute_process( @@ -58,4 +71,7 @@ if(BUILD_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/pyproject.toml b/pyproject.toml index c64e51f0..7b86ffe5 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -69,17 +69,4 @@ provider = "scikit_build_core.metadata.setuptools_scm" [tool.scikit-build.cmake.define] -CMAKE_POSITION_INDEPENDENT_CODE = "ON" -BUILD_FORTRAN = {env="BUILD_FORTRAN", default="OFF"} -USE_OPENMP = {env="USE_OPENMP", default="ON"} -USE_INTEL = {env="USE_INTEL", default="ON"} -# FORTRAN_COMPILER = {env="FORTRAN_COMPILER", default=""} -CMAKE_C_COMPILER = "icx" -CMAKE_Fortran_COMPILER = "ifx" -#CMAKE_C_FLAGS = "-xHost" -#CMAKE_Fortran_FLAGS = "-fpp -xHost -qopenmp -ipo -heap-arrays 32" -CMAKE_C_FLAGS = "-O3 -xHost -ipo -fma -align" -CMAKE_Fortran_FLAGS = "-fpp -O3 -xHost -qopenmp -ipo -heap-arrays 32 -unroll -fma -align" - -# command to build the fortran code -# BUILD_FORTRAN=ON USE_INTEL=ON USE_OPENMP=ON pip install . +BUILD_FORTRAN = {env="BUILD_FORTRAN", default="OFF"} \ No newline at end of file From 4071bcec85941e8db4ae4ebc98a3b609ad985de1 Mon Sep 17 00:00:00 2001 From: QG-phy Date: Thu, 25 Jul 2024 17:06:53 +0800 Subject: [PATCH 12/16] update readme --- README.md | 6 +++--- dptb/postprocess/optical/optical_cond.py | 10 ++++++++-- 2 files changed, 11 insertions(+), 5 deletions(-) diff --git a/README.md b/README.md index 84f45881..024cc273 100644 --- a/README.md +++ b/README.md @@ -56,10 +56,10 @@ Installing **DeePTB** is straightforward. We recommend using a virtual environme 3. Navigate to the root directory and install DeePTB: ```bash cd DeePTB - BUILD_FORTRAN=ON USE_INTEL=ON USE_OPENMP=ON pip install . + 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. The `USE_INTEL=ON` flag is used to enable the Intel MKL library, and the `USE_OPENMP=ON` flag is used to enable OpenMP parallelization. Both `USE_INTEL` and `USE_OPENMP` are `ON` by default. - + 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/postprocess/optical/optical_cond.py b/dptb/postprocess/optical/optical_cond.py index 9744d963..b2efa0bb 100644 --- a/dptb/postprocess/optical/optical_cond.py +++ b/dptb/postprocess/optical/optical_cond.py @@ -2,11 +2,10 @@ 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 -from dptb.postprocess.fortran import ac_cond as acdf2py import math from dptb.utils.make_kpoints import kmesh_sampling_negf import time @@ -15,6 +14,10 @@ from ase.io import read import matplotlib.pyplot as plt +try: + from dptb.postprocess.fortran import ac_cond as acdf2py +except ImportError: + acdf2py = None log = logging.getLogger(__name__) @@ -32,6 +35,9 @@ def __init__(self, model:torch.nn.Module, results_path: str=None, use_gui: bool= 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, From 74e7672f5f53643069abea0792e9c1be553be91c Mon Sep 17 00:00:00 2001 From: QG-phy Date: Fri, 26 Jul 2024 15:59:03 +0800 Subject: [PATCH 13/16] update the compile Flags in Cmakelist. --- CMakeLists.txt | 6 +++++- dptb/postprocess/optical/optical_cond.py | 8 +++++--- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index c5423a38..51ecf645 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -4,6 +4,10 @@ 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") @@ -15,7 +19,7 @@ if(BUILD_FORTRAN) project(${SKBUILD_PROJECT_NAME} LANGUAGES C Fortran) # 设置编译器标志 - set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -O3 -fPIC") + #set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -O3 -fPIC") if(USE_OPENMP) find_package(OpenMP REQUIRED) diff --git a/dptb/postprocess/optical/optical_cond.py b/dptb/postprocess/optical/optical_cond.py index b2efa0bb..53d5c5bc 100644 --- a/dptb/postprocess/optical/optical_cond.py +++ b/dptb/postprocess/optical/optical_cond.py @@ -54,6 +54,8 @@ def get_accond(self, direction='xx', g_s=2): + self.direction = direction + log.info('<><><><>'*5) # 调用from_ase方法，生成一个硅的AtomicData类型数据 dataset = AtomicData.from_ase(atoms=read(struct),**AtomicData_options) @@ -74,7 +76,7 @@ def get_accond(self, 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_cond_sig_{delta}.npy", {'energy':self.omegas, 'ac_cond_g': self.ac_cond_gauss, 'ac_cond_l': self.ac_cond_linhard}) + 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): @@ -84,8 +86,8 @@ def accond_plot(self): plt.plot(self.omegas, self.ac_cond_linhard.imag, label='Linhard:real') plt.legend() plt.xlabel("Energy (eV)") - plt.ylabel("sigma_xx") - plt.savefig("ACxx.png") + plt.ylabel(f"sigma_{self.direction}") + plt.savefig(f"{self.results_path}/AC_{self.direction}.png") if self.use_gui: plt.show() plt.close() From d576975653f0104407852e3a7c6c862833e5acd6 Mon Sep 17 00:00:00 2001 From: QG-phy Date: Fri, 26 Jul 2024 17:22:31 +0800 Subject: [PATCH 14/16] fix bugs in argcheck of valence_e and the bug for add scissor correction. --- dptb/postprocess/optical/optical_cond.py | 8 ++++---- dptb/utils/argcheck.py | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/dptb/postprocess/optical/optical_cond.py b/dptb/postprocess/optical/optical_cond.py index 53d5c5bc..2599478d 100644 --- a/dptb/postprocess/optical/optical_cond.py +++ b/dptb/postprocess/optical/optical_cond.py @@ -13,7 +13,7 @@ import os from ase.io import read import matplotlib.pyplot as plt - +from dptb.utils.constants import atomic_num_dict_r try: from dptb.postprocess.fortran import ac_cond as acdf2py except ImportError: @@ -60,11 +60,11 @@ def get_accond(self, # 调用from_ase方法，生成一个硅的AtomicData类型数据 dataset = AtomicData.from_ase(atoms=read(struct),**AtomicData_options) data = AtomicData.to_AtomicDataDict(dataset) - if valence_e is None and abs(gap_corr) > 1e-3: - uniq_type, counts = np.unique(data['atom_types'].numpy(), return_counts=True) + 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 = self.model.idp.type_to_chemical_symbol[uniq_type[i]] + symbol = atomic_num_dict_r[uniq_type[i]] assert symbol in valence_e tot_num_e += counts[i] * valence_e[symbol] diff --git a/dptb/utils/argcheck.py b/dptb/utils/argcheck.py index 6da6ae6d..575a4068 100644 --- a/dptb/utils/argcheck.py +++ b/dptb/utils/argcheck.py @@ -1136,7 +1136,7 @@ def ac_cond(): 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", [float, int, None], optional=True, default=None, doc=doc_valence_e), + 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), From dc1cb1f451eb8eb76d3ba0d36b0bc8c1b4e444f1 Mon Sep 17 00:00:00 2001 From: QG-phy Date: Fri, 26 Jul 2024 23:50:04 +0800 Subject: [PATCH 15/16] add scipy to diag. for the internal error of Intel MKL ERROR: Parameter 10 was incorrect on entry to CSTEDC. --- dptb/postprocess/optical/optical_cond.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/dptb/postprocess/optical/optical_cond.py b/dptb/postprocess/optical/optical_cond.py index 2599478d..3ad21827 100644 --- a/dptb/postprocess/optical/optical_cond.py +++ b/dptb/postprocess/optical/optical_cond.py @@ -14,6 +14,8 @@ 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: @@ -148,8 +150,14 @@ def cal_cond(model, data, e_fermi, mesh_grid, emax, num_val=None, gap_corr=0, nu # 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']) + if data['hamiltonian'].shape[0] == 1: + 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}') From 511e18068ef20ad7fbc51b18528ba8100178902a Mon Sep 17 00:00:00 2001 From: QG-phy Date: Fri, 26 Jul 2024 23:54:47 +0800 Subject: [PATCH 16/16] update diag --- dptb/postprocess/optical/optical_cond.py | 1 + 1 file changed, 1 insertion(+) diff --git a/dptb/postprocess/optical/optical_cond.py b/dptb/postprocess/optical/optical_cond.py index 3ad21827..f495a667 100644 --- a/dptb/postprocess/optical/optical_cond.py +++ b/dptb/postprocess/optical/optical_cond.py @@ -151,6 +151,7 @@ def cal_cond(model, data, e_fermi, mesh_grid, emax, num_val=None, gap_corr=0, nu 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,:,:]