Add finite-volume Kubo transport property (E. Prodan) (#5209)

Co-authored-by: Thomas D. Kuehne <tkuehne@cp2k.org>
This commit is contained in:
Dynamics of Condensed Matter 2026-05-15 20:38:26 +02:00 committed by GitHub
parent 7f57f22fd5
commit d43dca4c75
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
11 changed files with 983 additions and 1 deletions

View file

@ -627,6 +627,7 @@ list(
qs_ks_types.F
qs_ks_utils.F
qs_ks_reference.F
qs_kubo_transport.F
qs_linres_atom_current.F
qs_linres_current.F
qs_linres_current_utils.F

View file

@ -94,7 +94,7 @@ MODULE bibliography
VazdaCruz2021, Chen2025, Hernandez2025, Marek2025, Hehn2022, Hehn2024, Pasquier2025, &
Hanasaki2025, Tan2025, Broyden1965, Johnson1988, Kerker1981, &
FuHo1983, MethfesselPaxton1989, Marzari1999, dosSantos2023, Mermin1965, &
Schreder2021, Schreder2024_1, Schreder2024_2
KuhneHeskeProdan2020, Schreder2021, Schreder2024_1, Schreder2024_2
CONTAINS
@ -1660,6 +1660,12 @@ CONTAINS
source="J. Chem. Phys.", volume="152", pages="194103", &
year=2020, doi="10.1063/5.0007045", citation_key="Kühne2020")
CALL add_reference(key=KuhneHeskeProdan2020, &
authors=s2a("T. D. Kühne", "J. Heske", "E. Prodan"), &
title="Disordered crystals from first principles. II: Transport coefficients", &
source="Ann. Phys.", volume="421", pages="168290", &
year=2020, doi="10.1016/j.aop.2020.168290", citation_key="Kühne2020a")
CALL add_reference(key=Iannuzzi2026, &
authors=s2a("M. Iannuzzi", "J. Wilhelm", "F. Stein", "A. Bussy", &
"H. Elgabarty", "D. Golze", "A. Hehn", "M. Graml", &

View file

@ -17,6 +17,7 @@ MODULE input_cp2k_properties_dft
Hernandez2025, &
Iannuzzi2005, &
Kondov2007, &
KuhneHeskeProdan2020, &
Luber2014, &
Putrino2000, &
Putrino2002, &
@ -138,6 +139,10 @@ CONTAINS
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
CALL create_kubo_transport_section(subsection)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
CALL create_bandstructure_section(subsection)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
@ -148,6 +153,90 @@ CONTAINS
END SUBROUTINE create_properties_section
! **************************************************************************************************
!> \brief creates the input structure used to activate finite-volume Kubo transport
! **************************************************************************************************
SUBROUTINE create_kubo_transport_section(section)
TYPE(section_type), POINTER :: section
TYPE(keyword_type), POINTER :: keyword
CPASSERT(.NOT. ASSOCIATED(section))
NULLIFY (keyword)
CALL section_create(section, __LOCATION__, name="KUBO_TRANSPORT", &
description="Finite-volume Kubo-Greenwood transport coefficients from the "// &
"converged Quickstep Hamiltonian, overlap matrix, and atomic geometry.", &
n_keywords=8, n_subsections=0, repeats=.FALSE., &
citations=[KuhneHeskeProdan2020])
CALL keyword_create(keyword, __LOCATION__, name="_SECTION_PARAMETERS_", &
description="Controls the activation of the Kubo transport calculation.", &
default_l_val=.FALSE., &
lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="METHOD", &
description="Transport algorithm. DIAGONALIZATION evaluates the finite-volume "// &
"Kubo-Greenwood expression from the explicit spectrum. TD and CHEBYSHEV are "// &
"reserved for future diagonalization-free implementations.", &
usage="METHOD DIAGONALIZATION", &
n_var=1, type_of_var=char_t, default_c_val="DIAGONALIZATION")
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="TEMPERATURE", &
description="Electronic temperature used in the Fermi operator.", &
usage="TEMPERATURE 1.0", &
default_r_val=cp_unit_to_cp2k(value=1.0_dp, unit_str="K"), &
n_var=1, type_of_var=real_t, unit_str="K")
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="DISSIPATION", &
description="Dissipation/broadening parameter in the finite-temperature "// &
"Kubo formula.", &
usage="DISSIPATION 300.0", &
default_r_val=cp_unit_to_cp2k(value=300.0_dp, unit_str="K"), &
n_var=1, type_of_var=real_t, unit_str="K")
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="ENERGY_RANGE", &
description="Absolute chemical-potential range. The default 0 0 uses the "// &
"full eigenvalue range of the finite-volume Hamiltonian.", &
usage="ENERGY_RANGE -0.5 0.5", &
default_r_vals=[0.0_dp, 0.0_dp], n_var=2, type_of_var=real_t, &
unit_str="hartree")
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="NEUTRAL_MU", &
description="Optional fixed neutral chemical potential. If omitted, the "// &
"neutral point is found from the eigenvalue spectrum and total electron count.", &
usage="NEUTRAL_MU 0.245467619658419", &
n_var=1, type_of_var=real_t, unit_str="hartree")
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="N_MU", &
description="Number of chemical-potential grid points.", &
usage="N_MU 200", &
default_i_val=200, n_var=1, type_of_var=integer_t)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="NEUTRAL_GRID", &
description="Number of grid points used to locate the neutral chemical potential.", &
usage="NEUTRAL_GRID 10000", &
default_i_val=10000, n_var=1, type_of_var=integer_t)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
END SUBROUTINE create_kubo_transport_section
! **************************************************************************************************
!> \brief creates the input structure used to activate
!> a resonant inelastic xray scattering (RIXS) calculation

687
src/qs_kubo_transport.F Normal file
View file

@ -0,0 +1,687 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
! !
! SPDX-License-Identifier: GPL-2.0-or-later !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
!> \brief Finite-volume Kubo-Greenwood transport from converged Quickstep matrices.
! **************************************************************************************************
MODULE qs_kubo_transport
USE bibliography, ONLY: KuhneHeskeProdan2020,&
cite_reference
USE cell_types, ONLY: cell_type,&
pbc
USE cp_blacs_env, ONLY: cp_blacs_env_type
USE cp_dbcsr_api, ONLY: dbcsr_get_info,&
dbcsr_p_type,&
dbcsr_type
USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm
USE cp_fm_struct, ONLY: cp_fm_struct_create,&
cp_fm_struct_release,&
cp_fm_struct_type
USE cp_fm_types, ONLY: cp_fm_create,&
cp_fm_get_submatrix,&
cp_fm_release,&
cp_fm_type
USE cp_log_handling, ONLY: cp_get_default_logger,&
cp_logger_get_default_io_unit,&
cp_logger_type
USE input_section_types, ONLY: section_vals_get_subs_vals,&
section_vals_type,&
section_vals_val_get
USE kinds, ONLY: dp
USE mathconstants, ONLY: pi
USE message_passing, ONLY: mp_para_env_type
USE particle_types, ONLY: particle_type
USE physcon, ONLY: a_bohr,&
angstrom,&
e_charge,&
evolt,&
h_bar,&
kelvin
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_mo_types, ONLY: get_mo_set,&
mo_set_type
USE string_utilities, ONLY: uppercase
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_kubo_transport'
PUBLIC :: qs_scf_post_kubo_transport
CONTAINS
! **************************************************************************************************
!> \brief Compute and print the finite-volume Kubo-Greenwood conductivity tensor.
!> \param qs_env Quickstep environment
! **************************************************************************************************
SUBROUTINE qs_scf_post_kubo_transport(qs_env)
TYPE(qs_environment_type), POINTER :: qs_env
CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_kubo_transport'
CHARACTER(LEN=32) :: method
INTEGER :: handle, imu, ispin, nao, natom, &
nelectron_total, neutral_grid, nmu, &
nspins, output_unit
INTEGER, ALLOCATABLE, DIMENSION(:) :: ao_atom
INTEGER, DIMENSION(:), POINTER :: row_blk_sizes
LOGICAL :: allow_mo_reuse, do_kpoints, kubo_active, &
neutral_mu_explicit, ot_active
LOGICAL, ALLOCATABLE, DIMENSION(:) :: reused_mos
REAL(KIND=dp) :: density_factor, dissipation, emax, emin, &
mu, mu0, mu_step, ne, nh, sigma_iso, &
temperature, volume, volume_m3
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: maxocc
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: eig, s_dense
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dHH
REAL(KIND=dp), DIMENSION(3, 3) :: sigma
REAL(KIND=dp), DIMENSION(:), POINTER :: energy_range
TYPE(cell_type), POINTER :: cell
TYPE(cp_blacs_env_type), POINTER :: blacs_env
TYPE(cp_logger_type), POINTER :: logger
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
TYPE(mp_para_env_type), POINTER :: para_env
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(section_vals_type), POINTER :: input_section, kubo_section, ot_section
CALL timeset(routineN, handle)
NULLIFY (blacs_env, cell, energy_range, input_section, kubo_section, logger, matrix_ks, matrix_s, &
mos, ot_section, para_env, particle_set, row_blk_sizes)
CALL get_qs_env(qs_env, input=input_section)
kubo_section => section_vals_get_subs_vals(input_section, "PROPERTIES%KUBO_TRANSPORT")
CALL section_vals_val_get(kubo_section, "_SECTION_PARAMETERS_", l_val=kubo_active)
IF (.NOT. kubo_active) THEN
CALL timestop(handle)
RETURN
END IF
CALL cite_reference(KuhneHeskeProdan2020)
CALL section_vals_val_get(kubo_section, "METHOD", c_val=method)
CALL uppercase(method)
SELECT CASE (TRIM(method))
CASE ("DIAGONALIZATION")
CASE ("TD", "CHEBYSHEV")
CPABORT("KUBO_TRANSPORT METHOD "//TRIM(method)//" is reserved but not implemented yet.")
CASE DEFAULT
CPABORT("Unknown KUBO_TRANSPORT METHOD: "//TRIM(method))
END SELECT
ot_section => section_vals_get_subs_vals(input_section, "DFT%SCF%OT")
CALL section_vals_val_get(ot_section, "_SECTION_PARAMETERS_", l_val=ot_active)
allow_mo_reuse = .NOT. ot_active
CALL get_qs_env(qs_env, blacs_env=blacs_env, cell=cell, do_kpoints=do_kpoints, &
matrix_ks=matrix_ks, matrix_s=matrix_s, mos=mos, natom=natom, &
nelectron_total=nelectron_total, para_env=para_env, particle_set=particle_set)
IF (do_kpoints) CPABORT("KUBO_TRANSPORT currently expects a finite Gamma-point supercell.")
CPASSERT(ASSOCIATED(matrix_s))
CPASSERT(ASSOCIATED(matrix_ks))
CPASSERT(ASSOCIATED(cell))
CPASSERT(ASSOCIATED(particle_set))
CALL section_vals_val_get(kubo_section, "TEMPERATURE", r_val=temperature)
CALL section_vals_val_get(kubo_section, "DISSIPATION", r_val=dissipation)
CALL section_vals_val_get(kubo_section, "ENERGY_RANGE", r_vals=energy_range)
CALL section_vals_val_get(kubo_section, "NEUTRAL_MU", explicit=neutral_mu_explicit)
IF (neutral_mu_explicit) CALL section_vals_val_get(kubo_section, "NEUTRAL_MU", r_val=mu0)
CALL section_vals_val_get(kubo_section, "N_MU", i_val=nmu)
CALL section_vals_val_get(kubo_section, "NEUTRAL_GRID", i_val=neutral_grid)
IF (temperature <= 0.0_dp) CPABORT("KUBO_TRANSPORT TEMPERATURE must be positive.")
IF (dissipation <= 0.0_dp) CPABORT("KUBO_TRANSPORT DISSIPATION must be positive.")
IF (nmu < 1) CPABORT("KUBO_TRANSPORT N_MU must be at least one.")
IF (neutral_grid < 10) CPABORT("KUBO_TRANSPORT NEUTRAL_GRID must be at least ten.")
nspins = SIZE(matrix_ks)
CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nao, row_blk_size=row_blk_sizes)
CALL build_ao_atom_map(row_blk_sizes, natom, ao_atom)
CALL dbcsr_to_dense(matrix_s(1)%matrix, blacs_env, para_env, s_dense)
ALLOCATE (eig(nao, nspins), dHH(nao, nao, 3, nspins), maxocc(nspins), reused_mos(nspins))
DO ispin = 1, nspins
maxocc(ispin) = mos(ispin)%maxocc
CALL setup_spin_transport(matrix_ks(ispin)%matrix, s_dense, blacs_env, para_env, &
particle_set, cell, ao_atom, mos(ispin), allow_mo_reuse, &
eig(:, ispin), dHH(:, :, :, ispin), reused_mos(ispin))
END DO
emin = MINVAL(eig)
emax = MAXVAL(eig)
IF (.NOT. neutral_mu_explicit) THEN
CALL find_neutral_mu(eig, maxocc, temperature, REAL(nelectron_total, KIND=dp), neutral_grid, &
emin, emax, mu0)
END IF
IF (energy_range(1) < energy_range(2)) THEN
emin = energy_range(1)
emax = energy_range(2)
END IF
volume = cell%deth
volume_m3 = volume*a_bohr**3
density_factor = 1.0_dp/volume_m3
logger => cp_get_default_logger()
output_unit = cp_logger_get_default_io_unit(logger)
IF (output_unit > 0) THEN
WRITE (output_unit, '(/,T2,A)') "Kubo-Greenwood finite-volume transport"
WRITE (output_unit, '(T3,A,T38,A)') "Method:", TRIM(method)
IF (ALL(reused_mos)) THEN
WRITE (output_unit, '(T3,A,T38,A)') "Eigensystem source:", "CP2K canonical MOs"
ELSEIF (ANY(reused_mos)) THEN
WRITE (output_unit, '(T3,A,T38,A)') "Eigensystem source:", "mixed canonical/internal"
ELSE
WRITE (output_unit, '(T3,A,T38,A)') "Eigensystem source:", "internal diagonalization"
END IF
WRITE (output_unit, '(T3,A,T38,I12)') "Number of atoms:", natom
WRITE (output_unit, '(T3,A,T38,I12)') "Number of AO functions:", nao
WRITE (output_unit, '(T3,A,T38,F12.4)') "Temperature [K]:", temperature*kelvin
WRITE (output_unit, '(T3,A,T38,F12.4)') "Dissipation [K]:", dissipation*kelvin
WRITE (output_unit, '(T3,A,T38,F12.6)') "Dissipation [eV]:", dissipation*evolt
WRITE (output_unit, '(T3,A,T38,F12.6)') "Neutral chemical potential [eV]:", mu0*evolt
WRITE (output_unit, '(/,T3,A)') "mu[Ha] mu-mu0[eV] Ne[m^-3] Nh[m^-3] "// &
"sigma_iso[S/cm] sigma_xx[S/m] sigma_xy[S/m] sigma_xz[S/m] sigma_yx[S/m] "// &
"sigma_yy[S/m] sigma_yz[S/m] sigma_zx[S/m] sigma_zy[S/m] sigma_zz[S/m]"
END IF
IF (nmu == 1) THEN
mu_step = 0.0_dp
ELSE
mu_step = (emax - emin)/REAL(nmu - 1, KIND=dp)
END IF
DO imu = 1, nmu
mu = emin + REAL(imu - 1, KIND=dp)*mu_step
CALL electron_hole_density(eig, maxocc, mu, mu0, temperature, ne, nh)
CALL conductivity_at_mu(eig, dHH, maxocc, mu, temperature, dissipation, volume, sigma)
IF (output_unit > 0) THEN
sigma_iso = (sigma(1, 1) + sigma(2, 2) + sigma(3, 3))/3.0_dp*1.0E8_dp
WRITE (output_unit, '(T3,F12.6,F14.6,2ES16.7,10ES17.8)') mu, (mu - mu0)*evolt, &
ne*density_factor, nh*density_factor, &
sigma_iso, &
sigma(1, 1)*1.0E10_dp, sigma(1, 2)*1.0E10_dp, sigma(1, 3)*1.0E10_dp, &
sigma(2, 1)*1.0E10_dp, sigma(2, 2)*1.0E10_dp, sigma(2, 3)*1.0E10_dp, &
sigma(3, 1)*1.0E10_dp, sigma(3, 2)*1.0E10_dp, sigma(3, 3)*1.0E10_dp
WRITE (output_unit, '(T3,A,1X,ES17.8)') "KUBO_TRANSPORT| sigma_iso[S/cm]", sigma_iso
END IF
END DO
DEALLOCATE (ao_atom, dHH, eig, maxocc, reused_mos, s_dense)
CALL timestop(handle)
END SUBROUTINE qs_scf_post_kubo_transport
! **************************************************************************************************
!> \brief Build a global AO -> atom map from DBCSR atomic block sizes.
!> \param row_blk_sizes ...
!> \param natom ...
!> \param ao_atom ...
! **************************************************************************************************
SUBROUTINE build_ao_atom_map(row_blk_sizes, natom, ao_atom)
INTEGER, DIMENSION(:), INTENT(IN) :: row_blk_sizes
INTEGER, INTENT(IN) :: natom
INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: ao_atom
INTEGER :: iao, iatom, nao, u
nao = SUM(row_blk_sizes(1:natom))
ALLOCATE (ao_atom(nao))
iao = 0
DO iatom = 1, natom
DO u = 1, row_blk_sizes(iatom)
iao = iao + 1
ao_atom(iao) = iatom
END DO
END DO
END SUBROUTINE build_ao_atom_map
! **************************************************************************************************
!> \brief Copy a DBCSR matrix to a replicated dense matrix.
!> \param matrix ...
!> \param blacs_env ...
!> \param para_env ...
!> \param dense ...
! **************************************************************************************************
SUBROUTINE dbcsr_to_dense(matrix, blacs_env, para_env, dense)
TYPE(dbcsr_type), INTENT(IN) :: matrix
TYPE(cp_blacs_env_type), POINTER :: blacs_env
TYPE(mp_para_env_type), POINTER :: para_env
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
INTENT(OUT) :: dense
INTEGER :: ncol, nrow
TYPE(cp_fm_struct_type), POINTER :: fm_struct
TYPE(cp_fm_type) :: fm
NULLIFY (fm_struct)
CALL dbcsr_get_info(matrix, nfullrows_total=nrow, nfullcols_total=ncol)
ALLOCATE (dense(nrow, ncol))
CALL cp_fm_struct_create(fm_struct, context=blacs_env, para_env=para_env, &
nrow_global=nrow, ncol_global=ncol)
CALL cp_fm_create(fm, fm_struct)
CALL copy_dbcsr_to_fm(matrix, fm)
CALL cp_fm_get_submatrix(fm, dense, n_rows=nrow, n_cols=ncol)
CALL cp_fm_release(fm)
CALL cp_fm_struct_release(fm_struct)
END SUBROUTINE dbcsr_to_dense
! **************************************************************************************************
!> \brief Prepare eigenvalues and the transformed [X,H] matrices for one spin channel.
!> \param ks_matrix ...
!> \param s_dense ...
!> \param blacs_env ...
!> \param para_env ...
!> \param particle_set ...
!> \param cell ...
!> \param ao_atom ...
!> \param mo_set ...
!> \param allow_mo_reuse ...
!> \param eig ...
!> \param dHH ...
!> \param reused_mos ...
! **************************************************************************************************
SUBROUTINE setup_spin_transport(ks_matrix, s_dense, blacs_env, para_env, particle_set, cell, &
ao_atom, mo_set, allow_mo_reuse, eig, dHH, reused_mos)
TYPE(dbcsr_type), INTENT(IN) :: ks_matrix
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: s_dense
TYPE(cp_blacs_env_type), POINTER :: blacs_env
TYPE(mp_para_env_type), POINTER :: para_env
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(cell_type), POINTER :: cell
INTEGER, DIMENSION(:), INTENT(IN) :: ao_atom
TYPE(mo_set_type), INTENT(IN) :: mo_set
LOGICAL, INTENT(IN) :: allow_mo_reuse
REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eig
REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: dHH
LOGICAL, INTENT(OUT) :: reused_mos
INTEGER :: idir, mo_nao, nao, nmo
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: coeff_dense, gmat, h_dense, horth, op, &
op_orth, uvec, work
REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues
TYPE(cp_fm_type), POINTER :: mo_coeff
NULLIFY (mo_coeff, mo_eigenvalues)
nao = SIZE(s_dense, 1)
reused_mos = .FALSE.
CALL dbcsr_to_dense(ks_matrix, blacs_env, para_env, h_dense)
ALLOCATE (op(nao, nao), work(nao, nao))
IF (allow_mo_reuse) THEN
CALL get_mo_set(mo_set=mo_set, nao=mo_nao, nmo=nmo, eigenvalues=mo_eigenvalues, &
mo_coeff=mo_coeff)
reused_mos = mo_nao == nao .AND. nmo >= nao .AND. ASSOCIATED(mo_eigenvalues) .AND. &
ASSOCIATED(mo_coeff)
IF (reused_mos) reused_mos = SIZE(mo_eigenvalues) >= nao
END IF
IF (reused_mos) THEN
ALLOCATE (coeff_dense(nao, nao))
eig(:) = mo_eigenvalues(1:nao)
CALL cp_fm_get_submatrix(mo_coeff, coeff_dense, n_rows=nao, n_cols=nao)
DO idir = 1, 3
CALL build_commutator_kernel(h_dense, particle_set, cell, ao_atom, idir, op)
CALL transform_to_eigenbasis(op, coeff_dense, dHH(:, :, idir), work)
END DO
DEALLOCATE (coeff_dense)
ELSE
ALLOCATE (gmat(nao, nao), horth(nao, nao), op_orth(nao, nao), uvec(nao, nao))
CALL inverse_sqrt_symmetric(s_dense, gmat)
work(:, :) = MATMUL(h_dense, gmat)
horth(:, :) = MATMUL(gmat, work)
CALL symmetrize(horth)
uvec(:, :) = horth
CALL diagonalize_symmetric(uvec, eig)
DO idir = 1, 3
CALL build_commutator_kernel(h_dense, particle_set, cell, ao_atom, idir, op)
work(:, :) = MATMUL(op, gmat)
op_orth(:, :) = MATMUL(gmat, work)
CALL transform_to_eigenbasis(op_orth, uvec, dHH(:, :, idir), work)
END DO
DEALLOCATE (gmat, horth, op_orth, uvec)
END IF
DEALLOCATE (h_dense, op, work)
END SUBROUTINE setup_spin_transport
! **************************************************************************************************
!> \brief Symmetric inverse square root via dense diagonalization.
!> \param matrix ...
!> \param inverse_sqrt ...
! **************************************************************************************************
SUBROUTINE inverse_sqrt_symmetric(matrix, inverse_sqrt)
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: matrix
REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: inverse_sqrt
INTEGER :: i, info, lwork, n
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eig, work
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: scaled, vectors
n = SIZE(matrix, 1)
lwork = MAX(1, 8*n)
ALLOCATE (eig(n), scaled(n, n), vectors(n, n), work(lwork))
vectors(:, :) = matrix
CALL dsyev("V", "L", n, vectors, n, eig, work, lwork, info)
IF (info /= 0) CPABORT("Diagonalization of the overlap matrix failed.")
IF (MINVAL(eig) <= 0.0_dp) CPABORT("Overlap matrix is not positive definite.")
scaled(:, :) = vectors
DO i = 1, n
scaled(:, i) = scaled(:, i)/SQRT(eig(i))
END DO
CALL dgemm("N", "T", n, n, n, 1.0_dp, scaled, n, vectors, n, 0.0_dp, inverse_sqrt, n)
DEALLOCATE (eig, scaled, vectors, work)
END SUBROUTINE inverse_sqrt_symmetric
! **************************************************************************************************
!> \brief Dense symmetric diagonalization.
!> \param matrix ...
!> \param eig ...
! **************************************************************************************************
SUBROUTINE diagonalize_symmetric(matrix, eig)
REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: matrix
REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eig
INTEGER :: info, lwork, n
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work
n = SIZE(matrix, 1)
lwork = MAX(1, 8*n)
ALLOCATE (work(lwork))
CALL dsyev("V", "L", n, matrix, n, eig, work, lwork, info)
IF (info /= 0) CPABORT("Diagonalization of the orthogonalized KS matrix failed.")
DEALLOCATE (work)
END SUBROUTINE diagonalize_symmetric
! **************************************************************************************************
!> \brief Enforce exact symmetry after dense products.
!> \param matrix ...
! **************************************************************************************************
SUBROUTINE symmetrize(matrix)
REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: matrix
INTEGER :: i, j, n
REAL(KIND=dp) :: value
n = SIZE(matrix, 1)
DO j = 1, n
DO i = j + 1, n
value = 0.5_dp*(matrix(i, j) + matrix(j, i))
matrix(i, j) = value
matrix(j, i) = value
END DO
END DO
END SUBROUTINE symmetrize
! **************************************************************************************************
!> \brief Build element-wise kernel (r_col-r_row)_idir * matrix(row,col).
!> \param matrix ...
!> \param particle_set ...
!> \param cell ...
!> \param ao_atom ...
!> \param idir ...
!> \param op ...
! **************************************************************************************************
SUBROUTINE build_commutator_kernel(matrix, particle_set, cell, ao_atom, idir, op)
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: matrix
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(cell_type), POINTER :: cell
INTEGER, DIMENSION(:), INTENT(IN) :: ao_atom
INTEGER, INTENT(IN) :: idir
REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: op
INTEGER :: i, j, n
REAL(KIND=dp), DIMENSION(3) :: dr
n = SIZE(matrix, 1)
DO j = 1, n
DO i = 1, n
dr(:) = pbc(particle_set(ao_atom(j))%r(:) - particle_set(ao_atom(i))%r(:), cell)
op(i, j) = dr(idir)*matrix(i, j)
END DO
END DO
END SUBROUTINE build_commutator_kernel
! **************************************************************************************************
!> \brief Transform an AO matrix to the eigenvector basis.
!> \param op ...
!> \param uvec ...
!> \param transformed ...
!> \param work ...
! **************************************************************************************************
SUBROUTINE transform_to_eigenbasis(op, uvec, transformed, work)
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: op, uvec
REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: transformed
REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: work
INTEGER :: n
n = SIZE(op, 1)
CALL dgemm("N", "N", n, n, n, 1.0_dp, op, n, uvec, n, 0.0_dp, work, n)
CALL dgemm("T", "N", n, n, n, 1.0_dp, uvec, n, work, n, 0.0_dp, transformed, n)
END SUBROUTINE transform_to_eigenbasis
! **************************************************************************************************
!> \brief Fermi occupation with overflow protection.
!> \param eig ...
!> \param mu ...
!> \param kT ...
!> \param maxocc ...
!> \return ...
! **************************************************************************************************
PURE FUNCTION fermi_occ(eig, mu, kT, maxocc) RESULT(occ)
REAL(KIND=dp), INTENT(IN) :: eig, mu, kT, maxocc
REAL(KIND=dp) :: occ
REAL(KIND=dp) :: x
x = (eig - mu)/kT
IF (x > 80.0_dp) THEN
occ = 0.0_dp
ELSEIF (x < -80.0_dp) THEN
occ = maxocc
ELSE
occ = maxocc/(1.0_dp + EXP(x))
END IF
END FUNCTION fermi_occ
! **************************************************************************************************
!> \brief Fermi vacancy below the neutral point.
!> \param eig ...
!> \param mu ...
!> \param kT ...
!> \param maxocc ...
!> \return ...
! **************************************************************************************************
PURE FUNCTION fermi_hole(eig, mu, kT, maxocc) RESULT(hole)
REAL(KIND=dp), INTENT(IN) :: eig, mu, kT, maxocc
REAL(KIND=dp) :: hole
hole = maxocc - fermi_occ(eig, mu, kT, maxocc)
END FUNCTION fermi_hole
! **************************************************************************************************
!> \brief Total electron count at a chemical potential.
!> \param eig ...
!> \param maxocc ...
!> \param mu ...
!> \param kT ...
!> \return ...
! **************************************************************************************************
PURE FUNCTION electron_count(eig, maxocc, mu, kT) RESULT(count)
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: eig
REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: maxocc
REAL(KIND=dp), INTENT(IN) :: mu, kT
REAL(KIND=dp) :: count
INTEGER :: ispin, n
count = 0.0_dp
DO ispin = 1, SIZE(eig, 2)
DO n = 1, SIZE(eig, 1)
count = count + fermi_occ(eig(n, ispin), mu, kT, maxocc(ispin))
END DO
END DO
END FUNCTION electron_count
! **************************************************************************************************
!> \brief Determine the neutral chemical potential using the legacy two-step definition.
!> \param eig ...
!> \param maxocc ...
!> \param kT ...
!> \param neutral_electrons ...
!> \param neutral_grid ...
!> \param emin ...
!> \param emax ...
!> \param mu0 ...
! **************************************************************************************************
SUBROUTINE find_neutral_mu(eig, maxocc, kT, neutral_electrons, neutral_grid, emin, emax, mu0)
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: eig
REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: maxocc
REAL(KIND=dp), INTENT(IN) :: kT, neutral_electrons
INTEGER, INTENT(IN) :: neutral_grid
REAL(KIND=dp), INTENT(IN) :: emin, emax
REAL(KIND=dp), INTENT(OUT) :: mu0
INTEGER :: i
REAL(KIND=dp) :: best, count, mu, mup, ne, nh, &
previous_count
mup = 0.5_dp*(emin + emax)
previous_count = electron_count(eig, maxocc, emin, kT)
DO i = 1, neutral_grid
mu = emin + (emax - emin)*REAL(i, KIND=dp)/REAL(neutral_grid, KIND=dp)
count = electron_count(eig, maxocc, mu, kT)
IF ((count - neutral_electrons)*(previous_count - neutral_electrons) <= 0.0_dp) mup = mu
previous_count = count
END DO
mu0 = mup
best = HUGE(1.0_dp)
DO i = 1, neutral_grid
mu = emin + (emax - emin)*REAL(i, KIND=dp)/REAL(neutral_grid, KIND=dp)
CALL electron_hole_density(eig, maxocc, mu, mup, kT, ne, nh)
IF (ABS(nh - ne) <= best) THEN
mu0 = mu
best = ABS(nh - ne)
END IF
END DO
END SUBROUTINE find_neutral_mu
! **************************************************************************************************
!> \brief Electron and hole counts relative to a neutral point.
!> \param eig ...
!> \param maxocc ...
!> \param mu ...
!> \param mu0 ...
!> \param kT ...
!> \param ne ...
!> \param nh ...
! **************************************************************************************************
SUBROUTINE electron_hole_density(eig, maxocc, mu, mu0, kT, ne, nh)
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: eig
REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: maxocc
REAL(KIND=dp), INTENT(IN) :: mu, mu0, kT
REAL(KIND=dp), INTENT(OUT) :: ne, nh
INTEGER :: ispin, n
ne = 0.0_dp
nh = 0.0_dp
DO ispin = 1, SIZE(eig, 2)
DO n = 1, SIZE(eig, 1)
IF (eig(n, ispin) <= mu0) nh = nh + fermi_hole(eig(n, ispin), mu, kT, maxocc(ispin))
IF (eig(n, ispin) >= mu0) ne = ne + fermi_occ(eig(n, ispin), mu, kT, maxocc(ispin))
END DO
END DO
END SUBROUTINE electron_hole_density
! **************************************************************************************************
!> \brief Conductivity tensor at one chemical potential.
!> \param eig ...
!> \param dHH ...
!> \param maxocc ...
!> \param mu ...
!> \param kT ...
!> \param dissipation ...
!> \param volume ...
!> \param sigma ...
! **************************************************************************************************
SUBROUTINE conductivity_at_mu(eig, dHH, maxocc, mu, kT, dissipation, volume, sigma)
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: eig
REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: dHH
REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: maxocc
REAL(KIND=dp), INTENT(IN) :: mu, kT, dissipation, volume
REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: sigma
INTEGER :: idir, ispin, jdir, m, n, nao
REAL(KIND=dp) :: delta, eps, factor, g0
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: occ
nao = SIZE(eig, 1)
eps = 1.0E-12_dp
g0 = e_charge**2/(pi*h_bar)
sigma = 0.0_dp
ALLOCATE (occ(nao))
DO ispin = 1, SIZE(eig, 2)
DO n = 1, nao
occ(n) = fermi_occ(eig(n, ispin), mu, kT, maxocc(ispin))
END DO
DO idir = 1, 3
DO jdir = 1, 3
DO n = 1, nao
DO m = 1, nao
delta = eig(n, ispin) - eig(m, ispin)
factor = (occ(n) - occ(m))/(delta + eps)*dissipation/(dissipation**2 + delta**2)
sigma(jdir, idir) = sigma(jdir, idir) + &
dHH(m, n, jdir, ispin)*dHH(n, m, idir, ispin)*factor
END DO
END DO
END DO
END DO
END DO
sigma = pi*g0*sigma/(volume*angstrom)
DEALLOCATE (occ)
END SUBROUTINE conductivity_at_mu
END MODULE qs_kubo_transport

View file

@ -151,6 +151,7 @@ MODULE qs_scf_post_gpw
USE qs_ks_methods, ONLY: calc_rho_tot_gspace,&
qs_ks_update_qs_env
USE qs_ks_types, ONLY: qs_ks_did_change
USE qs_kubo_transport, ONLY: qs_scf_post_kubo_transport
USE qs_loc_dipole, ONLY: loc_dipole
USE qs_loc_states, ONLY: get_localization_info
USE qs_loc_types, ONLY: qs_loc_env_create,&
@ -3146,6 +3147,9 @@ CONTAINS
! Write cube files from the implicit Poisson solver
CALL qs_scf_post_ps_implicit(input, logger, qs_env)
! post SCF finite-volume Kubo transport
CALL qs_scf_post_kubo_transport(qs_env)
! post SCF Transport
CALL qs_scf_post_transport(qs_env)

View file

@ -0,0 +1,64 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT H2-kubo-diag-reuse
RUN_TYPE ENERGY
&END GLOBAL
&FORCE_EVAL
METHOD QS
&DFT
BASIS_SET_FILE_NAME BASIS_SET
POTENTIAL_FILE_NAME POTENTIAL
&MGRID
CUTOFF 80
NGRIDS 4
&END MGRID
&QS
EPS_DEFAULT 1.0E-10
METHOD GPW
&END QS
&SCF
ADDED_MOS -1
EPS_SCF 1.0E-7
MAX_SCF 30
SCF_GUESS ATOMIC
&DIAGONALIZATION
ALGORITHM STANDARD
&END DIAGONALIZATION
&MIXING
ALPHA 0.4
METHOD DIRECT_P_MIXING
&END MIXING
&OT OFF
&END OT
&END SCF
&XC
&XC_FUNCTIONAL PADE
&END XC_FUNCTIONAL
&END XC
&END DFT
&PROPERTIES
&KUBO_TRANSPORT ON
DISSIPATION [K] 300
ENERGY_RANGE [hartree] -0.35 -0.34
NEUTRAL_GRID 20
NEUTRAL_MU [hartree] -0.35
N_MU 1
TEMPERATURE [K] 300
&END KUBO_TRANSPORT
&END PROPERTIES
&SUBSYS
&CELL
ABC 6.0 6.0 6.0
PERIODIC NONE
&END CELL
&COORD
H 0.000000 0.000000 0.000000
H 0.720000 0.000000 0.000000
&END COORD
&KIND H
BASIS_SET DZVP-GTH-PADE
POTENTIAL GTH-PADE-q1
&END KIND
&END SUBSYS
&END FORCE_EVAL

View file

@ -0,0 +1,57 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT H2-kubo-gapw
RUN_TYPE ENERGY
&END GLOBAL
&FORCE_EVAL
METHOD QS
&DFT
BASIS_SET_FILE_NAME BASIS_SET
POTENTIAL_FILE_NAME POTENTIAL
&MGRID
CUTOFF 80
NGRIDS 4
&END MGRID
&QS
EPS_DEFAULT 1.0E-10
METHOD GAPW
&END QS
&SCF
EPS_SCF 1.0E-7
MAX_SCF 30
SCF_GUESS ATOMIC
&OT ON
PRECONDITIONER FULL_ALL
&END OT
&END SCF
&XC
&XC_FUNCTIONAL PADE
&END XC_FUNCTIONAL
&END XC
&END DFT
&PROPERTIES
&KUBO_TRANSPORT ON
DISSIPATION [K] 300
ENERGY_RANGE [hartree] -0.35 -0.34
NEUTRAL_GRID 20
NEUTRAL_MU [hartree] -0.35
N_MU 1
TEMPERATURE [K] 300
&END KUBO_TRANSPORT
&END PROPERTIES
&SUBSYS
&CELL
ABC 6.0 6.0 6.0
PERIODIC NONE
&END CELL
&COORD
H 0.000000 0.000000 0.000000
H 0.720000 0.000000 0.000000
&END COORD
&KIND H
BASIS_SET DZVP-GTH-PADE
POTENTIAL GTH-PADE-q1
&END KIND
&END SUBSYS
&END FORCE_EVAL

View file

@ -0,0 +1,57 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT H2-kubo-ot
RUN_TYPE ENERGY
&END GLOBAL
&FORCE_EVAL
METHOD QS
&DFT
BASIS_SET_FILE_NAME BASIS_SET
POTENTIAL_FILE_NAME POTENTIAL
&MGRID
CUTOFF 80
NGRIDS 4
&END MGRID
&QS
EPS_DEFAULT 1.0E-10
METHOD GPW
&END QS
&SCF
EPS_SCF 1.0E-7
MAX_SCF 30
SCF_GUESS ATOMIC
&OT ON
PRECONDITIONER FULL_ALL
&END OT
&END SCF
&XC
&XC_FUNCTIONAL PADE
&END XC_FUNCTIONAL
&END XC
&END DFT
&PROPERTIES
&KUBO_TRANSPORT ON
DISSIPATION [K] 300
ENERGY_RANGE [hartree] -0.35 -0.34
NEUTRAL_GRID 20
NEUTRAL_MU [hartree] -0.35
N_MU 1
TEMPERATURE [K] 300
&END KUBO_TRANSPORT
&END PROPERTIES
&SUBSYS
&CELL
ABC 6.0 6.0 6.0
PERIODIC NONE
&END CELL
&COORD
H 0.000000 0.000000 0.000000
H 0.720000 0.000000 0.000000
&END COORD
&KIND H
BASIS_SET DZVP-GTH-PADE
POTENTIAL GTH-PADE-q1
&END KIND
&END SUBSYS
&END FORCE_EVAL

View file

@ -0,0 +1,10 @@
"H2-kubo-ot.inp" = [{matcher="Kubo_sigma_iso", tol=1e-7, ref=0.290326773},
{matcher="Kubo_internal_diag", tol=0.0, ref=0.0},
{matcher="Kubo_reference", tol=0.0, ref=0.0}]
"H2-kubo-diag-reuse.inp" = [{matcher="Kubo_sigma_iso", tol=1e-7, ref=0.290326774},
{matcher="Kubo_reused_mos", tol=0.0, ref=0.0},
{matcher="Kubo_reference", tol=0.0, ref=0.0}]
"H2-kubo-gapw.inp" = [{matcher="Kubo_sigma_iso", tol=1e-7, ref=0.290798827},
{matcher="Kubo_internal_diag", tol=0.0, ref=0.0},
{matcher="Kubo_reference", tol=0.0, ref=0.0}]
#EOF

View file

@ -3,6 +3,7 @@
# Directories have been reordered according the execution time needed for a gfortran pdbg run using 2 MPI tasks
# in case a new directory is added just add it at the top of the list..
# the order will be regularly checked and modified...
QS/regtest-kubo-transport
QS/regtest-gw-realspace libint !ifx
QS/regtest-moments-kpoints
MIMIC/regtest-qmmm mimic mpiranks==4

View file

@ -118,6 +118,12 @@ registry["M009"] = GenericMatcher(r"PINT| Total energy =", col=5)
registry["M010"] = GenericMatcher(r"BAND TOTAL ENERGY [au]", col=6)
registry["M011"] = GenericMatcher(r"ENERGY| Total FORCE_EVAL", col=9)
registry["N_special_kpoints"] = GenericMatcher(r"Number of Special K-points:", col=5)
registry["Kubo_sigma_iso"] = GenericMatcher(r"KUBO_TRANSPORT| sigma_iso[S/cm]", col=3)
registry["Kubo_internal_diag"] = TextPresenceMatcher("internal diagonalization")
registry["Kubo_reused_mos"] = TextPresenceMatcher("CP2K canonical MOs")
registry["Kubo_reference"] = TextPresenceMatcher(
"https://doi.org/10.1016/j.aop.2020.168290"
)
registry["M012"] = GenericMatcher(r"B2(T) =", col=4)
registry["M013"] = GenericMatcher(r"sparseness function f2 =", col=5)
registry["M014"] = GenericMatcher(r"CheckSum Shifts =", col=4)