LCOV - code coverage report
Current view: top level - src - soc_pseudopotential_methods.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 70.5 % 112 79
Test Date: 2025-07-25 12:55:17 Functions: 66.7 % 3 2

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : MODULE soc_pseudopotential_methods
       9              :    USE atomic_kind_types,               ONLY: atomic_kind_type
      10              :    USE core_ppnl,                       ONLY: build_core_ppnl
      11              :    USE cp_cfm_types,                    ONLY: cp_cfm_get_info,&
      12              :                                               cp_cfm_type
      13              :    USE cp_control_types,                ONLY: dft_control_type
      14              :    USE cp_dbcsr_api,                    ONLY: dbcsr_add,&
      15              :                                               dbcsr_copy,&
      16              :                                               dbcsr_create,&
      17              :                                               dbcsr_desymmetrize,&
      18              :                                               dbcsr_p_type,&
      19              :                                               dbcsr_set,&
      20              :                                               dbcsr_type_antisymmetric,&
      21              :                                               dbcsr_type_no_symmetry
      22              :    USE cp_dbcsr_cp2k_link,              ONLY: cp_dbcsr_alloc_block_from_nbl
      23              :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      24              :                                               copy_fm_to_dbcsr,&
      25              :                                               dbcsr_allocate_matrix_set,&
      26              :                                               dbcsr_deallocate_matrix_set
      27              :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      28              :                                               cp_fm_get_info,&
      29              :                                               cp_fm_release,&
      30              :                                               cp_fm_type
      31              :    USE kinds,                           ONLY: dp
      32              :    USE kpoint_types,                    ONLY: get_kpoint_info,&
      33              :                                               kpoint_type
      34              :    USE parallel_gemm_api,               ONLY: parallel_gemm
      35              :    USE particle_types,                  ONLY: particle_type
      36              :    USE qs_environment_types,            ONLY: get_qs_env,&
      37              :                                               qs_environment_type
      38              :    USE qs_force_types,                  ONLY: qs_force_type
      39              :    USE qs_kind_types,                   ONLY: qs_kind_type
      40              :    USE qs_neighbor_list_types,          ONLY: get_neighbor_list_set_p,&
      41              :                                               neighbor_list_set_p_type
      42              :    USE virial_types,                    ONLY: virial_type
      43              : #include "./base/base_uses.f90"
      44              : 
      45              :    IMPLICIT NONE
      46              : 
      47              :    PRIVATE
      48              : 
      49              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'soc_pseudopotential_methods'
      50              : 
      51              :    PUBLIC :: V_SOC_xyz_from_pseudopotential, &
      52              :              remove_soc_outside_energy_window_ao, &
      53              :              remove_soc_outside_energy_window_mo
      54              : 
      55              : CONTAINS
      56              : 
      57              : ! **************************************************************************************************
      58              : !> \brief V^SOC_µν^(α),R = ħ/2 < ϕ_µ cell O | sum_ℓ ΔV_ℓ^SO(r,r') L^(α) | ϕ_ν cell R>, α = x,y,z
      59              : !>        see Hartwigsen, Goedecker, Hutter, Eq.(18), (19) (doi.org/10.1103/PhysRevB.58.3641)
      60              : !>        Caution: V^SOC_µν^(α) is purely imaginary and Hermitian; V^SOC_µν^(α) is stored as real
      61              : !>                 dbcsr matrix mat_V_SOC_xyz without symmetry; V^SOC_µν^(α) is stored without
      62              : !>                 the imaginary unit, i.e. mat_V_SOC_xyz is real and antisymmetric
      63              : !> \param qs_env ...
      64              : !> \param mat_V_SOC_xyz ...
      65              : !> \par History
      66              : !>    * 09.2023 created
      67              : !> \author Jan Wilhelm
      68              : ! **************************************************************************************************
      69           28 :    SUBROUTINE V_SOC_xyz_from_pseudopotential(qs_env, mat_V_SOC_xyz)
      70              :       TYPE(qs_environment_type), POINTER                 :: qs_env
      71              :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: mat_V_SOC_xyz
      72              : 
      73              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'V_SOC_xyz_from_pseudopotential'
      74              : 
      75              :       INTEGER                                            :: handle, img, nder, nimages, xyz
      76           14 :       INTEGER, DIMENSION(:, :, :), POINTER               :: cell_to_index
      77              :       LOGICAL                                            :: calculate_forces, do_kp, do_symmetric, &
      78              :                                                             use_virial
      79              :       REAL(KIND=dp)                                      :: eps_ppnl
      80           14 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      81           14 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: mat_l, mat_l_nosym, mat_pot_dummy, &
      82           14 :                                                             matrix_dummy, matrix_s
      83              :       TYPE(dft_control_type), POINTER                    :: dft_control
      84              :       TYPE(kpoint_type), POINTER                         :: kpoints
      85              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
      86           14 :          POINTER                                         :: sab_orb, sap_ppnl
      87           14 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      88           14 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
      89           14 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
      90              :       TYPE(virial_type), POINTER                         :: virial
      91              : 
      92           14 :       CALL timeset(routineN, handle)
      93              : 
      94           14 :       NULLIFY (qs_kind_set, dft_control, sab_orb, sap_ppnl, particle_set, atomic_kind_set, &
      95           14 :                cell_to_index)
      96              :       CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, dft_control=dft_control, &
      97              :                       matrix_s_kp=matrix_s, kpoints=kpoints, atomic_kind_set=atomic_kind_set, &
      98           14 :                       particle_set=particle_set, sab_orb=sab_orb, sap_ppnl=sap_ppnl)
      99              : 
     100           14 :       eps_ppnl = dft_control%qs_control%eps_ppnl
     101           14 :       nimages = dft_control%nimages
     102           14 :       do_kp = (nimages > 1)
     103           14 :       CALL get_neighbor_list_set_p(neighbor_list_sets=sab_orb, symmetric=do_symmetric)
     104           14 :       CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
     105              : 
     106           14 :       NULLIFY (mat_l, mat_pot_dummy)
     107           14 :       CALL dbcsr_allocate_matrix_set(mat_l, 3, nimages)
     108           56 :       DO xyz = 1, 3
     109          614 :          DO img = 1, nimages
     110          558 :             ALLOCATE (mat_l(xyz, img)%matrix)
     111              :             CALL dbcsr_create(mat_l(xyz, img)%matrix, template=matrix_s(1, 1)%matrix, &
     112          558 :                               matrix_type=dbcsr_type_antisymmetric)
     113          558 :             CALL cp_dbcsr_alloc_block_from_nbl(mat_l(xyz, img)%matrix, sab_orb)
     114          600 :             CALL dbcsr_set(mat_l(xyz, img)%matrix, 0.0_dp)
     115              :          END DO
     116              :       END DO
     117              : 
     118              :       ! get mat_l; the next CPASSERT fails if the atoms do not have any SOC parameters, i.e.
     119              :       ! SOC is zero and one should not activate the SOC section
     120           14 :       CPASSERT(ASSOCIATED(sap_ppnl))
     121           14 :       nder = 0
     122           14 :       use_virial = .FALSE.
     123           14 :       calculate_forces = .FALSE.
     124              : 
     125           14 :       NULLIFY (mat_pot_dummy)
     126           14 :       CALL dbcsr_allocate_matrix_set(mat_pot_dummy, 1, nimages)
     127          200 :       DO img = 1, nimages
     128          186 :          ALLOCATE (mat_pot_dummy(1, img)%matrix)
     129          186 :          CALL dbcsr_create(mat_pot_dummy(1, img)%matrix, template=matrix_s(1, 1)%matrix)
     130          186 :          CALL cp_dbcsr_alloc_block_from_nbl(mat_pot_dummy(1, img)%matrix, sab_orb)
     131          200 :          CALL dbcsr_set(mat_pot_dummy(1, img)%matrix, 0.0_dp)
     132              :       END DO
     133              : 
     134              :       CALL build_core_ppnl(mat_pot_dummy, matrix_dummy, force, virial, &
     135              :                            calculate_forces, use_virial, nder, &
     136              :                            qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, &
     137              :                            eps_ppnl, nimages=nimages, cell_to_index=cell_to_index, &
     138           14 :                            basis_type="ORB", matrix_l=mat_l)
     139              : 
     140           14 :       NULLIFY (mat_l_nosym)
     141           14 :       CALL dbcsr_allocate_matrix_set(mat_l_nosym, 3, nimages)
     142           56 :       DO xyz = 1, 3
     143          614 :          DO img = 1, nimages
     144              : 
     145          558 :             ALLOCATE (mat_l_nosym(xyz, img)%matrix)
     146          600 :             IF (do_kp) THEN
     147              :                CALL dbcsr_create(mat_l_nosym(xyz, img)%matrix, template=matrix_s(1, 1)%matrix, &
     148          534 :                                  matrix_type=dbcsr_type_antisymmetric)
     149          534 :                CALL dbcsr_copy(mat_l_nosym(xyz, img)%matrix, mat_l(xyz, img)%matrix)
     150              :             ELSE
     151              :                CALL dbcsr_create(mat_l_nosym(xyz, img)%matrix, template=matrix_s(1, 1)%matrix, &
     152           24 :                                  matrix_type=dbcsr_type_no_symmetry)
     153           24 :                CALL dbcsr_desymmetrize(mat_l(xyz, img)%matrix, mat_l_nosym(xyz, img)%matrix)
     154              :             END IF
     155              : 
     156              :          END DO
     157              :       END DO
     158              : 
     159           14 :       NULLIFY (mat_V_SOC_xyz)
     160           14 :       CALL dbcsr_allocate_matrix_set(mat_V_SOC_xyz, 3, nimages)
     161           56 :       DO xyz = 1, 3
     162          614 :          DO img = 1, nimages
     163          558 :             ALLOCATE (mat_V_SOC_xyz(xyz, img)%matrix)
     164          558 :             IF (do_kp) THEN
     165              :                ! mat_V_SOC_xyz^R with neighbor cell R actually has no symmetry
     166              :                ! mat_V_SOC_xyz^R_µν = mat_V_SOC_xyz^R_νµ* (the actual symmetry is
     167              :                ! mat_V_SOC_xyz^R_µν = mat_V_SOC_xyz^-R_νµ* )  but rskp_transform
     168              :                ! for mat_V_SOC_xyz^R -> mat_V_SOC_xyz(k) requires symmetry...
     169              :                CALL dbcsr_create(mat_V_SOC_xyz(xyz, img)%matrix, template=matrix_s(1, 1)%matrix, &
     170          534 :                                  matrix_type=dbcsr_type_antisymmetric)
     171              :             ELSE
     172              :                CALL dbcsr_create(mat_V_SOC_xyz(xyz, img)%matrix, template=matrix_s(1, 1)%matrix, &
     173           24 :                                  matrix_type=dbcsr_type_no_symmetry)
     174              :             END IF
     175          558 :             CALL cp_dbcsr_alloc_block_from_nbl(mat_V_SOC_xyz(xyz, img)%matrix, sab_orb)
     176              :             ! factor 0.5 from ħ/2 prefactor
     177              :             CALL dbcsr_add(mat_V_SOC_xyz(xyz, img)%matrix, mat_l_nosym(xyz, img)%matrix, &
     178          600 :                            0.0_dp, 0.5_dp)
     179              :          END DO
     180              :       END DO
     181              : 
     182           14 :       CALL dbcsr_deallocate_matrix_set(mat_pot_dummy)
     183           14 :       CALL dbcsr_deallocate_matrix_set(mat_l_nosym)
     184           14 :       CALL dbcsr_deallocate_matrix_set(mat_l)
     185              : 
     186           14 :       CALL timestop(handle)
     187              : 
     188           14 :    END SUBROUTINE V_SOC_xyz_from_pseudopotential
     189              : 
     190              : ! **************************************************************************************************
     191              : !> \brief Remove SOC outside of energy window (otherwise, numerical problems arise
     192              : !>        because energetically low semicore states and energetically very high
     193              : !>        unbound states couple to the states around the Fermi level).
     194              : !>        This routine is for mat_V_SOC_xyz being in the atomic-orbital (ao) basis.
     195              : !> \param mat_V_SOC_xyz ...
     196              : !> \param e_win ...
     197              : !> \param fm_mo_coeff ...
     198              : !> \param homo ...
     199              : !> \param eigenval ...
     200              : !> \param fm_s ...
     201              : ! **************************************************************************************************
     202            0 :    SUBROUTINE remove_soc_outside_energy_window_ao(mat_V_SOC_xyz, e_win, fm_mo_coeff, homo, &
     203            0 :                                                   eigenval, fm_s)
     204              :       TYPE(dbcsr_p_type), DIMENSION(:)                   :: mat_V_SOC_xyz
     205              :       REAL(KIND=dp)                                      :: e_win
     206              :       TYPE(cp_fm_type)                                   :: fm_mo_coeff
     207              :       INTEGER                                            :: homo
     208              :       REAL(KIND=dp), DIMENSION(:)                        :: eigenval
     209              :       TYPE(cp_fm_type)                                   :: fm_s
     210              : 
     211              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'remove_soc_outside_energy_window_ao'
     212              : 
     213              :       INTEGER                                            :: handle, i_glob, iiB, j_glob, jjB, nao, &
     214              :                                                             ncol_local, nrow_local, xyz
     215            0 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
     216              :       REAL(KIND=dp)                                      :: E_HOMO, E_i, E_j, E_LUMO
     217              :       TYPE(cp_fm_type)                                   :: fm_V_ao, fm_V_mo, fm_work
     218              : 
     219            0 :       CALL timeset(routineN, handle)
     220              : 
     221            0 :       CALL cp_fm_create(fm_work, fm_s%matrix_struct)
     222            0 :       CALL cp_fm_create(fm_V_ao, fm_s%matrix_struct)
     223            0 :       CALL cp_fm_create(fm_V_mo, fm_s%matrix_struct)
     224              : 
     225              :       CALL cp_fm_get_info(matrix=fm_s, &
     226              :                           nrow_local=nrow_local, &
     227              :                           ncol_local=ncol_local, &
     228              :                           row_indices=row_indices, &
     229            0 :                           col_indices=col_indices)
     230              : 
     231            0 :       nao = SIZE(eigenval)
     232              : 
     233            0 :       E_HOMO = eigenval(homo)
     234            0 :       E_LUMO = eigenval(homo + 1)
     235              : 
     236            0 :       DO xyz = 1, 3
     237              : 
     238            0 :          CALL copy_dbcsr_to_fm(mat_V_SOC_xyz(xyz)%matrix, fm_V_ao)
     239              : 
     240              :          ! V_MO = C^T*V_AO*C
     241              :          CALL parallel_gemm(transa="N", transb="N", m=nao, n=nao, k=nao, alpha=1.0_dp, &
     242            0 :                             matrix_a=fm_V_ao, matrix_b=fm_mo_coeff, beta=0.0_dp, matrix_c=fm_work)
     243              : 
     244              :          CALL parallel_gemm(transa="T", transb="N", m=nao, n=nao, k=nao, alpha=1.0_dp, &
     245            0 :                             matrix_a=fm_mo_coeff, matrix_b=fm_work, beta=0.0_dp, matrix_c=fm_V_mo)
     246              : 
     247            0 :          DO jjB = 1, ncol_local
     248            0 :             j_glob = col_indices(jjB)
     249            0 :             DO iiB = 1, nrow_local
     250            0 :                i_glob = row_indices(iiB)
     251              : 
     252            0 :                E_i = eigenval(i_glob)
     253            0 :                E_j = eigenval(j_glob)
     254              : 
     255              :                IF (E_i < E_HOMO - 0.5_dp*e_win .OR. E_i > E_LUMO + 0.5_dp*e_win .OR. &
     256            0 :                    E_j < E_HOMO - 0.5_dp*e_win .OR. E_j > E_LUMO + 0.5_dp*e_win) THEN
     257            0 :                   fm_V_mo%local_data(iiB, jjB) = 0.0_dp
     258              :                END IF
     259              : 
     260              :             END DO
     261              :          END DO
     262              : 
     263              :          ! V_AO = S*C*V_MO*C^T*S
     264              :          CALL parallel_gemm(transa="N", transb="T", m=nao, n=nao, k=nao, alpha=1.0_dp, &
     265            0 :                             matrix_a=fm_V_mo, matrix_b=fm_mo_coeff, beta=0.0_dp, matrix_c=fm_work)
     266              : 
     267              :          CALL parallel_gemm(transa="N", transb="N", m=nao, n=nao, k=nao, alpha=1.0_dp, &
     268            0 :                             matrix_a=fm_mo_coeff, matrix_b=fm_work, beta=0.0_dp, matrix_c=fm_V_ao)
     269              : 
     270              :          CALL parallel_gemm(transa="N", transb="N", m=nao, n=nao, k=nao, alpha=1.0_dp, &
     271            0 :                             matrix_a=fm_s, matrix_b=fm_V_ao, beta=0.0_dp, matrix_c=fm_work)
     272              : 
     273              :          CALL parallel_gemm(transa="N", transb="N", m=nao, n=nao, k=nao, alpha=1.0_dp, &
     274            0 :                             matrix_a=fm_work, matrix_b=fm_s, beta=0.0_dp, matrix_c=fm_V_ao)
     275              : 
     276            0 :          CALL copy_fm_to_dbcsr(fm_V_ao, mat_V_SOC_xyz(xyz)%matrix)
     277              : 
     278              :       END DO
     279              : 
     280            0 :       CALL cp_fm_release(fm_work)
     281            0 :       CALL cp_fm_release(fm_V_ao)
     282            0 :       CALL cp_fm_release(fm_V_mo)
     283              : 
     284            0 :       CALL timestop(handle)
     285              : 
     286            0 :    END SUBROUTINE remove_soc_outside_energy_window_ao
     287              : 
     288              : ! **************************************************************************************************
     289              : !> \brief ...
     290              : !> \param cfm_ks_spinor ...
     291              : !> \param e_win ...
     292              : !> \param eigenval ...
     293              : !> \param E_HOMO ...
     294              : !> \param E_LUMO ...
     295              : ! **************************************************************************************************
     296          336 :    SUBROUTINE remove_soc_outside_energy_window_mo(cfm_ks_spinor, e_win, eigenval, E_HOMO, E_LUMO)
     297              :       TYPE(cp_cfm_type)                                  :: cfm_ks_spinor
     298              :       REAL(KIND=dp)                                      :: e_win
     299              :       REAL(KIND=dp), DIMENSION(:)                        :: eigenval
     300              :       REAL(KIND=dp)                                      :: E_HOMO, E_LUMO
     301              : 
     302              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'remove_soc_outside_energy_window_mo'
     303              : 
     304              :       INTEGER                                            :: handle, i_glob, iiB, j_glob, jjB, &
     305              :                                                             ncol_global, ncol_local, nrow_global, &
     306              :                                                             nrow_local
     307          336 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
     308              :       REAL(KIND=dp)                                      :: E_i, E_j
     309              : 
     310              :       ! Remove SOC outside of energy window (otherwise, numerical problems arise
     311              :       ! because energetically low semicore states and energetically very high
     312              :       ! unbound states couple to the states around the Fermi level).
     313              :       ! This routine is for cfm_ks_spinor being in the molecular-orbital (mo) with
     314              :       ! corresponding eigenvalues "eigenval".
     315              : 
     316          336 :       CALL timeset(routineN, handle)
     317              : 
     318              :       CALL cp_cfm_get_info(matrix=cfm_ks_spinor, &
     319              :                            nrow_global=nrow_global, &
     320              :                            ncol_global=ncol_global, &
     321              :                            nrow_local=nrow_local, &
     322              :                            ncol_local=ncol_local, &
     323              :                            row_indices=row_indices, &
     324          336 :                            col_indices=col_indices)
     325              : 
     326          336 :       CPASSERT(nrow_global == SIZE(eigenval))
     327          336 :       CPASSERT(ncol_global == SIZE(eigenval))
     328              : 
     329         7664 :       DO jjB = 1, ncol_local
     330         7328 :          j_glob = col_indices(jjB)
     331        91536 :          DO iiB = 1, nrow_local
     332        83872 :             i_glob = row_indices(iiB)
     333              : 
     334        83872 :             E_i = eigenval(i_glob)
     335        83872 :             E_j = eigenval(j_glob)
     336              : 
     337              :             IF (E_i < E_HOMO - 0.5_dp*e_win .OR. E_i > E_LUMO + 0.5_dp*e_win .OR. &
     338        91200 :                 E_j < E_HOMO - 0.5_dp*e_win .OR. E_j > E_LUMO + 0.5_dp*e_win) THEN
     339        61340 :                cfm_ks_spinor%local_data(iiB, jjB) = 0.0_dp
     340              :             END IF
     341              : 
     342              :          END DO
     343              :       END DO
     344              : 
     345          336 :       CALL timestop(handle)
     346              : 
     347          336 :    END SUBROUTINE remove_soc_outside_energy_window_mo
     348              : 
     349              : END MODULE soc_pseudopotential_methods
        

Generated by: LCOV version 2.0-1