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_create,&
12 : cp_cfm_get_info,&
13 : cp_cfm_set_all,&
14 : cp_cfm_to_cfm,&
15 : cp_cfm_type
16 : USE cp_control_types, ONLY: dft_control_type
17 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
18 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
19 : copy_fm_to_dbcsr,&
20 : dbcsr_allocate_matrix_set,&
21 : dbcsr_deallocate_matrix_set
22 : USE cp_fm_struct, ONLY: cp_fm_struct_type
23 : USE cp_fm_types, ONLY: cp_fm_create,&
24 : cp_fm_get_info,&
25 : cp_fm_release,&
26 : cp_fm_type
27 : USE dbcsr_api, ONLY: dbcsr_add,&
28 : dbcsr_create,&
29 : dbcsr_desymmetrize,&
30 : dbcsr_p_type,&
31 : dbcsr_set,&
32 : dbcsr_type_antisymmetric,&
33 : dbcsr_type_no_symmetry
34 : USE kinds, ONLY: dp
35 : USE mathconstants, ONLY: gaussi,&
36 : z_one,&
37 : z_zero
38 : USE parallel_gemm_api, ONLY: parallel_gemm
39 : USE particle_types, ONLY: particle_type
40 : USE qs_environment_types, ONLY: get_qs_env,&
41 : qs_environment_type
42 : USE qs_force_types, ONLY: qs_force_type
43 : USE qs_kind_types, ONLY: qs_kind_type
44 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
45 : USE soc_pseudopotential_utils, ONLY: add_dbcsr_submat,&
46 : add_fm_submat,&
47 : create_cfm_double
48 : USE virial_types, ONLY: virial_type
49 : #include "./base/base_uses.f90"
50 :
51 : IMPLICIT NONE
52 :
53 : PRIVATE
54 :
55 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'soc_pseudopotential_methods'
56 :
57 : PUBLIC :: V_SOC_xyz_from_pseudopotential, H_KS_spinor, remove_soc_outside_energy_window_ao, &
58 : remove_soc_outside_energy_window_mo
59 :
60 : CONTAINS
61 :
62 : ! **************************************************************************************************
63 : !> \brief Compute V^SOC_µν^(α) = ħ/2 < ϕ_µ | sum_ℓ ΔV_ℓ^SO(r,r') L^(α) | ϕ_ν >, α = x, y, z, see
64 : !> Hartwigsen, Goedecker, Hutter, Eq.(18), (19) (doi.org/10.1103/PhysRevB.58.3641)
65 : !> Caution: V^SOC_µν^(α) is purely imaginary and Hermitian; V^SOC_µν^(α) is stored as real
66 : !> dbcsr matrix mat_V_SOC_xyz without symmetry; V^SOC_µν^(α) is stored without
67 : !> the imaginary unit, i.e. mat_V_SOC_xyz is real and antisymmetric
68 : !> \param qs_env ...
69 : !> \param mat_V_SOC_xyz ...
70 : !> \par History
71 : !> * 09.2023 created
72 : ! **************************************************************************************************
73 18 : SUBROUTINE V_SOC_xyz_from_pseudopotential(qs_env, mat_V_SOC_xyz)
74 : TYPE(qs_environment_type), POINTER :: qs_env
75 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_V_SOC_xyz
76 :
77 : CHARACTER(LEN=*), PARAMETER :: routineN = 'V_SOC_xyz_from_pseudopotential'
78 :
79 : INTEGER :: handle, nder, xyz
80 : LOGICAL :: calculate_forces, use_virial
81 : REAL(KIND=dp) :: eps_ppnl
82 18 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
83 18 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
84 18 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_l, mat_l_nosym, mat_pot_dummy, &
85 18 : matrix_dummy
86 : TYPE(dft_control_type), POINTER :: dft_control
87 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
88 18 : POINTER :: sab_orb, sap_ppnl
89 18 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
90 18 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
91 18 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
92 : TYPE(virial_type), POINTER :: virial
93 :
94 18 : CALL timeset(routineN, handle)
95 :
96 18 : NULLIFY (qs_kind_set, dft_control, sab_orb, sap_ppnl, particle_set, atomic_kind_set)
97 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, dft_control=dft_control, &
98 : matrix_s=matrix_s, atomic_kind_set=atomic_kind_set, &
99 18 : particle_set=particle_set, sab_orb=sab_orb, sap_ppnl=sap_ppnl)
100 :
101 18 : eps_ppnl = dft_control%qs_control%eps_ppnl
102 :
103 18 : NULLIFY (mat_l)
104 18 : CALL dbcsr_allocate_matrix_set(mat_l, 3, 1)
105 72 : DO xyz = 1, 3
106 54 : ALLOCATE (mat_l(xyz, 1)%matrix)
107 : CALL dbcsr_create(mat_l(xyz, 1)%matrix, template=matrix_s(1)%matrix, &
108 54 : matrix_type=dbcsr_type_antisymmetric)
109 54 : CALL cp_dbcsr_alloc_block_from_nbl(mat_l(xyz, 1)%matrix, sab_orb)
110 72 : CALL dbcsr_set(mat_l(xyz, 1)%matrix, 0.0_dp)
111 : END DO
112 :
113 : ! get mat_l
114 18 : CPASSERT(ASSOCIATED(sap_ppnl))
115 18 : nder = 0
116 18 : use_virial = .FALSE.
117 18 : calculate_forces = .FALSE.
118 :
119 18 : NULLIFY (mat_pot_dummy)
120 18 : CALL dbcsr_allocate_matrix_set(mat_pot_dummy, 1, 1)
121 18 : ALLOCATE (mat_pot_dummy(1, 1)%matrix)
122 18 : CALL dbcsr_create(mat_pot_dummy(1, 1)%matrix, template=matrix_s(1)%matrix)
123 18 : CALL cp_dbcsr_alloc_block_from_nbl(mat_pot_dummy(1, 1)%matrix, sab_orb)
124 18 : CALL dbcsr_set(mat_pot_dummy(1, 1)%matrix, 0.0_dp)
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 18 : eps_ppnl, nimages=1, basis_type="ORB", matrix_l=mat_l)
130 :
131 18 : NULLIFY (mat_l_nosym)
132 18 : CALL dbcsr_allocate_matrix_set(mat_l_nosym, 3, 1)
133 72 : DO xyz = 1, 3
134 54 : ALLOCATE (mat_l_nosym(xyz, 1)%matrix)
135 : CALL dbcsr_create(mat_l_nosym(xyz, 1)%matrix, template=matrix_s(1)%matrix, &
136 54 : matrix_type=dbcsr_type_no_symmetry)
137 72 : CALL dbcsr_desymmetrize(mat_l(xyz, 1)%matrix, mat_l_nosym(xyz, 1)%matrix)
138 :
139 : END DO
140 :
141 18 : NULLIFY (mat_V_SOC_xyz)
142 18 : CALL dbcsr_allocate_matrix_set(mat_V_SOC_xyz, 3, 1)
143 72 : DO xyz = 1, 3
144 54 : ALLOCATE (mat_V_SOC_xyz(xyz, 1)%matrix)
145 : CALL dbcsr_create(mat_V_SOC_xyz(xyz, 1)%matrix, template=matrix_s(1)%matrix, &
146 54 : matrix_type=dbcsr_type_no_symmetry)
147 54 : CALL cp_dbcsr_alloc_block_from_nbl(mat_V_SOC_xyz(xyz, 1)%matrix, sab_orb)
148 : ! factor 0.5 from ħ/2 prefactor
149 72 : CALL dbcsr_add(mat_V_SOC_xyz(xyz, 1)%matrix, mat_l_nosym(xyz, 1)%matrix, 0.0_dp, 0.5_dp)
150 : END DO
151 :
152 18 : CALL dbcsr_deallocate_matrix_set(mat_pot_dummy)
153 18 : CALL dbcsr_deallocate_matrix_set(mat_l_nosym)
154 18 : CALL dbcsr_deallocate_matrix_set(mat_l)
155 :
156 18 : CALL timestop(handle)
157 :
158 18 : END SUBROUTINE V_SOC_xyz_from_pseudopotential
159 :
160 : ! **************************************************************************************************
161 : !> \brief Spinor KS-matrix H_µν,σσ' = h_µν,σ*δ_σσ' + sum_α V^SOC_µν^(α)*Pauli-matrix^(α)_σσ', see
162 : !> Hartwigsen, Goedecker, Hutter, Eq.(18) (doi.org/10.1103/PhysRevB.58.3641)
163 : !> \param cfm_ks_spinor_ao ...
164 : !> \param fm_ks ...
165 : !> \param n_spin ...
166 : !> \param mat_V_SOC_xyz ...
167 : !> \param cfm_s_double ...
168 : !> \param fm_s ...
169 : !> \param cfm_SOC_spinor_ao ...
170 : ! **************************************************************************************************
171 16 : SUBROUTINE H_KS_spinor(cfm_ks_spinor_ao, fm_ks, n_spin, mat_V_SOC_xyz, cfm_s_double, fm_s, &
172 : cfm_SOC_spinor_ao)
173 : TYPE(cp_cfm_type) :: cfm_ks_spinor_ao
174 : TYPE(cp_fm_type), DIMENSION(:) :: fm_ks
175 : INTEGER :: n_spin
176 : TYPE(dbcsr_p_type), DIMENSION(:) :: mat_V_SOC_xyz
177 : TYPE(cp_cfm_type), OPTIONAL :: cfm_s_double
178 : TYPE(cp_fm_type), OPTIONAL :: fm_s
179 : TYPE(cp_cfm_type), OPTIONAL :: cfm_SOC_spinor_ao
180 :
181 : CHARACTER(LEN=*), PARAMETER :: routineN = 'H_KS_spinor'
182 :
183 : INTEGER :: handle, nao, s
184 : TYPE(cp_fm_struct_type), POINTER :: str
185 :
186 16 : CALL timeset(routineN, handle)
187 :
188 16 : CALL cp_fm_get_info(fm_ks(1), nrow_global=nao)
189 :
190 16 : CALL create_cfm_double(fm_ks(1), cfm_ks_spinor_ao)
191 16 : CALL cp_cfm_set_all(cfm_ks_spinor_ao, z_zero)
192 :
193 16 : str => fm_ks(1)%matrix_struct
194 :
195 16 : s = nao + 1
196 :
197 16 : CALL add_dbcsr_submat(cfm_ks_spinor_ao, mat_V_SOC_xyz(1)%matrix, str, s, 1, z_one, .TRUE.)
198 16 : CALL add_dbcsr_submat(cfm_ks_spinor_ao, mat_V_SOC_xyz(2)%matrix, str, s, 1, gaussi, .TRUE.)
199 16 : CALL add_dbcsr_submat(cfm_ks_spinor_ao, mat_V_SOC_xyz(3)%matrix, str, 1, 1, z_one, .FALSE.)
200 16 : CALL add_dbcsr_submat(cfm_ks_spinor_ao, mat_V_SOC_xyz(3)%matrix, str, s, s, -z_one, .FALSE.)
201 :
202 16 : IF (PRESENT(cfm_SOC_spinor_ao)) THEN
203 16 : CALL cp_cfm_create(cfm_SOC_spinor_ao, cfm_ks_spinor_ao%matrix_struct)
204 16 : CALL cp_cfm_to_cfm(cfm_ks_spinor_ao, cfm_SOC_spinor_ao)
205 : END IF
206 :
207 16 : CALL add_fm_submat(cfm_ks_spinor_ao, fm_ks(1), 1, 1)
208 :
209 26 : SELECT CASE (n_spin)
210 : CASE (1)
211 10 : CALL add_fm_submat(cfm_ks_spinor_ao, fm_ks(1), s, s)
212 : CASE (2)
213 6 : CPASSERT(SIZE(fm_ks) == 2)
214 6 : CALL add_fm_submat(cfm_ks_spinor_ao, fm_ks(2), s, s)
215 : CASE DEFAULT
216 16 : CPABORT("We have either one or two spin channels.")
217 : END SELECT
218 :
219 16 : IF (PRESENT(cfm_s_double)) THEN
220 16 : CPASSERT(PRESENT(fm_s))
221 16 : CALL create_cfm_double(fm_s, cfm_s_double)
222 16 : CALL cp_cfm_set_all(cfm_s_double, z_zero)
223 16 : CALL add_fm_submat(cfm_s_double, fm_s, 1, 1)
224 16 : CALL add_fm_submat(cfm_s_double, fm_s, s, s)
225 : END IF
226 :
227 16 : CALL timestop(handle)
228 :
229 16 : END SUBROUTINE H_KS_spinor
230 :
231 : ! **************************************************************************************************
232 : !> \brief Remove SOC outside of energy window (otherwise, numerical problems arise
233 : !> because energetically low semicore states and energetically very high
234 : !> unbound states couple to the states around the Fermi level).
235 : !> This routine is for mat_V_SOC_xyz being in the atomic-orbital (ao) basis.
236 : !> \param mat_V_SOC_xyz ...
237 : !> \param e_win ...
238 : !> \param fm_mo_coeff ...
239 : !> \param homo ...
240 : !> \param eigenval ...
241 : !> \param fm_s ...
242 : ! **************************************************************************************************
243 0 : SUBROUTINE remove_soc_outside_energy_window_ao(mat_V_SOC_xyz, e_win, fm_mo_coeff, homo, &
244 0 : eigenval, fm_s)
245 : TYPE(dbcsr_p_type), DIMENSION(:) :: mat_V_SOC_xyz
246 : REAL(KIND=dp) :: e_win
247 : TYPE(cp_fm_type) :: fm_mo_coeff
248 : INTEGER :: homo
249 : REAL(KIND=dp), DIMENSION(:) :: eigenval
250 : TYPE(cp_fm_type) :: fm_s
251 :
252 : CHARACTER(LEN=*), PARAMETER :: routineN = 'remove_soc_outside_energy_window_ao'
253 :
254 : INTEGER :: handle, i_glob, iiB, j_glob, jjB, nao, &
255 : ncol_local, nrow_local, xyz
256 0 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
257 : REAL(KIND=dp) :: E_HOMO, E_i, E_j, E_LUMO
258 : TYPE(cp_fm_type) :: fm_V_ao, fm_V_mo, fm_work
259 :
260 0 : CALL timeset(routineN, handle)
261 :
262 0 : CALL cp_fm_create(fm_work, fm_s%matrix_struct)
263 0 : CALL cp_fm_create(fm_V_ao, fm_s%matrix_struct)
264 0 : CALL cp_fm_create(fm_V_mo, fm_s%matrix_struct)
265 :
266 : CALL cp_fm_get_info(matrix=fm_s, &
267 : nrow_local=nrow_local, &
268 : ncol_local=ncol_local, &
269 : row_indices=row_indices, &
270 0 : col_indices=col_indices)
271 :
272 0 : nao = SIZE(eigenval)
273 :
274 0 : E_HOMO = eigenval(homo)
275 0 : E_LUMO = eigenval(homo + 1)
276 :
277 0 : DO xyz = 1, 3
278 :
279 0 : CALL copy_dbcsr_to_fm(mat_V_SOC_xyz(xyz)%matrix, fm_V_ao)
280 :
281 : ! V_MO = C^T*V_AO*C
282 : CALL parallel_gemm(transa="N", transb="N", m=nao, n=nao, k=nao, alpha=1.0_dp, &
283 0 : matrix_a=fm_V_ao, matrix_b=fm_mo_coeff, beta=0.0_dp, matrix_c=fm_work)
284 :
285 : CALL parallel_gemm(transa="T", transb="N", m=nao, n=nao, k=nao, alpha=1.0_dp, &
286 0 : matrix_a=fm_mo_coeff, matrix_b=fm_work, beta=0.0_dp, matrix_c=fm_V_mo)
287 :
288 0 : DO jjB = 1, ncol_local
289 0 : j_glob = col_indices(jjB)
290 0 : DO iiB = 1, nrow_local
291 0 : i_glob = row_indices(iiB)
292 :
293 0 : E_i = eigenval(i_glob)
294 0 : E_j = eigenval(j_glob)
295 :
296 : IF (E_i < E_HOMO - 0.5_dp*e_win .OR. E_i > E_LUMO + 0.5_dp*e_win .OR. &
297 0 : E_j < E_HOMO - 0.5_dp*e_win .OR. E_j > E_LUMO + 0.5_dp*e_win) THEN
298 0 : fm_V_mo%local_data(iiB, jjB) = 0.0_dp
299 : END IF
300 :
301 : END DO
302 : END DO
303 :
304 : ! V_AO = S*C*V_MO*C^T*S
305 : CALL parallel_gemm(transa="N", transb="T", m=nao, n=nao, k=nao, alpha=1.0_dp, &
306 0 : matrix_a=fm_V_mo, matrix_b=fm_mo_coeff, beta=0.0_dp, matrix_c=fm_work)
307 :
308 : CALL parallel_gemm(transa="N", transb="N", m=nao, n=nao, k=nao, alpha=1.0_dp, &
309 0 : matrix_a=fm_mo_coeff, matrix_b=fm_work, beta=0.0_dp, matrix_c=fm_V_ao)
310 :
311 : CALL parallel_gemm(transa="N", transb="N", m=nao, n=nao, k=nao, alpha=1.0_dp, &
312 0 : matrix_a=fm_s, matrix_b=fm_V_ao, beta=0.0_dp, matrix_c=fm_work)
313 :
314 : CALL parallel_gemm(transa="N", transb="N", m=nao, n=nao, k=nao, alpha=1.0_dp, &
315 0 : matrix_a=fm_work, matrix_b=fm_s, beta=0.0_dp, matrix_c=fm_V_ao)
316 :
317 0 : CALL copy_fm_to_dbcsr(fm_V_ao, mat_V_SOC_xyz(xyz)%matrix)
318 :
319 : END DO
320 :
321 0 : CALL cp_fm_release(fm_work)
322 0 : CALL cp_fm_release(fm_V_ao)
323 0 : CALL cp_fm_release(fm_V_mo)
324 :
325 0 : CALL timestop(handle)
326 :
327 0 : END SUBROUTINE remove_soc_outside_energy_window_ao
328 :
329 : ! **************************************************************************************************
330 : !> \brief ...
331 : !> \param cfm_ks_spinor ...
332 : !> \param e_win ...
333 : !> \param eigenval ...
334 : !> \param E_HOMO ...
335 : !> \param E_LUMO ...
336 : ! **************************************************************************************************
337 56 : SUBROUTINE remove_soc_outside_energy_window_mo(cfm_ks_spinor, e_win, eigenval, E_HOMO, E_LUMO)
338 : TYPE(cp_cfm_type) :: cfm_ks_spinor
339 : REAL(KIND=dp) :: e_win
340 : REAL(KIND=dp), DIMENSION(:) :: eigenval
341 : REAL(KIND=dp) :: E_HOMO, E_LUMO
342 :
343 : CHARACTER(LEN=*), PARAMETER :: routineN = 'remove_soc_outside_energy_window_mo'
344 :
345 : INTEGER :: handle, i_glob, iiB, j_glob, jjB, &
346 : ncol_global, ncol_local, nrow_global, &
347 : nrow_local
348 56 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
349 : REAL(KIND=dp) :: E_i, E_j
350 :
351 : ! Remove SOC outside of energy window (otherwise, numerical problems arise
352 : ! because energetically low semicore states and energetically very high
353 : ! unbound states couple to the states around the Fermi level).
354 : ! This routine is for cfm_ks_spinor being in the molecular-orbital (mo) with
355 : ! corresponding eigenvalues "eigenval".
356 :
357 56 : CALL timeset(routineN, handle)
358 :
359 : CALL cp_cfm_get_info(matrix=cfm_ks_spinor, &
360 : nrow_global=nrow_global, &
361 : ncol_global=ncol_global, &
362 : nrow_local=nrow_local, &
363 : ncol_local=ncol_local, &
364 : row_indices=row_indices, &
365 56 : col_indices=col_indices)
366 :
367 56 : CPASSERT(nrow_global == SIZE(eigenval))
368 56 : CPASSERT(ncol_global == SIZE(eigenval))
369 :
370 2360 : DO jjB = 1, ncol_local
371 2304 : j_glob = col_indices(jjB)
372 54200 : DO iiB = 1, nrow_local
373 51840 : i_glob = row_indices(iiB)
374 :
375 51840 : E_i = eigenval(i_glob)
376 51840 : E_j = eigenval(j_glob)
377 :
378 : IF (E_i < E_HOMO - 0.5_dp*e_win .OR. E_i > E_LUMO + 0.5_dp*e_win .OR. &
379 54144 : E_j < E_HOMO - 0.5_dp*e_win .OR. E_j > E_LUMO + 0.5_dp*e_win) THEN
380 49848 : cfm_ks_spinor%local_data(iiB, jjB) = 0.0_dp
381 : END IF
382 :
383 : END DO
384 : END DO
385 :
386 56 : CALL timestop(handle)
387 :
388 56 : END SUBROUTINE remove_soc_outside_energy_window_mo
389 :
390 : END MODULE soc_pseudopotential_methods
|
