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 : !USE physcon, ONLY: evolt
37 : #include "./base/base_uses.f90"
38 :
39 : IMPLICIT NONE
40 :
41 : PRIVATE
42 :
43 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'soc_pseudopotential_methods'
44 :
45 : PUBLIC :: V_SOC_xyz_from_pseudopotential, &
46 : remove_soc_outside_energy_window_mo
47 :
48 : CONTAINS
49 :
50 : ! **************************************************************************************************
51 : !> \brief V^SOC_µν^(α),R = ħ/2 < ϕ_µ cell O | sum_ℓ ΔV_ℓ^SO(r,r') L^(α) | ϕ_ν cell R>, α = x,y,z
52 : !> see Hartwigsen, Goedecker, Hutter, Eq.(18), (19) (doi.org/10.1103/PhysRevB.58.3641)
53 : !> Caution: V^SOC_µν^(α) is purely imaginary and Hermitian; V^SOC_µν^(α) is stored as real
54 : !> dbcsr matrix mat_V_SOC_xyz without symmetry; V^SOC_µν^(α) is stored without
55 : !> the imaginary unit, i.e. mat_V_SOC_xyz is real and antisymmetric
56 : !> \param qs_env ...
57 : !> \param mat_V_SOC_xyz ...
58 : !> \par History
59 : !> * 09.2023 created
60 : !> \author Jan Wilhelm
61 : ! **************************************************************************************************
62 136 : SUBROUTINE V_SOC_xyz_from_pseudopotential(qs_env, mat_V_SOC_xyz)
63 : TYPE(qs_environment_type), POINTER :: qs_env
64 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_V_SOC_xyz
65 :
66 : CHARACTER(LEN=*), PARAMETER :: routineN = 'V_SOC_xyz_from_pseudopotential'
67 :
68 : INTEGER :: handle, img, nder, nimages, xyz
69 68 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
70 : LOGICAL :: calculate_forces, do_kp, do_symmetric, &
71 : use_virial
72 : REAL(KIND=dp) :: eps_ppnl
73 68 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
74 68 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_l, mat_l_nosym, mat_pot_dummy, &
75 68 : matrix_dummy, matrix_s
76 : TYPE(dft_control_type), POINTER :: dft_control
77 : TYPE(kpoint_type), POINTER :: kpoints
78 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
79 68 : POINTER :: sab_orb, sap_ppnl
80 68 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
81 68 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
82 68 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
83 : TYPE(virial_type), POINTER :: virial
84 :
85 68 : CALL timeset(routineN, handle)
86 :
87 68 : NULLIFY (qs_kind_set, dft_control, sab_orb, sap_ppnl, particle_set, atomic_kind_set, &
88 68 : cell_to_index)
89 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, dft_control=dft_control, &
90 : matrix_s_kp=matrix_s, kpoints=kpoints, atomic_kind_set=atomic_kind_set, &
91 68 : particle_set=particle_set, sab_orb=sab_orb, sap_ppnl=sap_ppnl, do_kpoints=do_kp)
92 :
93 68 : eps_ppnl = dft_control%qs_control%eps_ppnl
94 68 : nimages = dft_control%nimages
95 : ! Even a one-image k-point environment needs the triangular storage expected by
96 : ! rskp_transform. Desymmetrizing redistributes blocks away from their neighbour-list owners.
97 68 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_orb, symmetric=do_symmetric)
98 68 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
99 :
100 68 : NULLIFY (mat_l, mat_pot_dummy)
101 68 : CALL dbcsr_allocate_matrix_set(mat_l, 3, nimages)
102 272 : DO xyz = 1, 3
103 11342 : DO img = 1, nimages
104 11070 : ALLOCATE (mat_l(xyz, img)%matrix)
105 : CALL dbcsr_create(mat_l(xyz, img)%matrix, template=matrix_s(1, 1)%matrix, &
106 11070 : matrix_type=dbcsr_type_antisymmetric)
107 11070 : CALL cp_dbcsr_alloc_block_from_nbl(mat_l(xyz, img)%matrix, sab_orb)
108 11274 : CALL dbcsr_set(mat_l(xyz, img)%matrix, 0.0_dp)
109 : END DO
110 : END DO
111 :
112 : ! get mat_l; the next CPASSERT fails if the atoms do not have any SOC parameters, i.e.
113 : ! SOC is zero and one should not activate the SOC section
114 68 : CPASSERT(ASSOCIATED(sap_ppnl))
115 68 : nder = 0
116 68 : use_virial = .FALSE.
117 68 : calculate_forces = .FALSE.
118 :
119 68 : NULLIFY (mat_pot_dummy)
120 68 : CALL dbcsr_allocate_matrix_set(mat_pot_dummy, 1, nimages)
121 3758 : DO img = 1, nimages
122 3690 : ALLOCATE (mat_pot_dummy(1, img)%matrix)
123 3690 : CALL dbcsr_create(mat_pot_dummy(1, img)%matrix, template=matrix_s(1, 1)%matrix)
124 3690 : CALL cp_dbcsr_alloc_block_from_nbl(mat_pot_dummy(1, img)%matrix, sab_orb)
125 3758 : CALL dbcsr_set(mat_pot_dummy(1, img)%matrix, 0.0_dp)
126 : END DO
127 :
128 : CALL build_core_ppnl(mat_pot_dummy, matrix_dummy, force, virial, &
129 : calculate_forces, use_virial, nder, &
130 : qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, &
131 : eps_ppnl, nimages=nimages, cell_to_index=cell_to_index, &
132 68 : basis_type="ORB", matrix_l=mat_l)
133 :
134 68 : NULLIFY (mat_l_nosym)
135 68 : CALL dbcsr_allocate_matrix_set(mat_l_nosym, 3, nimages)
136 272 : DO xyz = 1, 3
137 11342 : DO img = 1, nimages
138 :
139 11070 : ALLOCATE (mat_l_nosym(xyz, img)%matrix)
140 11274 : IF (do_kp) THEN
141 : CALL dbcsr_create(mat_l_nosym(xyz, img)%matrix, template=matrix_s(1, 1)%matrix, &
142 11016 : matrix_type=dbcsr_type_antisymmetric)
143 11016 : CALL dbcsr_copy(mat_l_nosym(xyz, img)%matrix, mat_l(xyz, img)%matrix)
144 : ELSE
145 : CALL dbcsr_create(mat_l_nosym(xyz, img)%matrix, template=matrix_s(1, 1)%matrix, &
146 54 : matrix_type=dbcsr_type_no_symmetry)
147 54 : CALL dbcsr_desymmetrize(mat_l(xyz, img)%matrix, mat_l_nosym(xyz, img)%matrix)
148 : END IF
149 :
150 : END DO
151 : END DO
152 :
153 68 : NULLIFY (mat_V_SOC_xyz)
154 68 : CALL dbcsr_allocate_matrix_set(mat_V_SOC_xyz, 3, nimages)
155 272 : DO xyz = 1, 3
156 11342 : DO img = 1, nimages
157 11070 : ALLOCATE (mat_V_SOC_xyz(xyz, img)%matrix)
158 11070 : IF (do_kp) THEN
159 : ! mat_V_SOC_xyz^R with neighbor cell R actually has no symmetry
160 : ! mat_V_SOC_xyz^R_µν = mat_V_SOC_xyz^R_νµ* (the actual symmetry is
161 : ! mat_V_SOC_xyz^R_µν = mat_V_SOC_xyz^-R_νµ* ) but rskp_transform
162 : ! for mat_V_SOC_xyz^R -> mat_V_SOC_xyz(k) requires symmetry...
163 : CALL dbcsr_create(mat_V_SOC_xyz(xyz, img)%matrix, template=matrix_s(1, 1)%matrix, &
164 11016 : matrix_type=dbcsr_type_antisymmetric)
165 : ELSE
166 : CALL dbcsr_create(mat_V_SOC_xyz(xyz, img)%matrix, template=matrix_s(1, 1)%matrix, &
167 54 : matrix_type=dbcsr_type_no_symmetry)
168 : END IF
169 11070 : CALL cp_dbcsr_alloc_block_from_nbl(mat_V_SOC_xyz(xyz, img)%matrix, sab_orb)
170 : ! factor 0.5 from ħ/2 prefactor
171 : CALL dbcsr_add(mat_V_SOC_xyz(xyz, img)%matrix, mat_l_nosym(xyz, img)%matrix, &
172 11274 : 0.0_dp, 0.5_dp)
173 : END DO
174 : END DO
175 :
176 68 : CALL dbcsr_deallocate_matrix_set(mat_pot_dummy)
177 68 : CALL dbcsr_deallocate_matrix_set(mat_l_nosym)
178 68 : CALL dbcsr_deallocate_matrix_set(mat_l)
179 :
180 68 : CALL timestop(handle)
181 :
182 68 : END SUBROUTINE V_SOC_xyz_from_pseudopotential
183 :
184 : ! **************************************************************************************************
185 : !> \brief ...
186 : !> \param cfm_ks_spinor ...
187 : !> \param e_win_cbm ...
188 : !> \param temp_smear ...
189 : !> \param eigenval ...
190 : !> \param e_fermi ...
191 : ! **************************************************************************************************
192 768 : SUBROUTINE remove_soc_outside_energy_window_mo(cfm_ks_spinor, e_win_cbm, temp_smear, &
193 256 : eigenval, e_fermi)
194 : TYPE(cp_cfm_type) :: cfm_ks_spinor
195 : REAL(KIND=dp) :: e_win_cbm, temp_smear
196 : REAL(KIND=dp), DIMENSION(:) :: eigenval
197 : REAL(KIND=dp) :: e_fermi
198 :
199 : CHARACTER(LEN=*), PARAMETER :: routineN = 'remove_soc_outside_energy_window_mo'
200 :
201 : INTEGER :: handle, i_glob, iiB, j_glob, jjB, &
202 : ncol_global, ncol_local, nrow_global, &
203 : nrow_local
204 256 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
205 : REAL(KIND=dp) :: E_i, E_j, w_i, w_j, x_i, x_j
206 :
207 : !REAL(KIND=dp) :: E_HOMO, E_LUMO, e_fermi
208 :
209 : ! Remove SOC outside of energy window (otherwise, numerical problems arise
210 : ! because energetically low semicore states and energetically very high
211 : ! unbound states couple to the states around the Fermi level).
212 : ! This routine is for cfm_ks_spinor being in the molecular-orbital (mo) with
213 : ! corresponding eigenvalues "eigenval".
214 :
215 256 : CALL timeset(routineN, handle)
216 :
217 : CALL cp_cfm_get_info(matrix=cfm_ks_spinor, &
218 : nrow_global=nrow_global, &
219 : ncol_global=ncol_global, &
220 : nrow_local=nrow_local, &
221 : ncol_local=ncol_local, &
222 : row_indices=row_indices, &
223 256 : col_indices=col_indices)
224 :
225 256 : CPASSERT(nrow_global == SIZE(eigenval))
226 256 : CPASSERT(ncol_global == SIZE(eigenval))
227 :
228 : ! Apply a smooth energy window to the SOC coupling matrix. Each state gets a
229 : ! weight in [0,1]: 1 = full SOC, decaying to 0 outside the window via a
230 : ! Fermi-like function of width temp_smear. A matrix element coupling states i
231 : ! and j is scaled by SQRT(w_i*w_j), applying the damping symmetrically to both.
232 7424 : DO jjB = 1, ncol_local
233 7168 : j_glob = col_indices(jjB)
234 107776 : DO iiB = 1, nrow_local
235 100352 : i_glob = row_indices(iiB)
236 :
237 100352 : E_i = eigenval(i_glob)
238 100352 : E_j = eigenval(j_glob)
239 100352 : IF (E_i <= e_fermi) THEN
240 : !x_i = ABS(E_i - e_fermi) - 0.5_dp * e_win_vbm
241 : w_i = 1.0_dp
242 : ELSE
243 71680 : x_i = ABS(E_i - e_fermi) - e_win_cbm
244 71680 : w_i = 1.0_dp/(EXP(x_i/temp_smear) + 1.0_dp)
245 : END IF
246 :
247 100352 : IF (E_j <= e_fermi) THEN
248 : !x_j = ABS(E_j - e_fermi) - 0.5_dp * e_win_vbm
249 : w_j = 1.0_dp
250 : ELSE
251 71680 : x_j = ABS(E_j - e_fermi) - e_win_cbm
252 71680 : w_j = 1.0_dp/(EXP(x_j/temp_smear) + 1.0_dp)
253 : END IF
254 :
255 : cfm_ks_spinor%local_data(iiB, jjB) = &
256 107520 : cfm_ks_spinor%local_data(iiB, jjB)*SQRT(w_i*w_j)
257 :
258 : END DO
259 : END DO
260 :
261 256 : CALL timestop(handle)
262 :
263 256 : END SUBROUTINE remove_soc_outside_energy_window_mo
264 :
265 : END MODULE soc_pseudopotential_methods
|
