Skip to content

Commit

Permalink
Splitting of electronic states due to spin-orbit coupling
Browse files Browse the repository at this point in the history 
The SOC parameters from Hartwigsen, Goedecker, Hutter (https://doi.org/10.1103/PhysRevB.58.3641), Eq. 18, 19 + Table I, are used to compute the spin-orbit shift of electronic states (see Eq. 11 in https://iopscience.iop.org/article/10.1088/0953-8984/25/10/105503).
  • Loading branch information
@JWilhelm
JWilhelm committed Jun 23, 2023
1 parent 6ad6eaa commit 08caac2
Show file tree
Hide file tree
Showing 16 changed files with 3,161 additions and 108 deletions.
1 change: 1 addition & 0 deletions src/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1180,6 +1180,7 @@ list(
subsys/cp_subsys_types.F
subsys/damping_dipole_types.F
subsys/external_potential_types.F
subsys/external_potential_soc_parameters.F
subsys/force_field_kind_types.F
subsys/molecule_kind_list_types.F
subsys/molecule_kind_types.F
Expand Down
110 changes: 91 additions & 19 deletions src/core_ppnl.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@
!> - Extended by the derivatives for DFPT [Sandra Luber, Edward Ditler, 2021]
! **************************************************************************************************
MODULE core_ppnl
USE ai_angmom, ONLY: angmom
USE ai_overlap, ONLY: overlap
USE atomic_kind_types, ONLY: atomic_kind_type,&
get_atomic_kind_set
Expand Down Expand Up @@ -80,10 +81,11 @@ MODULE core_ppnl
!> \param cell_to_index ...
!> \param basis_type ...
!> \param deltaR Weighting factors of the derivatives wrt. nuclear positions
!> \param matrix_l ...
! **************************************************************************************************
SUBROUTINE build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, eps_ppnl, &
nimages, cell_to_index, basis_type, deltaR)
nimages, cell_to_index, basis_type, deltaR, matrix_l)

TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p
TYPE(qs_force_type), DIMENSION(:), POINTER :: force
Expand All @@ -102,24 +104,28 @@ SUBROUTINE build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces,
CHARACTER(LEN=*), INTENT(IN) :: basis_type
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
OPTIONAL :: deltaR
TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
POINTER :: matrix_l

CHARACTER(LEN=*), PARAMETER :: routineN = 'build_core_ppnl'

INTEGER :: atom_a, first_col, handle, i, iab, iac, iatom, ibc, icol, ikind, ilist, img, &
irow, iset, j, jatom, jkind, jneighbor, kac, katom, kbc, kkind, l, lc_max, lc_min, ldai, &
ldsab, lppnl, maxco, maxder, maxl, maxlgto, maxlppnl, maxppnl, maxsgf, na, natom, nb, &
ncoa, ncoc, nkind, nlist, nneighbor, nnl, np, nppnl, nprjc, nseta, nsgfa, prjc, sgfa, slot
INTEGER :: atom_a, first_col, handle, i, i_dim, iab, iac, iatom, ib, ibc, icol, ikind, &
ilist, img, irow, iset, j, jatom, jb, jkind, jneighbor, kac, katom, kbc, kkind, l, &
lc_max, lc_min, ldai, ldsab, lppnl, maxco, maxder, maxl, maxlgto, maxlppnl, maxppnl, &
maxsgf, na, natom, nb, ncoa, ncoc, nkind, nlist, nneighbor, nnl, np, nppnl, nprjc, nseta, &
nsgfa, prjc, sgfa, slot
INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
INTEGER, DIMENSION(3) :: cell_b, cell_c
INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nprj_ppnl, &
nsgf_seta
INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
LOGICAL :: do_dR, dogth, dokp, found, ppnl_present
LOGICAL :: do_dR, do_soc, dogth, dokp, found, &
ppnl_present
LOGICAL, DIMENSION(0:9) :: is_nonlocal
REAL(KIND=dp) :: dac, f0, ppnl_radius
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: radp
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: sab, work
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ai_work
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ai_work, lab, work_l
REAL(KIND=dp), DIMENSION(1) :: rprjc, zetc
REAL(KIND=dp), DIMENSION(3) :: fa, fb, rab, rac, rbc
REAL(KIND=dp), DIMENSION(3, 3) :: pv_thread
Expand All @@ -129,10 +135,10 @@ SUBROUTINE build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces,
TYPE(gth_potential_p_type), DIMENSION(:), POINTER :: gpotential
TYPE(clist_type), POINTER :: clist
TYPE(alist_type), POINTER :: alist_ac, alist_bc
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: achint, acint, bchint, bcint, c_nl
REAL(KIND=dp), DIMENSION(:, :), POINTER :: cprj, h_block, h_nl, p_block, r_2block, &
r_3block, rpgfa, sphi_a, vprj_ppnl, &
zeta
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: achint, acint, alkint, bchint, bcint, &
blkint, c_nl
REAL(KIND=dp), DIMENSION(:, :), POINTER :: cprj, h_block, h_nl, l_block_x, l_block_y, &
l_block_z, p_block, r_2block, r_3block, rpgfa, sphi_a, vprj_ppnl, wprj_ppnl, zeta
REAL(KIND=dp), DIMENSION(:), POINTER :: a_nl, alpha_ppnl, hprj, set_radius_a
REAL(KIND=dp), DIMENSION(3, SIZE(particle_set)) :: force_thread
TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
Expand All @@ -156,6 +162,8 @@ SUBROUTINE build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces,
CALL timeset(routineN, handle)
END IF

do_soc = PRESENT(matrix_l)

ppnl_present = ASSOCIATED(sap_ppnl)

IF (ppnl_present) THEN
Expand Down Expand Up @@ -262,19 +270,23 @@ SUBROUTINE build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces,
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP SHARED (basis_set, gpotential, spotential, maxder, ncoset, &
!$OMP sap_ppnl, sap_int, nkind, ldsab, ldai, nder, nco ) &
!$OMP sap_ppnl, sap_int, nkind, ldsab, ldai, nder, nco, do_soc ) &
!$OMP PRIVATE (ikind, kkind, iatom, katom, nlist, ilist, nneighbor, jneighbor, &
!$OMP cell_c, rac, iac, first_sgfa, la_max, la_min, npgfa, nseta, nsgfa, nsgf_seta, &
!$OMP slot, sphi_a, zeta, cprj, hprj, lppnl, nppnl, nprj_ppnl, &
!$OMP clist, iset, ncoa, sgfa, prjc, work, sab, ai_work, nprjc, ppnl_radius, &
!$OMP ncoc, rpgfa, first_col, vprj_ppnl, i, j, l, dogth, &
!$OMP clist, iset, ncoa, sgfa, prjc, work, work_l, sab, lab, ai_work, nprjc, &
!$OMP ppnl_radius, ncoc, rpgfa, first_col, vprj_ppnl, wprj_ppnl, i, j, l, dogth, &
!$OMP set_radius_a, rprjc, dac, lc_max, lc_min, zetc, alpha_ppnl, &
!$OMP na, nb, np, nnl, is_nonlocal, a_nl, h_nl, c_nl, radp)
!$OMP na, nb, np, nnl, is_nonlocal, a_nl, h_nl, c_nl, radp, i_dim, ib, jb)

ALLOCATE (sab(ldsab, ldsab*maxder), work(ldsab, ldsab*maxder))
sab = 0.0_dp
ALLOCATE (ai_work(ldai, ldai, ncoset(nder + 1)))
ai_work = 0.0_dp
IF (do_soc) THEN
ALLOCATE (lab(ldsab, ldsab, 3), work_l(ldsab, ldsab, 3))
lab = 0.0_dp
END IF

!$OMP DO SCHEDULE(GUIDED)
DO slot = 1, sap_ppnl(1)%nl_size
Expand Down Expand Up @@ -315,6 +327,7 @@ SUBROUTINE build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces,
nprj_ppnl => gpotential(kkind)%gth_potential%nprj_ppnl
ppnl_radius = gpotential(kkind)%gth_potential%ppnl_radius
vprj_ppnl => gpotential(kkind)%gth_potential%vprj_ppnl
wprj_ppnl => gpotential(kkind)%gth_potential%wprj_ppnl
ELSE IF (ASSOCIATED(spotential(kkind)%sgp_potential)) THEN
! SGP potential
dogth = .FALSE.
Expand Down Expand Up @@ -343,9 +356,14 @@ SUBROUTINE build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces,
clist%cell = cell_c
clist%rac = rac
ALLOCATE (clist%acint(nsgfa, nppnl, maxder), &
clist%achint(nsgfa, nppnl, maxder))
clist%achint(nsgfa, nppnl, maxder), &
clist%alint(nsgfa, nppnl, 3), &
clist%alkint(nsgfa, nppnl, 3))
clist%acint = 0._dp
clist%achint = 0._dp
clist%alint = 0._dp
clist%alkint = 0._dp

clist%nsgf_cnt = 0
NULLIFY (clist%sgf_list)
DO iset = 1, nseta
Expand Down Expand Up @@ -379,7 +397,20 @@ SUBROUTINE build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces,
work(1:na, first_col + prjc:first_col + prjc + nb - 1) = &
MATMUL(sab(1:na, first_col + 1:first_col + np), cprj(1:np, prjc:prjc + nb - 1))
END DO

IF (do_soc) THEN
! *** Calculate the primitive angular momentum integrals needed for spin-orbit coupling ***
lab = 0.0_dp
CALL angmom(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
lc_max, 1, zetc, rprjc, -rac, (/0._dp, 0._dp, 0._dp/), lab)
DO i_dim = 1, 3
work_l(1:na, prjc:prjc + nb - 1, i_dim) = &
MATMUL(lab(1:na, 1:np, i_dim), cprj(1:np, prjc:prjc + nb - 1))
END DO
END IF

prjc = prjc + nprjc

END DO
na = nsgf_seta(iset)
nb = nppnl
Expand All @@ -391,12 +422,22 @@ SUBROUTINE build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces,
! work(1, first_col), ldsab, 0.0_dp, clist%acint(sgfa, 1, i), nsgfa)
clist%acint(sgfa:sgfa + na - 1, 1:nb, i) = &
MATMUL(TRANSPOSE(sphi_a(1:np, sgfa:sgfa + na - 1)), work(1:np, first_col:first_col + nb - 1))

! *** Multiply with interaction matrix(h) ***
! CALL dgemm("N", "N", nsgf_seta(iset), nppnl, nppnl, 1.0_dp, clist%acint(sgfa, 1, i), nsgfa, &
! vprj_ppnl(1, 1), SIZE(vprj_ppnl, 1), 0.0_dp, clist%achint(sgfa, 1, i), nsgfa)
clist%achint(sgfa:sgfa + na - 1, 1:nb, i) = &
MATMUL(clist%acint(sgfa:sgfa + na - 1, 1:nb, i), vprj_ppnl(1:nb, 1:nb))

END DO
IF (do_soc) THEN
DO i_dim = 1, 3
clist%alint(sgfa:sgfa + na - 1, 1:nb, i_dim) = &
MATMUL(TRANSPOSE(sphi_a(1:np, sgfa:sgfa + na - 1)), work_l(1:np, 1:nb, i_dim))
clist%alkint(sgfa:sgfa + na - 1, 1:nb, i_dim) = &
MATMUL(clist%alint(sgfa:sgfa + na - 1, 1:nb, i_dim), wprj_ppnl(1:nb, 1:nb))
END DO
END IF
ELSE
! SGP potential
! *** Calculate the primitive overlap integrals ***
Expand Down Expand Up @@ -430,6 +471,7 @@ SUBROUTINE build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces,
END DO

DEALLOCATE (sab, ai_work, work)
IF (do_soc) DEALLOCATE (lab, work_l)
!$OMP END PARALLEL

! *** Set up a sorting index
Expand All @@ -442,15 +484,17 @@ SUBROUTINE build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces,

!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP SHARED (dokp, basis_set, matrix_h, cell_to_index,&
!$OMP SHARED (dokp, basis_set, matrix_h, matrix_l, cell_to_index,&
!$OMP sab_orb, matrix_p, sap_int, nkind, eps_ppnl, force, &
!$OMP do_dR, deltaR, maxder, nder, &
!$OMP locks, virial, use_virial, calculate_forces) &
!$OMP locks, virial, use_virial, calculate_forces, do_soc) &
!$OMP PRIVATE (ikind, jkind, iatom, jatom, cell_b, rab, &
!$OMP slot, iab, atom_a, f0, irow, icol, h_block, &
!$OMP l_block_x, l_block_y, l_block_z, &
!$OMP r_2block, r_3block, &
!$OMP found,p_block, iac, ibc, alist_ac, alist_bc, acint, bcint, &
!$OMP achint, bchint, na, np, nb, katom, j, fa, fb, rbc, rac, &
!$OMP achint, bchint, alkint, blkint, &
!$OMP na, np, nb, katom, j, fa, fb, rbc, rac, &
!$OMP kkind, kac, kbc, i, img, hash, iatom8, natom) &
!$OMP REDUCTION (+ : pv_thread, force_thread )

Expand Down Expand Up @@ -502,6 +546,12 @@ SUBROUTINE build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces,
END IF
NULLIFY (h_block)
CALL dbcsr_get_block_p(matrix_h(1, img)%matrix, irow, icol, h_block, found)
IF (do_soc) THEN
NULLIFY (l_block_x, l_block_y, l_block_z)
CALL dbcsr_get_block_p(matrix_l(1, img)%matrix, irow, icol, l_block_x, found)
CALL dbcsr_get_block_p(matrix_l(2, img)%matrix, irow, icol, l_block_y, found)
CALL dbcsr_get_block_p(matrix_l(3, img)%matrix, irow, icol, l_block_z, found)
END IF

IF (do_dR) THEN
NULLIFY (r_2block, r_3block)
Expand Down Expand Up @@ -534,6 +584,10 @@ SUBROUTINE build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces,
bcint => alist_bc%clist(kbc)%acint
achint => alist_ac%clist(kac)%achint
bchint => alist_bc%clist(kbc)%achint
IF (do_soc) THEN
alkint => alist_ac%clist(kac)%alkint
blkint => alist_bc%clist(kbc)%alkint
END IF
na = SIZE(acint, 1)
np = SIZE(acint, 2)
nb = SIZE(bcint, 1)
Expand All @@ -549,6 +603,24 @@ SUBROUTINE build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces,
MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 1)))
END IF
END IF
IF (do_soc) THEN
IF (iatom <= jatom) THEN
l_block_x(1:na, 1:nb) = l_block_x(1:na, 1:nb) + &
MATMUL(alkint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 1)))
l_block_y(1:na, 1:nb) = l_block_y(1:na, 1:nb) + &
MATMUL(alkint(1:na, 1:np, 2), TRANSPOSE(bcint(1:nb, 1:np, 1)))
l_block_z(1:na, 1:nb) = l_block_z(1:na, 1:nb) + &
MATMUL(alkint(1:na, 1:np, 3), TRANSPOSE(bcint(1:nb, 1:np, 1)))

ELSE
l_block_x(1:nb, 1:na) = l_block_x(1:nb, 1:na) - &
MATMUL(blkint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 1)))
l_block_y(1:nb, 1:na) = l_block_y(1:nb, 1:na) - &
MATMUL(blkint(1:nb, 1:np, 2), TRANSPOSE(acint(1:na, 1:np, 1)))
l_block_z(1:nb, 1:na) = l_block_z(1:nb, 1:na) - &
MATMUL(blkint(1:nb, 1:np, 3), TRANSPOSE(acint(1:na, 1:np, 1)))
END IF
END IF
!$ CALL omp_unset_lock(locks(hash))
IF (calculate_forces) THEN
IF (ASSOCIATED(p_block)) THEN
Expand Down
5 changes: 4 additions & 1 deletion src/input_constants.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1068,7 +1068,10 @@ MODULE input_constants
kp_weights_W_uniform = 12, &
kp_weights_W_tailored = 13, &
gw_gf_mic = 14, &
gw_gf_gamma = 15
gw_gf_gamma = 15, &
soc_none = 16, &
soc_lda = 17, &
soc_pbe = 18

! periodic RESP parameters
INTEGER, PARAMETER, PUBLIC :: do_resp_x_dir = 0, &
Expand Down
31 changes: 29 additions & 2 deletions src/input_cp2k_mp2.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,8 +33,9 @@ MODULE input_cp2k_mp2
mp2_method_none, numerical, ot_precond_full_all, ot_precond_full_kinetic, &
ot_precond_full_single, ot_precond_full_single_inverse, ot_precond_none, &
ot_precond_s_inverse, ri_default, ri_rpa_g0w0_crossing_bisection, &
ri_rpa_g0w0_crossing_newton, ri_rpa_g0w0_crossing_z_shot, wfc_mm_style_gemm, &
wfc_mm_style_syrk, z_solver_cg, z_solver_pople, z_solver_richardson, z_solver_sd
ri_rpa_g0w0_crossing_newton, ri_rpa_g0w0_crossing_z_shot, soc_lda, soc_none, soc_pbe, &
wfc_mm_style_gemm, wfc_mm_style_syrk, z_solver_cg, z_solver_pople, z_solver_richardson, &
z_solver_sd
USE input_cp2k_hfx, ONLY: create_hfx_section
USE input_cp2k_kpoints, ONLY: create_kpoint_set_section
USE input_keyword_types, ONLY: keyword_create,&
Expand Down Expand Up @@ -931,6 +932,32 @@ SUBROUTINE create_ri_g0w0(section)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

CALL keyword_create(keyword, __LOCATION__, name="SOC", &
description="Calculate the spin-orbit splitting of the eigenvalues/band structure "// &
"using the spin-orbit part of the GTH pseudos parametrized in Hartwigsen, Goedecker, "// &
"Hutter, Phys. Rev. B 58, 3641 (1998), Eq. 19, "// &
"parameters in Table I.", &
usage="SOC", &
enum_c_vals=s2a("NONE", "LDA", "PBE"), &
enum_i_vals=(/soc_none, soc_lda, soc_pbe/), &
enum_desc=s2a("No SOC.", &
"Use parameters from LDA (PADE) pseudopotential.", &
"Use parameters from PBE pseudopotential."), &
default_i_val=soc_none)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

CALL keyword_create(keyword, __LOCATION__, name="SOC_ENERGY_WINDOW", &
description="For perturbative SOC calculation, only "// &
"take frontier levels in an energy window "// &
"[E_HOMO - SOC_ENERGY_WINDOW/2 , E_LUMO + SOC_ENERGY_WINDOW/2 "// &
"into account for the diagonalization of H^GW,SOC.", &
usage="SOC_ENERGY_WINDOW 20.0_eV", &
default_r_val=cp_unit_to_cp2k(value=50.0_dp, unit_str="eV"), &
unit_str="eV")
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

! here we generate a subsection for the periodic GW correction
CALL create_periodic_gw_correction_section(subsection)
CALL section_add_subsection(section, subsection)
Expand Down
2 changes: 1 addition & 1 deletion src/mp2_integrals.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -379,7 +379,7 @@ SUBROUTINE mp2_ri_gpw_compute_in(BIb_C, BIb_C_gw, BIb_C_bse_ij, BIb_C_bse_ab, gd
WRITE (unit_nr, '(T3,A,T61,F20.10)') &
"RI_INFO| Omega: ", ri_metric%omega
CASE (do_potential_id)
WRITE (unit_nr, FMT="(T3,A,T73,A)") &
WRITE (unit_nr, FMT="(T3,A,T74,A)") &
"RI_INFO| RI metric: ", "OVERLAP"
CASE (do_potential_truncated)
WRITE (unit_nr, FMT="(T3,A,T64,A)") &
Expand Down

0 comments on commit 08caac2

Please sign in to comment.