(git:ed6f26b)
Loading...
Searching...
No Matches
soc_pseudopotential_methods.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
10 USE core_ppnl, ONLY: build_core_ppnl
11 USE cp_cfm_types, ONLY: cp_cfm_get_info,&
14 USE cp_dbcsr_api, ONLY: dbcsr_add,&
19 dbcsr_set,&
20 dbcsr_type_antisymmetric,&
21 dbcsr_type_no_symmetry
27 USE cp_fm_types, ONLY: cp_fm_create,&
31 USE kinds, ONLY: dp
32 USE kpoint_types, ONLY: get_kpoint_info,&
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
54
55CONTAINS
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 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 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
77 LOGICAL :: calculate_forces, do_kp, do_symmetric, &
78 use_virial
79 REAL(kind=dp) :: eps_ppnl
80 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
81 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_l, mat_l_nosym, mat_pot_dummy, &
82 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 POINTER :: sab_orb, sap_ppnl
87 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
88 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
89 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
90 TYPE(virial_type), POINTER :: virial
91
92 CALL timeset(routinen, handle)
93
94 NULLIFY (qs_kind_set, dft_control, sab_orb, sap_ppnl, particle_set, atomic_kind_set, &
95 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 particle_set=particle_set, sab_orb=sab_orb, sap_ppnl=sap_ppnl)
99
100 eps_ppnl = dft_control%qs_control%eps_ppnl
101 nimages = dft_control%nimages
102 do_kp = (nimages > 1)
103 CALL get_neighbor_list_set_p(neighbor_list_sets=sab_orb, symmetric=do_symmetric)
104 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
105
106 NULLIFY (mat_l, mat_pot_dummy)
107 CALL dbcsr_allocate_matrix_set(mat_l, 3, nimages)
108 DO xyz = 1, 3
109 DO img = 1, nimages
110 ALLOCATE (mat_l(xyz, img)%matrix)
111 CALL dbcsr_create(mat_l(xyz, img)%matrix, template=matrix_s(1, 1)%matrix, &
112 matrix_type=dbcsr_type_antisymmetric)
113 CALL cp_dbcsr_alloc_block_from_nbl(mat_l(xyz, img)%matrix, sab_orb)
114 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 cpassert(ASSOCIATED(sap_ppnl))
121 nder = 0
122 use_virial = .false.
123 calculate_forces = .false.
124
125 NULLIFY (mat_pot_dummy)
126 CALL dbcsr_allocate_matrix_set(mat_pot_dummy, 1, nimages)
127 DO img = 1, nimages
128 ALLOCATE (mat_pot_dummy(1, img)%matrix)
129 CALL dbcsr_create(mat_pot_dummy(1, img)%matrix, template=matrix_s(1, 1)%matrix)
130 CALL cp_dbcsr_alloc_block_from_nbl(mat_pot_dummy(1, img)%matrix, sab_orb)
131 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 basis_type="ORB", matrix_l=mat_l)
139
140 NULLIFY (mat_l_nosym)
141 CALL dbcsr_allocate_matrix_set(mat_l_nosym, 3, nimages)
142 DO xyz = 1, 3
143 DO img = 1, nimages
144
145 ALLOCATE (mat_l_nosym(xyz, img)%matrix)
146 IF (do_kp) THEN
147 CALL dbcsr_create(mat_l_nosym(xyz, img)%matrix, template=matrix_s(1, 1)%matrix, &
148 matrix_type=dbcsr_type_antisymmetric)
149 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 matrix_type=dbcsr_type_no_symmetry)
153 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 NULLIFY (mat_v_soc_xyz)
160 CALL dbcsr_allocate_matrix_set(mat_v_soc_xyz, 3, nimages)
161 DO xyz = 1, 3
162 DO img = 1, nimages
163 ALLOCATE (mat_v_soc_xyz(xyz, img)%matrix)
164 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 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 matrix_type=dbcsr_type_no_symmetry)
174 END IF
175 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 0.0_dp, 0.5_dp)
179 END DO
180 END DO
181
182 CALL dbcsr_deallocate_matrix_set(mat_pot_dummy)
183 CALL dbcsr_deallocate_matrix_set(mat_l_nosym)
185
186 CALL timestop(handle)
187
188 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 SUBROUTINE remove_soc_outside_energy_window_ao(mat_V_SOC_xyz, e_win, fm_mo_coeff, homo, &
203 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 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 CALL timeset(routinen, handle)
220
221 CALL cp_fm_create(fm_work, fm_s%matrix_struct)
222 CALL cp_fm_create(fm_v_ao, fm_s%matrix_struct)
223 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 col_indices=col_indices)
230
231 nao = SIZE(eigenval)
232
233 e_homo = eigenval(homo)
234 e_lumo = eigenval(homo + 1)
235
236 DO xyz = 1, 3
237
238 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 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 matrix_a=fm_mo_coeff, matrix_b=fm_work, beta=0.0_dp, matrix_c=fm_v_mo)
246
247 DO jjb = 1, ncol_local
248 j_glob = col_indices(jjb)
249 DO iib = 1, nrow_local
250 i_glob = row_indices(iib)
251
252 e_i = eigenval(i_glob)
253 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 e_j < e_homo - 0.5_dp*e_win .OR. e_j > e_lumo + 0.5_dp*e_win) THEN
257 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 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 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 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 matrix_a=fm_work, matrix_b=fm_s, beta=0.0_dp, matrix_c=fm_v_ao)
275
276 CALL copy_fm_to_dbcsr(fm_v_ao, mat_v_soc_xyz(xyz)%matrix)
277
278 END DO
279
280 CALL cp_fm_release(fm_work)
281 CALL cp_fm_release(fm_v_ao)
282 CALL cp_fm_release(fm_v_mo)
283
284 CALL timestop(handle)
285
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 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 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 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 col_indices=col_indices)
325
326 cpassert(nrow_global == SIZE(eigenval))
327 cpassert(ncol_global == SIZE(eigenval))
328
329 DO jjb = 1, ncol_local
330 j_glob = col_indices(jjb)
331 DO iib = 1, nrow_local
332 i_glob = row_indices(iib)
333
334 e_i = eigenval(i_glob)
335 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 e_j < e_homo - 0.5_dp*e_win .OR. e_j > e_lumo + 0.5_dp*e_win) THEN
339 cfm_ks_spinor%local_data(iib, jjb) = 0.0_dp
340 END IF
341
342 END DO
343 END DO
344
345 CALL timestop(handle)
346
348
Define the atomic kind types and their sub types.
Calculation of the non-local pseudopotential contribution to the core Hamiltonian <a|V(non-local)|b> ...
Definition core_ppnl.F:15
subroutine, public build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, eps_ppnl, nimages, cell_to_index, basis_type, deltar, matrix_l, atcore)
...
Definition core_ppnl.F:90
Represents a complex full matrix distributed on many processors.
subroutine, public cp_cfm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, matrix_struct, para_env)
Returns information about a full matrix.
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
subroutine, public dbcsr_desymmetrize(matrix_a, matrix_b)
...
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
subroutine, public dbcsr_set(matrix, alpha)
...
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
DBCSR operations in CP2K.
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Types and basic routines needed for a kpoint calculation.
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym)
Retrieve information from a kpoint environment.
basic linear algebra operations for full matrixes
Define the data structure for the particle information.
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs)
Get the QUICKSTEP environment.
Define the quickstep kind type and their sub types.
Define the neighbor list data types and the corresponding functionality.
subroutine, public get_neighbor_list_set_p(neighbor_list_sets, nlist, symmetric)
Return the components of the first neighbor list set.
subroutine, public v_soc_xyz_from_pseudopotential(qs_env, mat_v_soc_xyz)
V^SOC_µν^(α),R = ħ/2 < ϕ_µ cell O | sum_ℓ ΔV_ℓ^SO(r,r') L^(α) | ϕ_ν cell R>, α = x,...
subroutine, public remove_soc_outside_energy_window_ao(mat_v_soc_xyz, e_win, fm_mo_coeff, homo, eigenval, fm_s)
Remove SOC outside of energy window (otherwise, numerical problems arise because energetically low se...
subroutine, public remove_soc_outside_energy_window_mo(cfm_ks_spinor, e_win, eigenval, e_homo, e_lumo)
...
Provides all information about an atomic kind.
Represent a complex full matrix.
represent a full matrix
Contains information about kpoints.
Provides all information about a quickstep kind.