1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief from the response current density calculates the shift tensor 8!> and the susceptibility 9!> \par History 10!> created 02-2006 [MI] 11!> \author MI 12! ************************************************************************************************** 13MODULE qs_linres_nmr_shift 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_get_default_io_unit,& 21 cp_logger_type 22 USE cp_output_handling, ONLY: cp_p_file,& 23 cp_print_key_finished_output,& 24 cp_print_key_should_output,& 25 cp_print_key_unit_nr 26 USE cp_para_types, ONLY: cp_para_env_type 27 USE input_section_types, ONLY: section_vals_get_subs_vals,& 28 section_vals_type,& 29 section_vals_val_get 30 USE kinds, ONLY: default_string_length,& 31 dp 32 USE mathconstants, ONLY: gaussi,& 33 twopi 34 USE mathlib, ONLY: diamat_all 35 USE message_passing, ONLY: mp_sum 36 USE particle_types, ONLY: particle_type 37 USE pw_env_types, ONLY: pw_env_get,& 38 pw_env_type 39 USE pw_grid_types, ONLY: pw_grid_type 40 USE pw_methods, ONLY: pw_axpy,& 41 pw_transfer,& 42 pw_zero 43 USE pw_pool_types, ONLY: pw_pool_create_pw,& 44 pw_pool_give_back_pw,& 45 pw_pool_p_type,& 46 pw_pool_type 47 USE pw_spline_utils, ONLY: Eval_Interp_Spl3_pbc,& 48 find_coeffs,& 49 pw_spline_do_precond,& 50 pw_spline_precond_create,& 51 pw_spline_precond_release,& 52 pw_spline_precond_set_kind,& 53 pw_spline_precond_type,& 54 spl3_pbc 55 USE pw_types, ONLY: COMPLEXDATA1D,& 56 REALDATA3D,& 57 REALSPACE,& 58 RECIPROCALSPACE,& 59 pw_p_type,& 60 pw_type 61 USE qs_environment_types, ONLY: get_qs_env,& 62 qs_environment_type 63 USE qs_grid_atom, ONLY: grid_atom_type 64 USE qs_harmonics_atom, ONLY: harmonics_atom_type 65 USE qs_kind_types, ONLY: get_qs_kind,& 66 qs_kind_type 67 USE qs_linres_nmr_epr_common_utils, ONLY: mult_G_ov_G2_grid 68 USE qs_linres_op, ONLY: fac_vecp,& 69 set_vecp,& 70 set_vecp_rev 71 USE qs_linres_types, ONLY: current_env_type,& 72 get_current_env,& 73 get_nmr_env,& 74 jrho_atom_type,& 75 nmr_env_type 76 USE qs_rho_types, ONLY: qs_rho_get 77 USE realspace_grid_types, ONLY: realspace_grid_desc_type 78 USE util, ONLY: get_limit 79#include "./base/base_uses.f90" 80 81 IMPLICIT NONE 82 83 PRIVATE 84 85 ! *** Public subroutines *** 86 PUBLIC :: nmr_shift_print, & 87 nmr_shift 88 89 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_nmr_shift' 90 91! ************** 92 93CONTAINS 94 95! ************************************************************************************************** 96!> \brief ... 97!> \param nmr_env ... 98!> \param current_env ... 99!> \param qs_env ... 100!> \param iB ... 101! ************************************************************************************************** 102 SUBROUTINE nmr_shift(nmr_env, current_env, qs_env, iB) 103 104 TYPE(nmr_env_type) :: nmr_env 105 TYPE(current_env_type) :: current_env 106 TYPE(qs_environment_type), POINTER :: qs_env 107 INTEGER, INTENT(IN) :: iB 108 109 CHARACTER(LEN=*), PARAMETER :: routineN = 'nmr_shift', routineP = moduleN//':'//routineN 110 111 INTEGER :: handle, idir, idir2, idir3, iiB, iiiB, & 112 ispin, natom, nspins 113 LOGICAL :: gapw, interpolate_shift 114 REAL(dp) :: scale_fac 115 REAL(dp), DIMENSION(:, :, :), POINTER :: chemical_shift, chemical_shift_loc, & 116 chemical_shift_nics, & 117 chemical_shift_nics_loc 118 TYPE(cell_type), POINTER :: cell 119 TYPE(dft_control_type), POINTER :: dft_control 120 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 121 TYPE(pw_env_type), POINTER :: pw_env 122 TYPE(pw_p_type) :: pw_gspace_work, shift_pw_rspace 123 TYPE(pw_p_type), DIMENSION(:), POINTER :: jrho1_g 124 TYPE(pw_p_type), DIMENSION(:, :), POINTER :: shift_pw_gspace 125 TYPE(pw_p_type), POINTER :: jrho_gspace 126 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools 127 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 128 TYPE(realspace_grid_desc_type), POINTER :: auxbas_rs_desc 129 TYPE(section_vals_type), POINTER :: nmr_section 130 131 CALL timeset(routineN, handle) 132 133 NULLIFY (chemical_shift, chemical_shift_loc, chemical_shift_nics, chemical_shift_nics_loc, & 134 cell, dft_control, pw_env, auxbas_rs_desc, auxbas_pw_pool, jrho_gspace, & 135 pw_pools, particle_set, jrho1_g) 136 137 CALL get_qs_env(qs_env=qs_env, cell=cell, dft_control=dft_control, & 138 particle_set=particle_set) 139 140 gapw = dft_control%qs_control%gapw 141 natom = SIZE(particle_set, 1) 142 nspins = dft_control%nspins 143 144 CALL get_nmr_env(nmr_env=nmr_env, chemical_shift=chemical_shift, & 145 chemical_shift_loc=chemical_shift_loc, & 146 chemical_shift_nics=chemical_shift_nics, & 147 chemical_shift_nics_loc=chemical_shift_nics_loc, & 148 interpolate_shift=interpolate_shift) 149 150 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env) 151 CALL pw_env_get(pw_env, auxbas_rs_desc=auxbas_rs_desc, & 152 auxbas_pw_pool=auxbas_pw_pool, pw_pools=pw_pools) 153 ! 154 ! 155 nmr_section => section_vals_get_subs_vals(qs_env%input, & 156 & "PROPERTIES%LINRES%NMR") 157 ! 158 ! Initialize 159 ! Allocate grids for the calculation of jrho and the shift 160 ALLOCATE (shift_pw_gspace(3, nspins)) 161 DO ispin = 1, nspins 162 DO idir = 1, 3 163 CALL pw_pool_create_pw(auxbas_pw_pool, shift_pw_gspace(idir, ispin)%pw, & 164 use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE) 165 CALL pw_zero(shift_pw_gspace(idir, ispin)%pw) 166 ENDDO 167 ENDDO 168 ! 169 ! 170 CALL set_vecp(iB, iiB, iiiB) 171 ! 172 CALL pw_pool_create_pw(auxbas_pw_pool, pw_gspace_work%pw, & 173 use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE) 174 CALL pw_zero(pw_gspace_work%pw) 175 DO ispin = 1, nspins 176 ! 177 DO idir = 1, 3 178 CALL qs_rho_get(current_env%jrho1_set(idir)%rho, rho_g=jrho1_g) 179 jrho_gspace => jrho1_g(ispin) 180 ! Field gradient 181 ! loop over the Gvec components: x,y,z 182 DO idir2 = 1, 3 183 IF (idir /= idir2) THEN 184 ! in reciprocal space multiply (G_idir2(i)/G(i)^2)J_(idir)(G(i)) 185 CALL mult_G_ov_G2_grid(cell, auxbas_pw_pool, jrho_gspace, & 186 pw_gspace_work, idir2, 0.0_dp) 187 ! 188 ! scale and add to the correct component of the shift column 189 CALL set_vecp_rev(idir, idir2, idir3) 190 scale_fac = fac_vecp(idir3, idir2, idir) 191 CALL pw_axpy(pw_gspace_work%pw, shift_pw_gspace(idir3, ispin)%pw, scale_fac) 192 ENDIF 193 ENDDO 194 ! 195 ENDDO ! idir 196 ENDDO ! ispin 197 ! 198 CALL pw_pool_give_back_pw(auxbas_pw_pool, pw_gspace_work%pw) 199 ! 200 ! compute shildings 201 IF (interpolate_shift) THEN 202 CALL pw_pool_create_pw(auxbas_pw_pool, shift_pw_rspace%pw, & 203 use_data=REALDATA3D, in_space=REALSPACE) 204 DO ispin = 1, nspins 205 DO idir = 1, 3 206 ! Here first G->R and then interpolation to get the shifts. 207 ! The interpolation doesnt work in parallel yet. 208 ! The choice between both methods should be left to the user. 209 CALL pw_transfer(shift_pw_gspace(idir, ispin)%pw, shift_pw_rspace%pw) 210 CALL interpolate_shift_pwgrid(nmr_env, pw_env, particle_set, cell, shift_pw_rspace, & 211 iB, idir, nmr_section) 212 ENDDO 213 ENDDO 214 CALL pw_pool_give_back_pw(auxbas_pw_pool, shift_pw_rspace%pw) 215 ELSE 216 DO ispin = 1, nspins 217 DO idir = 1, 3 218 ! Here the shifts are computed from summation of the coeff on the G-grip . 219 CALL gsum_shift_pwgrid(nmr_env, particle_set, cell, & 220 shift_pw_gspace(idir, ispin), iB, idir) 221 ENDDO 222 ENDDO 223 ENDIF 224 ! 225 IF (gapw) THEN 226 DO idir = 1, 3 227 ! Finally the radial functions are multiplied by the YLM and properly summed 228 ! The resulting array is J on the local grid. One array per atom. 229 ! Local contributions by numerical integration over the spherical grids 230 CALL nmr_shift_gapw(nmr_env, current_env, qs_env, iB, idir) 231 ENDDO ! idir 232 ENDIF 233 ! 234 ! Dellocate grids for the calculation of jrho and the shift 235 DO ispin = 1, nspins 236 DO idir = 1, 3 237 CALL pw_pool_give_back_pw(auxbas_pw_pool, shift_pw_gspace(idir, ispin)%pw) 238 ENDDO 239 ENDDO 240 DEALLOCATE (shift_pw_gspace) 241 ! 242 ! Finalize 243 CALL timestop(handle) 244 ! 245 END SUBROUTINE nmr_shift 246 247! ************************************************************************************************** 248!> \brief ... 249!> \param nmr_env ... 250!> \param current_env ... 251!> \param qs_env ... 252!> \param iB ... 253!> \param idir ... 254! ************************************************************************************************** 255 SUBROUTINE nmr_shift_gapw(nmr_env, current_env, qs_env, iB, idir) 256 257 TYPE(nmr_env_type) :: nmr_env 258 TYPE(current_env_type) :: current_env 259 TYPE(qs_environment_type), POINTER :: qs_env 260 INTEGER, INTENT(IN) :: IB, idir 261 262 CHARACTER(len=*), PARAMETER :: routineN = 'nmr_shift_gapw', routineP = moduleN//':'//routineN 263 264 INTEGER :: handle, ia, iat, iatom, idir2_1, idir3_1, ikind, ir, ira, ispin, j, jatom, mepos, & 265 n_nics, na, natom, natom_local, natom_tot, nkind, nnics_local, nr, nra, nspins, num_pe, & 266 output_unit 267 INTEGER, ALLOCATABLE, DIMENSION(:) :: list_j, list_nics_j 268 INTEGER, DIMENSION(2) :: bo 269 INTEGER, DIMENSION(:), POINTER :: atom_list 270 LOGICAL :: do_nics, paw_atom 271 REAL(dp) :: ddiff, dist, dum, itegrated_jrho, & 272 r_jatom(3), rdiff(3), rij(3), s_1, & 273 s_2, scale_fac_1, shift_gapw_radius 274 REAL(dp), ALLOCATABLE, DIMENSION(:) :: j_grid 275 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: cs_loc_tmp, cs_nics_loc_tmp, dist_ij, & 276 dist_nics_ij, r_grid 277 REAL(dp), DIMENSION(:, :), POINTER :: jrho_h_grid, jrho_s_grid, r_nics 278 REAL(dp), DIMENSION(:, :, :), POINTER :: chemical_shift_loc, & 279 chemical_shift_nics_loc 280 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 281 TYPE(cell_type), POINTER :: cell 282 TYPE(cp_logger_type), POINTER :: logger 283 TYPE(cp_para_env_type), POINTER :: para_env 284 TYPE(dft_control_type), POINTER :: dft_control 285 TYPE(grid_atom_type), POINTER :: grid_atom 286 TYPE(harmonics_atom_type), POINTER :: harmonics 287 TYPE(jrho_atom_type), DIMENSION(:), POINTER :: jrho1_atom_set 288 TYPE(jrho_atom_type), POINTER :: jrho1_atom 289 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 290 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 291 292 CALL timeset(routineN, handle) 293 ! 294 NULLIFY (atomic_kind_set, qs_kind_set, cell, dft_control, para_env, particle_set, & 295 chemical_shift_loc, chemical_shift_nics_loc, jrho1_atom_set, & 296 jrho1_atom, r_nics, jrho_h_grid, jrho_s_grid, & 297 atom_list, grid_atom, harmonics, logger) 298 ! 299 logger => cp_get_default_logger() 300 output_unit = cp_logger_get_default_io_unit(logger) 301 ! 302 CALL get_qs_env(qs_env=qs_env, & 303 atomic_kind_set=atomic_kind_set, & 304 qs_kind_set=qs_kind_set, & 305 cell=cell, & 306 dft_control=dft_control, & 307 para_env=para_env, & 308 particle_set=particle_set) 309 310 CALL get_nmr_env(nmr_env=nmr_env, & 311 chemical_shift_loc=chemical_shift_loc, & 312 chemical_shift_nics_loc=chemical_shift_nics_loc, & 313 shift_gapw_radius=shift_gapw_radius, & 314 n_nics=n_nics, & 315 r_nics=r_nics, & 316 do_nics=do_nics) 317 318 CALL get_current_env(current_env=current_env, & 319 jrho1_atom_set=jrho1_atom_set) 320 ! 321 nkind = SIZE(atomic_kind_set, 1) 322 natom_tot = SIZE(particle_set, 1) 323 nspins = dft_control%nspins 324 itegrated_jrho = 0.0_dp 325 ! 326 idir2_1 = MODULO(idir, 3) + 1 327 idir3_1 = MODULO(idir + 1, 3) + 1 328 scale_fac_1 = fac_vecp(idir3_1, idir2_1, idir) 329 ! 330 ALLOCATE (cs_loc_tmp(3, natom_tot), list_j(natom_tot), & 331 dist_ij(3, natom_tot)) 332 cs_loc_tmp = 0.0_dp 333 IF (do_nics) THEN 334 ALLOCATE (cs_nics_loc_tmp(3, n_nics), list_nics_j(n_nics), & 335 dist_nics_ij(3, n_nics)) 336 cs_nics_loc_tmp = 0.0_dp 337 ENDIF 338 ! 339 ! Loop over atoms to collocate the current on each atomic grid, JA 340 ! Per each JA, loop over the points where the shift needs to be computed 341 DO ikind = 1, nkind 342 343 NULLIFY (atom_list, grid_atom, harmonics) 344 CALL get_atomic_kind(atomic_kind_set(ikind), & 345 atom_list=atom_list, & 346 natom=natom) 347 348 CALL get_qs_kind(qs_kind_set(ikind), & 349 paw_atom=paw_atom, & 350 harmonics=harmonics, & 351 grid_atom=grid_atom) 352 ! 353 na = grid_atom%ng_sphere 354 nr = grid_atom%nr 355 nra = nr*na 356 ALLOCATE (r_grid(3, nra), j_grid(nra)) 357 ira = 1 358 DO ia = 1, na 359 DO ir = 1, nr 360 r_grid(:, ira) = grid_atom%rad(ir)*harmonics%a(:, ia) 361 ira = ira + 1 362 ENDDO 363 ENDDO 364 ! 365 ! Quick cycle if needed 366 IF (paw_atom) THEN 367 ! 368 ! Distribute the atoms of this kind 369 num_pe = para_env%num_pe 370 mepos = para_env%mepos 371 bo = get_limit(natom, num_pe, mepos) 372 ! 373 DO iat = bo(1), bo(2) 374 iatom = atom_list(iat) 375 ! 376 ! find all the atoms within the radius 377 natom_local = 0 378 DO jatom = 1, natom_tot 379 rij(:) = pbc(particle_set(iatom)%r, particle_set(jatom)%r, cell) 380 dist = SQRT(rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3)) 381 IF (dist .LE. shift_gapw_radius) THEN 382 natom_local = natom_local + 1 383 list_j(natom_local) = jatom 384 dist_ij(:, natom_local) = rij(:) 385 ENDIF 386 ENDDO 387 ! 388 ! ... also for nics 389 IF (do_nics) THEN 390 nnics_local = 0 391 DO jatom = 1, n_nics 392 r_jatom(:) = r_nics(:, jatom) 393 rij(:) = pbc(particle_set(iatom)%r, r_jatom, cell) 394 dist = SQRT(rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3)) 395 IF (dist .LE. shift_gapw_radius) THEN 396 nnics_local = nnics_local + 1 397 list_nics_j(nnics_local) = jatom 398 dist_nics_ij(:, nnics_local) = rij(:) 399 ENDIF 400 ENDDO 401 ENDIF 402 ! 403 NULLIFY (jrho1_atom, jrho_h_grid, jrho_s_grid) 404 jrho1_atom => jrho1_atom_set(iatom) 405 ! 406 DO ispin = 1, nspins 407 jrho_h_grid => jrho1_atom%jrho_vec_rad_h(idir, ispin)%r_coef 408 jrho_s_grid => jrho1_atom%jrho_vec_rad_s(idir, ispin)%r_coef 409 ! 410 ! loop over the atoms neighbors of iatom in terms of the current density 411 ! for each compute the contribution to the shift coming from the 412 ! local current density at iatom 413 ira = 1 414 DO ia = 1, na 415 DO ir = 1, nr 416 j_grid(ira) = (jrho_h_grid(ir, ia) - jrho_s_grid(ir, ia))*grid_atom%weight(ia, ir) 417 itegrated_jrho = itegrated_jrho + j_grid(ira) 418 ira = ira + 1 419 ENDDO 420 ENDDO 421 ! 422 DO j = 1, natom_local 423 jatom = list_j(j) 424 rij(:) = dist_ij(:, j) 425 ! 426 s_1 = 0.0_dp 427 s_2 = 0.0_dp 428 DO ira = 1, nra 429 ! 430 rdiff(:) = rij(:) - r_grid(:, ira) 431 ddiff = SQRT(rdiff(1)*rdiff(1) + rdiff(2)*rdiff(2) + rdiff(3)*rdiff(3)) 432 IF (ddiff .GT. 1.0E-12_dp) THEN 433 dum = scale_fac_1*j_grid(ira)/(ddiff*ddiff*ddiff) 434 s_1 = s_1 + rdiff(idir2_1)*dum 435 s_2 = s_2 + rdiff(idir3_1)*dum 436 ENDIF ! ddiff 437 ENDDO ! ira 438 cs_loc_tmp(idir3_1, jatom) = cs_loc_tmp(idir3_1, jatom) + s_1 439 cs_loc_tmp(idir2_1, jatom) = cs_loc_tmp(idir2_1, jatom) - s_2 440 ENDDO ! j 441 ! 442 IF (do_nics) THEN 443 DO j = 1, nnics_local 444 jatom = list_nics_j(j) 445 rij(:) = dist_nics_ij(:, j) 446 ! 447 s_1 = 0.0_dp 448 s_2 = 0.0_dp 449 DO ira = 1, nra 450 ! 451 rdiff(:) = rij(:) - r_grid(:, ira) 452 ddiff = SQRT(rdiff(1)*rdiff(1) + rdiff(2)*rdiff(2) + rdiff(3)*rdiff(3)) 453 IF (ddiff .GT. 1.0E-12_dp) THEN 454 dum = scale_fac_1*j_grid(ira)/(ddiff*ddiff*ddiff) 455 s_1 = s_1 + rdiff(idir2_1)*dum 456 s_2 = s_2 + rdiff(idir3_1)*dum 457 ENDIF ! ddiff 458 ENDDO ! ira 459 cs_nics_loc_tmp(idir3_1, jatom) = cs_nics_loc_tmp(idir3_1, jatom) + s_1 460 cs_nics_loc_tmp(idir2_1, jatom) = cs_nics_loc_tmp(idir2_1, jatom) - s_2 461 ENDDO ! j 462 ENDIF ! do_nics 463 ENDDO ! ispin 464 ENDDO ! iat 465 ENDIF 466 DEALLOCATE (r_grid, j_grid) 467 ENDDO ! ikind 468 ! 469 ! 470 CALL mp_sum(itegrated_jrho, para_env%group) 471 IF (output_unit > 0) THEN 472 WRITE (output_unit, '(T2,A,E24.16)') 'Integrated local j_'& 473 &//ACHAR(idir + 119)//ACHAR(iB + 119)//'(r)=', itegrated_jrho 474 ENDIF 475 ! 476 CALL mp_sum(cs_loc_tmp, para_env%group) 477 chemical_shift_loc(:, iB, :) = chemical_shift_loc(:, iB, :) & 478 & - nmr_env%shift_factor_gapw*cs_loc_tmp(:, :)/2.0_dp 479 ! 480 DEALLOCATE (cs_loc_tmp, list_j, dist_ij) 481 ! 482 IF (do_nics) THEN 483 CALL mp_sum(cs_nics_loc_tmp, para_env%group) 484 chemical_shift_nics_loc(:, iB, :) = chemical_shift_nics_loc(:, iB, :) & 485 & - nmr_env%shift_factor_gapw*cs_nics_loc_tmp(:, :)/2.0_dp 486 ! 487 DEALLOCATE (cs_nics_loc_tmp, list_nics_j, dist_nics_ij) 488 ENDIF 489 ! 490 CALL timestop(handle) 491 ! 492 END SUBROUTINE nmr_shift_gapw 493 494! ************************************************************************************************** 495!> \brief interpolate the shift calculated on the PW grid in order to ger 496!> the value on arbitrary points in real space 497!> \param nmr_env to get the shift tensor and the list of additional points 498!> \param pw_env ... 499!> \param particle_set for the atomic position 500!> \param cell to take into account the pbs, and to have the volume 501!> \param shift_pw_rspace specific component of the shift tensor on the pw grid 502!> \param i_B component of the magnetic field for which the shift is calculated (row) 503!> \param idir component of the vector \int_{r}[ ((r-r') x j(r))/|r-r'|^3 ] = Bind(r') 504!> \param nmr_section ... 505!> \author MI 506! ************************************************************************************************** 507 SUBROUTINE interpolate_shift_pwgrid(nmr_env, pw_env, particle_set, cell, shift_pw_rspace, & 508 i_B, idir, nmr_section) 509 510 TYPE(nmr_env_type) :: nmr_env 511 TYPE(pw_env_type), POINTER :: pw_env 512 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 513 TYPE(cell_type), POINTER :: cell 514 TYPE(pw_p_type) :: shift_pw_rspace 515 INTEGER, INTENT(IN) :: i_B, idir 516 TYPE(section_vals_type), POINTER :: nmr_section 517 518 CHARACTER(LEN=*), PARAMETER :: routineN = 'interpolate_shift_pwgrid', & 519 routineP = moduleN//':'//routineN 520 521 INTEGER :: aint_precond, handle, iat, iatom, & 522 max_iter, n_nics, natom, precond_kind 523 INTEGER, DIMENSION(:), POINTER :: cs_atom_list 524 LOGICAL :: do_nics, success 525 REAL(dp) :: eps_r, eps_x, R_iatom(3), ra(3), & 526 shift_val 527 REAL(dp), DIMENSION(:, :), POINTER :: r_nics 528 REAL(dp), DIMENSION(:, :, :), POINTER :: chemical_shift, chemical_shift_nics 529 TYPE(pw_p_type) :: shiftspl 530 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 531 TYPE(pw_spline_precond_type), POINTER :: precond 532 TYPE(section_vals_type), POINTER :: interp_section 533 534! 535 536 CALL timeset(routineN, handle) 537 ! 538 NULLIFY (interp_section) 539 NULLIFY (auxbas_pw_pool, precond) 540 NULLIFY (cs_atom_list, chemical_shift, chemical_shift_nics, r_nics) 541 542 CPASSERT(ASSOCIATED(shift_pw_rspace%pw)) 543 544 interp_section => section_vals_get_subs_vals(nmr_section, & 545 "INTERPOLATOR") 546 CALL section_vals_val_get(interp_section, "aint_precond", & 547 i_val=aint_precond) 548 CALL section_vals_val_get(interp_section, "precond", i_val=precond_kind) 549 CALL section_vals_val_get(interp_section, "max_iter", i_val=max_iter) 550 CALL section_vals_val_get(interp_section, "eps_r", r_val=eps_r) 551 CALL section_vals_val_get(interp_section, "eps_x", r_val=eps_x) 552 553 ! calculate spline coefficients 554 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool) 555 CALL pw_pool_create_pw(auxbas_pw_pool, shiftspl%pw, & 556 use_data=REALDATA3D, in_space=REALSPACE) 557 558 CALL pw_spline_precond_create(precond, precond_kind=aint_precond, & 559 pool=auxbas_pw_pool, pbc=.TRUE., transpose=.FALSE.) 560 CALL pw_spline_do_precond(precond, shift_pw_rspace%pw, shiftspl%pw) 561 CALL pw_spline_precond_set_kind(precond, precond_kind) 562 success = find_coeffs(values=shift_pw_rspace%pw, coeffs=shiftspl%pw, & 563 linOp=spl3_pbc, preconditioner=precond, pool=auxbas_pw_pool, & 564 eps_r=eps_r, eps_x=eps_x, max_iter=max_iter) 565 CPASSERT(success) 566 CALL pw_spline_precond_release(precond) 567 568 CALL get_nmr_env(nmr_env=nmr_env, cs_atom_list=cs_atom_list, & 569 chemical_shift=chemical_shift, & 570 chemical_shift_nics=chemical_shift_nics, & 571 n_nics=n_nics, r_nics=r_nics, & 572 do_nics=do_nics) 573 574 IF (ASSOCIATED(cs_atom_list)) THEN 575 natom = SIZE(cs_atom_list, 1) 576 ELSE 577 natom = -1 578 ENDIF 579 580 DO iat = 1, natom 581 iatom = cs_atom_list(iat) 582 R_iatom = pbc(particle_set(iatom)%r, cell) 583 shift_val = Eval_Interp_Spl3_pbc(R_iatom, shiftspl%pw) 584 chemical_shift(idir, i_B, iatom) = chemical_shift(idir, i_B, iatom) + & 585 nmr_env%shift_factor*twopi**2*shift_val 586 END DO 587 588 IF (do_nics) THEN 589 DO iatom = 1, n_nics 590 ra(:) = r_nics(:, iatom) 591 R_iatom = pbc(ra, cell) 592 shift_val = Eval_Interp_Spl3_pbc(R_iatom, shiftspl%pw) 593 chemical_shift_nics(idir, i_B, iatom) = chemical_shift_nics(idir, i_B, iatom) + & 594 nmr_env%shift_factor*twopi**2*shift_val 595 END DO 596 END IF 597 598 CALL pw_pool_give_back_pw(auxbas_pw_pool, shiftspl%pw) 599 600 ! 601 CALL timestop(handle) 602 ! 603 END SUBROUTINE interpolate_shift_pwgrid 604 605! ************************************************************************************************** 606!> \brief ... 607!> \param nmr_env ... 608!> \param particle_set ... 609!> \param cell ... 610!> \param shift_pw_gspace ... 611!> \param i_B ... 612!> \param idir ... 613! ************************************************************************************************** 614 SUBROUTINE gsum_shift_pwgrid(nmr_env, particle_set, cell, shift_pw_gspace, & 615 i_B, idir) 616 TYPE(nmr_env_type) :: nmr_env 617 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 618 TYPE(cell_type), POINTER :: cell 619 TYPE(pw_p_type) :: shift_pw_gspace 620 INTEGER, INTENT(IN) :: i_B, idir 621 622 CHARACTER(LEN=*), PARAMETER :: routineN = 'gsum_shift_pwgrid', & 623 routineP = moduleN//':'//routineN 624 625 COMPLEX(dp) :: cplx 626 INTEGER :: handle, iat, iatom, n_nics, natom 627 INTEGER, DIMENSION(:), POINTER :: cs_atom_list 628 LOGICAL :: do_nics 629 REAL(dp) :: R_iatom(3), ra(3) 630 REAL(dp), DIMENSION(:, :), POINTER :: r_nics 631 REAL(dp), DIMENSION(:, :, :), POINTER :: chemical_shift, chemical_shift_nics 632 633! 634! 635 636 CALL timeset(routineN, handle) 637 ! 638 NULLIFY (cs_atom_list, chemical_shift, chemical_shift_nics, r_nics) 639 CPASSERT(ASSOCIATED(shift_pw_gspace%pw)) 640 ! 641 CALL get_nmr_env(nmr_env=nmr_env, cs_atom_list=cs_atom_list, & 642 chemical_shift=chemical_shift, & 643 chemical_shift_nics=chemical_shift_nics, & 644 n_nics=n_nics, r_nics=r_nics, do_nics=do_nics) 645 ! 646 IF (ASSOCIATED(cs_atom_list)) THEN 647 natom = SIZE(cs_atom_list, 1) 648 ELSE 649 natom = -1 650 ENDIF 651 ! 652 ! compute the chemical shift 653 DO iat = 1, natom 654 iatom = cs_atom_list(iat) 655 R_iatom = pbc(particle_set(iatom)%r, cell) 656 CALL gsumr(R_iatom, shift_pw_gspace%pw, cplx) 657 chemical_shift(idir, i_B, iatom) = chemical_shift(idir, i_B, iatom) + & 658 nmr_env%shift_factor*twopi**2*REAL(cplx, dp) 659 ENDDO 660 ! 661 ! compute nics 662 IF (do_nics) THEN 663 DO iat = 1, n_nics 664 ra = pbc(r_nics(:, iat), cell) 665 CALL gsumr(ra, shift_pw_gspace%pw, cplx) 666 chemical_shift_nics(idir, i_B, iat) = chemical_shift_nics(idir, i_B, iat) + & 667 nmr_env%shift_factor*twopi**2*REAL(cplx, dp) 668 ENDDO 669 ENDIF 670 671 CALL timestop(handle) 672 673 END SUBROUTINE gsum_shift_pwgrid 674 675! ************************************************************************************************** 676!> \brief ... 677!> \param r ... 678!> \param pw ... 679!> \param cplx ... 680! ************************************************************************************************** 681 SUBROUTINE gsumr(r, pw, cplx) 682 REAL(dp), INTENT(IN) :: r(3) 683 TYPE(pw_type), POINTER :: pw 684 COMPLEX(dp) :: cplx 685 686 COMPLEX(dp) :: rg 687 INTEGER :: ig 688 TYPE(pw_grid_type), POINTER :: grid 689 690 grid => pw%pw_grid 691 cplx = CMPLX(0.0_dp, 0.0_dp, KIND=dp) 692 DO ig = grid%first_gne0, grid%ngpts_cut_local 693 rg = (grid%g(1, ig)*r(1) + grid%g(2, ig)*r(2) + grid%g(3, ig)*r(3))*gaussi 694 cplx = cplx + pw%cc(ig)*EXP(rg) 695 ENDDO 696 IF (grid%have_g0) cplx = cplx + pw%cc(1) 697 CALL mp_sum(cplx, grid%para%group) 698 END SUBROUTINE gsumr 699 700! ************************************************************************************************** 701!> \brief Shielding tensor and Chi are printed into a file 702!> if required from input 703!> It is possible to print only for a subset of atoms or 704!> or points in non-ionic positions 705!> \param nmr_env ... 706!> \param current_env ... 707!> \param qs_env ... 708!> \author MI 709! ************************************************************************************************** 710 SUBROUTINE nmr_shift_print(nmr_env, current_env, qs_env) 711 TYPE(nmr_env_type) :: nmr_env 712 TYPE(current_env_type) :: current_env 713 TYPE(qs_environment_type), POINTER :: qs_env 714 715 CHARACTER(len=*), PARAMETER :: routineN = 'nmr_shift_print', & 716 routineP = moduleN//':'//routineN 717 718 CHARACTER(LEN=2) :: element_symbol 719 CHARACTER(LEN=default_string_length) :: name, title 720 INTEGER :: iatom, ir, n_nics, nat_print, natom, & 721 output_unit, unit_atoms, unit_nics 722 INTEGER, DIMENSION(:), POINTER :: cs_atom_list 723 LOGICAL :: do_nics, gapw 724 REAL(dp) :: chi_aniso, chi_iso, chi_sym_tot(3, 3), chi_tensor(3, 3, 2), & 725 chi_tensor_loc(3, 3, 2), chi_tensor_loc_tmp(3, 3), chi_tensor_tmp(3, 3), chi_tmp(3, 3), & 726 eig(3), rpos(3), shift_aniso, shift_iso, shift_sym_tot(3, 3) 727 REAL(dp), DIMENSION(:, :), POINTER :: r_nics 728 REAL(dp), DIMENSION(:, :, :), POINTER :: cs, cs_loc, cs_nics, cs_nics_loc, & 729 cs_nics_tot, cs_tot 730 REAL(dp), EXTERNAL :: DDOT 731 TYPE(atomic_kind_type), POINTER :: atom_kind 732 TYPE(cp_logger_type), POINTER :: logger 733 TYPE(dft_control_type), POINTER :: dft_control 734 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 735 TYPE(section_vals_type), POINTER :: nmr_section 736 737 NULLIFY (cs, cs_nics, r_nics, cs_loc, cs_nics_loc, logger, particle_set, atom_kind, dft_control) 738 739 logger => cp_get_default_logger() 740 output_unit = cp_logger_get_default_io_unit(logger) 741 742 nmr_section => section_vals_get_subs_vals(qs_env%input, & 743 "PROPERTIES%LINRES%NMR") 744 745 CALL get_nmr_env(nmr_env=nmr_env, & 746 chemical_shift=cs, & 747 chemical_shift_nics=cs_nics, & 748 chemical_shift_loc=cs_loc, & 749 chemical_shift_nics_loc=cs_nics_loc, & 750 cs_atom_list=cs_atom_list, & 751 n_nics=n_nics, & 752 r_nics=r_nics, & 753 do_nics=do_nics) 754 ! 755 CALL get_current_env(current_env=current_env, & 756 chi_tensor=chi_tensor, & 757 chi_tensor_loc=chi_tensor_loc) 758 ! 759 ! multiply by the appropriate factor 760 chi_tensor_tmp(:, :) = 0.0_dp 761 chi_tensor_loc_tmp(:, :) = 0.0_dp 762 chi_tensor_tmp(:, :) = (chi_tensor(:, :, 1) + chi_tensor(:, :, 2))*nmr_env%chi_factor 763 !chi_tensor_loc_tmp(:,:) = (chi_tensor_loc(:,:,1) + chi_tensor_loc(:,:,2)) * here there is another factor 764 ! 765 CALL get_qs_env(qs_env=qs_env, & 766 dft_control=dft_control, & 767 particle_set=particle_set) 768 769 natom = SIZE(particle_set, 1) 770 gapw = dft_control%qs_control%gapw 771 nat_print = SIZE(cs_atom_list, 1) 772 773 ALLOCATE (cs_tot(3, 3, nat_print)) 774 IF (do_nics) THEN 775 ALLOCATE (cs_nics_tot(3, 3, n_nics)) 776 ENDIF 777 ! Finalize Chi calculation 778 ! Symmetrize 779 chi_sym_tot(:, :) = (chi_tensor_tmp(:, :) + TRANSPOSE(chi_tensor_tmp(:, :)))/2.0_dp 780 IF (gapw) THEN 781 chi_sym_tot(:, :) = chi_sym_tot(:, :) & 782 & + (chi_tensor_loc_tmp(:, :) + TRANSPOSE(chi_tensor_loc_tmp(:, :)))/2.0_dp 783 ENDIF 784 chi_tmp(:, :) = chi_sym_tot(:, :) 785 CALL diamat_all(chi_tmp, eig) 786 chi_iso = (eig(1) + eig(2) + eig(3))/3.0_dp 787 chi_aniso = eig(3) - (eig(2) + eig(1))/2.0_dp 788 ! 789 IF (output_unit > 0) THEN 790 WRITE (output_unit, '(T2,A,E14.6)') 'CheckSum Chi =', & 791 SQRT(DDOT(9, chi_tensor_tmp(1, 1), 1, chi_tensor_tmp(1, 1), 1)) 792 ENDIF 793 ! 794 IF (BTEST(cp_print_key_should_output(logger%iter_info, nmr_section, & 795 "PRINT%CHI_TENSOR"), cp_p_file)) THEN 796 797 unit_atoms = cp_print_key_unit_nr(logger, nmr_section, "PRINT%CHI_TENSOR", & 798 extension=".data", middle_name="CHI", log_filename=.FALSE.) 799 800 WRITE (title, '(A)') "Magnetic Susceptibility Tensor " 801 IF (unit_atoms > 0) THEN 802 WRITE (unit_atoms, '(T2,A)') title 803 WRITE (unit_atoms, '(T1,A)') " CHI from SOFT J in 10^-30 J/T^2 units" 804 WRITE (unit_atoms, '(3(A,f10.4))') ' XX = ', chi_tensor_tmp(1, 1),& 805 & ' XY = ', chi_tensor_tmp(1, 2),& 806 & ' XZ = ', chi_tensor_tmp(1, 3) 807 WRITE (unit_atoms, '(3(A,f10.4))') ' YX = ', chi_tensor_tmp(2, 1),& 808 & ' YY = ', chi_tensor_tmp(2, 2),& 809 & ' YZ = ', chi_tensor_tmp(2, 3) 810 WRITE (unit_atoms, '(3(A,f10.4))') ' ZX = ', chi_tensor_tmp(3, 1),& 811 & ' ZY = ', chi_tensor_tmp(3, 2),& 812 & ' ZZ = ', chi_tensor_tmp(3, 3) 813 IF (gapw) THEN 814 WRITE (unit_atoms, '(T1,A)') " CHI from LOCAL J in 10^-30 J/T^2 units" 815 WRITE (unit_atoms, '(3(A,f10.4))') ' XX = ', chi_tensor_loc_tmp(1, 1),& 816 & ' XY = ', chi_tensor_loc_tmp(1, 2),& 817 & ' XZ = ', chi_tensor_loc_tmp(1, 3) 818 WRITE (unit_atoms, '(3(A,f10.4))') ' YX = ', chi_tensor_loc_tmp(2, 1),& 819 & ' YY = ', chi_tensor_loc_tmp(2, 2),& 820 & ' YZ = ', chi_tensor_loc_tmp(2, 3) 821 WRITE (unit_atoms, '(3(A,f10.4))') ' ZX = ', chi_tensor_loc_tmp(3, 1),& 822 & ' ZY = ', chi_tensor_loc_tmp(3, 2),& 823 & ' ZZ = ', chi_tensor_loc_tmp(3, 3) 824 ENDIF 825 WRITE (unit_atoms, '(T1,A)') " Total CHI in 10^-30 J/T^2 units" 826 WRITE (unit_atoms, '(3(A,f10.4))') ' XX = ', chi_sym_tot(1, 1),& 827 & ' XY = ', chi_sym_tot(1, 2),& 828 & ' XZ = ', chi_sym_tot(1, 3) 829 WRITE (unit_atoms, '(3(A,f10.4))') ' YX = ', chi_sym_tot(2, 1),& 830 & ' YY = ', chi_sym_tot(2, 2),& 831 & ' YZ = ', chi_sym_tot(2, 3) 832 WRITE (unit_atoms, '(3(A,f10.4))') ' ZX = ', chi_sym_tot(3, 1),& 833 & ' ZY = ', chi_sym_tot(3, 2),& 834 & ' ZZ = ', chi_sym_tot(3, 3) 835 chi_sym_tot(:, :) = chi_sym_tot(:, :)*nmr_env%chi_SI2ppmcgs 836 WRITE (unit_atoms, '(T1,A)') " Total CHI in ppm-cgs units" 837 WRITE (unit_atoms, '(3(A,f10.4))') ' XX = ', chi_sym_tot(1, 1),& 838 & ' XY = ', chi_sym_tot(1, 2),& 839 & ' XZ = ', chi_sym_tot(1, 3) 840 WRITE (unit_atoms, '(3(A,f10.4))') ' YX = ', chi_sym_tot(2, 1),& 841 & ' YY = ', chi_sym_tot(2, 2),& 842 & ' YZ = ', chi_sym_tot(2, 3) 843 WRITE (unit_atoms, '(3(A,f10.4))') ' ZX = ', chi_sym_tot(3, 1),& 844 & ' ZY = ', chi_sym_tot(3, 2),& 845 & ' ZZ = ', chi_sym_tot(3, 3) 846 WRITE (unit_atoms, '(/T1,3(A,f10.4))') & 847 ' PV1=', nmr_env%chi_SI2ppmcgs*eig(1), & 848 ' PV2=', nmr_env%chi_SI2ppmcgs*eig(2), & 849 ' PV3=', nmr_env%chi_SI2ppmcgs*eig(3) 850 WRITE (unit_atoms, '(T1,A,F10.4,10X,A,F10.4)') & 851 ' ISO=', nmr_env%chi_SI2ppmcgs*chi_iso, & 852 'ANISO=', nmr_env%chi_SI2ppmcgs*chi_aniso 853 ENDIF 854 855 CALL cp_print_key_finished_output(unit_atoms, logger, nmr_section,& 856 & "PRINT%CHI_TENSOR") 857 ENDIF ! print chi 858 ! 859 ! Add the chi part to the shifts 860 cs_tot(:, :, :) = 0.0_dp 861 DO ir = 1, nat_print 862 iatom = cs_atom_list(ir) 863 rpos(1:3) = particle_set(iatom)%r(1:3) 864 atom_kind => particle_set(iatom)%atomic_kind 865 CALL get_atomic_kind(atom_kind, name=name, element_symbol=element_symbol) 866 cs_tot(:, :, ir) = chi_tensor_tmp(:, :)*nmr_env%chi_SI2shiftppm + cs(:, :, iatom) 867 IF (gapw) cs_tot(:, :, ir) = cs_tot(:, :, ir) + cs_loc(:, :, iatom) 868 ENDDO ! ir 869 IF (output_unit > 0) THEN 870 WRITE (output_unit, '(T2,A,E14.6)') 'CheckSum Shifts =', & 871 SQRT(DDOT(9*SIZE(cs_tot, 3), cs_tot(1, 1, 1), 1, cs_tot(1, 1, 1), 1)) 872 ENDIF 873 ! 874 ! print shifts 875 IF (BTEST(cp_print_key_should_output(logger%iter_info, nmr_section, & 876 "PRINT%SHIELDING_TENSOR"), cp_p_file)) THEN 877 878 unit_atoms = cp_print_key_unit_nr(logger, nmr_section, "PRINT%SHIELDING_TENSOR", & 879 extension=".data", middle_name="SHIFT", & 880 log_filename=.FALSE.) 881 882 nat_print = SIZE(cs_atom_list, 1) 883 IF (unit_atoms > 0) THEN 884 WRITE (title, '(A,1X,I5)') "Shielding atom at atomic positions. # tensors printed ", nat_print 885 WRITE (unit_atoms, '(T2,A)') title 886 DO ir = 1, nat_print 887 iatom = cs_atom_list(ir) 888 rpos(1:3) = particle_set(iatom)%r(1:3) 889 atom_kind => particle_set(iatom)%atomic_kind 890 CALL get_atomic_kind(atom_kind, name=name, element_symbol=element_symbol) 891 shift_sym_tot(:, :) = 0.5_dp*(cs_tot(:, :, ir) + TRANSPOSE(cs_tot(:, :, ir))) 892 CALL diamat_all(shift_sym_tot, eig) 893 shift_iso = (eig(1) + eig(2) + eig(3))/3.0_dp 894 shift_aniso = eig(3) - (eig(2) + eig(1))/2.0_dp 895 ! 896 WRITE (unit_atoms, '(T2,I5,A,2X,A2,2X,3f15.6)') iatom, TRIM(name), element_symbol, rpos(1:3) 897 ! 898 IF (gapw) THEN 899 WRITE (unit_atoms, '(T1,A)') " SIGMA from SOFT J" 900 WRITE (unit_atoms, '(3(A,f10.4))') ' XX = ', cs(1, 1, iatom),& 901 & ' XY = ', cs(1, 2, iatom),& 902 & ' XZ = ', cs(1, 3, iatom) 903 WRITE (unit_atoms, '(3(A,f10.4))') ' YX = ', cs(2, 1, iatom),& 904 & ' YY = ', cs(2, 2, iatom),& 905 & ' YZ = ', cs(2, 3, iatom) 906 WRITE (unit_atoms, '(3(A,f10.4))') ' ZX = ', cs(3, 1, iatom),& 907 & ' ZY = ', cs(3, 2, iatom),& 908 & ' ZZ = ', cs(3, 3, iatom) 909 WRITE (unit_atoms, '(T1,A)') " SIGMA from LOCAL J" 910 WRITE (unit_atoms, '(3(A,f10.4))') ' XX = ', cs_loc(1, 1, iatom),& 911 & ' XY = ', cs_loc(1, 2, iatom),& 912 & ' XZ = ', cs_loc(1, 3, iatom) 913 WRITE (unit_atoms, '(3(A,f10.4))') ' YX = ', cs_loc(2, 1, iatom),& 914 & ' YY = ', cs_loc(2, 2, iatom),& 915 & ' YZ = ', cs_loc(2, 3, iatom) 916 WRITE (unit_atoms, '(3(A,f10.4))') ' ZX = ', cs_loc(3, 1, iatom),& 917 & ' ZY = ', cs_loc(3, 2, iatom),& 918 & ' ZZ = ', cs_loc(3, 3, iatom) 919 ENDIF 920 WRITE (unit_atoms, '(T1,A)') " SIGMA TOTAL" 921 WRITE (unit_atoms, '(3(A,f10.4))') ' XX = ', cs_tot(1, 1, ir),& 922 & ' XY = ', cs_tot(1, 2, ir),& 923 & ' XZ = ', cs_tot(1, 3, ir) 924 WRITE (unit_atoms, '(3(A,f10.4))') ' YX = ', cs_tot(2, 1, ir),& 925 & ' YY = ', cs_tot(2, 2, ir),& 926 & ' YZ = ', cs_tot(2, 3, ir) 927 WRITE (unit_atoms, '(3(A,f10.4))') ' ZX = ', cs_tot(3, 1, ir),& 928 & ' ZY = ', cs_tot(3, 2, ir),& 929 & ' ZZ = ', cs_tot(3, 3, ir) 930 WRITE (unit_atoms, '(T1,2(A,f12.4))') ' ISOTROPY = ', shift_iso,& 931 & ' ANISOTROPY = ', shift_aniso 932 ENDDO ! ir 933 ENDIF 934 CALL cp_print_key_finished_output(unit_atoms, logger, nmr_section,& 935 & "PRINT%SHIELDING_TENSOR") 936 937 IF (do_nics) THEN 938 ! 939 ! Add the chi part to the nics 940 cs_nics_tot(:, :, :) = 0.0_dp 941 DO ir = 1, n_nics 942 cs_nics_tot(:, :, ir) = chi_tensor_tmp(:, :)*nmr_env%chi_SI2shiftppm + cs_nics(:, :, ir) 943 IF (gapw) cs_nics_tot(:, :, ir) = cs_nics_tot(:, :, ir) + cs_nics_loc(:, :, ir) 944 ENDDO ! ir 945 IF (output_unit > 0) THEN 946 WRITE (output_unit, '(T2,A,E14.6)') 'CheckSum NICS =', & 947 SQRT(DDOT(9*SIZE(cs_nics_tot, 3), cs_nics_tot(1, 1, 1), 1, cs_nics_tot(1, 1, 1), 1)) 948 ENDIF 949 ! 950 unit_nics = cp_print_key_unit_nr(logger, nmr_section, "PRINT%SHIELDING_TENSOR", & 951 extension=".data", middle_name="NICS", & 952 log_filename=.FALSE.) 953 IF (unit_nics > 0) THEN 954 WRITE (title, '(A,1X,I5)') "Shielding at nics positions. # tensors printed ", n_nics 955 WRITE (unit_nics, '(T2,A)') title 956 DO ir = 1, n_nics 957 shift_sym_tot(:, :) = 0.5_dp*(cs_nics_tot(:, :, ir) + TRANSPOSE(cs_nics_tot(:, :, ir))) 958 CALL diamat_all(shift_sym_tot, eig) 959 shift_iso = (eig(1) + eig(2) + eig(3))/3.0_dp 960 shift_aniso = eig(3) - (eig(2) + eig(1))/2.0_dp 961 ! 962 WRITE (unit_nics, '(T2,I5,2X,3f15.6)') ir, r_nics(1:3, ir) 963 ! 964 IF (gapw) THEN 965 WRITE (unit_nics, '(T1,A)') " SIGMA from SOFT J" 966 WRITE (unit_nics, '(3(A,f10.4))') ' XX = ', cs_nics(1, 1, ir),& 967 & ' XY = ', cs_nics(1, 2, ir),& 968 & ' XZ = ', cs_nics(1, 3, ir) 969 WRITE (unit_nics, '(3(A,f10.4))') ' YX = ', cs_nics(2, 1, ir),& 970 & ' YY = ', cs_nics(2, 2, ir),& 971 & ' YZ = ', cs_nics(2, 3, ir) 972 WRITE (unit_nics, '(3(A,f10.4))') ' ZX = ', cs_nics(3, 1, ir),& 973 & ' ZY = ', cs_nics(3, 2, ir),& 974 & ' ZZ = ', cs_nics(3, 3, ir) 975 WRITE (unit_nics, '(T1,A)') " SIGMA from LOCAL J" 976 WRITE (unit_nics, '(3(A,f10.4))') ' XX = ', cs_nics_loc(1, 1, ir),& 977 & ' XY = ', cs_nics_loc(1, 2, ir),& 978 & ' XZ = ', cs_nics_loc(1, 3, ir) 979 WRITE (unit_nics, '(3(A,f10.4))') ' YX = ', cs_nics_loc(2, 1, ir),& 980 & ' YY = ', cs_nics_loc(2, 2, ir),& 981 & ' YZ = ', cs_nics_loc(2, 3, ir) 982 WRITE (unit_nics, '(3(A,f10.4))') ' ZX = ', cs_nics_loc(3, 1, ir),& 983 & ' ZY = ', cs_nics_loc(3, 2, ir),& 984 & ' ZZ = ', cs_nics_loc(3, 3, ir) 985 ENDIF 986 WRITE (unit_nics, '(T1,A)') " SIGMA TOTAL" 987 WRITE (unit_nics, '(3(A,f10.4))') ' XX = ', cs_nics_tot(1, 1, ir),& 988 & ' XY = ', cs_nics_tot(1, 2, ir),& 989 & ' XZ = ', cs_nics_tot(1, 3, ir) 990 WRITE (unit_nics, '(3(A,f10.4))') ' YX = ', cs_nics_tot(2, 1, ir),& 991 & ' YY = ', cs_nics_tot(2, 2, ir),& 992 & ' YZ = ', cs_nics_tot(2, 3, ir) 993 WRITE (unit_nics, '(3(A,f10.4))') ' ZX = ', cs_nics_tot(3, 1, ir),& 994 & ' ZY = ', cs_nics_tot(3, 2, ir),& 995 & ' ZZ = ', cs_nics_tot(3, 3, ir) 996 WRITE (unit_nics, '(T1,2(A,f12.4))') ' ISOTROPY = ', shift_iso,& 997 & ' ANISOTROPY = ', shift_aniso 998 ENDDO 999 ENDIF 1000 CALL cp_print_key_finished_output(unit_nics, logger, nmr_section,& 1001 & "PRINT%SHIELDING_TENSOR") 1002 ENDIF 1003 ENDIF ! print shift 1004 ! 1005 ! clean up 1006 DEALLOCATE (cs_tot) 1007 IF (do_nics) THEN 1008 DEALLOCATE (cs_nics_tot) 1009 ENDIF 1010 ! 1011!100 FORMAT(A,1X,I5) 1012!101 FORMAT(T2,A) 1013!102 FORMAT(T2,I5,A,2X,A2,2X,3f15.6) 1014!103 FORMAT(T2,I5,2X,3f15.6) 1015!104 FORMAT(T1,A) 1016!105 FORMAT(3(A,f10.4)) 1017!106 FORMAT(T1,2(A,f12.4)) 1018 END SUBROUTINE nmr_shift_print 1019 1020END MODULE qs_linres_nmr_shift 1021 1022