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