1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief given the response wavefunctions obtained by the application 8!> of the (rxp), p, and ((dk-dl)xp) operators, 9!> here the current density vector (jx, jy, jz) 10!> is computed for the 3 directions of the magnetic field (Bx, By, Bz) 11!> \par History 12!> created 02-2006 [MI] 13!> \author MI 14! ************************************************************************************************** 15MODULE qs_linres_atom_current 16 USE atomic_kind_types, ONLY: atomic_kind_type,& 17 get_atomic_kind 18 USE basis_set_types, ONLY: get_gto_basis_set,& 19 gto_basis_set_p_type,& 20 gto_basis_set_type 21 USE cell_types, ONLY: cell_type,& 22 pbc 23 USE cp_control_types, ONLY: dft_control_type 24 USE cp_log_handling, ONLY: cp_logger_get_default_io_unit 25 USE cp_para_types, ONLY: cp_para_env_type 26 USE dbcsr_api, ONLY: dbcsr_get_block_p,& 27 dbcsr_p_type 28 USE input_constants, ONLY: current_gauge_atom,& 29 current_gauge_r,& 30 current_gauge_r_and_step_func 31 USE kinds, ONLY: dp 32 USE message_passing, ONLY: mp_sum 33 USE orbital_pointers, ONLY: indso,& 34 nsoset 35 USE particle_types, ONLY: particle_type 36 USE paw_proj_set_types, ONLY: get_paw_proj_set,& 37 paw_proj_set_type 38 USE qs_environment_types, ONLY: get_qs_env,& 39 qs_environment_type 40 USE qs_grid_atom, ONLY: grid_atom_type 41 USE qs_harmonics_atom, ONLY: get_none0_cg_list,& 42 harmonics_atom_type 43 USE qs_kind_types, ONLY: get_qs_kind,& 44 get_qs_kind_set,& 45 qs_kind_type 46 USE qs_linres_op, ONLY: fac_vecp,& 47 set_vecp,& 48 set_vecp_rev 49 USE qs_linres_types, ONLY: allocate_jrho_atom_rad,& 50 allocate_jrho_coeff,& 51 current_env_type,& 52 get_current_env,& 53 jrho_atom_type,& 54 set2zero_jrho_atom_rad 55 USE qs_neighbor_list_types, ONLY: get_iterator_info,& 56 neighbor_list_iterate,& 57 neighbor_list_iterator_create,& 58 neighbor_list_iterator_p_type,& 59 neighbor_list_iterator_release,& 60 neighbor_list_set_p_type 61 USE qs_oce_methods, ONLY: proj_blk 62 USE qs_oce_types, ONLY: oce_matrix_type 63 USE qs_rho_atom_types, ONLY: rho_atom_coeff 64 USE sap_kind_types, ONLY: alist_pre_align_blk,& 65 alist_type,& 66 get_alist 67 USE util, ONLY: get_limit 68 69!$ USE OMP_LIB, ONLY: omp_get_max_threads, & 70!$ omp_get_thread_num, & 71!$ omp_lock_kind, & 72!$ omp_init_lock, omp_set_lock, & 73!$ omp_unset_lock, omp_destroy_lock 74 75#include "./base/base_uses.f90" 76 77 IMPLICIT NONE 78 79 PRIVATE 80 81 ! *** Public subroutines *** 82 PUBLIC :: calculate_jrho_atom_rad, calculate_jrho_atom, calculate_jrho_atom_coeff 83 84 ! *** Global parameters (only in this module) 85 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_atom_current' 86 87CONTAINS 88 89! ************************************************************************************************** 90!> \brief Calculate the expansion coefficients for the atomic terms 91!> of the current densitiy in GAPW 92!> \param qs_env ... 93!> \param current_env ... 94!> \param mat_d0 ... 95!> \param mat_jp ... 96!> \param mat_jp_rii ... 97!> \param mat_jp_riii ... 98!> \param iB ... 99!> \param idir ... 100!> \par History 101!> 07.2006 created [MI] 102!> 02.2009 using new setup of projector-basis overlap [jgh] 103!> 08.2016 add OpenMP [EPCC] 104!> 09.2016 use automatic arrays [M Tucker] 105! ************************************************************************************************** 106 SUBROUTINE calculate_jrho_atom_coeff(qs_env, current_env, mat_d0, mat_jp, mat_jp_rii, & 107 mat_jp_riii, iB, idir) 108 109 TYPE(qs_environment_type), POINTER :: qs_env 110 TYPE(current_env_type) :: current_env 111 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_d0, mat_jp, mat_jp_rii, mat_jp_riii 112 INTEGER, INTENT(IN) :: iB, idir 113 114 CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_jrho_atom_coeff', & 115 routineP = moduleN//':'//routineN 116 117 INTEGER :: bo(2), handle, iac, iat, iatom, ibc, idir2, ii, iii, ikind, ispin, istat, jatom, & 118 jkind, kac, katom, kbc, kkind, len_CPC, len_PC1, max_gau, max_nsgf, mepos, n_cont_a, & 119 n_cont_b, nat, natom, nkind, nsatbas, nsgfa, nsgfb, nso, nsoctot, nspins, num_pe, & 120 output_unit 121 INTEGER, DIMENSION(3) :: cell_b 122 INTEGER, DIMENSION(:), POINTER :: atom_list, list_a, list_b 123 LOGICAL :: den_found, dista, distab, distb, & 124 is_not_associated, paw_atom, & 125 sgf_soft_only_a, sgf_soft_only_b 126 REAL(dp) :: eps_cpc, hard_radius_c, jmax, nbr_dbl, & 127 rab(3), rbc(3) 128 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: a_matrix, b_matrix, c_matrix, d_matrix 129 REAL(KIND=dp), DIMENSION(:, :), POINTER :: C_coeff_hh_a, C_coeff_hh_b, & 130 C_coeff_ss_a, C_coeff_ss_b, r_coef_h, & 131 r_coef_s, tmp_coeff, zero_coeff 132 TYPE(alist_type), POINTER :: alist_ac, alist_bc 133 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 134 TYPE(cp_para_env_type), POINTER :: para_env 135 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_a, mat_b, mat_c, mat_d 136 TYPE(dft_control_type), POINTER :: dft_control 137 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list 138 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b, orb_basis_set 139 TYPE(jrho_atom_type), DIMENSION(:), POINTER :: jrho1_atom_set 140 TYPE(neighbor_list_iterator_p_type), & 141 DIMENSION(:), POINTER :: nl_iterator 142 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 143 POINTER :: sab_all 144 TYPE(oce_matrix_type), POINTER :: oce 145 TYPE(paw_proj_set_type), POINTER :: paw_proj 146 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 147 TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: a_block, b_block, c_block, d_block, & 148 jp2_RARnu, jp_RARnu 149 150!$ INTEGER(kind=omp_lock_kind), ALLOCATABLE, DIMENSION(:) :: proj_blk_lock, alloc_lock 151!$ INTEGER :: lock 152 153 CALL timeset(routineN, handle) 154 155 NULLIFY (atomic_kind_set, qs_kind_set, dft_control, sab_all, jrho1_atom_set, oce, & 156 para_env, zero_coeff, tmp_coeff) 157 158 CALL get_qs_env(qs_env=qs_env, & 159 atomic_kind_set=atomic_kind_set, & 160 qs_kind_set=qs_kind_set, & 161 dft_control=dft_control, & 162 oce=oce, & 163 sab_all=sab_all, & 164 para_env=para_env) 165 CPASSERT(ASSOCIATED(oce)) 166 167 CALL get_current_env(current_env=current_env, jrho1_atom_set=jrho1_atom_set) 168 CPASSERT(ASSOCIATED(jrho1_atom_set)) 169 170 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, & 171 maxsgf=max_nsgf, & 172 maxgtops=max_gau) 173 174 eps_cpc = dft_control%qs_control%gapw_control%eps_cpc 175 176 idir2 = 1 177 IF (idir .NE. iB) THEN 178 CALL set_vecp_rev(idir, iB, idir2) 179 ENDIF 180 CALL set_vecp(iB, ii, iii) 181 182 ! Set pointers for the different gauge 183 mat_a => mat_d0 184 mat_b => mat_jp 185 mat_c => mat_jp_rii 186 mat_d => mat_jp_riii 187 188 ! Density-like matrices 189 nkind = SIZE(qs_kind_set) 190 natom = SIZE(jrho1_atom_set) 191 nspins = dft_control%nspins 192 193 ! Reset CJC coefficients and local density arrays 194 DO ikind = 1, nkind 195 NULLIFY (atom_list) 196 CALL get_atomic_kind(atomic_kind_set(ikind), & 197 atom_list=atom_list, & 198 natom=nat) 199 CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom) 200 201 ! Quick cycle if needed. 202 IF (.NOT. paw_atom) CYCLE 203 204 ! Initialize the density matrix-like arrays. 205 DO iat = 1, nat 206 iatom = atom_list(iat) 207 DO ispin = 1, nspins 208 IF (ASSOCIATED(jrho1_atom_set(iatom)%cjc0_h(1)%r_coef)) THEN 209 jrho1_atom_set(iatom)%cjc0_h(ispin)%r_coef = 0.0_dp 210 jrho1_atom_set(iatom)%cjc0_s(ispin)%r_coef = 0.0_dp 211 jrho1_atom_set(iatom)%cjc_h(ispin)%r_coef = 0.0_dp 212 jrho1_atom_set(iatom)%cjc_s(ispin)%r_coef = 0.0_dp 213 jrho1_atom_set(iatom)%cjc_ii_h(ispin)%r_coef = 0.0_dp 214 jrho1_atom_set(iatom)%cjc_ii_s(ispin)%r_coef = 0.0_dp 215 jrho1_atom_set(iatom)%cjc_iii_h(ispin)%r_coef = 0.0_dp 216 jrho1_atom_set(iatom)%cjc_iii_s(ispin)%r_coef = 0.0_dp 217 ENDIF 218 ENDDO ! ispin 219 ENDDO ! iat 220 ENDDO ! ikind 221 222 ! Three centers 223 ALLOCATE (basis_set_list(nkind)) 224 DO ikind = 1, nkind 225 CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set_a) 226 IF (ASSOCIATED(basis_set_a)) THEN 227 basis_set_list(ikind)%gto_basis_set => basis_set_a 228 ELSE 229 NULLIFY (basis_set_list(ikind)%gto_basis_set) 230 END IF 231 END DO 232 233 len_PC1 = max_nsgf*max_gau 234 len_CPC = max_gau*max_gau 235 236 num_pe = 1 237!$ num_pe = omp_get_max_threads() 238 CALL neighbor_list_iterator_create(nl_iterator, sab_all, nthread=num_pe) 239 240!$OMP PARALLEL DEFAULT( NONE ) & 241!$OMP SHARED( nspins, max_nsgf, max_gau & 242!$OMP , len_PC1, len_CPC & 243!$OMP , nl_iterator, basis_set_list & 244!$OMP , mat_a, mat_b, mat_c, mat_d & 245!$OMP , nkind, qs_kind_set, eps_cpc, oce & 246!$OMP , ii, iii, jrho1_atom_set & 247!$OMP , natom, proj_blk_lock, alloc_lock & 248!$OMP ) & 249!$OMP PRIVATE( a_block, b_block, c_block, d_block & 250!$OMP , jp_RARnu, jp2_RARnu & 251!$OMP , a_matrix, b_matrix, c_matrix, d_matrix, istat & 252!$OMP , mepos & 253!$OMP , ikind, jkind, kkind, iatom, jatom, katom & 254!$OMP , cell_b, rab, rbc & 255!$OMP , basis_set_a, nsgfa & 256!$OMP , basis_set_b, nsgfb & 257!$OMP , orb_basis_set, jmax, den_found & 258!$OMP , hard_radius_c, paw_proj & 259!$OMP , nsatbas, nsoctot, nso, paw_atom & 260!$OMP , iac , alist_ac, kac, n_cont_a, list_a, sgf_soft_only_a & 261!$OMP , ibc , alist_bc, kbc, n_cont_b, list_b, sgf_soft_only_b & 262!$OMP , C_coeff_hh_a, C_coeff_ss_a, dista & 263!$OMP , C_coeff_hh_b, C_coeff_ss_b, distb & 264!$OMP , distab & 265!$OMP , r_coef_s, r_coef_h & 266!$OMP ) 267 268 NULLIFY (a_block, b_block, c_block, d_block) 269 NULLIFY (orb_basis_set, jp_RARnu, jp2_RARnu) 270 271 ALLOCATE (a_block(nspins), b_block(nspins), c_block(nspins), d_block(nspins), & 272 jp_RARnu(nspins), jp2_RARnu(nspins), & 273 STAT=istat) 274 CPASSERT(istat == 0) 275 276 ALLOCATE (a_matrix(max_nsgf, max_nsgf), b_matrix(max_nsgf, max_nsgf), & 277 c_matrix(max_nsgf, max_nsgf), d_matrix(max_nsgf, max_nsgf), & 278 STAT=istat) 279 CPASSERT(istat == 0) 280 281!$OMP SINGLE 282!$ ALLOCATE (alloc_lock(natom)) 283!$ ALLOCATE (proj_blk_lock(nspins*natom)) 284!$OMP END SINGLE 285 286!$OMP DO 287!$ DO lock = 1, natom 288!$ call omp_init_lock(alloc_lock(lock)) 289!$ END DO 290!$OMP END DO 291 292!$OMP DO 293!$ DO lock = 1, nspins*natom 294!$ call omp_init_lock(proj_blk_lock(lock)) 295!$ END DO 296!$OMP END DO 297 298 mepos = 0 299!$ mepos = omp_get_thread_num() 300 301 DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0) 302 303 CALL get_iterator_info(nl_iterator, mepos=mepos, & 304 ikind=ikind, jkind=jkind, & 305 iatom=iatom, jatom=jatom, cell=cell_b, r=rab) 306 basis_set_a => basis_set_list(ikind)%gto_basis_set 307 IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE 308 basis_set_b => basis_set_list(jkind)%gto_basis_set 309 IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE 310 nsgfa = basis_set_a%nsgf 311 nsgfb = basis_set_b%nsgf 312 DO ispin = 1, nspins 313 NULLIFY (jp_RARnu(ispin)%r_coef, jp2_RARnu(ispin)%r_coef) 314 ALLOCATE (jp_RARnu(ispin)%r_coef(nsgfa, nsgfb), & 315 jp2_RARnu(ispin)%r_coef(nsgfa, nsgfb)) 316 ENDDO 317 318 ! Take the block \mu\nu of jpab, jpab_ii and jpab_iii 319 jmax = 0._dp 320 DO ispin = 1, nspins 321 NULLIFY (a_block(ispin)%r_coef) 322 NULLIFY (b_block(ispin)%r_coef) 323 NULLIFY (c_block(ispin)%r_coef) 324 NULLIFY (d_block(ispin)%r_coef) 325 CALL dbcsr_get_block_p(matrix=mat_a(ispin)%matrix, & 326 row=iatom, col=jatom, block=a_block(ispin)%r_coef, & 327 found=den_found) 328 jmax = jmax + MAXVAL(ABS(a_block(ispin)%r_coef)) 329 CALL dbcsr_get_block_p(matrix=mat_b(ispin)%matrix, & 330 row=iatom, col=jatom, block=b_block(ispin)%r_coef, & 331 found=den_found) 332 jmax = jmax + MAXVAL(ABS(b_block(ispin)%r_coef)) 333 CALL dbcsr_get_block_p(matrix=mat_c(ispin)%matrix, & 334 row=iatom, col=jatom, block=c_block(ispin)%r_coef, & 335 found=den_found) 336 jmax = jmax + MAXVAL(ABS(c_block(ispin)%r_coef)) 337 CALL dbcsr_get_block_p(matrix=mat_d(ispin)%matrix, & 338 row=iatom, col=jatom, block=d_block(ispin)%r_coef, & 339 found=den_found) 340 jmax = jmax + MAXVAL(ABS(d_block(ispin)%r_coef)) 341 ENDDO 342 343 ! Loop over atoms 344 DO kkind = 1, nkind 345 NULLIFY (paw_proj) 346 CALL get_qs_kind(qs_kind_set(kkind), & 347 basis_set=orb_basis_set, & 348 hard_radius=hard_radius_c, & 349 paw_proj_set=paw_proj, & 350 paw_atom=paw_atom) 351 352 ! Quick cycle if needed. 353 IF (.NOT. paw_atom) CYCLE 354 355 CALL get_paw_proj_set(paw_proj_set=paw_proj, nsatbas=nsatbas) 356 nsoctot = nsatbas 357 358 iac = ikind + nkind*(kkind - 1) 359 ibc = jkind + nkind*(kkind - 1) 360 IF (.NOT. ASSOCIATED(oce%intac(iac)%alist)) CYCLE 361 IF (.NOT. ASSOCIATED(oce%intac(ibc)%alist)) CYCLE 362 363 CALL get_alist(oce%intac(iac), alist_ac, iatom) 364 CALL get_alist(oce%intac(ibc), alist_bc, jatom) 365 IF (.NOT. ASSOCIATED(alist_ac)) CYCLE 366 IF (.NOT. ASSOCIATED(alist_bc)) CYCLE 367 368 DO kac = 1, alist_ac%nclist 369 DO kbc = 1, alist_bc%nclist 370 IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) CYCLE 371 372 IF (ALL(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN 373 IF (jmax*alist_bc%clist(kbc)%maxac*alist_ac%clist(kac)%maxac < eps_cpc) CYCLE 374 375 n_cont_a = alist_ac%clist(kac)%nsgf_cnt 376 n_cont_b = alist_bc%clist(kbc)%nsgf_cnt 377 sgf_soft_only_a = alist_ac%clist(kac)%sgf_soft_only 378 sgf_soft_only_b = alist_bc%clist(kbc)%sgf_soft_only 379 IF (n_cont_a .EQ. 0 .OR. n_cont_b .EQ. 0) CYCLE 380 381 ! thanks to the linearity of the response, we 382 ! can avoid computing soft-soft interations. 383 ! those terms are already included in the 384 ! regular grid. 385 IF (sgf_soft_only_a .AND. sgf_soft_only_b) CYCLE 386 387 list_a => alist_ac%clist(kac)%sgf_list 388 list_b => alist_bc%clist(kbc)%sgf_list 389 390 katom = alist_ac%clist(kac)%catom 391!$ CALL omp_set_lock(alloc_lock(katom)) 392 IF (.NOT. ASSOCIATED(jrho1_atom_set(katom)%cjc0_h(1)%r_coef)) THEN 393 CALL allocate_jrho_coeff(jrho1_atom_set, katom, nsoctot) 394 ENDIF 395!$ CALL omp_unset_lock(alloc_lock(katom)) 396 397 ! Compute the modified Qai matrix as 398 ! mQai_\mu\nu = Qai_\mu\nu - Qbi_\mu\nu * (R_A-R_\nu)_ii 399 ! + Qci_\mu\nu * (R_A-R_\nu)_iii 400 rbc = alist_bc%clist(kbc)%rac 401 DO ispin = 1, nspins 402 CALL DCOPY(nsgfa*nsgfb, b_block(ispin)%r_coef(1, 1), 1, & 403 jp_RARnu(ispin)%r_coef(1, 1), 1) 404 CALL DAXPY(nsgfa*nsgfb, -rbc(ii), d_block(ispin)%r_coef(1, 1), 1, & 405 jp_RARnu(ispin)%r_coef(1, 1), 1) 406 CALL DAXPY(nsgfa*nsgfb, rbc(iii), c_block(ispin)%r_coef(1, 1), 1, & 407 jp_RARnu(ispin)%r_coef(1, 1), 1) 408 ENDDO 409 410 ! Get the d_A's for the hard and soft densities. 411 IF (iatom == katom .AND. ALL(alist_ac%clist(kac)%cell == 0)) THEN 412 C_coeff_hh_a => alist_ac%clist(kac)%achint(:, :, 1) 413 C_coeff_ss_a => alist_ac%clist(kac)%acint(:, :, 1) 414 dista = .FALSE. 415 ELSE 416 C_coeff_hh_a => alist_ac%clist(kac)%acint(:, :, 1) 417 C_coeff_ss_a => alist_ac%clist(kac)%acint(:, :, 1) 418 dista = .TRUE. 419 END IF 420 ! Get the d_B's for the hard and soft densities. 421 IF (jatom == katom .AND. ALL(alist_bc%clist(kbc)%cell == 0)) THEN 422 C_coeff_hh_b => alist_bc%clist(kbc)%achint(:, :, 1) 423 C_coeff_ss_b => alist_bc%clist(kbc)%acint(:, :, 1) 424 distb = .FALSE. 425 ELSE 426 C_coeff_hh_b => alist_bc%clist(kbc)%acint(:, :, 1) 427 C_coeff_ss_b => alist_bc%clist(kbc)%acint(:, :, 1) 428 distb = .TRUE. 429 END IF 430 431 distab = dista .AND. distb 432 433 nso = nsoctot 434 435 DO ispin = 1, nspins 436 437 ! align the blocks 438 CALL alist_pre_align_blk(a_block(ispin)%r_coef, SIZE(a_block(ispin)%r_coef, 1), & 439 a_matrix, SIZE(a_matrix, 1), & 440 list_a, n_cont_a, list_b, n_cont_b) 441 442 CALL alist_pre_align_blk(jp_RARnu(ispin)%r_coef, SIZE(jp_RARnu(ispin)%r_coef, 1), & 443 b_matrix, SIZE(b_matrix, 1), & 444 list_a, n_cont_a, list_b, n_cont_b) 445 446 CALL alist_pre_align_blk(c_block(ispin)%r_coef, SIZE(c_block(ispin)%r_coef, 1), & 447 c_matrix, SIZE(c_matrix, 1), & 448 list_a, n_cont_a, list_b, n_cont_b) 449 CALL alist_pre_align_blk(d_block(ispin)%r_coef, SIZE(d_block(ispin)%r_coef, 1), & 450 d_matrix, SIZE(d_matrix, 1), & 451 list_a, n_cont_a, list_b, n_cont_b) 452 453!$ CALL omp_set_lock(proj_blk_lock((katom - 1)*nspins + ispin)) 454 !------------------------------------------------------------------ 455 ! P_\alpha\alpha' 456 r_coef_h => jrho1_atom_set(katom)%cjc0_h(ispin)%r_coef 457 r_coef_s => jrho1_atom_set(katom)%cjc0_s(ispin)%r_coef 458 CALL proj_blk(C_coeff_hh_a, C_coeff_ss_a, n_cont_a, & 459 C_coeff_hh_b, C_coeff_ss_b, n_cont_b, & 460 a_matrix, max_nsgf, r_coef_h, r_coef_s, nso, & 461 len_PC1, len_CPC, 1.0_dp, distab) 462 !------------------------------------------------------------------ 463 ! mQai_\alpha\alpha' 464 r_coef_h => jrho1_atom_set(katom)%cjc_h(ispin)%r_coef 465 r_coef_s => jrho1_atom_set(katom)%cjc_s(ispin)%r_coef 466 CALL proj_blk(C_coeff_hh_a, C_coeff_ss_a, n_cont_a, & 467 C_coeff_hh_b, C_coeff_ss_b, n_cont_b, & 468 b_matrix, max_nsgf, r_coef_h, r_coef_s, nso, & 469 len_PC1, len_CPC, 1.0_dp, distab) 470 !------------------------------------------------------------------ 471 ! Qci_\alpha\alpha' 472 r_coef_h => jrho1_atom_set(katom)%cjc_ii_h(ispin)%r_coef 473 r_coef_s => jrho1_atom_set(katom)%cjc_ii_s(ispin)%r_coef 474 CALL proj_blk(C_coeff_hh_a, C_coeff_ss_a, n_cont_a, & 475 C_coeff_hh_b, C_coeff_ss_b, n_cont_b, & 476 c_matrix, max_nsgf, r_coef_h, r_coef_s, nso, & 477 len_PC1, len_CPC, 1.0_dp, distab) 478 !------------------------------------------------------------------ 479 ! Qbi_\alpha\alpha' 480 r_coef_h => jrho1_atom_set(katom)%cjc_iii_h(ispin)%r_coef 481 r_coef_s => jrho1_atom_set(katom)%cjc_iii_s(ispin)%r_coef 482 CALL proj_blk(C_coeff_hh_a, C_coeff_ss_a, n_cont_a, & 483 C_coeff_hh_b, C_coeff_ss_b, n_cont_b, & 484 d_matrix, max_nsgf, r_coef_h, r_coef_s, nso, & 485 len_PC1, len_CPC, 1.0_dp, distab) 486 !------------------------------------------------------------------ 487!$ CALL omp_unset_lock(proj_blk_lock((katom - 1)*nspins + ispin)) 488 489 ENDDO ! ispin 490 491 EXIT !search loop over jatom-katom list 492 END IF 493 END DO 494 END DO 495 496 ENDDO ! kkind 497 DO ispin = 1, nspins 498 DEALLOCATE (jp_RARnu(ispin)%r_coef, jp2_RARnu(ispin)%r_coef) 499 ENDDO 500 END DO 501 ! Wait for all threads to finish the loop before locks can be freed 502!$OMP BARRIER 503 504!$OMP DO 505!$ DO lock = 1, natom 506!$ call omp_destroy_lock(alloc_lock(lock)) 507!$ END DO 508!$OMP END DO 509 510!$OMP DO 511!$ DO lock = 1, nspins*natom 512!$ call omp_destroy_lock(proj_blk_lock(lock)) 513!$ END DO 514!$OMP END DO 515 516!$OMP SINGLE 517!$ DEALLOCATE (alloc_lock) 518!$ DEALLOCATE (proj_blk_lock) 519!$OMP END SINGLE NOWAIT 520 521 DEALLOCATE (a_matrix, b_matrix, c_matrix, d_matrix, & 522 a_block, b_block, c_block, d_block, & 523 jp_RARnu, jp2_RARnu & 524 ) 525 526!$OMP END PARALLEL 527 528 CALL neighbor_list_iterator_release(nl_iterator) 529 DEALLOCATE (basis_set_list) 530 531 ! parallel sum up 532 nbr_dbl = 0.0_dp 533 DO ikind = 1, nkind 534 CALL get_atomic_kind(atomic_kind_set(ikind), & 535 atom_list=atom_list, & 536 natom=nat) 537 CALL get_qs_kind(qs_kind_set(ikind), & 538 basis_set=orb_basis_set, & 539 paw_proj_set=paw_proj, & 540 paw_atom=paw_atom) 541 542 IF (.NOT. paw_atom) CYCLE 543 544 CALL get_paw_proj_set(paw_proj_set=paw_proj, nsatbas=nsatbas) 545 nsoctot = nsatbas 546 547 num_pe = para_env%num_pe 548 mepos = para_env%mepos 549 bo = get_limit(nat, num_pe, mepos) 550 551 ALLOCATE (zero_coeff(nsoctot, nsoctot)) 552 DO iat = 1, nat 553 iatom = atom_list(iat) 554 is_not_associated = .NOT. ASSOCIATED(jrho1_atom_set(iatom)%cjc0_h(1)%r_coef) 555 556 IF (iat .GE. bo(1) .AND. iat .LE. bo(2) .AND. is_not_associated) THEN 557 CALL allocate_jrho_coeff(jrho1_atom_set, iatom, nsoctot) 558 ENDIF 559 560 DO ispin = 1, nspins 561 562 tmp_coeff => jrho1_atom_set(iatom)%cjc0_h(ispin)%r_coef 563 IF (is_not_associated) THEN 564 zero_coeff = 0.0_dp; tmp_coeff => zero_coeff 565 ENDIF 566 CALL mp_sum(tmp_coeff, para_env%group) 567 568 tmp_coeff => jrho1_atom_set(iatom)%cjc0_s(ispin)%r_coef 569 IF (is_not_associated) THEN 570 zero_coeff = 0.0_dp; tmp_coeff => zero_coeff 571 ENDIF 572 CALL mp_sum(tmp_coeff, para_env%group) 573 574 tmp_coeff => jrho1_atom_set(iatom)%cjc_h(ispin)%r_coef 575 IF (is_not_associated) THEN 576 zero_coeff = 0.0_dp; tmp_coeff => zero_coeff 577 ENDIF 578 579 CALL mp_sum(tmp_coeff, para_env%group) 580 tmp_coeff => jrho1_atom_set(iatom)%cjc_s(ispin)%r_coef 581 IF (is_not_associated) THEN 582 zero_coeff = 0.0_dp; tmp_coeff => zero_coeff 583 ENDIF 584 CALL mp_sum(tmp_coeff, para_env%group) 585 586 tmp_coeff => jrho1_atom_set(iatom)%cjc_ii_h(ispin)%r_coef 587 IF (is_not_associated) THEN 588 zero_coeff = 0.0_dp; tmp_coeff => zero_coeff 589 ENDIF 590 CALL mp_sum(tmp_coeff, para_env%group) 591 592 tmp_coeff => jrho1_atom_set(iatom)%cjc_ii_s(ispin)%r_coef 593 IF (is_not_associated) THEN 594 zero_coeff = 0.0_dp; tmp_coeff => zero_coeff 595 ENDIF 596 CALL mp_sum(tmp_coeff, para_env%group) 597 598 tmp_coeff => jrho1_atom_set(iatom)%cjc_iii_h(ispin)%r_coef 599 IF (is_not_associated) THEN 600 zero_coeff = 0.0_dp; tmp_coeff => zero_coeff 601 ENDIF 602 CALL mp_sum(tmp_coeff, para_env%group) 603 604 tmp_coeff => jrho1_atom_set(iatom)%cjc_iii_s(ispin)%r_coef 605 IF (is_not_associated) THEN 606 zero_coeff = 0.0_dp; tmp_coeff => zero_coeff 607 ENDIF 608 CALL mp_sum(tmp_coeff, para_env%group) 609 610 IF (ASSOCIATED(jrho1_atom_set(iatom)%cjc0_h(ispin)%r_coef)) & 611 nbr_dbl = nbr_dbl + 8.0_dp*REAL(SIZE(jrho1_atom_set(iatom)%cjc0_h(ispin)%r_coef), dp) 612 ENDDO ! ispin 613 ENDDO ! iat 614 615 DEALLOCATE (zero_coeff) 616 617 ENDDO ! ikind 618 619 output_unit = cp_logger_get_default_io_unit() 620 IF (output_unit > 0) THEN 621 WRITE (output_unit, '(A,E8.2)') 'calculate_jrho_atom_coeff: nbr_dbl=', nbr_dbl 622 ENDIF 623 624 CALL timestop(handle) 625 626 END SUBROUTINE calculate_jrho_atom_coeff 627 628! ************************************************************************************************** 629!> \brief ... 630!> \param qs_env ... 631!> \param current_env ... 632!> \param idir ... 633! ************************************************************************************************** 634 SUBROUTINE calculate_jrho_atom_rad(qs_env, current_env, idir) 635 ! 636 TYPE(qs_environment_type), POINTER :: qs_env 637 TYPE(current_env_type) :: current_env 638 INTEGER, INTENT(IN) :: idir 639 640 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_jrho_atom_rad', & 641 routineP = moduleN//':'//routineN 642 643 INTEGER :: damax_iso_not0, damax_iso_not0_local, handle, i1, i2, iat, iatom, icg, ikind, & 644 ipgf1, ipgf2, ir, iset1, iset2, iso, iso1, iso1_first, iso1_last, iso2, iso2_first, & 645 iso2_last, ispin, l, l_iso, llmax, lmax12, lmax_expansion, lmin12, m1s, m2s, m_iso, & 646 max_iso_not0, max_iso_not0_local, max_max_iso_not0, max_nso, max_s_harm, maxl, maxlgto, & 647 maxso, mepos, n1s, n2s, na, natom, natom_tot, nkind, nr, nset, nspins, num_pe, size1, & 648 size2 649 INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list, dacg_n_list 650 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list, dacg_list 651 INTEGER, DIMENSION(2) :: bo 652 INTEGER, DIMENSION(:), POINTER :: atom_list, lmax, lmin, npgf, o2nindex 653 LOGICAL :: paw_atom 654 LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: is_set_to_0 655 REAL(dp) :: hard_radius 656 REAL(dp), ALLOCATABLE, DIMENSION(:) :: g1, g2, gauge_h, gauge_s 657 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: cjc0_h_block, cjc0_s_block, cjc_h_block, & 658 cjc_ii_h_block, cjc_ii_s_block, cjc_iii_h_block, cjc_iii_s_block, cjc_s_block, dgg_1, gg, & 659 gg_lm1 660 REAL(dp), DIMENSION(:, :), POINTER :: coeff, Fr_a_h, Fr_a_h_ii, Fr_a_h_iii, Fr_a_s, & 661 Fr_a_s_ii, Fr_a_s_iii, Fr_b_h, Fr_b_h_ii, Fr_b_h_iii, Fr_b_s, Fr_b_s_ii, Fr_b_s_iii, & 662 Fr_h, Fr_s, zet 663 REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG 664 REAL(dp), DIMENSION(:, :, :, :), POINTER :: my_CG_dxyz_asym 665 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 666 TYPE(cp_para_env_type), POINTER :: para_env 667 TYPE(dft_control_type), POINTER :: dft_control 668 TYPE(grid_atom_type), POINTER :: grid_atom 669 TYPE(gto_basis_set_type), POINTER :: orb_basis_set 670 TYPE(harmonics_atom_type), POINTER :: harmonics 671 TYPE(jrho_atom_type), DIMENSION(:), POINTER :: jrho1_atom_set 672 TYPE(jrho_atom_type), POINTER :: jrho1_atom 673 TYPE(paw_proj_set_type), POINTER :: paw_proj 674 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 675 676 CALL timeset(routineN, handle) 677 ! 678 NULLIFY (atomic_kind_set, qs_kind_set, dft_control, para_env, & 679 coeff, Fr_h, Fr_s, Fr_a_h, Fr_a_s, Fr_a_h_ii, Fr_a_s_ii, & 680 Fr_a_h_iii, Fr_a_s_iii, Fr_b_h, Fr_b_s, Fr_b_h_ii, & 681 Fr_b_s_ii, Fr_b_h_iii, Fr_b_s_iii, jrho1_atom_set, & 682 jrho1_atom) 683 ! 684 CALL get_qs_env(qs_env=qs_env, & 685 atomic_kind_set=atomic_kind_set, & 686 qs_kind_set=qs_kind_set, & 687 dft_control=dft_control, & 688 para_env=para_env) 689 690 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxlgto=maxlgto) 691 692 ! 693 CALL get_current_env(current_env=current_env, & 694 jrho1_atom_set=jrho1_atom_set) 695 ! 696 697 nkind = SIZE(qs_kind_set) 698 nspins = dft_control%nspins 699 ! 700 natom_tot = SIZE(jrho1_atom_set, 1) 701 ALLOCATE (is_set_to_0(natom_tot, nspins)) 702 is_set_to_0(:, :) = .FALSE. 703 704 ! 705 DO ikind = 1, nkind 706 NULLIFY (atom_list, grid_atom, harmonics, orb_basis_set, & 707 lmax, lmin, npgf, zet, grid_atom, harmonics, my_CG, my_CG_dxyz_asym) 708 ! 709 CALL get_atomic_kind(atomic_kind_set(ikind), & 710 atom_list=atom_list, & 711 natom=natom) 712 CALL get_qs_kind(qs_kind_set(ikind), & 713 grid_atom=grid_atom, & 714 paw_proj_set=paw_proj, & 715 paw_atom=paw_atom, & 716 harmonics=harmonics, & 717 hard_radius=hard_radius, & 718 basis_set=orb_basis_set) 719 ! 720 ! Quick cycle if needed. 721 IF (.NOT. paw_atom) CYCLE 722 ! 723 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, & 724 lmax=lmax, lmin=lmin, & 725 maxl=maxl, npgf=npgf, & 726 nset=nset, zet=zet, & 727 maxso=maxso) 728 CALL get_paw_proj_set(paw_proj_set=paw_proj, o2nindex=o2nindex) 729 ! 730 nr = grid_atom%nr 731 na = grid_atom%ng_sphere 732 max_iso_not0 = harmonics%max_iso_not0 733 damax_iso_not0 = harmonics%damax_iso_not0 734 max_max_iso_not0 = MAX(max_iso_not0, damax_iso_not0) 735 lmax_expansion = indso(1, max_max_iso_not0) 736 max_s_harm = harmonics%max_s_harm 737 llmax = harmonics%llmax 738 ! 739 ! Distribute the atoms of this kind 740 num_pe = para_env%num_pe 741 mepos = para_env%mepos 742 bo = get_limit(natom, num_pe, mepos) 743 ! 744 my_CG => harmonics%my_CG 745 my_CG_dxyz_asym => harmonics%my_CG_dxyz_asym 746 ! 747 ! Allocate some arrays. 748 max_nso = nsoset(maxl) 749 ALLOCATE (g1(nr), g2(nr), gg(nr, 0:2*maxl), gg_lm1(nr, 0:2*maxl), dgg_1(nr, 0:2*maxl), & 750 cjc0_h_block(max_nso, max_nso), cjc0_s_block(max_nso, max_nso), & 751 cjc_h_block(max_nso, max_nso), cjc_s_block(max_nso, max_nso), & 752 cjc_ii_h_block(max_nso, max_nso), cjc_ii_s_block(max_nso, max_nso), & 753 cjc_iii_h_block(max_nso, max_nso), cjc_iii_s_block(max_nso, max_nso), & 754 cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm), & 755 dacg_list(2, nsoset(maxl)**2, max_s_harm), dacg_n_list(max_s_harm), & 756 gauge_h(nr), gauge_s(nr)) 757 ! 758 ! Compute the gauge 759 SELECT CASE (current_env%gauge) 760 CASE (current_gauge_r) 761 ! d(r)=r 762 gauge_h(1:nr) = grid_atom%rad(1:nr) 763 gauge_s(1:nr) = grid_atom%rad(1:nr) 764 CASE (current_gauge_r_and_step_func) 765 ! step function 766 gauge_h(1:nr) = 0e0_dp 767 DO ir = 1, nr 768 IF (grid_atom%rad(ir) .LE. hard_radius) THEN 769 gauge_s(ir) = grid_atom%rad(ir) 770 ELSE 771 gauge_s(ir) = gauge_h(ir) 772 ENDIF 773 ENDDO 774 CASE (current_gauge_atom) 775 ! d(r)=A 776 gauge_h(1:nr) = HUGE(0e0_dp) !0e0_dp 777 gauge_s(1:nr) = HUGE(0e0_dp) !0e0_dp 778 CASE DEFAULT 779 CPABORT("Unknown gauge, try again...") 780 END SELECT 781 ! 782 ! 783 m1s = 0 784 DO iset1 = 1, nset 785 m2s = 0 786 DO iset2 = 1, nset 787 CALL get_none0_cg_list(my_CG, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), & 788 max_s_harm, lmax_expansion, cg_list, cg_n_list, max_iso_not0_local) 789 CPASSERT(max_iso_not0_local .LE. max_iso_not0) 790 CALL get_none0_cg_list(my_CG_dxyz_asym, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), & 791 max_s_harm, lmax_expansion, dacg_list, dacg_n_list, damax_iso_not0_local) 792 CPASSERT(damax_iso_not0_local .LE. damax_iso_not0) 793 794 n1s = nsoset(lmax(iset1)) 795 DO ipgf1 = 1, npgf(iset1) 796 iso1_first = nsoset(lmin(iset1) - 1) + 1 + n1s*(ipgf1 - 1) + m1s 797 iso1_last = nsoset(lmax(iset1)) + n1s*(ipgf1 - 1) + m1s 798 size1 = iso1_last - iso1_first + 1 799 iso1_first = o2nindex(iso1_first) 800 iso1_last = o2nindex(iso1_last) 801 i1 = iso1_last - iso1_first + 1 802 CPASSERT(size1 == i1) 803 i1 = nsoset(lmin(iset1) - 1) + 1 804 ! 805 g1(1:nr) = EXP(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr)) 806 ! 807 n2s = nsoset(lmax(iset2)) 808 DO ipgf2 = 1, npgf(iset2) 809 iso2_first = nsoset(lmin(iset2) - 1) + 1 + n2s*(ipgf2 - 1) + m2s 810 iso2_last = nsoset(lmax(iset2)) + n2s*(ipgf2 - 1) + m2s 811 size2 = iso2_last - iso2_first + 1 812 iso2_first = o2nindex(iso2_first) 813 iso2_last = o2nindex(iso2_last) 814 i2 = iso2_last - iso2_first + 1 815 CPASSERT(size2 == i2) 816 i2 = nsoset(lmin(iset2) - 1) + 1 817 ! 818 g2(1:nr) = EXP(-zet(ipgf2, iset2)*grid_atom%rad2(1:nr)) 819 ! 820 lmin12 = lmin(iset1) + lmin(iset2) 821 lmax12 = lmax(iset1) + lmax(iset2) 822 ! 823 gg = 0.0_dp 824 gg_lm1 = 0.0_dp 825 dgg_1 = 0.0_dp 826 ! 827 ! Take only the terms of expansion for L < lmax_expansion 828 IF (lmin12 .LE. lmax_expansion) THEN 829 ! 830 IF (lmax12 .GT. lmax_expansion) lmax12 = lmax_expansion 831 ! 832 IF (lmin12 == 0) THEN 833 gg(1:nr, lmin12) = g1(1:nr)*g2(1:nr) 834 gg_lm1(1:nr, lmin12) = 0.0_dp 835 ELSE 836 gg(1:nr, lmin12) = grid_atom%rad2l(1:nr, lmin12)*g1(1:nr)*g2(1:nr) 837 gg_lm1(1:nr, lmin12) = grid_atom%rad2l(1:nr, lmin12 - 1)*g1(1:nr)*g2(1:nr) 838 ENDIF 839 ! 840 DO l = lmin12 + 1, lmax12 841 gg(1:nr, l) = grid_atom%rad(1:nr)*gg(1:nr, l - 1) 842 gg_lm1(1:nr, l) = gg(1:nr, l - 1) 843 ENDDO 844 ! 845 DO l = lmin12, lmax12 846 dgg_1(1:nr, l) = 2.0_dp*(zet(ipgf1, iset1) - zet(ipgf2, iset2))& 847 & *gg(1:nr, l)*grid_atom%rad(1:nr) 848 ENDDO 849 ELSE 850 CYCLE 851 ENDIF ! lmin12 852 ! 853 DO iat = bo(1), bo(2) 854 iatom = atom_list(iat) 855 ! 856 DO ispin = 1, nspins 857 !------------------------------------------------------------------ 858 ! P_\alpha\alpha' 859 cjc0_h_block = HUGE(1.0_dp) 860 cjc0_s_block = HUGE(1.0_dp) 861 ! 862 ! Hard term 863 coeff => jrho1_atom_set(iatom)%cjc0_h(ispin)%r_coef 864 cjc0_h_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = & 865 coeff(iso1_first:iso1_last, iso2_first:iso2_last) 866 ! 867 ! Soft term 868 coeff => jrho1_atom_set(iatom)%cjc0_s(ispin)%r_coef 869 cjc0_s_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = & 870 coeff(iso1_first:iso1_last, iso2_first:iso2_last) 871 !------------------------------------------------------------------ 872 ! mQai_\alpha\alpha' 873 cjc_h_block = HUGE(1.0_dp) 874 cjc_s_block = HUGE(1.0_dp) 875 ! 876 ! Hard term 877 coeff => jrho1_atom_set(iatom)%cjc_h(ispin)%r_coef 878 cjc_h_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = & 879 coeff(iso1_first:iso1_last, iso2_first:iso2_last) 880 ! 881 ! Soft term 882 coeff => jrho1_atom_set(iatom)%cjc_s(ispin)%r_coef 883 cjc_s_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = & 884 coeff(iso1_first:iso1_last, iso2_first:iso2_last) 885 !------------------------------------------------------------------ 886 ! Qci_\alpha\alpha' 887 cjc_ii_h_block = HUGE(1.0_dp) 888 cjc_ii_s_block = HUGE(1.0_dp) 889 ! 890 ! Hard term 891 coeff => jrho1_atom_set(iatom)%cjc_ii_h(ispin)%r_coef 892 cjc_ii_h_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = & 893 coeff(iso1_first:iso1_last, iso2_first:iso2_last) 894 ! 895 ! Soft term 896 coeff => jrho1_atom_set(iatom)%cjc_ii_s(ispin)%r_coef 897 cjc_ii_s_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = & 898 coeff(iso1_first:iso1_last, iso2_first:iso2_last) 899 !------------------------------------------------------------------ 900 ! Qbi_\alpha\alpha' 901 cjc_iii_h_block = HUGE(1.0_dp) 902 cjc_iii_s_block = HUGE(1.0_dp) 903 ! 904 ! 905 ! Hard term 906 coeff => jrho1_atom_set(iatom)%cjc_iii_h(ispin)%r_coef 907 cjc_iii_h_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = & 908 coeff(iso1_first:iso1_last, iso2_first:iso2_last) 909 ! 910 ! Soft term 911 coeff => jrho1_atom_set(iatom)%cjc_iii_s(ispin)%r_coef 912 cjc_iii_s_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = & 913 coeff(iso1_first:iso1_last, iso2_first:iso2_last) 914 !------------------------------------------------------------------ 915 ! 916 ! Allocation radial functions 917 jrho1_atom => jrho1_atom_set(iatom) 918 IF (.NOT. ASSOCIATED(jrho1_atom%jrho_a_h(ispin)%r_coef)) THEN 919 CALL allocate_jrho_atom_rad(jrho1_atom, ispin, nr, na, & 920 max_max_iso_not0) 921 is_set_to_0(iatom, ispin) = .TRUE. 922 ELSE 923 IF (.NOT. is_set_to_0(iatom, ispin)) THEN 924 CALL set2zero_jrho_atom_rad(jrho1_atom, ispin) 925 is_set_to_0(iatom, ispin) = .TRUE. 926 ENDIF 927 ENDIF 928 !------------------------------------------------------------------ 929 ! 930 Fr_h => jrho1_atom%jrho_h(ispin)%r_coef 931 Fr_s => jrho1_atom%jrho_s(ispin)%r_coef 932 !------------------------------------------------------------------ 933 ! 934 Fr_a_h => jrho1_atom%jrho_a_h(ispin)%r_coef 935 Fr_a_s => jrho1_atom%jrho_a_s(ispin)%r_coef 936 Fr_b_h => jrho1_atom%jrho_b_h(ispin)%r_coef 937 Fr_b_s => jrho1_atom%jrho_b_s(ispin)%r_coef 938 !------------------------------------------------------------------ 939 ! 940 Fr_a_h_ii => jrho1_atom%jrho_a_h_ii(ispin)%r_coef 941 Fr_a_s_ii => jrho1_atom%jrho_a_s_ii(ispin)%r_coef 942 Fr_b_h_ii => jrho1_atom%jrho_b_h_ii(ispin)%r_coef 943 Fr_b_s_ii => jrho1_atom%jrho_b_s_ii(ispin)%r_coef 944 !------------------------------------------------------------------ 945 ! 946 Fr_a_h_iii => jrho1_atom%jrho_a_h_iii(ispin)%r_coef 947 Fr_a_s_iii => jrho1_atom%jrho_a_s_iii(ispin)%r_coef 948 Fr_b_h_iii => jrho1_atom%jrho_b_h_iii(ispin)%r_coef 949 Fr_b_s_iii => jrho1_atom%jrho_b_s_iii(ispin)%r_coef 950 !------------------------------------------------------------------ 951 ! 952 DO iso = 1, max_iso_not0_local 953 l_iso = indso(1, iso) ! not needed 954 m_iso = indso(2, iso) ! not needed 955 ! 956 DO icg = 1, cg_n_list(iso) 957 ! 958 iso1 = cg_list(1, icg, iso) 959 iso2 = cg_list(2, icg, iso) 960 ! 961 IF (.NOT. (iso2 > 0 .AND. iso1 > 0)) THEN 962 WRITE (*, *) 'iso1=', iso1, ' iso2=', iso2, ' iso=', iso, ' icg=', icg 963 WRITE (*, *) '.... will stop!' 964 ENDIF 965 CPASSERT(iso2 > 0 .AND. iso1 > 0) 966 ! 967 l = indso(1, iso1) + indso(1, iso2) 968 IF (l .GT. lmax_expansion .OR. l .LT. .0) THEN 969 WRITE (*, *) 'calculate_jrho_atom_rad: 1 l', l 970 WRITE (*, *) 'calculate_jrho_atom_rad: 1 lmax_expansion', lmax_expansion 971 WRITE (*, *) '.... will stop!' 972 ENDIF 973 CPASSERT(l <= lmax_expansion) 974 !------------------------------------------------------------------ 975 ! P0 976 ! 977 IF (current_env%gauge .EQ. current_gauge_atom) THEN 978 ! Hard term 979 Fr_h(1:nr, iso) = Fr_h(1:nr, iso) + & 980 gg(1:nr, l)*cjc0_h_block(iso1, iso2)* & 981 my_CG(iso1, iso2, iso) 982 ! Soft term 983 Fr_s(1:nr, iso) = Fr_s(1:nr, iso) + & 984 gg(1:nr, l)*cjc0_s_block(iso1, iso2)* & 985 my_CG(iso1, iso2, iso) 986 ELSE 987 ! Hard term 988 Fr_h(1:nr, iso) = Fr_h(1:nr, iso) + & 989 gg(1:nr, l)*cjc0_h_block(iso1, iso2)* & 990 my_CG(iso1, iso2, iso)*(grid_atom%rad(1:nr) - gauge_h(1:nr)) 991 ! Soft term 992 Fr_s(1:nr, iso) = Fr_s(1:nr, iso) + & 993 gg(1:nr, l)*cjc0_s_block(iso1, iso2)* & 994 my_CG(iso1, iso2, iso)*(grid_atom%rad(1:nr) - gauge_s(1:nr)) 995 ENDIF 996 !------------------------------------------------------------------ 997 ! Rai 998 ! 999 ! Hard term 1000 Fr_a_h(1:nr, iso) = Fr_a_h(1:nr, iso) + & 1001 dgg_1(1:nr, l)*cjc_h_block(iso1, iso2)* & 1002 my_CG(iso1, iso2, iso) 1003 ! 1004 ! Soft term 1005 Fr_a_s(1:nr, iso) = Fr_a_s(1:nr, iso) + & 1006 dgg_1(1:nr, l)*cjc_s_block(iso1, iso2)* & 1007 my_CG(iso1, iso2, iso) 1008 !------------------------------------------------------------------ 1009 ! Rci 1010 ! 1011 IF (current_env%gauge .EQ. current_gauge_atom) THEN 1012 ! Hard term 1013 Fr_a_h_ii(1:nr, iso) = Fr_a_h_ii(1:nr, iso) + & 1014 dgg_1(1:nr, l)* & 1015 cjc_ii_h_block(iso1, iso2)* & 1016 my_CG(iso1, iso2, iso) 1017 ! Soft term 1018 Fr_a_s_ii(1:nr, iso) = Fr_a_s_ii(1:nr, iso) + & 1019 dgg_1(1:nr, l)* & 1020 cjc_ii_s_block(iso1, iso2)* & 1021 my_CG(iso1, iso2, iso) 1022 ELSE 1023 ! Hard term 1024 Fr_a_h_ii(1:nr, iso) = Fr_a_h_ii(1:nr, iso) + & 1025 dgg_1(1:nr, l)*gauge_h(1:nr)* & 1026 cjc_ii_h_block(iso1, iso2)* & 1027 my_CG(iso1, iso2, iso) 1028 ! Soft term 1029 Fr_a_s_ii(1:nr, iso) = Fr_a_s_ii(1:nr, iso) + & 1030 dgg_1(1:nr, l)*gauge_s(1:nr)* & 1031 cjc_ii_s_block(iso1, iso2)* & 1032 my_CG(iso1, iso2, iso) 1033 ENDIF 1034 !------------------------------------------------------------------ 1035 ! Rbi 1036 ! 1037 IF (current_env%gauge .EQ. current_gauge_atom) THEN 1038 ! Hard term 1039 Fr_a_h_iii(1:nr, iso) = Fr_a_h_iii(1:nr, iso) + & 1040 dgg_1(1:nr, l)* & 1041 cjc_iii_h_block(iso1, iso2)* & 1042 my_CG(iso1, iso2, iso) 1043 ! Soft term 1044 Fr_a_s_iii(1:nr, iso) = Fr_a_s_iii(1:nr, iso) + & 1045 dgg_1(1:nr, l)* & 1046 cjc_iii_s_block(iso1, iso2)* & 1047 my_CG(iso1, iso2, iso) 1048 ELSE 1049 ! Hard term 1050 Fr_a_h_iii(1:nr, iso) = Fr_a_h_iii(1:nr, iso) + & 1051 dgg_1(1:nr, l)*gauge_h(1:nr)* & 1052 cjc_iii_h_block(iso1, iso2)* & 1053 my_CG(iso1, iso2, iso) 1054 ! Soft term 1055 Fr_a_s_iii(1:nr, iso) = Fr_a_s_iii(1:nr, iso) + & 1056 dgg_1(1:nr, l)*gauge_s(1:nr)* & 1057 cjc_iii_s_block(iso1, iso2)* & 1058 my_CG(iso1, iso2, iso) 1059 ENDIF 1060 !------------------------------------------------------------------ 1061 ENDDO !icg 1062 ! 1063 ENDDO ! iso 1064 ! 1065 ! 1066 DO iso = 1, damax_iso_not0_local 1067 l_iso = indso(1, iso) ! not needed 1068 m_iso = indso(2, iso) ! not needed 1069 ! 1070 DO icg = 1, dacg_n_list(iso) 1071 ! 1072 iso1 = dacg_list(1, icg, iso) 1073 iso2 = dacg_list(2, icg, iso) 1074 ! 1075 IF (.NOT. (iso2 > 0 .AND. iso1 > 0)) THEN 1076 WRITE (*, *) 'iso1=', iso1, ' iso2=', iso2, ' iso=', iso, ' icg=', icg 1077 WRITE (*, *) '.... will stop!' 1078 ENDIF 1079 CPASSERT(iso2 > 0 .AND. iso1 > 0) 1080 ! 1081 l = indso(1, iso1) + indso(1, iso2) 1082 IF (l .GT. lmax_expansion) THEN 1083 WRITE (*, *) 'calculate_jrho_atom_rad: 1 l', l 1084 WRITE (*, *) 'calculate_jrho_atom_rad: 1 lmax_expansion', lmax_expansion 1085 WRITE (*, *) '.... will stop!' 1086 ENDIF 1087 CPASSERT(l <= lmax_expansion) 1088 !------------------------------------------------------------------ 1089 ! Daij 1090 ! 1091 ! Hard term 1092 Fr_b_h(1:nr, iso) = Fr_b_h(1:nr, iso) + & 1093 gg_lm1(1:nr, l)*cjc_h_block(iso1, iso2)* & 1094 my_CG_dxyz_asym(idir, iso1, iso2, iso) 1095 ! 1096 ! Soft term 1097 Fr_b_s(1:nr, iso) = Fr_b_s(1:nr, iso) + & 1098 gg_lm1(1:nr, l)*cjc_s_block(iso1, iso2)* & 1099 my_CG_dxyz_asym(idir, iso1, iso2, iso) 1100 ! 1101 !------------------------------------------------------------------ 1102 ! Dcij 1103 ! 1104 IF (current_env%gauge .EQ. current_gauge_atom) THEN 1105 ! Hard term 1106 Fr_b_h_ii(1:nr, iso) = Fr_b_h_ii(1:nr, iso) + & 1107 gg_lm1(1:nr, l)* & 1108 cjc_ii_h_block(iso1, iso2)* & 1109 my_CG_dxyz_asym(idir, iso1, iso2, iso) 1110 ! Soft term 1111 Fr_b_s_ii(1:nr, iso) = Fr_b_s_ii(1:nr, iso) + & 1112 gg_lm1(1:nr, l)* & 1113 cjc_ii_s_block(iso1, iso2)* & 1114 my_CG_dxyz_asym(idir, iso1, iso2, iso) 1115 ELSE 1116 ! Hard term 1117 Fr_b_h_ii(1:nr, iso) = Fr_b_h_ii(1:nr, iso) + & 1118 gg_lm1(1:nr, l)*gauge_h(1:nr)* & 1119 cjc_ii_h_block(iso1, iso2)* & 1120 my_CG_dxyz_asym(idir, iso1, iso2, iso) 1121 ! Soft term 1122 Fr_b_s_ii(1:nr, iso) = Fr_b_s_ii(1:nr, iso) + & 1123 gg_lm1(1:nr, l)*gauge_s(1:nr)* & 1124 cjc_ii_s_block(iso1, iso2)* & 1125 my_CG_dxyz_asym(idir, iso1, iso2, iso) 1126 ENDIF 1127 !------------------------------------------------------------------ 1128 ! Dbij 1129 ! 1130 IF (current_env%gauge .EQ. current_gauge_atom) THEN 1131 ! Hard term 1132 Fr_b_h_iii(1:nr, iso) = Fr_b_h_iii(1:nr, iso) + & 1133 gg_lm1(1:nr, l)* & 1134 cjc_iii_h_block(iso1, iso2)* & 1135 my_CG_dxyz_asym(idir, iso1, iso2, iso) 1136 ! Soft term 1137 Fr_b_s_iii(1:nr, iso) = Fr_b_s_iii(1:nr, iso) + & 1138 gg_lm1(1:nr, l)* & 1139 cjc_iii_s_block(iso1, iso2)* & 1140 my_CG_dxyz_asym(idir, iso1, iso2, iso) 1141 ELSE 1142 ! Hard term 1143 Fr_b_h_iii(1:nr, iso) = Fr_b_h_iii(1:nr, iso) + & 1144 gg_lm1(1:nr, l)*gauge_h(1:nr)* & 1145 cjc_iii_h_block(iso1, iso2)* & 1146 my_CG_dxyz_asym(idir, iso1, iso2, iso) 1147 ! Soft term 1148 Fr_b_s_iii(1:nr, iso) = Fr_b_s_iii(1:nr, iso) + & 1149 gg_lm1(1:nr, l)*gauge_s(1:nr)* & 1150 cjc_iii_s_block(iso1, iso2)* & 1151 my_CG_dxyz_asym(idir, iso1, iso2, iso) 1152 ENDIF 1153 !------------------------------------------------------------------ 1154 ENDDO ! icg 1155 ENDDO ! iso 1156 ! 1157 ENDDO ! ispin 1158 ENDDO ! iat 1159 ! 1160 !------------------------------------------------------------------ 1161 ! 1162 ENDDO !ipgf2 1163 ENDDO ! ipgf1 1164 m2s = m2s + maxso 1165 ENDDO ! iset2 1166 m1s = m1s + maxso 1167 ENDDO ! iset1 1168 ! 1169 DEALLOCATE (cjc0_h_block, cjc0_s_block, cjc_h_block, cjc_s_block, cjc_ii_h_block, cjc_ii_s_block, & 1170 cjc_iii_h_block, cjc_iii_s_block, g1, g2, gg, gg_lm1, dgg_1, gauge_h, gauge_s, & 1171 cg_list, cg_n_list, dacg_list, dacg_n_list) 1172 ENDDO ! ikind 1173 ! 1174 ! 1175 DEALLOCATE (is_set_to_0) 1176 ! 1177 CALL timestop(handle) 1178 ! 1179 END SUBROUTINE calculate_jrho_atom_rad 1180 1181! ************************************************************************************************** 1182!> \brief ... 1183!> \param jrho1_atom ... 1184!> \param jrho_h ... 1185!> \param jrho_s ... 1186!> \param grid_atom ... 1187!> \param harmonics ... 1188!> \param do_igaim ... 1189!> \param ratom ... 1190!> \param natm_gauge ... 1191!> \param iB ... 1192!> \param idir ... 1193!> \param ispin ... 1194! ************************************************************************************************** 1195 SUBROUTINE calculate_jrho_atom_ang(jrho1_atom, jrho_h, jrho_s, grid_atom, & 1196 harmonics, do_igaim, ratom, natm_gauge, & 1197 iB, idir, ispin) 1198 ! 1199 TYPE(jrho_atom_type), POINTER :: jrho1_atom 1200 REAL(dp), DIMENSION(:, :), POINTER :: jrho_h, jrho_s 1201 TYPE(grid_atom_type), POINTER :: grid_atom 1202 TYPE(harmonics_atom_type), POINTER :: harmonics 1203 LOGICAL, INTENT(IN) :: do_igaim 1204 INTEGER, INTENT(IN) :: iB, idir, ispin, natm_gauge 1205 REAL(dp), INTENT(IN) :: ratom(:, :) 1206 1207 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_jrho_atom_ang', & 1208 routineP = moduleN//':'//routineN 1209 1210 INTEGER :: ia, idir2, iiB, iiiB, ir, & 1211 iso, max_max_iso_not0, na, nr 1212 REAL(dp) :: rad_part, scale 1213 REAL(dp), DIMENSION(:, :), POINTER :: a, Fr_a_h, Fr_a_h_ii, Fr_a_h_iii, & 1214 Fr_a_s, Fr_a_s_ii, Fr_a_s_iii, Fr_b_h, Fr_b_h_ii, Fr_b_h_iii, Fr_b_s, & 1215 Fr_b_s_ii, Fr_b_s_iii, Fr_h, Fr_s, slm 1216 REAL(dp), DIMENSION(:), POINTER :: r 1217 REAL(dp), DIMENSION(:, :, :), ALLOCATABLE :: g 1218! 1219! 1220 NULLIFY (Fr_h, Fr_s, Fr_a_h, Fr_a_s, Fr_a_h_ii, Fr_a_s_ii, Fr_a_h_iii, Fr_a_s_iii, & 1221 Fr_b_h, Fr_b_s, Fr_b_h_ii, Fr_b_s_ii, Fr_b_h_iii, Fr_b_s_iii, & 1222 a, slm) 1223 ! 1224 CPASSERT(ASSOCIATED(jrho_h)) 1225 CPASSERT(ASSOCIATED(jrho_s)) 1226 CPASSERT(ASSOCIATED(jrho1_atom)) 1227 ! just to be sure... 1228 CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_h)) 1229 CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_s)) 1230 CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_h)) 1231 CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_s)) 1232 CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_h(ispin)%r_coef)) 1233 CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_s(ispin)%r_coef)) 1234 CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_h(ispin)%r_coef)) 1235 CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_s(ispin)%r_coef)) 1236 CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_h_ii)) 1237 CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_s_ii)) 1238 CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_h_ii)) 1239 CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_s_ii)) 1240 CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_h_ii(ispin)%r_coef)) 1241 CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_s_ii(ispin)%r_coef)) 1242 CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_h_ii(ispin)%r_coef)) 1243 CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_s_ii(ispin)%r_coef)) 1244 CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_h_iii)) 1245 CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_s_iii)) 1246 CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_h_iii)) 1247 CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_s_iii)) 1248 CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_h_iii(ispin)%r_coef)) 1249 CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_s_iii(ispin)%r_coef)) 1250 CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_h_iii(ispin)%r_coef)) 1251 CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_s_iii(ispin)%r_coef)) 1252 ! 1253 ! 1254 nr = grid_atom%nr 1255 na = grid_atom%ng_sphere 1256 max_max_iso_not0 = MAX(harmonics%max_iso_not0, harmonics%damax_iso_not0) 1257 ALLOCATE (g(3, nr, na)) 1258 !------------------------------------------------------------------ 1259 ! 1260 Fr_h => jrho1_atom%jrho_h(ispin)%r_coef 1261 Fr_s => jrho1_atom%jrho_s(ispin)%r_coef 1262 !------------------------------------------------------------------ 1263 ! 1264 Fr_a_h => jrho1_atom%jrho_a_h(ispin)%r_coef !Rai 1265 Fr_a_s => jrho1_atom%jrho_a_s(ispin)%r_coef 1266 Fr_b_h => jrho1_atom%jrho_b_h(ispin)%r_coef !Daij 1267 Fr_b_s => jrho1_atom%jrho_b_s(ispin)%r_coef 1268 !------------------------------------------------------------------ 1269 ! 1270 Fr_a_h_ii => jrho1_atom%jrho_a_h_ii(ispin)%r_coef !Rci 1271 Fr_a_s_ii => jrho1_atom%jrho_a_s_ii(ispin)%r_coef 1272 Fr_b_h_ii => jrho1_atom%jrho_b_h_ii(ispin)%r_coef !Dcij 1273 Fr_b_s_ii => jrho1_atom%jrho_b_s_ii(ispin)%r_coef 1274 !------------------------------------------------------------------ 1275 ! 1276 Fr_a_h_iii => jrho1_atom%jrho_a_h_iii(ispin)%r_coef !Rbi 1277 Fr_a_s_iii => jrho1_atom%jrho_a_s_iii(ispin)%r_coef 1278 Fr_b_h_iii => jrho1_atom%jrho_b_h_iii(ispin)%r_coef !Dbij 1279 Fr_b_s_iii => jrho1_atom%jrho_b_s_iii(ispin)%r_coef 1280 !------------------------------------------------------------------ 1281 ! 1282 a => harmonics%a 1283 slm => harmonics%slm 1284 r => grid_atom%rad 1285 ! 1286 CALL set_vecp(iB, iiB, iiiB) 1287 ! 1288 scale = 0.0_dp 1289 idir2 = 1 1290 IF (idir .NE. iB) THEN 1291 CALL set_vecp_rev(idir, iB, idir2) 1292 scale = fac_vecp(idir, iB, idir2) 1293 ENDIF 1294 ! 1295 ! Set the gauge 1296 CALL get_gauge() 1297 ! 1298 DO ir = 1, nr 1299 DO iso = 1, max_max_iso_not0 1300 DO ia = 1, na 1301 IF (do_igaim) THEN 1302 !------------------------------------------------------------------ 1303 ! Hard current density response 1304 ! radial(ia,ir) = ( aj(ia) * Rai(ir,iso) + Daij 1305 ! - aii(ia) * ( aj(ia) * Rbi(ir,iso) + Dbij ) 1306 ! + aiii(ia) * ( aj(ia) * Rci(ir,iso) + Dcij ) 1307 ! ) * Ylm(ia) 1308 rad_part = a(idir, ia)*Fr_a_h(ir, iso) + Fr_b_h(ir, iso) & 1309 & - g(iiB, ir, ia)*(a(idir, ia)*Fr_a_h_iii(ir, iso) + Fr_b_h_iii(ir, iso))& 1310 & + g(iiiB, ir, ia)*(a(idir, ia)*Fr_a_h_ii(ir, iso) + Fr_b_h_ii(ir, iso))& 1311 & + scale*(a(idir2, ia)*r(ir) - g(idir2, ir, ia))*Fr_h(ir, iso) 1312 ! 1313 jrho_h(ir, ia) = jrho_h(ir, ia) + rad_part*slm(ia, iso) 1314 !------------------------------------------------------------------ 1315 ! Soft current density response 1316 rad_part = a(idir, ia)*Fr_a_s(ir, iso) + Fr_b_s(ir, iso) & 1317 & - g(iiB, ir, ia)*(a(idir, ia)*Fr_a_s_iii(ir, iso) + Fr_b_s_iii(ir, iso))& 1318 & + g(iiiB, ir, ia)*(a(idir, ia)*Fr_a_s_ii(ir, iso) + Fr_b_s_ii(ir, iso))& 1319 & + scale*(a(idir2, ia)*r(ir) - g(idir2, ir, ia))*Fr_s(ir, iso) 1320 ! 1321 jrho_s(ir, ia) = jrho_s(ir, ia) + rad_part*slm(ia, iso) 1322 !------------------------------------------------------------------ 1323 ELSE 1324 !------------------------------------------------------------------ 1325 ! Hard current density response 1326 ! radial(ia,ir) = ( aj(ia) * Rai(ir,iso) + Daij 1327 ! - aii(ia) * ( aj(ia) * Rbi(ir,iso) + Dbij ) 1328 ! + aiii(ia) * ( aj(ia) * Rci(ir,iso) + Dcij ) 1329 ! ) * Ylm(ia) 1330 rad_part = a(idir, ia)*Fr_a_h(ir, iso) + Fr_b_h(ir, iso) & 1331 & - a(iiB, ia)*(a(idir, ia)*Fr_a_h_iii(ir, iso) + Fr_b_h_iii(ir, iso))& 1332 & + a(iiiB, ia)*(a(idir, ia)*Fr_a_h_ii(ir, iso) + Fr_b_h_ii(ir, iso))& 1333 & + scale*a(idir2, ia)*Fr_h(ir, iso) 1334 ! 1335 jrho_h(ir, ia) = jrho_h(ir, ia) + rad_part*slm(ia, iso) 1336 !------------------------------------------------------------------ 1337 ! Soft current density response 1338 rad_part = a(idir, ia)*Fr_a_s(ir, iso) + Fr_b_s(ir, iso) & 1339 & - a(iiB, ia)*(a(idir, ia)*Fr_a_s_iii(ir, iso) + Fr_b_s_iii(ir, iso))& 1340 & + a(iiiB, ia)*(a(idir, ia)*Fr_a_s_ii(ir, iso) + Fr_b_s_ii(ir, iso))& 1341 & + scale*a(idir2, ia)*Fr_s(ir, iso) 1342 ! 1343 jrho_s(ir, ia) = jrho_s(ir, ia) + rad_part*slm(ia, iso) 1344 !------------------------------------------------------------------ 1345 ENDIF 1346 ENDDO ! ia 1347 ENDDO ! iso 1348 ENDDO ! ir 1349 ! 1350 DEALLOCATE (g) 1351 ! 1352 CONTAINS 1353 ! 1354! ************************************************************************************************** 1355!> \brief ... 1356! ************************************************************************************************** 1357 SUBROUTINE get_gauge() 1358 INTEGER :: iatom, ixyz, jatom 1359 REAL(dp) :: ab, pa, pb, point(3), pra(3), prb(3), & 1360 res, tmp 1361 REAL(dp), ALLOCATABLE, DIMENSION(:) :: buf 1362 1363 ALLOCATE (buf(natm_gauge)) 1364 DO ir = 1, nr 1365 DO ia = 1, na 1366 DO ixyz = 1, 3 1367 g(ixyz, ir, ia) = 0.0_dp 1368 ENDDO 1369 point(1) = r(ir)*a(1, ia) 1370 point(2) = r(ir)*a(2, ia) 1371 point(3) = r(ir)*a(3, ia) 1372 DO iatom = 1, natm_gauge 1373 buf(iatom) = 1.0_dp 1374 pra = point - ratom(:, iatom) 1375 pa = SQRT(pra(1)**2 + pra(2)**2 + pra(3)**2) 1376 DO jatom = 1, natm_gauge 1377 IF (iatom .EQ. jatom) CYCLE 1378 prb = point - ratom(:, jatom) 1379 pb = SQRT(prb(1)**2 + prb(2)**2 + prb(3)**2) 1380 ab = SQRT((pra(1) - prb(1))**2 + (pra(2) - prb(2))**2 + (pra(3) - prb(3))**2) 1381 ! 1382 tmp = (pa - pb)/ab 1383 tmp = 0.5_dp*(3.0_dp - tmp**2)*tmp 1384 tmp = 0.5_dp*(3.0_dp - tmp**2)*tmp 1385 tmp = 0.5_dp*(3.0_dp - tmp**2)*tmp 1386 buf(iatom) = buf(iatom)*0.5_dp*(1.0_dp - tmp) 1387 ENDDO 1388 ENDDO 1389 DO ixyz = 1, 3 1390 res = 0.0_dp 1391 DO iatom = 1, natm_gauge 1392 res = res + ratom(ixyz, iatom)*buf(iatom) 1393 ENDDO 1394 res = res/SUM(buf(1:natm_gauge)) 1395 ! 1396 g(ixyz, ir, ia) = res 1397 ENDDO 1398 ENDDO 1399 ENDDO 1400 DEALLOCATE (buf) 1401 END SUBROUTINE get_gauge 1402 END SUBROUTINE calculate_jrho_atom_ang 1403 1404! ************************************************************************************************** 1405!> \brief ... 1406!> \param current_env ... 1407!> \param qs_env ... 1408!> \param iB ... 1409!> \param idir ... 1410! ************************************************************************************************** 1411 SUBROUTINE calculate_jrho_atom(current_env, qs_env, iB, idir) 1412 TYPE(current_env_type) :: current_env 1413 TYPE(qs_environment_type), POINTER :: qs_env 1414 INTEGER, INTENT(IN) :: iB, idir 1415 1416 CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_jrho_atom', & 1417 routineP = moduleN//':'//routineN 1418 1419 INTEGER :: iat, iatom, ikind, ispin, jatom, & 1420 natm_gauge, natm_tot, natom, nkind, & 1421 nspins 1422 INTEGER, DIMENSION(2) :: bo 1423 INTEGER, DIMENSION(:), POINTER :: atom_list 1424 LOGICAL :: do_igaim, gapw, paw_atom 1425 REAL(dp) :: hard_radius, r(3) 1426 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: ratom 1427 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 1428 TYPE(cell_type), POINTER :: cell 1429 TYPE(cp_para_env_type), POINTER :: para_env 1430 TYPE(dft_control_type), POINTER :: dft_control 1431 TYPE(grid_atom_type), POINTER :: grid_atom 1432 TYPE(harmonics_atom_type), POINTER :: harmonics 1433 TYPE(jrho_atom_type), DIMENSION(:), POINTER :: jrho1_atom_set 1434 TYPE(jrho_atom_type), POINTER :: jrho1_atom 1435 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 1436 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 1437 1438 NULLIFY (para_env, dft_control) 1439 NULLIFY (jrho1_atom_set, grid_atom, harmonics) 1440 NULLIFY (atomic_kind_set, qs_kind_set, atom_list) 1441 1442 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, & 1443 atomic_kind_set=atomic_kind_set, & 1444 qs_kind_set=qs_kind_set, & 1445 particle_set=particle_set, & 1446 cell=cell, & 1447 para_env=para_env) 1448 1449 CALL get_current_env(current_env=current_env, & 1450 jrho1_atom_set=jrho1_atom_set) 1451 1452 do_igaim = .FALSE. 1453 IF (current_env%gauge .EQ. current_gauge_atom) do_igaim = .TRUE. 1454 1455 gapw = dft_control%qs_control%gapw 1456 nkind = SIZE(qs_kind_set, 1) 1457 nspins = dft_control%nspins 1458 1459 natm_tot = SIZE(particle_set) 1460 ALLOCATE (ratom(3, natm_tot)) 1461 1462 IF (gapw) THEN 1463 DO ikind = 1, nkind 1464 NULLIFY (atom_list, grid_atom, harmonics) 1465 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom) 1466 CALL get_qs_kind(qs_kind_set(ikind), & 1467 grid_atom=grid_atom, & 1468 harmonics=harmonics, & 1469 hard_radius=hard_radius, & 1470 paw_atom=paw_atom) 1471 1472 IF (.NOT. paw_atom) CYCLE 1473 1474 ! Distribute the atoms of this kind 1475 1476 bo = get_limit(natom, para_env%num_pe, para_env%mepos) 1477 1478 DO iat = bo(1), bo(2) 1479 iatom = atom_list(iat) 1480 NULLIFY (jrho1_atom) 1481 jrho1_atom => jrho1_atom_set(iatom) 1482 1483 natm_gauge = 0 1484 DO jatom = 1, natm_tot 1485 r(:) = pbc(particle_set(jatom)%r(:) - particle_set(iatom)%r(:), cell) 1486 ! SQRT(SUM(r(:)**2)) .LE. 2.0_dp*hard_radius 1487 IF (SUM(r(:)**2) .LE. (4.0_dp*hard_radius*hard_radius)) THEN 1488 natm_gauge = natm_gauge + 1 1489 ratom(:, natm_gauge) = r(:) 1490 ENDIF 1491 ENDDO 1492 1493 DO ispin = 1, nspins 1494 jrho1_atom%jrho_vec_rad_h(idir, ispin)%r_coef = 0.0_dp 1495 jrho1_atom%jrho_vec_rad_s(idir, ispin)%r_coef = 0.0_dp 1496 CALL calculate_jrho_atom_ang(jrho1_atom, & 1497 jrho1_atom%jrho_vec_rad_h(idir, ispin)%r_coef, & 1498 jrho1_atom%jrho_vec_rad_s(idir, ispin)%r_coef, & 1499 grid_atom, harmonics, & 1500 do_igaim, & 1501 ratom, natm_gauge, iB, idir, ispin) 1502 END DO !ispin 1503 END DO !iat 1504 END DO !ikind 1505 END IF 1506 1507 DEALLOCATE (ratom) 1508 1509 END SUBROUTINE calculate_jrho_atom 1510 1511END MODULE qs_linres_atom_current 1512