diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index e71cc59127..eae6055e39 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -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 diff --git a/src/common/bibliography.F b/src/common/bibliography.F index 95db19cf45..0373ca92da 100644 --- a/src/common/bibliography.F +++ b/src/common/bibliography.F @@ -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", & diff --git a/src/input_cp2k_properties_dft.F b/src/input_cp2k_properties_dft.F index f9d4bac1aa..37bf74cfae 100644 --- a/src/input_cp2k_properties_dft.F +++ b/src/input_cp2k_properties_dft.F @@ -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 diff --git a/src/qs_kubo_transport.F b/src/qs_kubo_transport.F new file mode 100644 index 0000000000..cd0ff470ec --- /dev/null +++ b/src/qs_kubo_transport.F @@ -0,0 +1,687 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2026 CP2K developers group ! +! ! +! 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 diff --git a/src/qs_scf_post_gpw.F b/src/qs_scf_post_gpw.F index 9474446608..b41212ca07 100644 --- a/src/qs_scf_post_gpw.F +++ b/src/qs_scf_post_gpw.F @@ -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) diff --git a/tests/QS/regtest-kubo-transport/H2-kubo-diag-reuse.inp b/tests/QS/regtest-kubo-transport/H2-kubo-diag-reuse.inp new file mode 100644 index 0000000000..4ac35ec69c --- /dev/null +++ b/tests/QS/regtest-kubo-transport/H2-kubo-diag-reuse.inp @@ -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 diff --git a/tests/QS/regtest-kubo-transport/H2-kubo-gapw.inp b/tests/QS/regtest-kubo-transport/H2-kubo-gapw.inp new file mode 100644 index 0000000000..e7ab368765 --- /dev/null +++ b/tests/QS/regtest-kubo-transport/H2-kubo-gapw.inp @@ -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 diff --git a/tests/QS/regtest-kubo-transport/H2-kubo-ot.inp b/tests/QS/regtest-kubo-transport/H2-kubo-ot.inp new file mode 100644 index 0000000000..856e74e5b8 --- /dev/null +++ b/tests/QS/regtest-kubo-transport/H2-kubo-ot.inp @@ -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 diff --git a/tests/QS/regtest-kubo-transport/TEST_FILES.toml b/tests/QS/regtest-kubo-transport/TEST_FILES.toml new file mode 100644 index 0000000000..da34fd8faa --- /dev/null +++ b/tests/QS/regtest-kubo-transport/TEST_FILES.toml @@ -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 diff --git a/tests/TEST_DIRS b/tests/TEST_DIRS index 3b0d39d2e4..4e4d14d277 100644 --- a/tests/TEST_DIRS +++ b/tests/TEST_DIRS @@ -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 diff --git a/tests/matchers.py b/tests/matchers.py index 9d6785cd69..5f5ca39627 100644 --- a/tests/matchers.py +++ b/tests/matchers.py @@ -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)