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 : ! **************************************************************************************************
9 : !> \brief Calculates hyperfine values
10 : !> \par History
11 : !> created 04-2006 [RD]
12 : !> adapted 02-2007 [JGH]
13 : !> \author R. Declerck (Reinout.Declerck@UGent.be)
14 : ! **************************************************************************************************
15 : MODULE qs_epr_hyp
16 : USE atomic_kind_types, ONLY: atomic_kind_type,&
17 : get_atomic_kind
18 : USE cell_types, ONLY: cell_type,&
19 : pbc
20 : USE cp_control_types, ONLY: dft_control_type
21 : USE cp_log_handling, ONLY: cp_get_default_logger,&
22 : cp_logger_type
23 : USE cp_output_handling, ONLY: cp_print_key_unit_nr
24 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
25 : section_vals_type,&
26 : section_vals_val_get
27 : USE kinds, ONLY: dp
28 : USE mathconstants, ONLY: fourpi
29 : USE message_passing, ONLY: mp_para_env_type
30 : USE particle_types, ONLY: particle_type
31 : USE periodic_table, ONLY: ptable
32 : USE physcon, ONLY: a_bohr,&
33 : a_fine,&
34 : e_charge,&
35 : e_gfactor,&
36 : e_mass,&
37 : h_bar,&
38 : mu_perm
39 : USE pw_env_types, ONLY: pw_env_get,&
40 : pw_env_type
41 : USE pw_grid_types, ONLY: pw_grid_type
42 : USE pw_methods, ONLY: pw_axpy,&
43 : pw_dr2_gg,&
44 : pw_zero
45 : USE pw_pool_types, ONLY: pw_pool_type
46 : USE pw_types, ONLY: pw_c1d_gs_type
47 : USE qs_environment_types, ONLY: get_qs_env,&
48 : qs_environment_type
49 : USE qs_grid_atom, ONLY: grid_atom_type
50 : USE qs_harmonics_atom, ONLY: harmonics_atom_type
51 : USE qs_kind_types, ONLY: get_qs_kind,&
52 : qs_kind_type
53 : USE qs_rho_atom_types, ONLY: get_rho_atom,&
54 : rho_atom_coeff,&
55 : rho_atom_type
56 : USE qs_rho_types, ONLY: qs_rho_get,&
57 : qs_rho_type
58 : USE util, ONLY: get_limit
59 : #include "./base/base_uses.f90"
60 :
61 : IMPLICIT NONE
62 :
63 : PRIVATE
64 : PUBLIC :: qs_epr_hyp_calc
65 :
66 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_epr_hyp'
67 :
68 : CONTAINS
69 :
70 : ! **************************************************************************************************
71 : !> \brief ...
72 : !> \param qs_env ...
73 : ! **************************************************************************************************
74 30 : SUBROUTINE qs_epr_hyp_calc(qs_env)
75 :
76 : TYPE(qs_environment_type), POINTER :: qs_env
77 :
78 : CHARACTER(LEN=2) :: element_symbol
79 : INTEGER :: bo(2), ia, iat, iatom, idir1, idir2, ig, &
80 : ikind, ir, iso, jatom, mepos, natom, &
81 : natomkind, nkind, num_pe, output_unit, &
82 : z
83 30 : INTEGER, DIMENSION(:), POINTER :: atom_list
84 : LOGICAL :: lsd, paw_atom
85 : REAL(dp) :: arg, esum, hard_radius, hypanisotemp, &
86 : hypfactor, int_radius, rab2, rtemp
87 30 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: hypiso, hypiso_one
88 30 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: hypaniso
89 : REAL(KIND=dp), DIMENSION(3) :: ra, rab
90 30 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
91 : TYPE(cell_type), POINTER :: cell
92 : TYPE(cp_logger_type), POINTER :: logger
93 : TYPE(dft_control_type), POINTER :: dft_control
94 : TYPE(grid_atom_type), POINTER :: grid_atom
95 : TYPE(harmonics_atom_type), POINTER :: harmonics
96 : TYPE(mp_para_env_type), POINTER :: para_env
97 30 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
98 : TYPE(pw_c1d_gs_type) :: hypaniso_gspace, rhototspin_elec_gspace
99 30 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
100 : TYPE(pw_env_type), POINTER :: pw_env
101 : TYPE(pw_grid_type), POINTER :: pw_grid
102 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
103 30 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
104 : TYPE(qs_rho_type), POINTER :: rho
105 30 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: rho_rad_h, rho_rad_s
106 30 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
107 : TYPE(rho_atom_type), POINTER :: rho_atom
108 : TYPE(section_vals_type), POINTER :: dft_section
109 :
110 30 : NULLIFY (pw_env, cell, atomic_kind_set, qs_kind_set, auxbas_pw_pool, dft_control, &
111 30 : logger, dft_section, para_env, particle_set, rho, rho_atom, &
112 30 : rho_atom_set, rho_g)
113 :
114 60 : logger => cp_get_default_logger()
115 30 : dft_section => section_vals_get_subs_vals(qs_env%input, "DFT")
116 : output_unit = cp_print_key_unit_nr(logger, dft_section, &
117 : "PRINT%HYPERFINE_COUPLING_TENSOR", &
118 30 : extension=".eprhyp", log_filename=.FALSE.)
119 : CALL section_vals_val_get(dft_section, &
120 : "PRINT%HYPERFINE_COUPLING_TENSOR%INTERACTION_RADIUS", &
121 30 : r_val=int_radius)
122 30 : CALL section_vals_val_get(dft_section, "LSD", l_val=lsd)
123 :
124 30 : IF (.NOT. lsd) THEN
125 : ! EPR calculation only for LSD
126 28 : IF (output_unit > 0) THEN
127 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
128 14 : "Calculation of EPR hyperfine coupling tensors only for LSD"
129 : END IF
130 28 : NULLIFY (logger, dft_section)
131 28 : RETURN
132 : END IF
133 :
134 2 : hypfactor = -1.0_dp*mu_perm*e_charge*h_bar*e_gfactor/(2.0_dp*e_mass*a_bohr**3)
135 :
136 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, cell=cell, &
137 : rho=rho, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
138 : rho_atom_set=rho_atom_set, pw_env=pw_env, &
139 2 : particle_set=particle_set, para_env=para_env)
140 :
141 2 : IF (output_unit > 0) THEN
142 : WRITE (UNIT=output_unit, FMT="(/,T2,A,/,T2,A)") &
143 1 : "Calculation of EPR hyperfine coupling tensors", &
144 2 : REPEAT("-", 79)
145 : END IF
146 :
147 : ! allocate hyperfine matrices
148 2 : natom = SIZE(particle_set, 1)
149 6 : ALLOCATE (hypaniso(3, 3, natom))
150 6 : ALLOCATE (hypiso(natom))
151 6 : ALLOCATE (hypiso_one(natom))
152 :
153 : ! set the matrices to zero
154 8 : hypiso = 0.0_dp
155 8 : hypiso_one = 0.0_dp
156 80 : hypaniso = 0.0_dp
157 :
158 2 : nkind = SIZE(atomic_kind_set) ! nkind = number of atom types
159 :
160 6 : DO ikind = 1, nkind ! loop over atom types
161 4 : NULLIFY (atom_list, grid_atom, harmonics)
162 : CALL get_atomic_kind(atomic_kind_set(ikind), &
163 4 : atom_list=atom_list, natom=natomkind, z=z)
164 :
165 : CALL get_qs_kind(qs_kind_set(ikind), harmonics=harmonics, &
166 4 : grid_atom=grid_atom, paw_atom=paw_atom, hard_radius=hard_radius)
167 :
168 4 : IF (.NOT. paw_atom) CYCLE ! skip the rest and go to next atom type
169 :
170 4 : num_pe = para_env%num_pe
171 4 : mepos = para_env%mepos
172 4 : bo = get_limit(natomkind, num_pe, mepos)
173 :
174 13 : DO iat = bo(1), bo(2) ! natomkind = # atoms for ikind
175 3 : iatom = atom_list(iat)
176 3 : rho_atom => rho_atom_set(iatom)
177 3 : NULLIFY (rho_rad_h, rho_rad_s)
178 : CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=rho_rad_h, &
179 3 : rho_rad_s=rho_rad_s)
180 : ! Non-relativistic isotropic hyperfine value (hypiso_one)
181 153 : DO ia = 1, grid_atom%ng_sphere
182 2303 : DO iso = 1, harmonics%max_iso_not0
183 : hypiso_one(iatom) = hypiso_one(iatom) + &
184 : (rho_rad_h(1)%r_coef(grid_atom%nr, iso) - &
185 : rho_rad_h(2)%r_coef(grid_atom%nr, iso))* &
186 2300 : harmonics%slm(ia, iso)*grid_atom%wa(ia)/fourpi
187 : END DO
188 : END DO
189 : ! First calculate hard-soft contributions for the own nucleus
190 : ! + scalar relativistic isotropic hyperfine value (hypiso)
191 153 : DO ir = 1, grid_atom%nr
192 153 : IF (grid_atom%rad(ir) <= hard_radius) THEN
193 4080 : DO ia = 1, grid_atom%ng_sphere
194 4000 : hypanisotemp = 0.0_dp
195 62400 : DO iso = 1, harmonics%max_iso_not0
196 : hypiso(iatom) = hypiso(iatom) + &
197 : (rho_rad_h(1)%r_coef(ir, iso) - rho_rad_h(2)%r_coef(ir, iso))* &
198 : harmonics%slm(ia, iso)*grid_atom%wr(ir)*grid_atom%wa(ia)* &
199 : 2._dp/(REAL(z, KIND=dp)*a_fine**2* &
200 : (1._dp + 2._dp*grid_atom%rad(ir)/(REAL(z, KIND=dp)*a_fine**2))**2* &
201 58400 : fourpi*grid_atom%rad(ir)**2)
202 : hypanisotemp = hypanisotemp + &
203 : (rho_rad_h(1)%r_coef(ir, iso) - rho_rad_h(2)%r_coef(ir, iso) &
204 : - (rho_rad_s(1)%r_coef(ir, iso) - rho_rad_s(2)%r_coef(ir, iso)))* &
205 : harmonics%slm(ia, iso)*grid_atom%wr(ir)*grid_atom%wa(ia)/ &
206 62400 : grid_atom%rad(ir)**3
207 : END DO ! iso
208 : hypaniso(1, 1, iatom) = hypaniso(1, 1, iatom) + hypanisotemp* &
209 : (3._dp*grid_atom%sin_pol(ia)*grid_atom%cos_azi(ia)* &
210 4000 : grid_atom%sin_pol(ia)*grid_atom%cos_azi(ia) - 1._dp)
211 : hypaniso(1, 2, iatom) = hypaniso(1, 2, iatom) + hypanisotemp* &
212 : (3._dp*grid_atom%sin_pol(ia)*grid_atom%cos_azi(ia)* &
213 4000 : grid_atom%sin_pol(ia)*grid_atom%sin_azi(ia) - 0._dp)
214 : hypaniso(1, 3, iatom) = hypaniso(1, 3, iatom) + hypanisotemp* &
215 : (3._dp*grid_atom%sin_pol(ia)*grid_atom%cos_azi(ia)* &
216 4000 : grid_atom%cos_pol(ia) - 0._dp)
217 : hypaniso(2, 2, iatom) = hypaniso(2, 2, iatom) + hypanisotemp* &
218 : (3._dp*grid_atom%sin_pol(ia)*grid_atom%sin_azi(ia)* &
219 4000 : grid_atom%sin_pol(ia)*grid_atom%sin_azi(ia) - 1._dp)
220 : hypaniso(2, 3, iatom) = hypaniso(2, 3, iatom) + hypanisotemp* &
221 : (3._dp*grid_atom%sin_pol(ia)*grid_atom%sin_azi(ia)* &
222 4000 : grid_atom%cos_pol(ia) - 0._dp)
223 : hypaniso(3, 3, iatom) = hypaniso(3, 3, iatom) + hypanisotemp* &
224 : (3._dp*grid_atom%cos_pol(ia)* &
225 4080 : grid_atom%cos_pol(ia) - 1._dp)
226 : END DO ! ia
227 : END IF ! hard_radius
228 : END DO ! ir
229 :
230 : ! Now calculate hard-soft anisotropic contributions for the other nuclei
231 16 : DO jatom = 1, natom
232 9 : IF (jatom .EQ. iatom) CYCLE ! iatom already done
233 6 : rab = pbc(particle_set(iatom)%r, particle_set(jatom)%r, cell)
234 24 : rab2 = DOT_PRODUCT(rab, rab)
235 : ! SQRT(rab2) <= int_radius
236 9 : IF (rab2 <= (int_radius*int_radius)) THEN
237 306 : DO ir = 1, grid_atom%nr
238 306 : IF (grid_atom%rad(ir) <= hard_radius) THEN
239 8160 : DO ia = 1, grid_atom%ng_sphere
240 8000 : hypanisotemp = 0.0_dp
241 : rtemp = SQRT(rab2 + grid_atom%rad(ir)**2 + 2.0_dp*grid_atom%rad(ir)* &
242 : (rab(1)*grid_atom%sin_pol(ia)*grid_atom%cos_azi(ia) + &
243 : rab(2)*grid_atom%sin_pol(ia)*grid_atom%sin_azi(ia) + &
244 8000 : rab(3)*grid_atom%cos_pol(ia)))
245 124800 : DO iso = 1, harmonics%max_iso_not0
246 : hypanisotemp = hypanisotemp + &
247 : (rho_rad_h(1)%r_coef(ir, iso) - rho_rad_h(2)%r_coef(ir, iso) &
248 : - (rho_rad_s(1)%r_coef(ir, iso) - rho_rad_s(2)%r_coef(ir, iso)))* &
249 : harmonics%slm(ia, iso)*grid_atom%wr(ir)*grid_atom%wa(ia)/ &
250 124800 : rtemp**5
251 : END DO ! iso
252 : hypaniso(1, 1, jatom) = hypaniso(1, 1, jatom) + hypanisotemp* &
253 : (3._dp*(rab(1) + grid_atom%rad(ir)*grid_atom%sin_pol(ia)*grid_atom%cos_azi(ia))* &
254 8000 : (rab(1) + grid_atom%rad(ir)*grid_atom%sin_pol(ia)*grid_atom%cos_azi(ia)) - rtemp**2)
255 : hypaniso(1, 2, jatom) = hypaniso(1, 2, jatom) + hypanisotemp* &
256 : (3._dp*(rab(1) + grid_atom%rad(ir)*grid_atom%sin_pol(ia)*grid_atom%cos_azi(ia))* &
257 8000 : (rab(2) + grid_atom%rad(ir)*grid_atom%sin_pol(ia)*grid_atom%sin_azi(ia)) - 0._dp)
258 : hypaniso(1, 3, jatom) = hypaniso(1, 3, jatom) + hypanisotemp* &
259 : (3._dp*(rab(1) + grid_atom%rad(ir)*grid_atom%sin_pol(ia)*grid_atom%cos_azi(ia))* &
260 8000 : (rab(3) + grid_atom%rad(ir)*grid_atom%cos_pol(ia)) - 0._dp)
261 : hypaniso(2, 2, jatom) = hypaniso(2, 2, jatom) + hypanisotemp* &
262 : (3._dp*(rab(2) + grid_atom%rad(ir)*grid_atom%sin_pol(ia)*grid_atom%sin_azi(ia))* &
263 8000 : (rab(2) + grid_atom%rad(ir)*grid_atom%sin_pol(ia)*grid_atom%sin_azi(ia)) - rtemp**2)
264 : hypaniso(2, 3, jatom) = hypaniso(2, 3, jatom) + hypanisotemp* &
265 : (3._dp*(rab(2) + grid_atom%rad(ir)*grid_atom%sin_pol(ia)*grid_atom%sin_azi(ia))* &
266 8000 : (rab(3) + grid_atom%rad(ir)*grid_atom%cos_pol(ia)) - 0._dp)
267 : hypaniso(3, 3, jatom) = hypaniso(3, 3, jatom) + hypanisotemp* &
268 : (3._dp*(rab(3) + grid_atom%rad(ir)*grid_atom%cos_pol(ia))* &
269 8160 : (rab(3) + grid_atom%rad(ir)*grid_atom%cos_pol(ia)) - rtemp**2)
270 : END DO ! ia
271 : END IF ! hard_radius
272 : END DO ! ir
273 : END IF ! rab2
274 : END DO ! jatom
275 : END DO ! iat
276 : END DO ! ikind
277 :
278 : ! Now calculate the soft electronic spin density in reciprocal space (g-space)
279 : ! Plane waves grid to assemble the soft electronic spin density
280 : CALL pw_env_get(pw_env=pw_env, &
281 2 : auxbas_pw_pool=auxbas_pw_pool)
282 :
283 2 : CALL auxbas_pw_pool%create_pw(rhototspin_elec_gspace)
284 2 : CALL pw_zero(rhototspin_elec_gspace)
285 :
286 2 : pw_grid => rhototspin_elec_gspace%pw_grid
287 :
288 : ! Load the contribution of the soft electronic density
289 2 : CALL qs_rho_get(rho, rho_g=rho_g)
290 2 : CPASSERT(SIZE(rho_g) > 1)
291 2 : CALL pw_axpy(rho_g(1), rhototspin_elec_gspace)
292 2 : CALL pw_axpy(rho_g(2), rhototspin_elec_gspace, alpha=-1._dp)
293 : ! grid to assemble anisotropic hyperfine terms
294 2 : CALL auxbas_pw_pool%create_pw(hypaniso_gspace)
295 :
296 8 : DO idir1 = 1, 3
297 20 : DO idir2 = idir1, 3 ! tensor symmetry
298 12 : CALL pw_zero(hypaniso_gspace)
299 : CALL pw_dr2_gg(rhototspin_elec_gspace, hypaniso_gspace, &
300 12 : idir1, idir2)
301 54 : DO iatom = 1, natom
302 36 : esum = 0.0_dp
303 36 : ra(:) = pbc(particle_set(iatom)%r, cell)
304 4718628 : DO ig = 1, SIZE(hypaniso_gspace%array)
305 18874368 : arg = DOT_PRODUCT(pw_grid%g(:, ig), ra)
306 : esum = esum + COS(arg)*REAL(hypaniso_gspace%array(ig), dp) &
307 4718628 : - SIN(arg)*AIMAG(hypaniso_gspace%array(ig))
308 : END DO
309 : ! Actually, we need -1.0 * fourpi * hypaniso_gspace
310 36 : esum = esum*fourpi*(-1.0_dp)
311 48 : hypaniso(idir1, idir2, iatom) = hypaniso(idir1, idir2, iatom) + esum
312 : END DO
313 : END DO ! idir2
314 : END DO ! idir1
315 :
316 2 : CALL auxbas_pw_pool%give_back_pw(rhototspin_elec_gspace)
317 2 : CALL auxbas_pw_pool%give_back_pw(hypaniso_gspace)
318 :
319 : ! Multiply hyperfine matrices with constant*gyromagnetic ratio's
320 : ! to have it in units of Mhz.
321 :
322 8 : DO iatom = 1, natom
323 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
324 6 : z=z)
325 : hypiso(iatom) = hypiso(iatom)* &
326 6 : 2.0_dp/3.0_dp*hypfactor*ptable(z)%gyrom_ratio
327 : hypiso_one(iatom) = hypiso_one(iatom)* &
328 6 : 2.0_dp/3.0_dp*hypfactor*ptable(z)%gyrom_ratio
329 26 : DO idir1 = 1, 3
330 60 : DO idir2 = idir1, 3
331 : hypaniso(idir1, idir2, iatom) = hypaniso(idir1, idir2, iatom)* &
332 36 : hypfactor/fourpi*ptable(z)%gyrom_ratio
333 54 : IF (idir1 /= idir2) THEN
334 18 : hypaniso(idir2, idir1, iatom) = hypaniso(idir1, idir2, iatom)
335 : END IF
336 : END DO
337 : END DO
338 : END DO
339 :
340 : ! Global sum
341 2 : CALL para_env%sync()
342 2 : CALL para_env%sum(hypaniso)
343 2 : CALL para_env%sum(hypiso)
344 2 : CALL para_env%sum(hypiso_one)
345 :
346 : ! Print hyperfine matrices
347 2 : IF (output_unit > 0) THEN
348 4 : DO iatom = 1, natom
349 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
350 3 : element_symbol=element_symbol, z=z)
351 : WRITE (UNIT=output_unit, FMT="(T1,I5,T7,A,T10,I3,T14,F16.10,T31,A,T60,F20.10)") &
352 3 : iatom, element_symbol, ptable(z)%gyrom_ratio_isotope, ptable(z)%gyrom_ratio, &
353 6 : "[Mhz/T] Sca-Rel A_iso [Mhz]", hypiso(iatom)
354 : WRITE (UNIT=output_unit, FMT="(T31,A,T60,F20.10)") &
355 3 : " Non-Rel A_iso [Mhz]", hypiso_one(iatom)
356 : WRITE (UNIT=output_unit, FMT="(T4,A,T18,F20.10,1X,F20.10,1X,F20.10)") &
357 3 : " ", hypaniso(1, 1, iatom), hypaniso(1, 2, iatom), hypaniso(1, 3, iatom), &
358 3 : " A_ani [Mhz]", hypaniso(2, 1, iatom), hypaniso(2, 2, iatom), hypaniso(2, 3, iatom), &
359 10 : " ", hypaniso(3, 1, iatom), hypaniso(3, 2, iatom), hypaniso(3, 3, iatom)
360 : END DO
361 : END IF
362 :
363 : ! Deallocate the remainings ...
364 2 : DEALLOCATE (hypiso)
365 2 : DEALLOCATE (hypiso_one)
366 2 : DEALLOCATE (hypaniso)
367 :
368 90 : END SUBROUTINE qs_epr_hyp_calc
369 :
370 : END MODULE qs_epr_hyp
371 :
|
