1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \par History 8!> created 06-2006 [RD] 9!> \author RD 10! ************************************************************************************************** 11MODULE qs_linres_epr_ownutils 12 USE atomic_kind_types, ONLY: atomic_kind_type,& 13 get_atomic_kind 14 USE cell_types, ONLY: cell_type 15 USE cp_control_types, ONLY: dft_control_type 16 USE cp_log_handling, ONLY: cp_get_default_logger,& 17 cp_logger_type 18 USE cp_output_handling, ONLY: cp_p_file,& 19 cp_print_key_finished_output,& 20 cp_print_key_should_output,& 21 cp_print_key_unit_nr 22 USE cp_para_types, ONLY: cp_para_env_type 23 USE dbcsr_api, ONLY: dbcsr_p_type 24 USE input_section_types, ONLY: section_vals_get_subs_vals,& 25 section_vals_type,& 26 section_vals_val_get 27 USE kinds, ONLY: default_string_length,& 28 dp 29 USE mathlib, ONLY: diamat_all 30 USE message_passing, ONLY: mp_sum 31 USE particle_types, ONLY: particle_type 32 USE pw_env_types, ONLY: pw_env_get,& 33 pw_env_type 34 USE pw_methods, ONLY: pw_axpy,& 35 pw_integral_ab,& 36 pw_scale,& 37 pw_transfer,& 38 pw_zero 39 USE pw_pool_types, ONLY: pw_pool_create_pw,& 40 pw_pool_give_back_pw,& 41 pw_pool_p_type,& 42 pw_pool_type 43 USE pw_spline_utils, ONLY: Eval_Interp_Spl3_pbc,& 44 find_coeffs,& 45 pw_spline_do_precond,& 46 pw_spline_precond_create,& 47 pw_spline_precond_release,& 48 pw_spline_precond_set_kind,& 49 pw_spline_precond_type,& 50 spl3_pbc 51 USE pw_types, ONLY: COMPLEXDATA1D,& 52 REALDATA3D,& 53 REALSPACE,& 54 RECIPROCALSPACE,& 55 pw_p_type 56 USE qs_core_energies, ONLY: calculate_ptrace 57 USE qs_environment_types, ONLY: get_qs_env,& 58 qs_environment_type 59 USE qs_grid_atom, ONLY: grid_atom_type 60 USE qs_harmonics_atom, ONLY: harmonics_atom_type 61 USE qs_kind_types, ONLY: get_qs_kind,& 62 qs_kind_type 63 USE qs_linres_nmr_epr_common_utils, ONLY: mult_G_ov_G2_grid 64 USE qs_linres_op, ONLY: fac_vecp,& 65 set_vecp,& 66 set_vecp_rev 67 USE qs_linres_types, ONLY: current_env_type,& 68 epr_env_type,& 69 get_current_env,& 70 get_epr_env,& 71 jrho_atom_type,& 72 nablavks_atom_type 73 USE qs_rho_atom_types, ONLY: get_rho_atom,& 74 rho_atom_coeff,& 75 rho_atom_type 76 USE qs_rho_types, ONLY: qs_rho_get,& 77 qs_rho_p_type,& 78 qs_rho_type 79 USE realspace_grid_types, ONLY: realspace_grid_desc_type 80 USE util, ONLY: get_limit 81#include "./base/base_uses.f90" 82 83 IMPLICIT NONE 84 85 PRIVATE 86 PUBLIC :: epr_g_print, epr_g_zke, epr_g_so, epr_g_soo, epr_ind_magnetic_field 87 88 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_epr_ownutils' 89 90CONTAINS 91 92! ************************************************************************************************** 93!> \brief Prints the g tensor 94!> \param epr_env ... 95!> \param qs_env ... 96!> \par History 97!> 06.2006 created [RD] 98!> \author RD 99! ************************************************************************************************** 100 SUBROUTINE epr_g_print(epr_env, qs_env) 101 102 TYPE(epr_env_type) :: epr_env 103 TYPE(qs_environment_type), POINTER :: qs_env 104 105 CHARACTER(LEN=*), PARAMETER :: routineN = 'epr_g_print', routineP = moduleN//':'//routineN 106 107 CHARACTER(LEN=default_string_length) :: title 108 INTEGER :: idir1, idir2, output_unit, unit_nr 109 REAL(KIND=dp) :: eigenv_g(3), g_sym(3, 3), gsum 110 TYPE(cp_logger_type), POINTER :: logger 111 TYPE(section_vals_type), POINTER :: lr_section 112 113 NULLIFY (logger, lr_section) 114 115 logger => cp_get_default_logger() 116 lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES") 117 118 output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", & 119 extension=".linresLog") 120 121 gsum = 0.0_dp 122 123 DO idir1 = 1, 3 124 DO idir2 = 1, 3 125 gsum = gsum + epr_env%g_total(idir1, idir2) 126 END DO 127 END DO 128 129 IF (output_unit > 0) THEN 130 WRITE (UNIT=output_unit, FMT="(T2,A,T56,E14.6)") & 131 "epr|TOT:checksum", gsum 132 END IF 133 134 CALL cp_print_key_finished_output(output_unit, logger, lr_section, & 135 "PRINT%PROGRAM_RUN_INFO") 136 137 IF (BTEST(cp_print_key_should_output(logger%iter_info, lr_section, & 138 "EPR%PRINT%G_TENSOR"), cp_p_file)) THEN 139 140 unit_nr = cp_print_key_unit_nr(logger, lr_section, "EPR%PRINT%G_TENSOR", & 141 extension=".data", middle_name="GTENSOR", & 142 log_filename=.FALSE.) 143 144 IF (unit_nr > 0) THEN 145 146 WRITE (title, "(A)") "G tensor " 147 WRITE (unit_nr, "(T2,A)") title 148 149 WRITE (unit_nr, "(T2,A)") "gmatrix_zke" 150 WRITE (unit_nr, "(3(A,f15.10))") " XX=", epr_env%g_zke, & 151 " XY=", 0.0_dp, " XZ=", 0.0_dp 152 WRITE (unit_nr, "(3(A,f15.10))") " YX=", 0.0_dp, & 153 " YY=", epr_env%g_zke, " YZ=", 0.0_dp 154 WRITE (unit_nr, "(3(A,f15.10))") " ZX=", 0.0_dp, & 155 " ZY=", 0.0_dp, " ZZ=", epr_env%g_zke 156 157 WRITE (unit_nr, "(T2,A)") "gmatrix_so" 158 WRITE (unit_nr, "(3(A,f15.10))") " XX=", epr_env%g_so(1, 1), & 159 " XY=", epr_env%g_so(1, 2), " XZ=", epr_env%g_so(1, 3) 160 WRITE (unit_nr, "(3(A,f15.10))") " YX=", epr_env%g_so(2, 1), & 161 " YY=", epr_env%g_so(2, 2), " YZ=", epr_env%g_so(2, 3) 162 WRITE (unit_nr, "(3(A,f15.10))") " ZX=", epr_env%g_so(3, 1), & 163 " ZY=", epr_env%g_so(3, 2), " ZZ=", epr_env%g_so(3, 3) 164 165 WRITE (unit_nr, "(T2,A)") "gmatrix_soo" 166 WRITE (unit_nr, "(3(A,f15.10))") " XX=", epr_env%g_soo(1, 1), & 167 " XY=", epr_env%g_soo(1, 2), " XZ=", epr_env%g_soo(1, 3) 168 WRITE (unit_nr, "(3(A,f15.10))") " YX=", epr_env%g_soo(2, 1), & 169 " YY=", epr_env%g_soo(2, 2), " YZ=", epr_env%g_soo(2, 3) 170 WRITE (unit_nr, "(3(A,f15.10))") " ZX=", epr_env%g_soo(3, 1), & 171 " ZY=", epr_env%g_soo(3, 2), " ZZ=", epr_env%g_soo(3, 3) 172 173 WRITE (unit_nr, "(T2,A)") "gmatrix_total" 174 WRITE (unit_nr, "(3(A,f15.10))") " XX=", epr_env%g_total(1, 1) + epr_env%g_free_factor, & 175 " XY=", epr_env%g_total(1, 2), " XZ=", epr_env%g_total(1, 3) 176 WRITE (unit_nr, "(3(A,f15.10))") " YX=", epr_env%g_total(2, 1), & 177 " YY=", epr_env%g_total(2, 2) + epr_env%g_free_factor, " YZ=", epr_env%g_total(2, 3) 178 WRITE (unit_nr, "(3(A,f15.10))") " ZX=", epr_env%g_total(3, 1), & 179 " ZY=", epr_env%g_total(3, 2), " ZZ=", epr_env%g_total(3, 3) + epr_env%g_free_factor 180 181 DO idir1 = 1, 3 182 DO idir2 = 1, 3 183 g_sym(idir1, idir2) = (epr_env%g_total(idir1, idir2) + & 184 epr_env%g_total(idir2, idir1))/2.0_dp 185 END DO 186 END DO 187 188 WRITE (unit_nr, "(T2,A)") "gtensor_total" 189 WRITE (unit_nr, "(3(A,f15.10))") " XX=", g_sym(1, 1) + epr_env%g_free_factor, & 190 " XY=", g_sym(1, 2), " XZ=", g_sym(1, 3) 191 WRITE (unit_nr, "(3(A,f15.10))") " YX=", g_sym(2, 1), & 192 " YY=", g_sym(2, 2) + epr_env%g_free_factor, " YZ=", g_sym(2, 3) 193 WRITE (unit_nr, "(3(A,f15.10))") " ZX=", g_sym(3, 1), & 194 " ZY=", g_sym(3, 2), " ZZ=", g_sym(3, 3) + epr_env%g_free_factor 195 196 CALL diamat_all(g_sym, eigenv_g) 197 eigenv_g(:) = eigenv_g(:)*1.0e6_dp 198 199 WRITE (unit_nr, "(T2,A)") "delta_g principal values in ppm" 200 WRITE (unit_nr, "(f15.3,3(A,f15.10))") eigenv_g(1), " X=", g_sym(1, 1), & 201 " Y=", g_sym(2, 1), " Z=", g_sym(3, 1) 202 WRITE (unit_nr, "(f15.3,3(A,f15.10))") eigenv_g(2), " X=", g_sym(1, 2), & 203 " Y=", g_sym(2, 2), " Z=", g_sym(3, 2) 204 WRITE (unit_nr, "(f15.3,3(A,f15.10))") eigenv_g(3), " X=", g_sym(1, 3), & 205 " Y=", g_sym(2, 3), " Z=", g_sym(3, 3) 206 207 END IF 208 209 CALL cp_print_key_finished_output(unit_nr, logger, lr_section,& 210 & "EPR%PRINT%G_TENSOR") 211 212 ENDIF 213 214 END SUBROUTINE epr_g_print 215 216! ************************************************************************************************** 217!> \brief Calculate zke part of the g tensor 218!> \param epr_env ... 219!> \param qs_env ... 220!> \par History 221!> 06.2006 created [RD] 222!> \author RD 223! ************************************************************************************************** 224 SUBROUTINE epr_g_zke(epr_env, qs_env) 225 226 TYPE(epr_env_type) :: epr_env 227 TYPE(qs_environment_type), POINTER :: qs_env 228 229 CHARACTER(LEN=*), PARAMETER :: routineN = 'epr_g_zke', routineP = moduleN//':'//routineN 230 231 INTEGER :: i1, ispin, output_unit 232 REAL(KIND=dp) :: epr_g_zke_temp(2) 233 TYPE(cp_logger_type), POINTER :: logger 234 TYPE(cp_para_env_type), POINTER :: para_env 235 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: kinetic, rho_ao 236 TYPE(dft_control_type), POINTER :: dft_control 237 TYPE(qs_rho_type), POINTER :: rho 238 TYPE(section_vals_type), POINTER :: lr_section 239 240 NULLIFY (dft_control, logger, lr_section, rho, kinetic, para_env, rho_ao) 241 242 logger => cp_get_default_logger() 243 lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES") 244 245 output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", & 246 extension=".linresLog") 247 248 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, & 249 kinetic=kinetic, rho=rho, para_env=para_env) 250 251 CALL qs_rho_get(rho, rho_ao=rho_ao) 252 253 DO ispin = 1, dft_control%nspins 254 CALL calculate_ptrace(kinetic(1)%matrix, rho_ao(ispin)%matrix, & 255 ecore=epr_g_zke_temp(ispin)) 256 END DO 257 258 epr_env%g_zke = epr_env%g_zke_factor*(epr_g_zke_temp(1) - epr_g_zke_temp(2)) 259 DO i1 = 1, 3 260 epr_env%g_total(i1, i1) = epr_env%g_total(i1, i1) + epr_env%g_zke 261 END DO 262 263 IF (output_unit > 0) THEN 264 WRITE (UNIT=output_unit, FMT="(T2,A,T56,E24.16)") & 265 "epr|ZKE:g_zke", epr_env%g_zke 266 END IF 267 268 CALL cp_print_key_finished_output(output_unit, logger, lr_section, & 269 "PRINT%PROGRAM_RUN_INFO") 270 271 END SUBROUTINE epr_g_zke 272 273! ************************************************************************************************** 274!> \brief Calculates g_so 275!> \param epr_env ... 276!> \param current_env ... 277!> \param qs_env ... 278!> \param iB ... 279!> \par History 280!> 06.2006 created [RD] 281!> \author RD 282! ************************************************************************************************** 283 SUBROUTINE epr_g_so(epr_env, current_env, qs_env, iB) 284 285 TYPE(epr_env_type) :: epr_env 286 TYPE(current_env_type) :: current_env 287 TYPE(qs_environment_type), POINTER :: qs_env 288 INTEGER, INTENT(IN) :: iB 289 290 CHARACTER(LEN=*), PARAMETER :: routineN = 'epr_g_so', routineP = moduleN//':'//routineN 291 292 INTEGER :: aint_precond, ia, iat, iatom, idir1, & 293 idir2, idir3, ikind, ir, ispin, & 294 max_iter, natom, nkind, nspins, & 295 output_unit, precond_kind 296 INTEGER, DIMENSION(2) :: bo 297 INTEGER, DIMENSION(:), POINTER :: atom_list 298 LOGICAL :: gapw, paw_atom, success 299 REAL(dp) :: eps_r, eps_x, hard_radius, temp_so_soft, & 300 vks_ra_idir2, vks_ra_idir3 301 REAL(dp), DIMENSION(3, 3) :: temp_so_gapw 302 REAL(dp), DIMENSION(:, :), POINTER :: g_so, g_total 303 REAL(KIND=dp), DIMENSION(3) :: ra 304 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 305 TYPE(cp_logger_type), POINTER :: logger 306 TYPE(cp_para_env_type), POINTER :: para_env 307 TYPE(dft_control_type), POINTER :: dft_control 308 TYPE(grid_atom_type), POINTER :: grid_atom 309 TYPE(harmonics_atom_type), POINTER :: harmonics 310 TYPE(jrho_atom_type), DIMENSION(:), POINTER :: jrho1_atom_set 311 TYPE(jrho_atom_type), POINTER :: jrho1_atom 312 TYPE(nablavks_atom_type), DIMENSION(:), POINTER :: nablavks_atom_set 313 TYPE(nablavks_atom_type), POINTER :: nablavks_atom 314 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 315 TYPE(pw_env_type), POINTER :: pw_env 316 TYPE(pw_p_type), DIMENSION(:), POINTER :: jrho2_r, jrho3_r, nrho1_r, nrho2_r, & 317 nrho3_r 318 TYPE(pw_p_type), DIMENSION(:, :), POINTER :: vks_pw_spline 319 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 320 TYPE(pw_spline_precond_type), POINTER :: precond 321 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 322 TYPE(qs_rho_p_type), DIMENSION(:), POINTER :: jrho1_set 323 TYPE(qs_rho_p_type), DIMENSION(:, :), POINTER :: nablavks_set 324 TYPE(section_vals_type), POINTER :: interp_section, lr_section 325 326 NULLIFY (atomic_kind_set, qs_kind_set, atom_list, dft_control, & 327 grid_atom, g_so, g_total, harmonics, interp_section, jrho1_atom_set, & 328 jrho1_set, logger, lr_section, nablavks_atom, nablavks_atom_set, & 329 nablavks_set, para_env, particle_set, jrho2_r, jrho3_r, nrho1_r, nrho2_r, nrho3_r) 330 331 logger => cp_get_default_logger() 332 lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES") 333 334 output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", & 335 extension=".linresLog") 336 337 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, & 338 atomic_kind_set=atomic_kind_set, & 339 qs_kind_set=qs_kind_set, & 340 para_env=para_env, pw_env=pw_env, & 341 particle_set=particle_set) 342 343 CALL get_epr_env(epr_env=epr_env, & 344 nablavks_set=nablavks_set, & 345 nablavks_atom_set=nablavks_atom_set, & 346 g_total=g_total, g_so=g_so) 347 348 CALL get_current_env(current_env=current_env, & 349 jrho1_set=jrho1_set, jrho1_atom_set=jrho1_atom_set) 350 351 gapw = dft_control%qs_control%gapw 352 nkind = SIZE(qs_kind_set, 1) 353 nspins = dft_control%nspins 354 355 DO idir1 = 1, 3 356 CALL set_vecp(idir1, idir2, idir3) 357 ! j_pw x nabla_vks_pw 358 temp_so_soft = 0.0_dp 359 DO ispin = 1, nspins 360 CALL qs_rho_get(jrho1_set(idir2)%rho, rho_r=jrho2_r) 361 CALL qs_rho_get(jrho1_set(idir3)%rho, rho_r=jrho3_r) 362 CALL qs_rho_get(nablavks_set(idir2, ispin)%rho, rho_r=nrho2_r) 363 CALL qs_rho_get(nablavks_set(idir3, ispin)%rho, rho_r=nrho3_r) 364 temp_so_soft = temp_so_soft + (-1.0_dp)**(1 + ispin)*( & 365 pw_integral_ab(jrho2_r(ispin)%pw, nrho3_r(1)%pw) - & 366 pw_integral_ab(jrho3_r(ispin)%pw, nrho2_r(1)%pw)) 367 END DO 368 temp_so_soft = -1.0_dp*epr_env%g_so_factor*temp_so_soft 369 IF (output_unit > 0) THEN 370 WRITE (UNIT=output_unit, FMT="(T2,A,T18,I1,I1,T56,E24.16)") & 371 "epr|SOX:soft", iB, idir1, temp_so_soft 372 END IF 373 g_so(iB, idir1) = temp_so_soft 374 END DO !idir1 375 376 IF (gapw) THEN 377 378 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool) 379 ALLOCATE (vks_pw_spline(3, nspins)) 380 381 interp_section => section_vals_get_subs_vals(lr_section, & 382 "EPR%INTERPOLATOR") 383 CALL section_vals_val_get(interp_section, "aint_precond", & 384 i_val=aint_precond) 385 CALL section_vals_val_get(interp_section, "precond", i_val=precond_kind) 386 CALL section_vals_val_get(interp_section, "max_iter", i_val=max_iter) 387 CALL section_vals_val_get(interp_section, "eps_r", r_val=eps_r) 388 CALL section_vals_val_get(interp_section, "eps_x", r_val=eps_x) 389 390 DO ispin = 1, nspins 391 DO idir1 = 1, 3 392 CALL pw_pool_create_pw(auxbas_pw_pool, & 393 vks_pw_spline(idir1, ispin)%pw, & 394 use_data=REALDATA3D, in_space=REALSPACE) 395 ! calculate spline coefficients 396 CALL pw_spline_precond_create(precond, precond_kind=aint_precond, & 397 pool=auxbas_pw_pool, pbc=.TRUE., transpose=.FALSE.) 398 CALL qs_rho_get(nablavks_set(idir1, ispin)%rho, rho_r=nrho1_r) 399 CALL pw_spline_do_precond(precond, nrho1_r(1)%pw, & 400 vks_pw_spline(idir1, ispin)%pw) 401 CALL pw_spline_precond_set_kind(precond, precond_kind) 402 success = find_coeffs(values=nrho1_r(1)%pw, & 403 coeffs=vks_pw_spline(idir1, ispin)%pw, linOp=spl3_pbc, & 404 preconditioner=precond, pool=auxbas_pw_pool, & 405 eps_r=eps_r, eps_x=eps_x, max_iter=max_iter) 406 CPASSERT(success) 407 CALL pw_spline_precond_release(precond) 408 END DO ! idir1 409 END DO ! ispin 410 411 temp_so_gapw = 0.0_dp 412 413 DO ikind = 1, nkind 414 NULLIFY (atom_list, grid_atom, harmonics) 415 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom) 416 CALL get_qs_kind(qs_kind_set(ikind), & 417 hard_radius=hard_radius, & 418 grid_atom=grid_atom, & 419 harmonics=harmonics, & 420 paw_atom=paw_atom) 421 422 IF (.NOT. paw_atom) CYCLE 423 424 ! Distribute the atoms of this kind 425 426 bo = get_limit(natom, para_env%num_pe, para_env%mepos) 427 428 DO iat = 1, natom !bo(1),bo(2)! this partitioning blocks the interpolation 429 ! ! routines (i.e. waiting for mp_sum) 430 iatom = atom_list(iat) 431 NULLIFY (jrho1_atom, nablavks_atom) 432 jrho1_atom => jrho1_atom_set(iatom) 433 nablavks_atom => nablavks_atom_set(iatom) 434 DO idir1 = 1, 3 435 CALL set_vecp(idir1, idir2, idir3) 436 DO ispin = 1, nspins 437 DO ir = 1, grid_atom%nr 438 439 IF (grid_atom%rad(ir) >= hard_radius) CYCLE 440 441 DO ia = 1, grid_atom%ng_sphere 442 443 ra = particle_set(iatom)%r 444 ra(:) = ra(:) + grid_atom%rad(ir)*harmonics%a(:, ia) 445 vks_ra_idir2 = Eval_Interp_Spl3_pbc(ra, & 446 vks_pw_spline(idir2, ispin)%pw) 447 vks_ra_idir3 = Eval_Interp_Spl3_pbc(ra, & 448 vks_pw_spline(idir3, ispin)%pw) 449 450 IF (iat .LT. bo(1) .OR. iat .GT. bo(2)) CYCLE !quick and dirty: 451 ! !here take care of the partition 452 453 ! + sum_A j_loc_h_A x nabla_vks_s_A 454 temp_so_gapw(iB, idir1) = temp_so_gapw(iB, idir1) + & 455 (-1.0_dp)**(1 + ispin)*( & 456 jrho1_atom%jrho_vec_rad_h(idir2, ispin)%r_coef(ir, ia)* & 457 vks_ra_idir3 - & 458 jrho1_atom%jrho_vec_rad_h(idir3, ispin)%r_coef(ir, ia)* & 459 vks_ra_idir2 & 460 )*grid_atom%wr(ir)*grid_atom%wa(ia) 461 462 ! - sum_A j_loc_s_A x nabla_vks_s_A 463 temp_so_gapw(iB, idir1) = temp_so_gapw(iB, idir1) - & 464 (-1.0_dp)**(1 + ispin)*( & 465 jrho1_atom%jrho_vec_rad_s(idir2, ispin)%r_coef(ir, ia)* & 466 vks_ra_idir3 - & 467 jrho1_atom%jrho_vec_rad_s(idir3, ispin)%r_coef(ir, ia)* & 468 vks_ra_idir2 & 469 )*grid_atom%wr(ir)*grid_atom%wa(ia) 470 471 ! + sum_A j_loc_h_A x nabla_vks_loc_h_A 472 temp_so_gapw(iB, idir1) = temp_so_gapw(iB, idir1) + & 473 (-1.0_dp)**(1 + ispin)*( & 474 jrho1_atom%jrho_vec_rad_h(idir2, ispin)%r_coef(ir, ia)* & 475 nablavks_atom%nablavks_vec_rad_h(idir3, ispin)%r_coef(ir, ia) - & 476 jrho1_atom%jrho_vec_rad_h(idir3, ispin)%r_coef(ir, ia)* & 477 nablavks_atom%nablavks_vec_rad_h(idir2, ispin)%r_coef(ir, ia) & 478 )*grid_atom%wr(ir)*grid_atom%wa(ia) 479 480 ! - sum_A j_loc_h_A x nabla_vks_loc_s_A 481 temp_so_gapw(iB, idir1) = temp_so_gapw(iB, idir1) - & 482 (-1.0_dp)**(1 + ispin)*( & 483 jrho1_atom%jrho_vec_rad_h(idir2, ispin)%r_coef(ir, ia)* & 484 nablavks_atom%nablavks_vec_rad_s(idir3, ispin)%r_coef(ir, ia) - & 485 jrho1_atom%jrho_vec_rad_h(idir3, ispin)%r_coef(ir, ia)* & 486 nablavks_atom%nablavks_vec_rad_s(idir2, ispin)%r_coef(ir, ia) & 487 )*grid_atom%wr(ir)*grid_atom%wa(ia) 488 489! ORIGINAL 490! ! + sum_A j_loc_h_A x nabla_vks_loc_h_A 491! temp_so_gapw(iB,idir1) = temp_so_gapw(iB,idir1) + & 492! (-1.0_dp)**(1.0_dp + REAL(ispin,KIND=dp)) * ( & 493! jrho1_atom%jrho_vec_rad_h(idir2,iB,ispin)%r_coef(ir,ia) * & 494! nablavks_atom%nablavks_vec_rad_h(idir3,ispin)%r_coef(ir,ia) - & 495! jrho1_atom%jrho_vec_rad_h(idir3,iB,ispin)%r_coef(ir,ia) * & 496! nablavks_atom%nablavks_vec_rad_h(idir2,ispin)%r_coef(ir,ia) & 497! ) * grid_atom%wr(ir)*grid_atom%wa(ia) 498! ! - sum_A j_loc_s_A x nabla_vks_loc_s_A 499! temp_so_gapw(iB,idir1) = temp_so_gapw(iB,idir1) - & 500! (-1.0_dp)**(1.0_dp + REAL(ispin,KIND=dp)) * ( & 501! jrho1_atom%jrho_vec_rad_s(idir2,iB,ispin)%r_coef(ir,ia) * & 502! nablavks_atom%nablavks_vec_rad_s(idir3,ispin)%r_coef(ir,ia) - & 503! jrho1_atom%jrho_vec_rad_s(idir3,iB,ispin)%r_coef(ir,ia) * & 504! nablavks_atom%nablavks_vec_rad_s(idir2,ispin)%r_coef(ir,ia) & 505! ) * grid_atom%wr(ir)*grid_atom%wa(ia) 506 END DO !ia 507 END DO !ir 508 END DO !ispin 509 END DO !idir1 510 END DO !iat 511 END DO !ikind 512 513 CALL mp_sum(temp_so_gapw, para_env%group) 514 temp_so_gapw(:, :) = -1.0_dp*epr_env%g_so_factor_gapw*temp_so_gapw(:, :) 515 516 IF (output_unit > 0) THEN 517 DO idir1 = 1, 3 518 WRITE (UNIT=output_unit, FMT="(T2,A,T18,I1,I1,T56,E24.16)") & 519 "epr|SOX:gapw", iB, idir1, temp_so_gapw(iB, idir1) 520 END DO 521 END IF 522 523 g_so(iB, :) = g_so(iB, :) + temp_so_gapw(iB, :) 524 525 DO ispin = 1, nspins 526 DO idir1 = 1, 3 527 CALL pw_pool_give_back_pw(auxbas_pw_pool, vks_pw_spline(idir1, ispin)%pw) 528 END DO 529 END DO 530 531 DEALLOCATE (vks_pw_spline) 532 533 END IF ! gapw 534 535 g_total(iB, :) = g_total(iB, :) + g_so(iB, :) 536 537 CALL cp_print_key_finished_output(output_unit, logger, lr_section, & 538 "PRINT%PROGRAM_RUN_INFO") 539 540 END SUBROUTINE epr_g_so 541 542! ************************************************************************************************** 543!> \brief Calculates g_soo (soft part only for now) 544!> \param epr_env ... 545!> \param current_env ... 546!> \param qs_env ... 547!> \param iB ... 548!> \par History 549!> 06.2006 created [RD] 550!> \author RD 551! ************************************************************************************************** 552 SUBROUTINE epr_g_soo(epr_env, current_env, qs_env, iB) 553 554 TYPE(epr_env_type) :: epr_env 555 TYPE(current_env_type) :: current_env 556 TYPE(qs_environment_type), POINTER :: qs_env 557 INTEGER, INTENT(IN) :: iB 558 559 CHARACTER(LEN=*), PARAMETER :: routineN = 'epr_g_soo', routineP = moduleN//':'//routineN 560 561 INTEGER :: aint_precond, ia, iat, iatom, idir1, & 562 ikind, ir, iso, ispin, max_iter, & 563 natom, nkind, nspins, output_unit, & 564 precond_kind 565 INTEGER, DIMENSION(2) :: bo 566 INTEGER, DIMENSION(:), POINTER :: atom_list 567 LOGICAL :: gapw, paw_atom, soo_rho_hard, success 568 REAL(dp) :: bind_ra_idir1, chi_tensor(3, 3, 2), & 569 eps_r, eps_x, hard_radius, rho_spin, & 570 temp_soo_soft 571 REAL(dp), DIMENSION(3, 3) :: temp_soo_gapw 572 REAL(dp), DIMENSION(:, :), POINTER :: g_soo, g_total 573 REAL(KIND=dp), DIMENSION(3) :: ra 574 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 575 TYPE(cp_logger_type), POINTER :: logger 576 TYPE(cp_para_env_type), POINTER :: para_env 577 TYPE(dft_control_type), POINTER :: dft_control 578 TYPE(grid_atom_type), POINTER :: grid_atom 579 TYPE(harmonics_atom_type), POINTER :: harmonics 580 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 581 TYPE(pw_env_type), POINTER :: pw_env 582 TYPE(pw_p_type), DIMENSION(:), POINTER :: brho1_r, rho_r 583 TYPE(pw_p_type), DIMENSION(:, :), POINTER :: bind_pw_spline 584 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 585 TYPE(pw_spline_precond_type), POINTER :: precond 586 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 587 TYPE(qs_rho_p_type), DIMENSION(:, :), POINTER :: bind_set 588 TYPE(qs_rho_type), POINTER :: rho 589 TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: rho_rad_h, rho_rad_s 590 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set 591 TYPE(rho_atom_type), POINTER :: rho_atom 592 TYPE(section_vals_type), POINTER :: g_section, interp_section, lr_section 593 594 NULLIFY (atomic_kind_set, qs_kind_set, atom_list, bind_set, dft_control, & 595 grid_atom, g_section, g_soo, g_total, harmonics, interp_section, & 596 logger, lr_section, para_env, particle_set, rho, rho_atom, & 597 rho_atom_set, rho_r, brho1_r) 598 599 logger => cp_get_default_logger() 600 lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES") 601 602 output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", & 603 extension=".linresLog") 604 605 g_section => section_vals_get_subs_vals(lr_section, & 606 "EPR%PRINT%G_TENSOR") 607 608 CALL section_vals_val_get(g_section, "soo_rho_hard", l_val=soo_rho_hard) 609 610 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, & 611 dft_control=dft_control, para_env=para_env, particle_set=particle_set, & 612 pw_env=pw_env, rho=rho, rho_atom_set=rho_atom_set) 613 614 CALL get_epr_env(epr_env=epr_env, bind_set=bind_set, & 615 g_soo=g_soo, g_total=g_total) 616 617 CALL get_current_env(current_env=current_env, & 618 chi_tensor=chi_tensor) 619 CALL qs_rho_get(rho, rho_r=rho_r) 620 621 gapw = dft_control%qs_control%gapw 622 nkind = SIZE(qs_kind_set, 1) 623 nspins = dft_control%nspins 624 625 DO idir1 = 1, 3 626 temp_soo_soft = 0.0_dp 627 DO ispin = 1, nspins 628 CALL qs_rho_get(bind_set(idir1, iB)%rho, rho_r=brho1_r) 629 temp_soo_soft = temp_soo_soft + (-1.0_dp)**(1 + ispin)* & 630 pw_integral_ab(brho1_r(1)%pw, rho_r(ispin)%pw) 631 END DO 632 temp_soo_soft = 1.0_dp*epr_env%g_soo_factor*temp_soo_soft 633 IF (output_unit > 0) THEN 634 WRITE (UNIT=output_unit, FMT="(T2,A,T18,i1,i1,T56,E24.16)") & 635 "epr|SOO:soft", iB, idir1, temp_soo_soft 636 END IF 637 g_soo(iB, idir1) = temp_soo_soft 638 END DO 639 640 DO idir1 = 1, 3 641 temp_soo_soft = 1.0_dp*epr_env%g_soo_chicorr_factor*chi_tensor(idir1, iB, 2)* & 642 (REAL(dft_control%multiplicity, KIND=dp) - 1.0_dp) 643 IF (output_unit > 0) THEN 644 WRITE (UNIT=output_unit, FMT="(T2,A,T18,i1,i1,T56,E24.16)") & 645 "epr|SOO:soft_g0", iB, idir1, temp_soo_soft 646 END IF 647 g_soo(iB, idir1) = g_soo(iB, idir1) + temp_soo_soft 648 END DO 649 650 IF (gapw .AND. soo_rho_hard) THEN 651 652 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool) 653 ALLOCATE (bind_pw_spline(3, 3)) 654 655 interp_section => section_vals_get_subs_vals(lr_section, & 656 "EPR%INTERPOLATOR") 657 CALL section_vals_val_get(interp_section, "aint_precond", & 658 i_val=aint_precond) 659 CALL section_vals_val_get(interp_section, "precond", i_val=precond_kind) 660 CALL section_vals_val_get(interp_section, "max_iter", i_val=max_iter) 661 CALL section_vals_val_get(interp_section, "eps_r", r_val=eps_r) 662 CALL section_vals_val_get(interp_section, "eps_x", r_val=eps_x) 663 664 DO idir1 = 1, 3 665 CALL pw_pool_create_pw(auxbas_pw_pool, & 666 bind_pw_spline(idir1, iB)%pw, & 667 use_data=REALDATA3D, in_space=REALSPACE) 668 ! calculate spline coefficients 669 CALL pw_spline_precond_create(precond, precond_kind=aint_precond, & 670 pool=auxbas_pw_pool, pbc=.TRUE., transpose=.FALSE.) 671 CALL qs_rho_get(bind_set(idir1, iB)%rho, rho_r=brho1_r) 672 CALL pw_spline_do_precond(precond, brho1_r(1)%pw, & 673 bind_pw_spline(idir1, iB)%pw) 674 CALL pw_spline_precond_set_kind(precond, precond_kind) 675 success = find_coeffs(values=brho1_r(1)%pw, & 676 coeffs=bind_pw_spline(idir1, iB)%pw, linOp=spl3_pbc, & 677 preconditioner=precond, pool=auxbas_pw_pool, & 678 eps_r=eps_r, eps_x=eps_x, max_iter=max_iter) 679 CPASSERT(success) 680 CALL pw_spline_precond_release(precond) 681 END DO ! idir1 682 683 temp_soo_gapw = 0.0_dp 684 685 DO ikind = 1, nkind 686 NULLIFY (atom_list, grid_atom, harmonics) 687 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom) 688 CALL get_qs_kind(qs_kind_set(ikind), & 689 hard_radius=hard_radius, & 690 grid_atom=grid_atom, & 691 harmonics=harmonics, & 692 paw_atom=paw_atom) 693 694 IF (.NOT. paw_atom) CYCLE 695 696 ! Distribute the atoms of this kind 697 698 bo = get_limit(natom, para_env%num_pe, para_env%mepos) 699 700 DO iat = 1, natom !bo(1),bo(2)! this partitioning blocks the interpolation 701 ! ! routines (i.e. waiting for mp_sum) 702 iatom = atom_list(iat) 703 rho_atom => rho_atom_set(iatom) 704 NULLIFY (rho_rad_h, rho_rad_s) 705 CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=rho_rad_h, & 706 rho_rad_s=rho_rad_s) 707 DO idir1 = 1, 3 708 DO ispin = 1, nspins 709 DO ir = 1, grid_atom%nr 710 711 IF (grid_atom%rad(ir) >= hard_radius) CYCLE 712 713 DO ia = 1, grid_atom%ng_sphere 714 715 ra = particle_set(iatom)%r 716 ra(:) = ra(:) + grid_atom%rad(ir)*harmonics%a(:, ia) 717 bind_ra_idir1 = Eval_Interp_Spl3_pbc(ra, & 718 bind_pw_spline(idir1, iB)%pw) 719 720 IF (iat .LT. bo(1) .OR. iat .GT. bo(2)) CYCLE !quick and dirty: 721 ! !here take care of the partition 722 723 rho_spin = 0.0_dp 724 725 DO iso = 1, harmonics%max_iso_not0 726 rho_spin = rho_spin + & 727 (rho_rad_h(ispin)%r_coef(ir, iso) - & 728 rho_rad_s(ispin)%r_coef(ir, iso))* & 729 harmonics%slm(ia, iso) 730 END DO 731 732 temp_soo_gapw(iB, idir1) = temp_soo_gapw(iB, idir1) + & 733 (-1.0_dp)**(1 + ispin)*( & 734 bind_ra_idir1*rho_spin & 735 )*grid_atom%wr(ir)*grid_atom%wa(ia) 736 737 END DO !ia 738 END DO !ir 739 END DO ! ispin 740 END DO !idir1 741 END DO !iat 742 END DO !ikind 743 744 CALL mp_sum(temp_soo_gapw, para_env%group) 745 temp_soo_gapw(:, :) = 1.0_dp*epr_env%g_soo_factor*temp_soo_gapw(:, :) 746 747 IF (output_unit > 0) THEN 748 DO idir1 = 1, 3 749 WRITE (UNIT=output_unit, FMT="(T2,A,T18,I1,I1,T56,E24.16)") & 750 "epr|SOO:gapw", iB, idir1, temp_soo_gapw(iB, idir1) 751 END DO 752 END IF 753 754 g_soo(iB, :) = g_soo(iB, :) + temp_soo_gapw(iB, :) 755 756 DO idir1 = 1, 3 757 CALL pw_pool_give_back_pw(auxbas_pw_pool, bind_pw_spline(idir1, iB)%pw) 758 END DO 759 760 DEALLOCATE (bind_pw_spline) 761 762 END IF ! gapw 763 764 g_total(iB, :) = g_total(iB, :) + g_soo(iB, :) 765 766 CALL cp_print_key_finished_output(output_unit, logger, lr_section, & 767 "PRINT%PROGRAM_RUN_INFO") 768 769 END SUBROUTINE epr_g_soo 770 771! ************************************************************************************************** 772!> \brief ... 773!> \param epr_env ... 774!> \param current_env ... 775!> \param qs_env ... 776!> \param iB ... 777! ************************************************************************************************** 778 SUBROUTINE epr_ind_magnetic_field(epr_env, current_env, qs_env, iB) 779 780 TYPE(epr_env_type) :: epr_env 781 TYPE(current_env_type) :: current_env 782 TYPE(qs_environment_type), POINTER :: qs_env 783 INTEGER, INTENT(IN) :: iB 784 785 CHARACTER(LEN=*), PARAMETER :: routineN = 'epr_ind_magnetic_field', & 786 routineP = moduleN//':'//routineN 787 788 INTEGER :: handle, idir, idir2, idir3, iiB, iiiB, & 789 ispin, natom, nspins 790 LOGICAL :: gapw 791 REAL(dp) :: scale_fac 792 TYPE(cell_type), POINTER :: cell 793 TYPE(dft_control_type), POINTER :: dft_control 794 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 795 TYPE(pw_env_type), POINTER :: pw_env 796 TYPE(pw_p_type) :: pw_gspace_work, shift_pw_rspace 797 TYPE(pw_p_type), DIMENSION(:), POINTER :: epr_rho_r, jrho1_g 798 TYPE(pw_p_type), DIMENSION(:, :), POINTER :: shift_pw_gspace 799 TYPE(pw_p_type), POINTER :: rho_gspace 800 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools 801 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 802 TYPE(realspace_grid_desc_type), POINTER :: auxbas_rs_desc 803 804 CALL timeset(routineN, handle) 805 806 NULLIFY (cell, dft_control, pw_env, auxbas_rs_desc, auxbas_pw_pool, & 807 rho_gspace, pw_pools, particle_set, jrho1_g, epr_rho_r) 808 809 CALL get_qs_env(qs_env=qs_env, cell=cell, dft_control=dft_control, & 810 particle_set=particle_set) 811 812 gapw = dft_control%qs_control%gapw 813 natom = SIZE(particle_set, 1) 814 nspins = dft_control%nspins 815 816 CALL get_epr_env(epr_env=epr_env) 817 818 CALL get_current_env(current_env=current_env) 819 820 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env) 821 CALL pw_env_get(pw_env, auxbas_rs_desc=auxbas_rs_desc, & 822 auxbas_pw_pool=auxbas_pw_pool, pw_pools=pw_pools) 823 ! 824 ! Initialize 825 ! Allocate grids for the calculation of jrho and the shift 826 ALLOCATE (shift_pw_gspace(3, nspins)) 827 DO ispin = 1, nspins 828 DO idir = 1, 3 829 CALL pw_pool_create_pw(auxbas_pw_pool, shift_pw_gspace(idir, ispin)%pw, & 830 use_data=COMPLEXDATA1D, & 831 in_space=RECIPROCALSPACE) 832 CALL pw_zero(shift_pw_gspace(idir, ispin)%pw) 833 ENDDO 834 ENDDO 835 CALL pw_pool_create_pw(auxbas_pw_pool, shift_pw_rspace%pw, & 836 use_data=REALDATA3D, in_space=REALSPACE) 837 CALL pw_zero(shift_pw_rspace%pw) 838 CALL pw_pool_create_pw(auxbas_pw_pool, pw_gspace_work%pw, & 839 use_data=COMPLEXDATA1D, & 840 in_space=RECIPROCALSPACE) 841 842 CALL pw_zero(pw_gspace_work%pw) 843 ! 844 CALL set_vecp(iB, iiB, iiiB) 845 ! 846 DO ispin = 1, nspins 847 ! 848 DO idir3 = 1, 3 849 ! set to zero for the calculation of the shift 850 CALL pw_zero(shift_pw_gspace(idir3, ispin)%pw) 851 ENDDO 852 DO idir = 1, 3 853 CALL qs_rho_get(current_env%jrho1_set(idir)%rho, rho_g=jrho1_g) 854 rho_gspace => jrho1_g(ispin) 855 ! Field gradient 856 ! loop over the Gvec components: x,y,z 857 DO idir2 = 1, 3 858 IF (idir /= idir2) THEN 859 ! in reciprocal space multiply (G_idir2(i)/G(i)^2)J_(idir)(G(i)) 860 CALL mult_G_ov_G2_grid(cell, auxbas_pw_pool, rho_gspace, & 861 pw_gspace_work, idir2, 0.0_dp) 862 ! 863 ! scale and add to the correct component of the shift column 864 CALL set_vecp_rev(idir, idir2, idir3) 865 scale_fac = fac_vecp(idir3, idir2, idir) 866 CALL pw_scale(pw_gspace_work%pw, scale_fac) 867 CALL pw_axpy(pw_gspace_work%pw, shift_pw_gspace(idir3, ispin)%pw) 868 ENDIF 869 ENDDO 870 ! 871 ENDDO ! idir 872 ENDDO ! ispin 873 874 ! Store the total soft induced magnetic field (corrected for sic) 875 IF (dft_control%nspins == 2) THEN 876 DO idir = 1, 3 877 CALL qs_rho_get(epr_env%bind_set(idir, iB)%rho, rho_r=epr_rho_r) 878 CALL pw_transfer(shift_pw_gspace(idir, 2)%pw, epr_rho_r(1)%pw) 879 ENDDO 880 END IF 881 ! 882 ! Dellocate grids for the calculation of jrho and the shift 883 CALL pw_pool_give_back_pw(auxbas_pw_pool, pw_gspace_work%pw) 884 DO ispin = 1, dft_control%nspins 885 DO idir = 1, 3 886 CALL pw_pool_give_back_pw(auxbas_pw_pool, shift_pw_gspace(idir, ispin)%pw) 887 ENDDO 888 ENDDO 889 DEALLOCATE (shift_pw_gspace) 890 CALL pw_pool_give_back_pw(auxbas_pw_pool, shift_pw_rspace%pw) 891 ! 892 ! Finalize 893 CALL timestop(handle) 894 ! 895 END SUBROUTINE epr_ind_magnetic_field 896 897END MODULE qs_linres_epr_ownutils 898 899