soc_pseudopotential_methods.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   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                                                      !
!--------------------------------------------------------------------------------------------------!
 
MODULE soc_pseudopotential_methods
   USE atomic_kind_types,               ONLY: atomic_kind_type
   USE core_ppnl,                       ONLY: build_core_ppnl
   USE cp_cfm_types,                    ONLY: cp_cfm_get_info,&
                                              cp_cfm_type
   USE cp_control_types,                ONLY: dft_control_type
   USE cp_dbcsr_api,                    ONLY: dbcsr_add,&
                                              dbcsr_copy,&
                                              dbcsr_create,&
                                              dbcsr_desymmetrize,&
                                              dbcsr_p_type,&
                                              dbcsr_set,&
                                              dbcsr_type_antisymmetric,&
                                              dbcsr_type_no_symmetry
   USE cp_dbcsr_cp2k_link,              ONLY: cp_dbcsr_alloc_block_from_nbl
   USE cp_dbcsr_operations,             ONLY: dbcsr_allocate_matrix_set,&
                                              dbcsr_deallocate_matrix_set
   USE kinds,                           ONLY: dp
   USE kpoint_types,                    ONLY: get_kpoint_info,&
                                              kpoint_type
   USE particle_types,                  ONLY: particle_type
   USE qs_environment_types,            ONLY: get_qs_env,&
                                              qs_environment_type
   USE qs_force_types,                  ONLY: qs_force_type
   USE qs_kind_types,                   ONLY: qs_kind_type
   USE qs_neighbor_list_types,          ONLY: get_neighbor_list_set_p,&
                                              neighbor_list_set_p_type
   USE virial_types,                    ONLY: virial_type
   !USE physcon,                         ONLY: evolt
#include "./base/base_uses.f90"
 
   IMPLICIT NONE
 
   PRIVATE
 
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'soc_pseudopotential_methods'
 
   PUBLIC :: v_soc_xyz_from_pseudopotential, &
             remove_soc_outside_energy_window_mo
 
CONTAINS
 
! **************************************************************************************************
!> \brief V^SOC_µν^(α),R = ħ/2 < ϕ_µ cell O | sum_ℓ ΔV_ℓ^SO(r,r') L^(α) | ϕ_ν cell R>, α = x,y,z
!>        see Hartwigsen, Goedecker, Hutter, Eq.(18), (19) (doi.org/10.1103/PhysRevB.58.3641)
!>        Caution: V^SOC_µν^(α) is purely imaginary and Hermitian; V^SOC_µν^(α) is stored as real
!>                 dbcsr matrix mat_V_SOC_xyz without symmetry; V^SOC_µν^(α) is stored without
!>                 the imaginary unit, i.e. mat_V_SOC_xyz is real and antisymmetric
!> \param qs_env ...
!> \param mat_V_SOC_xyz ...
!> \par History
!>    * 09.2023 created
!> \author Jan Wilhelm
! **************************************************************************************************
   SUBROUTINE v_soc_xyz_from_pseudopotential(qs_env, mat_V_SOC_xyz)
      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: mat_v_soc_xyz
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'V_SOC_xyz_from_pseudopotential'
 
      INTEGER                                            :: handle, img, nder, nimages, xyz
      INTEGER, DIMENSION(:, :, :), POINTER               :: cell_to_index
      LOGICAL                                            :: calculate_forces, do_kp, do_symmetric, &
                                                            use_virial
      REAL(kind=dp)                                      :: eps_ppnl
      TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: mat_l, mat_l_nosym, mat_pot_dummy, &
                                                            matrix_dummy, matrix_s
      TYPE(dft_control_type), POINTER                    :: dft_control
      TYPE(kpoint_type), POINTER                         :: kpoints
      TYPE(neighbor_list_set_p_type), DIMENSION(:), &
         POINTER                                         :: sab_orb, sap_ppnl
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
      TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
      TYPE(virial_type), POINTER                         :: virial
 
      CALL timeset(routinen, handle)
 
      NULLIFY (qs_kind_set, dft_control, sab_orb, sap_ppnl, particle_set, atomic_kind_set, &
               cell_to_index)
      CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, dft_control=dft_control, &
                      matrix_s_kp=matrix_s, kpoints=kpoints, atomic_kind_set=atomic_kind_set, &
                      particle_set=particle_set, sab_orb=sab_orb, sap_ppnl=sap_ppnl, do_kpoints=do_kp)
 
      eps_ppnl = dft_control%qs_control%eps_ppnl
      nimages = dft_control%nimages
      ! Even a one-image k-point environment needs the triangular storage expected by
      ! rskp_transform. Desymmetrizing redistributes blocks away from their neighbour-list owners.
      CALL get_neighbor_list_set_p(neighbor_list_sets=sab_orb, symmetric=do_symmetric)
      CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
 
      NULLIFY (mat_l, mat_pot_dummy)
      CALL dbcsr_allocate_matrix_set(mat_l, 3, nimages)
      DO xyz = 1, 3
         DO img = 1, nimages
            ALLOCATE (mat_l(xyz, img)%matrix)
            CALL dbcsr_create(mat_l(xyz, img)%matrix, template=matrix_s(1, 1)%matrix, &
                              matrix_type=dbcsr_type_antisymmetric)
            CALL cp_dbcsr_alloc_block_from_nbl(mat_l(xyz, img)%matrix, sab_orb)
            CALL dbcsr_set(mat_l(xyz, img)%matrix, 0.0_dp)
         END DO
      END DO
 
      ! get mat_l; the next CPASSERT fails if the atoms do not have any SOC parameters, i.e.
      ! SOC is zero and one should not activate the SOC section
      cpassert(ASSOCIATED(sap_ppnl))
      nder = 0
      use_virial = .false.
      calculate_forces = .false.
 
      NULLIFY (mat_pot_dummy)
      CALL dbcsr_allocate_matrix_set(mat_pot_dummy, 1, nimages)
      DO img = 1, nimages
         ALLOCATE (mat_pot_dummy(1, img)%matrix)
         CALL dbcsr_create(mat_pot_dummy(1, img)%matrix, template=matrix_s(1, 1)%matrix)
         CALL cp_dbcsr_alloc_block_from_nbl(mat_pot_dummy(1, img)%matrix, sab_orb)
         CALL dbcsr_set(mat_pot_dummy(1, img)%matrix, 0.0_dp)
      END DO
 
      CALL build_core_ppnl(mat_pot_dummy, matrix_dummy, force, virial, &
                           calculate_forces, use_virial, nder, &
                           qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, &
                           eps_ppnl, nimages=nimages, cell_to_index=cell_to_index, &
                           basis_type="ORB", matrix_l=mat_l)
 
      NULLIFY (mat_l_nosym)
      CALL dbcsr_allocate_matrix_set(mat_l_nosym, 3, nimages)
      DO xyz = 1, 3
         DO img = 1, nimages
 
            ALLOCATE (mat_l_nosym(xyz, img)%matrix)
            IF (do_kp) THEN
               CALL dbcsr_create(mat_l_nosym(xyz, img)%matrix, template=matrix_s(1, 1)%matrix, &
                                 matrix_type=dbcsr_type_antisymmetric)
               CALL dbcsr_copy(mat_l_nosym(xyz, img)%matrix, mat_l(xyz, img)%matrix)
            ELSE
               CALL dbcsr_create(mat_l_nosym(xyz, img)%matrix, template=matrix_s(1, 1)%matrix, &
                                 matrix_type=dbcsr_type_no_symmetry)
               CALL dbcsr_desymmetrize(mat_l(xyz, img)%matrix, mat_l_nosym(xyz, img)%matrix)
            END IF
 
         END DO
      END DO
 
      NULLIFY (mat_v_soc_xyz)
      CALL dbcsr_allocate_matrix_set(mat_v_soc_xyz, 3, nimages)
      DO xyz = 1, 3
         DO img = 1, nimages
            ALLOCATE (mat_v_soc_xyz(xyz, img)%matrix)
            IF (do_kp) THEN
               ! mat_V_SOC_xyz^R with neighbor cell R actually has no symmetry
               ! mat_V_SOC_xyz^R_µν = mat_V_SOC_xyz^R_νµ* (the actual symmetry is
               ! mat_V_SOC_xyz^R_µν = mat_V_SOC_xyz^-R_νµ* )  but rskp_transform
               ! for mat_V_SOC_xyz^R -> mat_V_SOC_xyz(k) requires symmetry...
               CALL dbcsr_create(mat_v_soc_xyz(xyz, img)%matrix, template=matrix_s(1, 1)%matrix, &
                                 matrix_type=dbcsr_type_antisymmetric)
            ELSE
               CALL dbcsr_create(mat_v_soc_xyz(xyz, img)%matrix, template=matrix_s(1, 1)%matrix, &
                                 matrix_type=dbcsr_type_no_symmetry)
            END IF
            CALL cp_dbcsr_alloc_block_from_nbl(mat_v_soc_xyz(xyz, img)%matrix, sab_orb)
            ! factor 0.5 from ħ/2 prefactor
            CALL dbcsr_add(mat_v_soc_xyz(xyz, img)%matrix, mat_l_nosym(xyz, img)%matrix, &
                           0.0_dp, 0.5_dp)
         END DO
      END DO
 
      CALL dbcsr_deallocate_matrix_set(mat_pot_dummy)
      CALL dbcsr_deallocate_matrix_set(mat_l_nosym)
      CALL dbcsr_deallocate_matrix_set(mat_l)
 
      CALL timestop(handle)
 
   END SUBROUTINE v_soc_xyz_from_pseudopotential
 
! **************************************************************************************************
!> \brief ...
!> \param cfm_ks_spinor ...
!> \param e_win_cbm ...
!> \param temp_smear ...
!> \param eigenval ...
!> \param e_fermi ...
! **************************************************************************************************
   SUBROUTINE remove_soc_outside_energy_window_mo(cfm_ks_spinor, e_win_cbm, temp_smear, &
                                                  eigenval, e_fermi)
      TYPE(cp_cfm_type)                                  :: cfm_ks_spinor
      REAL(kind=dp)                                      :: e_win_cbm, temp_smear
      REAL(kind=dp), DIMENSION(:)                        :: eigenval
      REAL(kind=dp)                                      :: e_fermi
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'remove_soc_outside_energy_window_mo'
 
      INTEGER                                            :: handle, i_glob, iib, j_glob, jjb, &
                                                            ncol_global, ncol_local, nrow_global, &
                                                            nrow_local
      INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
      REAL(kind=dp)                                      :: e_i, e_j, w_i, w_j, x_i, x_j
 
!REAL(KIND=dp)                                      :: E_HOMO, E_LUMO, e_fermi
 
      ! Remove SOC outside of energy window (otherwise, numerical problems arise
      ! because energetically low semicore states and energetically very high
      ! unbound states couple to the states around the Fermi level).
      ! This routine is for cfm_ks_spinor being in the molecular-orbital (mo) with
      ! corresponding eigenvalues "eigenval".
 
      CALL timeset(routinen, handle)
 
      CALL cp_cfm_get_info(matrix=cfm_ks_spinor, &
                           nrow_global=nrow_global, &
                           ncol_global=ncol_global, &
                           nrow_local=nrow_local, &
                           ncol_local=ncol_local, &
                           row_indices=row_indices, &
                           col_indices=col_indices)
 
      cpassert(nrow_global == SIZE(eigenval))
      cpassert(ncol_global == SIZE(eigenval))
 
      ! Apply a smooth energy window to the SOC coupling matrix. Each state gets a
      ! weight in [0,1]: 1 = full SOC, decaying to 0 outside the window via a
      ! Fermi-like function of width temp_smear. A matrix element coupling states i
      ! and j is scaled by SQRT(w_i*w_j), applying the damping symmetrically to both.
      DO jjb = 1, ncol_local
         j_glob = col_indices(jjb)
         DO iib = 1, nrow_local
            i_glob = row_indices(iib)
 
            e_i = eigenval(i_glob)
            e_j = eigenval(j_glob)
            IF (e_i <= e_fermi) THEN
               !x_i = ABS(E_i - e_fermi) - 0.5_dp * e_win_vbm
               w_i = 1.0_dp
            ELSE
               x_i = abs(e_i - e_fermi) - e_win_cbm
               w_i = 1.0_dp/(exp(x_i/temp_smear) + 1.0_dp)
            END IF
 
            IF (e_j <= e_fermi) THEN
               !x_j = ABS(E_j - e_fermi) - 0.5_dp * e_win_vbm
               w_j = 1.0_dp
            ELSE
               x_j = abs(e_j - e_fermi) - e_win_cbm
               w_j = 1.0_dp/(exp(x_j/temp_smear) + 1.0_dp)
            END IF
 
            cfm_ks_spinor%local_data(iib, jjb) = &
               cfm_ks_spinor%local_data(iib, jjb)*sqrt(w_i*w_j)
 
         END DO
      END DO
 
      CALL timestop(handle)
 
   END SUBROUTINE remove_soc_outside_energy_window_mo
 
END MODULE soc_pseudopotential_methods