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 48 : 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 24 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
70 : LOGICAL :: calculate_forces, do_kp, do_symmetric, &
71 : use_virial
72 : REAL(KIND=dp) :: eps_ppnl
73 24 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
74 24 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_l, mat_l_nosym, mat_pot_dummy, &
75 24 : 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 24 : POINTER :: sab_orb, sap_ppnl
80 24 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
81 24 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
82 24 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
83 : TYPE(virial_type), POINTER :: virial
84 :
85 24 : CALL timeset(routineN, handle)
86 :
87 24 : NULLIFY (qs_kind_set, dft_control, sab_orb, sap_ppnl, particle_set, atomic_kind_set, &
88 24 : 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 24 : particle_set=particle_set, sab_orb=sab_orb, sap_ppnl=sap_ppnl)
92 :
93 24 : eps_ppnl = dft_control%qs_control%eps_ppnl
94 24 : nimages = dft_control%nimages
95 24 : do_kp = (nimages > 1)
96 24 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_orb, symmetric=do_symmetric)
97 24 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
98 :
99 24 : NULLIFY (mat_l, mat_pot_dummy)
100 24 : CALL dbcsr_allocate_matrix_set(mat_l, 3, nimages)
101 96 : DO xyz = 1, 3
102 1704 : DO img = 1, nimages
103 1608 : ALLOCATE (mat_l(xyz, img)%matrix)
104 : CALL dbcsr_create(mat_l(xyz, img)%matrix, template=matrix_s(1, 1)%matrix, &
105 1608 : matrix_type=dbcsr_type_antisymmetric)
106 1608 : CALL cp_dbcsr_alloc_block_from_nbl(mat_l(xyz, img)%matrix, sab_orb)
107 1680 : CALL dbcsr_set(mat_l(xyz, img)%matrix, 0.0_dp)
108 : END DO
109 : END DO
110 :
111 : ! get mat_l; the next CPASSERT fails if the atoms do not have any SOC parameters, i.e.
112 : ! SOC is zero and one should not activate the SOC section
113 24 : CPASSERT(ASSOCIATED(sap_ppnl))
114 24 : nder = 0
115 24 : use_virial = .FALSE.
116 24 : calculate_forces = .FALSE.
117 :
118 24 : NULLIFY (mat_pot_dummy)
119 24 : CALL dbcsr_allocate_matrix_set(mat_pot_dummy, 1, nimages)
120 560 : DO img = 1, nimages
121 536 : ALLOCATE (mat_pot_dummy(1, img)%matrix)
122 536 : CALL dbcsr_create(mat_pot_dummy(1, img)%matrix, template=matrix_s(1, 1)%matrix)
123 536 : CALL cp_dbcsr_alloc_block_from_nbl(mat_pot_dummy(1, img)%matrix, sab_orb)
124 560 : CALL dbcsr_set(mat_pot_dummy(1, img)%matrix, 0.0_dp)
125 : END DO
126 :
127 : CALL build_core_ppnl(mat_pot_dummy, matrix_dummy, force, virial, &
128 : calculate_forces, use_virial, nder, &
129 : qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, &
130 : eps_ppnl, nimages=nimages, cell_to_index=cell_to_index, &
131 24 : basis_type="ORB", matrix_l=mat_l)
132 :
133 24 : NULLIFY (mat_l_nosym)
134 24 : CALL dbcsr_allocate_matrix_set(mat_l_nosym, 3, nimages)
135 96 : DO xyz = 1, 3
136 1704 : DO img = 1, nimages
137 :
138 1608 : ALLOCATE (mat_l_nosym(xyz, img)%matrix)
139 1680 : IF (do_kp) THEN
140 : CALL dbcsr_create(mat_l_nosym(xyz, img)%matrix, template=matrix_s(1, 1)%matrix, &
141 1584 : matrix_type=dbcsr_type_antisymmetric)
142 1584 : CALL dbcsr_copy(mat_l_nosym(xyz, img)%matrix, mat_l(xyz, img)%matrix)
143 : ELSE
144 : CALL dbcsr_create(mat_l_nosym(xyz, img)%matrix, template=matrix_s(1, 1)%matrix, &
145 24 : matrix_type=dbcsr_type_no_symmetry)
146 24 : CALL dbcsr_desymmetrize(mat_l(xyz, img)%matrix, mat_l_nosym(xyz, img)%matrix)
147 : END IF
148 :
149 : END DO
150 : END DO
151 :
152 24 : NULLIFY (mat_V_SOC_xyz)
153 24 : CALL dbcsr_allocate_matrix_set(mat_V_SOC_xyz, 3, nimages)
154 96 : DO xyz = 1, 3
155 1704 : DO img = 1, nimages
156 1608 : ALLOCATE (mat_V_SOC_xyz(xyz, img)%matrix)
157 1608 : IF (do_kp) THEN
158 : ! mat_V_SOC_xyz^R with neighbor cell R actually has no symmetry
159 : ! mat_V_SOC_xyz^R_µν = mat_V_SOC_xyz^R_νµ* (the actual symmetry is
160 : ! mat_V_SOC_xyz^R_µν = mat_V_SOC_xyz^-R_νµ* ) but rskp_transform
161 : ! for mat_V_SOC_xyz^R -> mat_V_SOC_xyz(k) requires symmetry...
162 : CALL dbcsr_create(mat_V_SOC_xyz(xyz, img)%matrix, template=matrix_s(1, 1)%matrix, &
163 1584 : matrix_type=dbcsr_type_antisymmetric)
164 : ELSE
165 : CALL dbcsr_create(mat_V_SOC_xyz(xyz, img)%matrix, template=matrix_s(1, 1)%matrix, &
166 24 : matrix_type=dbcsr_type_no_symmetry)
167 : END IF
168 1608 : CALL cp_dbcsr_alloc_block_from_nbl(mat_V_SOC_xyz(xyz, img)%matrix, sab_orb)
169 : ! factor 0.5 from ħ/2 prefactor
170 : CALL dbcsr_add(mat_V_SOC_xyz(xyz, img)%matrix, mat_l_nosym(xyz, img)%matrix, &
171 1680 : 0.0_dp, 0.5_dp)
172 : END DO
173 : END DO
174 :
175 24 : CALL dbcsr_deallocate_matrix_set(mat_pot_dummy)
176 24 : CALL dbcsr_deallocate_matrix_set(mat_l_nosym)
177 24 : CALL dbcsr_deallocate_matrix_set(mat_l)
178 :
179 24 : CALL timestop(handle)
180 :
181 24 : END SUBROUTINE V_SOC_xyz_from_pseudopotential
182 :
183 : ! **************************************************************************************************
184 : !> \brief ...
185 : !> \param cfm_ks_spinor ...
186 : !> \param e_win_cbm ...
187 : !> \param temp_smear ...
188 : !> \param eigenval ...
189 : !> \param e_fermi ...
190 : ! **************************************************************************************************
191 768 : SUBROUTINE remove_soc_outside_energy_window_mo(cfm_ks_spinor, e_win_cbm, temp_smear, &
192 256 : eigenval, e_fermi)
193 : TYPE(cp_cfm_type) :: cfm_ks_spinor
194 : REAL(KIND=dp) :: e_win_cbm, temp_smear
195 : REAL(KIND=dp), DIMENSION(:) :: eigenval
196 : REAL(KIND=dp) :: e_fermi
197 :
198 : CHARACTER(LEN=*), PARAMETER :: routineN = 'remove_soc_outside_energy_window_mo'
199 :
200 : INTEGER :: handle, i_glob, iiB, j_glob, jjB, &
201 : ncol_global, ncol_local, nrow_global, &
202 : nrow_local
203 256 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
204 : REAL(KIND=dp) :: E_i, E_j, w_i, w_j, x_i, x_j
205 :
206 : !REAL(KIND=dp) :: E_HOMO, E_LUMO, e_fermi
207 :
208 : ! Remove SOC outside of energy window (otherwise, numerical problems arise
209 : ! because energetically low semicore states and energetically very high
210 : ! unbound states couple to the states around the Fermi level).
211 : ! This routine is for cfm_ks_spinor being in the molecular-orbital (mo) with
212 : ! corresponding eigenvalues "eigenval".
213 :
214 256 : CALL timeset(routineN, handle)
215 :
216 : CALL cp_cfm_get_info(matrix=cfm_ks_spinor, &
217 : nrow_global=nrow_global, &
218 : ncol_global=ncol_global, &
219 : nrow_local=nrow_local, &
220 : ncol_local=ncol_local, &
221 : row_indices=row_indices, &
222 256 : col_indices=col_indices)
223 :
224 256 : CPASSERT(nrow_global == SIZE(eigenval))
225 256 : CPASSERT(ncol_global == SIZE(eigenval))
226 :
227 : ! Apply a smooth energy window to the SOC coupling matrix. Each state gets a
228 : ! weight in [0,1]: 1 = full SOC, decaying to 0 outside the window via a
229 : ! Fermi-like function of width temp_smear. A matrix element coupling states i
230 : ! and j is scaled by SQRT(w_i*w_j), applying the damping symmetrically to both.
231 7424 : DO jjB = 1, ncol_local
232 7168 : j_glob = col_indices(jjB)
233 107776 : DO iiB = 1, nrow_local
234 100352 : i_glob = row_indices(iiB)
235 :
236 100352 : E_i = eigenval(i_glob)
237 100352 : E_j = eigenval(j_glob)
238 100352 : IF (E_i <= e_fermi) THEN
239 : !x_i = ABS(E_i - e_fermi) - 0.5_dp * e_win_vbm
240 : w_i = 1.0_dp
241 : ELSE
242 71680 : x_i = ABS(E_i - e_fermi) - e_win_cbm
243 71680 : w_i = 1.0_dp/(EXP(x_i/temp_smear) + 1.0_dp)
244 : END IF
245 :
246 100352 : IF (E_j <= e_fermi) THEN
247 : !x_j = ABS(E_j - e_fermi) - 0.5_dp * e_win_vbm
248 : w_j = 1.0_dp
249 : ELSE
250 71680 : x_j = ABS(E_j - e_fermi) - e_win_cbm
251 71680 : w_j = 1.0_dp/(EXP(x_j/temp_smear) + 1.0_dp)
252 : END IF
253 :
254 : cfm_ks_spinor%local_data(iiB, jjB) = &
255 107520 : cfm_ks_spinor%local_data(iiB, jjB)*SQRT(w_i*w_j)
256 :
257 : END DO
258 : END DO
259 :
260 256 : CALL timestop(handle)
261 :
262 256 : END SUBROUTINE remove_soc_outside_energy_window_mo
263 :
264 : END MODULE soc_pseudopotential_methods
|
