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