Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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 40 : 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 20 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
77 : LOGICAL :: calculate_forces, do_kp, do_symmetric, &
78 : use_virial
79 : REAL(KIND=dp) :: eps_ppnl
80 20 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
81 20 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_l, mat_l_nosym, mat_pot_dummy, &
82 20 : 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 20 : POINTER :: sab_orb, sap_ppnl
87 20 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
88 20 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
89 20 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
90 : TYPE(virial_type), POINTER :: virial
91 :
92 20 : CALL timeset(routineN, handle)
93 :
94 20 : NULLIFY (qs_kind_set, dft_control, sab_orb, sap_ppnl, particle_set, atomic_kind_set, &
95 20 : 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 20 : particle_set=particle_set, sab_orb=sab_orb, sap_ppnl=sap_ppnl)
99 :
100 20 : eps_ppnl = dft_control%qs_control%eps_ppnl
101 20 : nimages = dft_control%nimages
102 20 : do_kp = (nimages > 1)
103 20 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_orb, symmetric=do_symmetric)
104 20 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
105 :
106 20 : NULLIFY (mat_l, mat_pot_dummy)
107 20 : CALL dbcsr_allocate_matrix_set(mat_l, 3, nimages)
108 80 : DO xyz = 1, 3
109 560 : DO img = 1, nimages
110 480 : ALLOCATE (mat_l(xyz, img)%matrix)
111 : CALL dbcsr_create(mat_l(xyz, img)%matrix, template=matrix_s(1, 1)%matrix, &
112 480 : matrix_type=dbcsr_type_antisymmetric)
113 480 : CALL cp_dbcsr_alloc_block_from_nbl(mat_l(xyz, img)%matrix, sab_orb)
114 540 : 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 20 : CPASSERT(ASSOCIATED(sap_ppnl))
121 20 : nder = 0
122 20 : use_virial = .FALSE.
123 20 : calculate_forces = .FALSE.
124 :
125 20 : NULLIFY (mat_pot_dummy)
126 20 : CALL dbcsr_allocate_matrix_set(mat_pot_dummy, 1, nimages)
127 180 : DO img = 1, nimages
128 160 : ALLOCATE (mat_pot_dummy(1, img)%matrix)
129 160 : CALL dbcsr_create(mat_pot_dummy(1, img)%matrix, template=matrix_s(1, 1)%matrix)
130 160 : CALL cp_dbcsr_alloc_block_from_nbl(mat_pot_dummy(1, img)%matrix, sab_orb)
131 180 : 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 20 : basis_type="ORB", matrix_l=mat_l)
139 :
140 20 : NULLIFY (mat_l_nosym)
141 20 : CALL dbcsr_allocate_matrix_set(mat_l_nosym, 3, nimages)
142 80 : DO xyz = 1, 3
143 560 : DO img = 1, nimages
144 :
145 480 : ALLOCATE (mat_l_nosym(xyz, img)%matrix)
146 540 : IF (do_kp) THEN
147 : CALL dbcsr_create(mat_l_nosym(xyz, img)%matrix, template=matrix_s(1, 1)%matrix, &
148 438 : matrix_type=dbcsr_type_antisymmetric)
149 438 : 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 42 : matrix_type=dbcsr_type_no_symmetry)
153 42 : 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 20 : NULLIFY (mat_V_SOC_xyz)
160 20 : CALL dbcsr_allocate_matrix_set(mat_V_SOC_xyz, 3, nimages)
161 80 : DO xyz = 1, 3
162 560 : DO img = 1, nimages
163 480 : ALLOCATE (mat_V_SOC_xyz(xyz, img)%matrix)
164 480 : 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 438 : 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 42 : matrix_type=dbcsr_type_no_symmetry)
174 : END IF
175 480 : 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 540 : 0.0_dp, 0.5_dp)
179 : END DO
180 : END DO
181 :
182 20 : CALL dbcsr_deallocate_matrix_set(mat_pot_dummy)
183 20 : CALL dbcsr_deallocate_matrix_set(mat_l_nosym)
184 20 : CALL dbcsr_deallocate_matrix_set(mat_l)
185 :
186 20 : CALL timestop(handle)
187 :
188 20 : 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 480 : 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 480 : 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 480 : 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 480 : col_indices=col_indices)
325 :
326 480 : CPASSERT(nrow_global == SIZE(eigenval))
327 480 : CPASSERT(ncol_global == SIZE(eigenval))
328 :
329 9760 : DO jjB = 1, ncol_local
330 9280 : j_glob = col_indices(jjB)
331 102240 : DO iiB = 1, nrow_local
332 92480 : i_glob = row_indices(iiB)
333 :
334 92480 : E_i = eigenval(i_glob)
335 92480 : 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 101760 : E_j < E_HOMO - 0.5_dp*e_win .OR. E_j > E_LUMO + 0.5_dp*e_win) THEN
339 71992 : cfm_ks_spinor%local_data(iiB, jjB) = 0.0_dp
340 : END IF
341 :
342 : END DO
343 : END DO
344 :
345 480 : CALL timestop(handle)
346 :
347 480 : END SUBROUTINE remove_soc_outside_energy_window_mo
348 :
349 : END MODULE soc_pseudopotential_methods
|
