Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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: dbcsr_allocate_matrix_set,&
24 : dbcsr_deallocate_matrix_set
25 : USE kinds, ONLY: dp
26 : USE kpoint_types, ONLY: get_kpoint_info,&
27 : kpoint_type
28 : USE particle_types, ONLY: particle_type
29 : USE qs_environment_types, ONLY: get_qs_env,&
30 : qs_environment_type
31 : USE qs_force_types, ONLY: qs_force_type
32 : USE qs_kind_types, ONLY: qs_kind_type
33 : USE qs_neighbor_list_types, ONLY: get_neighbor_list_set_p,&
34 : neighbor_list_set_p_type
35 : USE virial_types, ONLY: virial_type
36 : #include "./base/base_uses.f90"
37 :
38 : IMPLICIT NONE
39 :
40 : PRIVATE
41 :
42 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'soc_pseudopotential_methods'
43 :
44 : PUBLIC :: V_SOC_xyz_from_pseudopotential, &
45 : remove_soc_outside_energy_window_mo
46 :
47 : CONTAINS
48 :
49 : ! **************************************************************************************************
50 : !> \brief V^SOC_µν^(α),R = ħ/2 < ϕ_µ cell O | sum_ℓ ΔV_ℓ^SO(r,r') L^(α) | ϕ_ν cell R>, α = x,y,z
51 : !> see Hartwigsen, Goedecker, Hutter, Eq.(18), (19) (doi.org/10.1103/PhysRevB.58.3641)
52 : !> Caution: V^SOC_µν^(α) is purely imaginary and Hermitian; V^SOC_µν^(α) is stored as real
53 : !> dbcsr matrix mat_V_SOC_xyz without symmetry; V^SOC_µν^(α) is stored without
54 : !> the imaginary unit, i.e. mat_V_SOC_xyz is real and antisymmetric
55 : !> \param qs_env ...
56 : !> \param mat_V_SOC_xyz ...
57 : !> \par History
58 : !> * 09.2023 created
59 : !> \author Jan Wilhelm
60 : ! **************************************************************************************************
61 28 : SUBROUTINE V_SOC_xyz_from_pseudopotential(qs_env, mat_V_SOC_xyz)
62 : TYPE(qs_environment_type), POINTER :: qs_env
63 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_V_SOC_xyz
64 :
65 : CHARACTER(LEN=*), PARAMETER :: routineN = 'V_SOC_xyz_from_pseudopotential'
66 :
67 : INTEGER :: handle, img, nder, nimages, xyz
68 14 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
69 : LOGICAL :: calculate_forces, do_kp, do_symmetric, &
70 : use_virial
71 : REAL(KIND=dp) :: eps_ppnl
72 14 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
73 14 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_l, mat_l_nosym, mat_pot_dummy, &
74 14 : matrix_dummy, matrix_s
75 : TYPE(dft_control_type), POINTER :: dft_control
76 : TYPE(kpoint_type), POINTER :: kpoints
77 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
78 14 : POINTER :: sab_orb, sap_ppnl
79 14 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
80 14 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
81 14 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
82 : TYPE(virial_type), POINTER :: virial
83 :
84 14 : CALL timeset(routineN, handle)
85 :
86 14 : NULLIFY (qs_kind_set, dft_control, sab_orb, sap_ppnl, particle_set, atomic_kind_set, &
87 14 : cell_to_index)
88 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, dft_control=dft_control, &
89 : matrix_s_kp=matrix_s, kpoints=kpoints, atomic_kind_set=atomic_kind_set, &
90 14 : particle_set=particle_set, sab_orb=sab_orb, sap_ppnl=sap_ppnl)
91 :
92 14 : eps_ppnl = dft_control%qs_control%eps_ppnl
93 14 : nimages = dft_control%nimages
94 14 : do_kp = (nimages > 1)
95 14 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_orb, symmetric=do_symmetric)
96 14 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
97 :
98 14 : NULLIFY (mat_l, mat_pot_dummy)
99 14 : CALL dbcsr_allocate_matrix_set(mat_l, 3, nimages)
100 56 : DO xyz = 1, 3
101 614 : DO img = 1, nimages
102 558 : ALLOCATE (mat_l(xyz, img)%matrix)
103 : CALL dbcsr_create(mat_l(xyz, img)%matrix, template=matrix_s(1, 1)%matrix, &
104 558 : matrix_type=dbcsr_type_antisymmetric)
105 558 : CALL cp_dbcsr_alloc_block_from_nbl(mat_l(xyz, img)%matrix, sab_orb)
106 600 : CALL dbcsr_set(mat_l(xyz, img)%matrix, 0.0_dp)
107 : END DO
108 : END DO
109 :
110 : ! get mat_l; the next CPASSERT fails if the atoms do not have any SOC parameters, i.e.
111 : ! SOC is zero and one should not activate the SOC section
112 14 : CPASSERT(ASSOCIATED(sap_ppnl))
113 14 : nder = 0
114 14 : use_virial = .FALSE.
115 14 : calculate_forces = .FALSE.
116 :
117 14 : NULLIFY (mat_pot_dummy)
118 14 : CALL dbcsr_allocate_matrix_set(mat_pot_dummy, 1, nimages)
119 200 : DO img = 1, nimages
120 186 : ALLOCATE (mat_pot_dummy(1, img)%matrix)
121 186 : CALL dbcsr_create(mat_pot_dummy(1, img)%matrix, template=matrix_s(1, 1)%matrix)
122 186 : CALL cp_dbcsr_alloc_block_from_nbl(mat_pot_dummy(1, img)%matrix, sab_orb)
123 200 : CALL dbcsr_set(mat_pot_dummy(1, img)%matrix, 0.0_dp)
124 : END DO
125 :
126 : CALL build_core_ppnl(mat_pot_dummy, matrix_dummy, force, virial, &
127 : calculate_forces, use_virial, nder, &
128 : qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, &
129 : eps_ppnl, nimages=nimages, cell_to_index=cell_to_index, &
130 14 : basis_type="ORB", matrix_l=mat_l)
131 :
132 14 : NULLIFY (mat_l_nosym)
133 14 : CALL dbcsr_allocate_matrix_set(mat_l_nosym, 3, nimages)
134 56 : DO xyz = 1, 3
135 614 : DO img = 1, nimages
136 :
137 558 : ALLOCATE (mat_l_nosym(xyz, img)%matrix)
138 600 : IF (do_kp) THEN
139 : CALL dbcsr_create(mat_l_nosym(xyz, img)%matrix, template=matrix_s(1, 1)%matrix, &
140 534 : matrix_type=dbcsr_type_antisymmetric)
141 534 : CALL dbcsr_copy(mat_l_nosym(xyz, img)%matrix, mat_l(xyz, img)%matrix)
142 : ELSE
143 : CALL dbcsr_create(mat_l_nosym(xyz, img)%matrix, template=matrix_s(1, 1)%matrix, &
144 24 : matrix_type=dbcsr_type_no_symmetry)
145 24 : CALL dbcsr_desymmetrize(mat_l(xyz, img)%matrix, mat_l_nosym(xyz, img)%matrix)
146 : END IF
147 :
148 : END DO
149 : END DO
150 :
151 14 : NULLIFY (mat_V_SOC_xyz)
152 14 : CALL dbcsr_allocate_matrix_set(mat_V_SOC_xyz, 3, nimages)
153 56 : DO xyz = 1, 3
154 614 : DO img = 1, nimages
155 558 : ALLOCATE (mat_V_SOC_xyz(xyz, img)%matrix)
156 558 : IF (do_kp) THEN
157 : ! mat_V_SOC_xyz^R with neighbor cell R actually has no symmetry
158 : ! mat_V_SOC_xyz^R_µν = mat_V_SOC_xyz^R_νµ* (the actual symmetry is
159 : ! mat_V_SOC_xyz^R_µν = mat_V_SOC_xyz^-R_νµ* ) but rskp_transform
160 : ! for mat_V_SOC_xyz^R -> mat_V_SOC_xyz(k) requires symmetry...
161 : CALL dbcsr_create(mat_V_SOC_xyz(xyz, img)%matrix, template=matrix_s(1, 1)%matrix, &
162 534 : matrix_type=dbcsr_type_antisymmetric)
163 : ELSE
164 : CALL dbcsr_create(mat_V_SOC_xyz(xyz, img)%matrix, template=matrix_s(1, 1)%matrix, &
165 24 : matrix_type=dbcsr_type_no_symmetry)
166 : END IF
167 558 : CALL cp_dbcsr_alloc_block_from_nbl(mat_V_SOC_xyz(xyz, img)%matrix, sab_orb)
168 : ! factor 0.5 from ħ/2 prefactor
169 : CALL dbcsr_add(mat_V_SOC_xyz(xyz, img)%matrix, mat_l_nosym(xyz, img)%matrix, &
170 600 : 0.0_dp, 0.5_dp)
171 : END DO
172 : END DO
173 :
174 14 : CALL dbcsr_deallocate_matrix_set(mat_pot_dummy)
175 14 : CALL dbcsr_deallocate_matrix_set(mat_l_nosym)
176 14 : CALL dbcsr_deallocate_matrix_set(mat_l)
177 :
178 14 : CALL timestop(handle)
179 :
180 14 : END SUBROUTINE V_SOC_xyz_from_pseudopotential
181 :
182 : ! **************************************************************************************************
183 : !> \brief ...
184 : !> \param cfm_ks_spinor ...
185 : !> \param e_win ...
186 : !> \param eigenval ...
187 : !> \param E_HOMO ...
188 : !> \param E_LUMO ...
189 : ! **************************************************************************************************
190 336 : SUBROUTINE remove_soc_outside_energy_window_mo(cfm_ks_spinor, e_win, eigenval, E_HOMO, E_LUMO)
191 : TYPE(cp_cfm_type) :: cfm_ks_spinor
192 : REAL(KIND=dp) :: e_win
193 : REAL(KIND=dp), DIMENSION(:) :: eigenval
194 : REAL(KIND=dp) :: E_HOMO, E_LUMO
195 :
196 : CHARACTER(LEN=*), PARAMETER :: routineN = 'remove_soc_outside_energy_window_mo'
197 :
198 : INTEGER :: handle, i_glob, iiB, j_glob, jjB, &
199 : ncol_global, ncol_local, nrow_global, &
200 : nrow_local
201 336 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
202 : REAL(KIND=dp) :: E_i, E_j
203 :
204 : ! Remove SOC outside of energy window (otherwise, numerical problems arise
205 : ! because energetically low semicore states and energetically very high
206 : ! unbound states couple to the states around the Fermi level).
207 : ! This routine is for cfm_ks_spinor being in the molecular-orbital (mo) with
208 : ! corresponding eigenvalues "eigenval".
209 :
210 336 : CALL timeset(routineN, handle)
211 :
212 : CALL cp_cfm_get_info(matrix=cfm_ks_spinor, &
213 : nrow_global=nrow_global, &
214 : ncol_global=ncol_global, &
215 : nrow_local=nrow_local, &
216 : ncol_local=ncol_local, &
217 : row_indices=row_indices, &
218 336 : col_indices=col_indices)
219 :
220 336 : CPASSERT(nrow_global == SIZE(eigenval))
221 336 : CPASSERT(ncol_global == SIZE(eigenval))
222 :
223 7664 : DO jjB = 1, ncol_local
224 7328 : j_glob = col_indices(jjB)
225 91536 : DO iiB = 1, nrow_local
226 83872 : i_glob = row_indices(iiB)
227 :
228 83872 : E_i = eigenval(i_glob)
229 83872 : E_j = eigenval(j_glob)
230 :
231 : IF (E_i < E_HOMO - 0.5_dp*e_win .OR. E_i > E_LUMO + 0.5_dp*e_win .OR. &
232 91200 : E_j < E_HOMO - 0.5_dp*e_win .OR. E_j > E_LUMO + 0.5_dp*e_win) THEN
233 61340 : cfm_ks_spinor%local_data(iiB, jjB) = 0.0_dp
234 : END IF
235 :
236 : END DO
237 : END DO
238 :
239 336 : CALL timestop(handle)
240 :
241 336 : END SUBROUTINE remove_soc_outside_energy_window_mo
242 :
243 : END MODULE soc_pseudopotential_methods
|
