(git:ccc2433)
qs_linres_atom_current.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief given the response wavefunctions obtained by the application
10 !> of the (rxp), p, and ((dk-dl)xp) operators,
11 !> here the current density vector (jx, jy, jz)
12 !> is computed for the 3 directions of the magnetic field (Bx, By, Bz)
13 !> \par History
14 !> created 02-2006 [MI]
15 !> \author MI
16 ! **************************************************************************************************
17 MODULE qs_linres_atom_current
18  USE atomic_kind_types, ONLY: atomic_kind_type,&
19  get_atomic_kind
20  USE basis_set_types, ONLY: get_gto_basis_set,&
21  gto_basis_set_p_type,&
22  gto_basis_set_type
23  USE cell_types, ONLY: cell_type,&
24  pbc
25  USE cp_control_types, ONLY: dft_control_type
26  USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
27  USE dbcsr_api, ONLY: dbcsr_get_block_p,&
28  dbcsr_p_type
29  USE input_constants, ONLY: current_gauge_atom,&
30  current_gauge_r,&
31  current_gauge_r_and_step_func
32  USE kinds, ONLY: dp
33  USE message_passing, ONLY: mp_para_env_type
34  USE orbital_pointers, ONLY: indso,&
35  nsoset
36  USE particle_types, ONLY: particle_type
37  USE paw_basis_types, ONLY: get_paw_basis_info
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 
87 CONTAINS
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 
116  INTEGER :: bo(2), handle, iac, iat, iatom, ibc, idir2, ii, iii, ikind, ispin, istat, jatom, &
117  jkind, kac, katom, kbc, kkind, len_cpc, len_pc1, max_gau, max_nsgf, mepos, n_cont_a, &
118  n_cont_b, nat, natom, nkind, nsatbas, nsgfa, nsgfb, nso, nsoctot, nspins, num_pe, &
119  output_unit
120  INTEGER, DIMENSION(3) :: cell_b
121  INTEGER, DIMENSION(:), POINTER :: atom_list, list_a, list_b
122  LOGICAL :: den_found, dista, distab, distb, &
123  is_not_associated, paw_atom, &
124  sgf_soft_only_a, sgf_soft_only_b
125  REAL(dp) :: eps_cpc, jmax, nbr_dbl, rab(3), rbc(3)
126  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: a_matrix, b_matrix, c_matrix, d_matrix
127  REAL(kind=dp), DIMENSION(:, :), POINTER :: c_coeff_hh_a, c_coeff_hh_b, &
128  c_coeff_ss_a, c_coeff_ss_b, r_coef_h, &
129  r_coef_s, tmp_coeff, zero_coeff
130  TYPE(alist_type), POINTER :: alist_ac, alist_bc
131  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
132  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_a, mat_b, mat_c, mat_d
133  TYPE(dft_control_type), POINTER :: dft_control
134  TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
135  TYPE(gto_basis_set_type), POINTER :: basis_1c_set, basis_set_a, basis_set_b
136  TYPE(jrho_atom_type), DIMENSION(:), POINTER :: jrho1_atom_set
137  TYPE(mp_para_env_type), POINTER :: para_env
138  TYPE(neighbor_list_iterator_p_type), &
139  DIMENSION(:), POINTER :: nl_iterator
140  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
141  POINTER :: sab_all
142  TYPE(oce_matrix_type), POINTER :: oce
143  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
144  TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: a_block, b_block, c_block, d_block, &
145  jp2_rarnu, jp_rarnu
146 
147 !$ INTEGER(kind=omp_lock_kind), ALLOCATABLE, DIMENSION(:) :: proj_blk_lock, alloc_lock
148 !$ INTEGER :: lock
149 
150  CALL timeset(routinen, handle)
151 
152  NULLIFY (atomic_kind_set, qs_kind_set, dft_control, sab_all, jrho1_atom_set, oce, &
153  para_env, zero_coeff, tmp_coeff)
154 
155  CALL get_qs_env(qs_env=qs_env, &
156  atomic_kind_set=atomic_kind_set, &
157  qs_kind_set=qs_kind_set, &
158  dft_control=dft_control, &
159  oce=oce, &
160  sab_all=sab_all, &
161  para_env=para_env)
162  cpassert(ASSOCIATED(oce))
163 
164  CALL get_current_env(current_env=current_env, jrho1_atom_set=jrho1_atom_set)
165  cpassert(ASSOCIATED(jrho1_atom_set))
166 
167  CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
168  maxsgf=max_nsgf, maxgtops=max_gau, basis_type="GAPW_1C")
169 
170  eps_cpc = dft_control%qs_control%gapw_control%eps_cpc
171 
172  idir2 = 1
173  IF (idir .NE. ib) THEN
174  CALL set_vecp_rev(idir, ib, idir2)
175  END IF
176  CALL set_vecp(ib, ii, iii)
177 
178  ! Set pointers for the different gauge
179  mat_a => mat_d0
180  mat_b => mat_jp
181  mat_c => mat_jp_rii
182  mat_d => mat_jp_riii
183 
184  ! Density-like matrices
185  nkind = SIZE(qs_kind_set)
186  natom = SIZE(jrho1_atom_set)
187  nspins = dft_control%nspins
188 
189  ! Reset CJC coefficients and local density arrays
190  DO ikind = 1, nkind
191  NULLIFY (atom_list)
192  CALL get_atomic_kind(atomic_kind_set(ikind), &
193  atom_list=atom_list, &
194  natom=nat)
195  CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
196 
197  ! Quick cycle if needed.
198  IF (.NOT. paw_atom) cycle
199 
200  ! Initialize the density matrix-like arrays.
201  DO iat = 1, nat
202  iatom = atom_list(iat)
203  DO ispin = 1, nspins
204  IF (ASSOCIATED(jrho1_atom_set(iatom)%cjc0_h(1)%r_coef)) THEN
205  jrho1_atom_set(iatom)%cjc0_h(ispin)%r_coef = 0.0_dp
206  jrho1_atom_set(iatom)%cjc0_s(ispin)%r_coef = 0.0_dp
207  jrho1_atom_set(iatom)%cjc_h(ispin)%r_coef = 0.0_dp
208  jrho1_atom_set(iatom)%cjc_s(ispin)%r_coef = 0.0_dp
209  jrho1_atom_set(iatom)%cjc_ii_h(ispin)%r_coef = 0.0_dp
210  jrho1_atom_set(iatom)%cjc_ii_s(ispin)%r_coef = 0.0_dp
211  jrho1_atom_set(iatom)%cjc_iii_h(ispin)%r_coef = 0.0_dp
212  jrho1_atom_set(iatom)%cjc_iii_s(ispin)%r_coef = 0.0_dp
213  END IF
214  END DO ! ispin
215  END DO ! iat
216  END DO ! ikind
217 
218  ! Three centers
219  ALLOCATE (basis_set_list(nkind))
220  DO ikind = 1, nkind
221  CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set_a)
222  IF (ASSOCIATED(basis_set_a)) THEN
223  basis_set_list(ikind)%gto_basis_set => basis_set_a
224  ELSE
225  NULLIFY (basis_set_list(ikind)%gto_basis_set)
226  END IF
227  END DO
228 
229  len_pc1 = max_nsgf*max_gau
230  len_cpc = max_gau*max_gau
231 
232  num_pe = 1
233 !$ num_pe = omp_get_max_threads()
234  CALL neighbor_list_iterator_create(nl_iterator, sab_all, nthread=num_pe)
235 
236 !$OMP PARALLEL DEFAULT( NONE ) &
237 !$OMP SHARED( nspins, max_nsgf, max_gau &
238 !$OMP , len_PC1, len_CPC &
239 !$OMP , nl_iterator, basis_set_list &
240 !$OMP , mat_a, mat_b, mat_c, mat_d &
241 !$OMP , nkind, qs_kind_set, eps_cpc, oce &
242 !$OMP , ii, iii, jrho1_atom_set &
243 !$OMP , natom, proj_blk_lock, alloc_lock &
244 !$OMP ) &
245 !$OMP PRIVATE( a_block, b_block, c_block, d_block &
246 !$OMP , jp_RARnu, jp2_RARnu &
247 !$OMP , a_matrix, b_matrix, c_matrix, d_matrix, istat &
248 !$OMP , mepos &
249 !$OMP , ikind, jkind, kkind, iatom, jatom, katom &
250 !$OMP , cell_b, rab, rbc &
251 !$OMP , basis_set_a, nsgfa &
252 !$OMP , basis_set_b, nsgfb &
253 !$OMP , basis_1c_set, jmax, den_found &
254 !$OMP , nsatbas, nsoctot, nso, paw_atom &
255 !$OMP , iac , alist_ac, kac, n_cont_a, list_a, sgf_soft_only_a &
256 !$OMP , ibc , alist_bc, kbc, n_cont_b, list_b, sgf_soft_only_b &
257 !$OMP , C_coeff_hh_a, C_coeff_ss_a, dista &
258 !$OMP , C_coeff_hh_b, C_coeff_ss_b, distb &
259 !$OMP , distab &
260 !$OMP , r_coef_s, r_coef_h &
261 !$OMP )
262 
263  NULLIFY (a_block, b_block, c_block, d_block)
264  NULLIFY (basis_1c_set, jp_rarnu, jp2_rarnu)
265 
266  ALLOCATE (a_block(nspins), b_block(nspins), c_block(nspins), d_block(nspins), &
267  jp_rarnu(nspins), jp2_rarnu(nspins), &
268  stat=istat)
269  cpassert(istat == 0)
270 
271  ALLOCATE (a_matrix(max_nsgf, max_nsgf), b_matrix(max_nsgf, max_nsgf), &
272  c_matrix(max_nsgf, max_nsgf), d_matrix(max_nsgf, max_nsgf), &
273  stat=istat)
274  cpassert(istat == 0)
275 
276 !$OMP SINGLE
277 !$ ALLOCATE (alloc_lock(natom))
278 !$ ALLOCATE (proj_blk_lock(nspins*natom))
279 !$OMP END SINGLE
280 
281 !$OMP DO
282 !$ DO lock = 1, natom
283 !$ call omp_init_lock(alloc_lock(lock))
284 !$ END DO
285 !$OMP END DO
286 
287 !$OMP DO
288 !$ DO lock = 1, nspins*natom
289 !$ call omp_init_lock(proj_blk_lock(lock))
290 !$ END DO
291 !$OMP END DO
292 
293  mepos = 0
294 !$ mepos = omp_get_thread_num()
295 
296  DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
297 
298  CALL get_iterator_info(nl_iterator, mepos=mepos, &
299  ikind=ikind, jkind=jkind, &
300  iatom=iatom, jatom=jatom, cell=cell_b, r=rab)
301  basis_set_a => basis_set_list(ikind)%gto_basis_set
302  IF (.NOT. ASSOCIATED(basis_set_a)) cycle
303  basis_set_b => basis_set_list(jkind)%gto_basis_set
304  IF (.NOT. ASSOCIATED(basis_set_b)) cycle
305  nsgfa = basis_set_a%nsgf
306  nsgfb = basis_set_b%nsgf
307  DO ispin = 1, nspins
308  NULLIFY (jp_rarnu(ispin)%r_coef, jp2_rarnu(ispin)%r_coef)
309  ALLOCATE (jp_rarnu(ispin)%r_coef(nsgfa, nsgfb), &
310  jp2_rarnu(ispin)%r_coef(nsgfa, nsgfb))
311  END DO
312 
313  ! Take the block \mu\nu of jpab, jpab_ii and jpab_iii
314  jmax = 0._dp
315  DO ispin = 1, nspins
316  NULLIFY (a_block(ispin)%r_coef)
317  NULLIFY (b_block(ispin)%r_coef)
318  NULLIFY (c_block(ispin)%r_coef)
319  NULLIFY (d_block(ispin)%r_coef)
320  CALL dbcsr_get_block_p(matrix=mat_a(ispin)%matrix, &
321  row=iatom, col=jatom, block=a_block(ispin)%r_coef, &
322  found=den_found)
323  jmax = jmax + maxval(abs(a_block(ispin)%r_coef))
324  CALL dbcsr_get_block_p(matrix=mat_b(ispin)%matrix, &
325  row=iatom, col=jatom, block=b_block(ispin)%r_coef, &
326  found=den_found)
327  jmax = jmax + maxval(abs(b_block(ispin)%r_coef))
328  CALL dbcsr_get_block_p(matrix=mat_c(ispin)%matrix, &
329  row=iatom, col=jatom, block=c_block(ispin)%r_coef, &
330  found=den_found)
331  jmax = jmax + maxval(abs(c_block(ispin)%r_coef))
332  CALL dbcsr_get_block_p(matrix=mat_d(ispin)%matrix, &
333  row=iatom, col=jatom, block=d_block(ispin)%r_coef, &
334  found=den_found)
335  jmax = jmax + maxval(abs(d_block(ispin)%r_coef))
336  END DO
337 
338  ! Loop over atoms
339  DO kkind = 1, nkind
340  CALL get_qs_kind(qs_kind_set(kkind), &
341  basis_set=basis_1c_set, basis_type="GAPW_1C", &
342  paw_atom=paw_atom)
343 
344  ! Quick cycle if needed.
345  IF (.NOT. paw_atom) cycle
346 
347  CALL get_paw_basis_info(basis_1c_set, nsatbas=nsatbas)
348  nsoctot = nsatbas
349 
350  iac = ikind + nkind*(kkind - 1)
351  ibc = jkind + nkind*(kkind - 1)
352  IF (.NOT. ASSOCIATED(oce%intac(iac)%alist)) cycle
353  IF (.NOT. ASSOCIATED(oce%intac(ibc)%alist)) cycle
354 
355  CALL get_alist(oce%intac(iac), alist_ac, iatom)
356  CALL get_alist(oce%intac(ibc), alist_bc, jatom)
357  IF (.NOT. ASSOCIATED(alist_ac)) cycle
358  IF (.NOT. ASSOCIATED(alist_bc)) cycle
359 
360  DO kac = 1, alist_ac%nclist
361  DO kbc = 1, alist_bc%nclist
362  IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) cycle
363 
364  IF (all(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
365  IF (jmax*alist_bc%clist(kbc)%maxac*alist_ac%clist(kac)%maxac < eps_cpc) cycle
366 
367  n_cont_a = alist_ac%clist(kac)%nsgf_cnt
368  n_cont_b = alist_bc%clist(kbc)%nsgf_cnt
369  sgf_soft_only_a = alist_ac%clist(kac)%sgf_soft_only
370  sgf_soft_only_b = alist_bc%clist(kbc)%sgf_soft_only
371  IF (n_cont_a .EQ. 0 .OR. n_cont_b .EQ. 0) cycle
372 
373  ! thanks to the linearity of the response, we
374  ! can avoid computing soft-soft interactions.
375  ! those terms are already included in the
376  ! regular grid.
377  IF (sgf_soft_only_a .AND. sgf_soft_only_b) cycle
378 
379  list_a => alist_ac%clist(kac)%sgf_list
380  list_b => alist_bc%clist(kbc)%sgf_list
381 
382  katom = alist_ac%clist(kac)%catom
383 !$ CALL omp_set_lock(alloc_lock(katom))
384  IF (.NOT. ASSOCIATED(jrho1_atom_set(katom)%cjc0_h(1)%r_coef)) THEN
385  CALL allocate_jrho_coeff(jrho1_atom_set, katom, nsoctot)
386  END IF
387 !$ CALL omp_unset_lock(alloc_lock(katom))
388 
389  ! Compute the modified Qai matrix as
390  ! mQai_\mu\nu = Qai_\mu\nu - Qbi_\mu\nu * (R_A-R_\nu)_ii
391  ! + Qci_\mu\nu * (R_A-R_\nu)_iii
392  rbc = alist_bc%clist(kbc)%rac
393  DO ispin = 1, nspins
394  CALL dcopy(nsgfa*nsgfb, b_block(ispin)%r_coef(1, 1), 1, &
395  jp_rarnu(ispin)%r_coef(1, 1), 1)
396  CALL daxpy(nsgfa*nsgfb, -rbc(ii), d_block(ispin)%r_coef(1, 1), 1, &
397  jp_rarnu(ispin)%r_coef(1, 1), 1)
398  CALL daxpy(nsgfa*nsgfb, rbc(iii), c_block(ispin)%r_coef(1, 1), 1, &
399  jp_rarnu(ispin)%r_coef(1, 1), 1)
400  END DO
401 
402  ! Get the d_A's for the hard and soft densities.
403  IF (iatom == katom .AND. all(alist_ac%clist(kac)%cell == 0)) THEN
404  c_coeff_hh_a => alist_ac%clist(kac)%achint(:, :, 1)
405  c_coeff_ss_a => alist_ac%clist(kac)%acint(:, :, 1)
406  dista = .false.
407  ELSE
408  c_coeff_hh_a => alist_ac%clist(kac)%acint(:, :, 1)
409  c_coeff_ss_a => alist_ac%clist(kac)%acint(:, :, 1)
410  dista = .true.
411  END IF
412  ! Get the d_B's for the hard and soft densities.
413  IF (jatom == katom .AND. all(alist_bc%clist(kbc)%cell == 0)) THEN
414  c_coeff_hh_b => alist_bc%clist(kbc)%achint(:, :, 1)
415  c_coeff_ss_b => alist_bc%clist(kbc)%acint(:, :, 1)
416  distb = .false.
417  ELSE
418  c_coeff_hh_b => alist_bc%clist(kbc)%acint(:, :, 1)
419  c_coeff_ss_b => alist_bc%clist(kbc)%acint(:, :, 1)
420  distb = .true.
421  END IF
422 
423  distab = dista .AND. distb
424 
425  nso = nsoctot
426 
427  DO ispin = 1, nspins
428 
429  ! align the blocks
430  CALL alist_pre_align_blk(a_block(ispin)%r_coef, SIZE(a_block(ispin)%r_coef, 1), &
431  a_matrix, SIZE(a_matrix, 1), &
432  list_a, n_cont_a, list_b, n_cont_b)
433 
434  CALL alist_pre_align_blk(jp_rarnu(ispin)%r_coef, SIZE(jp_rarnu(ispin)%r_coef, 1), &
435  b_matrix, SIZE(b_matrix, 1), &
436  list_a, n_cont_a, list_b, n_cont_b)
437 
438  CALL alist_pre_align_blk(c_block(ispin)%r_coef, SIZE(c_block(ispin)%r_coef, 1), &
439  c_matrix, SIZE(c_matrix, 1), &
440  list_a, n_cont_a, list_b, n_cont_b)
441  CALL alist_pre_align_blk(d_block(ispin)%r_coef, SIZE(d_block(ispin)%r_coef, 1), &
442  d_matrix, SIZE(d_matrix, 1), &
443  list_a, n_cont_a, list_b, n_cont_b)
444 
445 !$ CALL omp_set_lock(proj_blk_lock((katom - 1)*nspins + ispin))
446  !------------------------------------------------------------------
447  ! P_\alpha\alpha'
448  r_coef_h => jrho1_atom_set(katom)%cjc0_h(ispin)%r_coef
449  r_coef_s => jrho1_atom_set(katom)%cjc0_s(ispin)%r_coef
450  CALL proj_blk(c_coeff_hh_a, c_coeff_ss_a, n_cont_a, &
451  c_coeff_hh_b, c_coeff_ss_b, n_cont_b, &
452  a_matrix, max_nsgf, r_coef_h, r_coef_s, nso, &
453  len_pc1, len_cpc, 1.0_dp, distab)
454  !------------------------------------------------------------------
455  ! mQai_\alpha\alpha'
456  r_coef_h => jrho1_atom_set(katom)%cjc_h(ispin)%r_coef
457  r_coef_s => jrho1_atom_set(katom)%cjc_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  b_matrix, max_nsgf, r_coef_h, r_coef_s, nso, &
461  len_pc1, len_cpc, 1.0_dp, distab)
462  !------------------------------------------------------------------
463  ! Qci_\alpha\alpha'
464  r_coef_h => jrho1_atom_set(katom)%cjc_ii_h(ispin)%r_coef
465  r_coef_s => jrho1_atom_set(katom)%cjc_ii_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  c_matrix, max_nsgf, r_coef_h, r_coef_s, nso, &
469  len_pc1, len_cpc, 1.0_dp, distab)
470  !------------------------------------------------------------------
471  ! Qbi_\alpha\alpha'
472  r_coef_h => jrho1_atom_set(katom)%cjc_iii_h(ispin)%r_coef
473  r_coef_s => jrho1_atom_set(katom)%cjc_iii_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  d_matrix, max_nsgf, r_coef_h, r_coef_s, nso, &
477  len_pc1, len_cpc, 1.0_dp, distab)
478  !------------------------------------------------------------------
479 !$ CALL omp_unset_lock(proj_blk_lock((katom - 1)*nspins + ispin))
480 
481  END DO ! ispin
482 
483  EXIT !search loop over jatom-katom list
484  END IF
485  END DO
486  END DO
487 
488  END DO ! kkind
489  DO ispin = 1, nspins
490  DEALLOCATE (jp_rarnu(ispin)%r_coef, jp2_rarnu(ispin)%r_coef)
491  END DO
492  END DO
493  ! Wait for all threads to finish the loop before locks can be freed
494 !$OMP BARRIER
495 
496 !$OMP DO
497 !$ DO lock = 1, natom
498 !$ call omp_destroy_lock(alloc_lock(lock))
499 !$ END DO
500 !$OMP END DO
501 
502 !$OMP DO
503 !$ DO lock = 1, nspins*natom
504 !$ call omp_destroy_lock(proj_blk_lock(lock))
505 !$ END DO
506 !$OMP END DO
507 
508 !$OMP SINGLE
509 !$ DEALLOCATE (alloc_lock)
510 !$ DEALLOCATE (proj_blk_lock)
511 !$OMP END SINGLE NOWAIT
512 
513  DEALLOCATE (a_matrix, b_matrix, c_matrix, d_matrix, &
514  a_block, b_block, c_block, d_block, &
515  jp_rarnu, jp2_rarnu &
516  )
517 
518 !$OMP END PARALLEL
519 
520  CALL neighbor_list_iterator_release(nl_iterator)
521  DEALLOCATE (basis_set_list)
522 
523  ! parallel sum up
524  nbr_dbl = 0.0_dp
525  DO ikind = 1, nkind
526  CALL get_atomic_kind(atomic_kind_set(ikind), &
527  atom_list=atom_list, &
528  natom=nat)
529  CALL get_qs_kind(qs_kind_set(ikind), &
530  basis_set=basis_1c_set, basis_type="GAPW_1C", &
531  paw_atom=paw_atom)
532 
533  IF (.NOT. paw_atom) cycle
534 
535  CALL get_paw_basis_info(basis_1c_set, nsatbas=nsatbas)
536  nsoctot = nsatbas
537 
538  num_pe = para_env%num_pe
539  mepos = para_env%mepos
540  bo = get_limit(nat, num_pe, mepos)
541 
542  ALLOCATE (zero_coeff(nsoctot, nsoctot))
543  DO iat = 1, nat
544  iatom = atom_list(iat)
545  is_not_associated = .NOT. ASSOCIATED(jrho1_atom_set(iatom)%cjc0_h(1)%r_coef)
546 
547  IF (iat .GE. bo(1) .AND. iat .LE. bo(2) .AND. is_not_associated) THEN
548  CALL allocate_jrho_coeff(jrho1_atom_set, iatom, nsoctot)
549  END IF
550 
551  DO ispin = 1, nspins
552 
553  tmp_coeff => jrho1_atom_set(iatom)%cjc0_h(ispin)%r_coef
554  IF (is_not_associated) THEN
555  zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
556  END IF
557  CALL para_env%sum(tmp_coeff)
558 
559  tmp_coeff => jrho1_atom_set(iatom)%cjc0_s(ispin)%r_coef
560  IF (is_not_associated) THEN
561  zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
562  END IF
563  CALL para_env%sum(tmp_coeff)
564 
565  tmp_coeff => jrho1_atom_set(iatom)%cjc_h(ispin)%r_coef
566  IF (is_not_associated) THEN
567  zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
568  END IF
569 
570  CALL para_env%sum(tmp_coeff)
571  tmp_coeff => jrho1_atom_set(iatom)%cjc_s(ispin)%r_coef
572  IF (is_not_associated) THEN
573  zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
574  END IF
575  CALL para_env%sum(tmp_coeff)
576 
577  tmp_coeff => jrho1_atom_set(iatom)%cjc_ii_h(ispin)%r_coef
578  IF (is_not_associated) THEN
579  zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
580  END IF
581  CALL para_env%sum(tmp_coeff)
582 
583  tmp_coeff => jrho1_atom_set(iatom)%cjc_ii_s(ispin)%r_coef
584  IF (is_not_associated) THEN
585  zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
586  END IF
587  CALL para_env%sum(tmp_coeff)
588 
589  tmp_coeff => jrho1_atom_set(iatom)%cjc_iii_h(ispin)%r_coef
590  IF (is_not_associated) THEN
591  zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
592  END IF
593  CALL para_env%sum(tmp_coeff)
594 
595  tmp_coeff => jrho1_atom_set(iatom)%cjc_iii_s(ispin)%r_coef
596  IF (is_not_associated) THEN
597  zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
598  END IF
599  CALL para_env%sum(tmp_coeff)
600 
601  IF (ASSOCIATED(jrho1_atom_set(iatom)%cjc0_h(ispin)%r_coef)) &
602  nbr_dbl = nbr_dbl + 8.0_dp*real(SIZE(jrho1_atom_set(iatom)%cjc0_h(ispin)%r_coef), dp)
603  END DO ! ispin
604  END DO ! iat
605 
606  DEALLOCATE (zero_coeff)
607 
608  END DO ! ikind
609 
610  output_unit = cp_logger_get_default_io_unit()
611  IF (output_unit > 0) THEN
612  WRITE (output_unit, '(A,E8.2)') 'calculate_jrho_atom_coeff: nbr_dbl=', nbr_dbl
613  END IF
614 
615  CALL timestop(handle)
616 
617  END SUBROUTINE calculate_jrho_atom_coeff
618 
619 ! **************************************************************************************************
620 !> \brief ...
621 !> \param qs_env ...
622 !> \param current_env ...
623 !> \param idir ...
624 ! **************************************************************************************************
625  SUBROUTINE calculate_jrho_atom_rad(qs_env, current_env, idir)
626  !
627  TYPE(qs_environment_type), POINTER :: qs_env
628  TYPE(current_env_type) :: current_env
629  INTEGER, INTENT(IN) :: idir
630 
631  CHARACTER(len=*), PARAMETER :: routinen = 'calculate_jrho_atom_rad'
632 
633  INTEGER :: damax_iso_not0, damax_iso_not0_local, handle, i1, i2, iat, iatom, icg, ikind, &
634  ipgf1, ipgf2, ir, iset1, iset2, iso, iso1, iso1_first, iso1_last, iso2, iso2_first, &
635  iso2_last, ispin, l, l_iso, llmax, lmax12, lmax_expansion, lmin12, m1s, m2s, m_iso, &
636  max_iso_not0, max_iso_not0_local, max_max_iso_not0, max_nso, max_s_harm, maxl, maxlgto, &
637  maxso, mepos, n1s, n2s, na, natom, natom_tot, nkind, nr, nset, nspins, num_pe, size1, &
638  size2
639  INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list, dacg_n_list
640  INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list, dacg_list
641  INTEGER, DIMENSION(2) :: bo
642  INTEGER, DIMENSION(:), POINTER :: atom_list, lmax, lmin, npgf, o2nindex
643  LOGICAL :: paw_atom
644  LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: is_set_to_0
645  REAL(dp) :: hard_radius
646  REAL(dp), ALLOCATABLE, DIMENSION(:) :: g1, g2, gauge_h, gauge_s
647  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: cjc0_h_block, cjc0_s_block, cjc_h_block, &
648  cjc_ii_h_block, cjc_ii_s_block, cjc_iii_h_block, cjc_iii_s_block, cjc_s_block, dgg_1, gg, &
649  gg_lm1
650  REAL(dp), DIMENSION(:, :), POINTER :: coeff, fr_a_h, fr_a_h_ii, fr_a_h_iii, fr_a_s, &
651  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, &
652  fr_h, fr_s, zet
653  REAL(dp), DIMENSION(:, :, :), POINTER :: my_cg
654  REAL(dp), DIMENSION(:, :, :, :), POINTER :: my_cg_dxyz_asym
655  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
656  TYPE(dft_control_type), POINTER :: dft_control
657  TYPE(grid_atom_type), POINTER :: grid_atom
658  TYPE(gto_basis_set_type), POINTER :: basis_1c_set
659  TYPE(harmonics_atom_type), POINTER :: harmonics
660  TYPE(jrho_atom_type), DIMENSION(:), POINTER :: jrho1_atom_set
661  TYPE(jrho_atom_type), POINTER :: jrho1_atom
662  TYPE(mp_para_env_type), POINTER :: para_env
663  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
664 
665  CALL timeset(routinen, handle)
666  !
667  NULLIFY (atomic_kind_set, qs_kind_set, dft_control, para_env, &
668  coeff, fr_h, fr_s, fr_a_h, fr_a_s, fr_a_h_ii, fr_a_s_ii, &
669  fr_a_h_iii, fr_a_s_iii, fr_b_h, fr_b_s, fr_b_h_ii, &
670  fr_b_s_ii, fr_b_h_iii, fr_b_s_iii, jrho1_atom_set, &
671  jrho1_atom)
672  !
673  CALL get_qs_env(qs_env=qs_env, &
674  atomic_kind_set=atomic_kind_set, &
675  qs_kind_set=qs_kind_set, &
676  dft_control=dft_control, &
677  para_env=para_env)
678 
679  CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxlgto=maxlgto)
680 
681  !
682  CALL get_current_env(current_env=current_env, &
683  jrho1_atom_set=jrho1_atom_set)
684  !
685 
686  nkind = SIZE(qs_kind_set)
687  nspins = dft_control%nspins
688  !
689  natom_tot = SIZE(jrho1_atom_set, 1)
690  ALLOCATE (is_set_to_0(natom_tot, nspins))
691  is_set_to_0(:, :) = .false.
692 
693  !
694  DO ikind = 1, nkind
695  NULLIFY (atom_list, grid_atom, harmonics, basis_1c_set, &
696  lmax, lmin, npgf, zet, grid_atom, harmonics, my_cg, my_cg_dxyz_asym)
697  !
698  CALL get_atomic_kind(atomic_kind_set(ikind), &
699  atom_list=atom_list, &
700  natom=natom)
701  CALL get_qs_kind(qs_kind_set(ikind), &
702  grid_atom=grid_atom, &
703  paw_atom=paw_atom, &
704  harmonics=harmonics, &
705  hard_radius=hard_radius, &
706  basis_set=basis_1c_set, &
707  basis_type="GAPW_1C")
708  !
709  ! Quick cycle if needed.
710  IF (.NOT. paw_atom) cycle
711  !
712  CALL get_gto_basis_set(gto_basis_set=basis_1c_set, &
713  lmax=lmax, lmin=lmin, &
714  maxl=maxl, npgf=npgf, &
715  nset=nset, zet=zet, &
716  maxso=maxso)
717  CALL get_paw_basis_info(basis_1c_set, o2nindex=o2nindex)
718  !
719  nr = grid_atom%nr
720  na = grid_atom%ng_sphere
721  max_iso_not0 = harmonics%max_iso_not0
722  damax_iso_not0 = harmonics%damax_iso_not0
723  max_max_iso_not0 = max(max_iso_not0, damax_iso_not0)
724  lmax_expansion = indso(1, max_max_iso_not0)
725  max_s_harm = harmonics%max_s_harm
726  llmax = harmonics%llmax
727  !
728  ! Distribute the atoms of this kind
729  num_pe = para_env%num_pe
730  mepos = para_env%mepos
731  bo = get_limit(natom, num_pe, mepos)
732  !
733  my_cg => harmonics%my_CG
734  my_cg_dxyz_asym => harmonics%my_CG_dxyz_asym
735  !
736  ! Allocate some arrays.
737  max_nso = nsoset(maxl)
738  ALLOCATE (g1(nr), g2(nr), gg(nr, 0:2*maxl), gg_lm1(nr, 0:2*maxl), dgg_1(nr, 0:2*maxl), &
739  cjc0_h_block(max_nso, max_nso), cjc0_s_block(max_nso, max_nso), &
740  cjc_h_block(max_nso, max_nso), cjc_s_block(max_nso, max_nso), &
741  cjc_ii_h_block(max_nso, max_nso), cjc_ii_s_block(max_nso, max_nso), &
742  cjc_iii_h_block(max_nso, max_nso), cjc_iii_s_block(max_nso, max_nso), &
743  cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm), &
744  dacg_list(2, nsoset(maxl)**2, max_s_harm), dacg_n_list(max_s_harm), &
745  gauge_h(nr), gauge_s(nr))
746  !
747  ! Compute the gauge
748  SELECT CASE (current_env%gauge)
749  CASE (current_gauge_r)
750  ! d(r)=r
751  gauge_h(1:nr) = grid_atom%rad(1:nr)
752  gauge_s(1:nr) = grid_atom%rad(1:nr)
753  CASE (current_gauge_r_and_step_func)
754  ! step function
755  gauge_h(1:nr) = 0e0_dp
756  DO ir = 1, nr
757  IF (grid_atom%rad(ir) .LE. hard_radius) THEN
758  gauge_s(ir) = grid_atom%rad(ir)
759  ELSE
760  gauge_s(ir) = gauge_h(ir)
761  END IF
762  END DO
763  CASE (current_gauge_atom)
764  ! d(r)=A
765  gauge_h(1:nr) = huge(0e0_dp) !0e0_dp
766  gauge_s(1:nr) = huge(0e0_dp) !0e0_dp
767  CASE DEFAULT
768  cpabort("Unknown gauge, try again...")
769  END SELECT
770  !
771  !
772  m1s = 0
773  DO iset1 = 1, nset
774  m2s = 0
775  DO iset2 = 1, nset
776  CALL get_none0_cg_list(my_cg, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
777  max_s_harm, lmax_expansion, cg_list, cg_n_list, max_iso_not0_local)
778  cpassert(max_iso_not0_local .LE. max_iso_not0)
779  CALL get_none0_cg_list(my_cg_dxyz_asym, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
780  max_s_harm, lmax_expansion, dacg_list, dacg_n_list, damax_iso_not0_local)
781  cpassert(damax_iso_not0_local .LE. damax_iso_not0)
782 
783  n1s = nsoset(lmax(iset1))
784  DO ipgf1 = 1, npgf(iset1)
785  iso1_first = nsoset(lmin(iset1) - 1) + 1 + n1s*(ipgf1 - 1) + m1s
786  iso1_last = nsoset(lmax(iset1)) + n1s*(ipgf1 - 1) + m1s
787  size1 = iso1_last - iso1_first + 1
788  iso1_first = o2nindex(iso1_first)
789  iso1_last = o2nindex(iso1_last)
790  i1 = iso1_last - iso1_first + 1
791  cpassert(size1 == i1)
792  i1 = nsoset(lmin(iset1) - 1) + 1
793  !
794  g1(1:nr) = exp(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))
795  !
796  n2s = nsoset(lmax(iset2))
797  DO ipgf2 = 1, npgf(iset2)
798  iso2_first = nsoset(lmin(iset2) - 1) + 1 + n2s*(ipgf2 - 1) + m2s
799  iso2_last = nsoset(lmax(iset2)) + n2s*(ipgf2 - 1) + m2s
800  size2 = iso2_last - iso2_first + 1
801  iso2_first = o2nindex(iso2_first)
802  iso2_last = o2nindex(iso2_last)
803  i2 = iso2_last - iso2_first + 1
804  cpassert(size2 == i2)
805  i2 = nsoset(lmin(iset2) - 1) + 1
806  !
807  g2(1:nr) = exp(-zet(ipgf2, iset2)*grid_atom%rad2(1:nr))
808  !
809  lmin12 = lmin(iset1) + lmin(iset2)
810  lmax12 = lmax(iset1) + lmax(iset2)
811  !
812  gg = 0.0_dp
813  gg_lm1 = 0.0_dp
814  dgg_1 = 0.0_dp
815  !
816  ! Take only the terms of expansion for L < lmax_expansion
817  IF (lmin12 .LE. lmax_expansion) THEN
818  !
819  IF (lmax12 .GT. lmax_expansion) lmax12 = lmax_expansion
820  !
821  IF (lmin12 == 0) THEN
822  gg(1:nr, lmin12) = g1(1:nr)*g2(1:nr)
823  gg_lm1(1:nr, lmin12) = 0.0_dp
824  ELSE
825  gg(1:nr, lmin12) = grid_atom%rad2l(1:nr, lmin12)*g1(1:nr)*g2(1:nr)
826  gg_lm1(1:nr, lmin12) = grid_atom%rad2l(1:nr, lmin12 - 1)*g1(1:nr)*g2(1:nr)
827  END IF
828  !
829  DO l = lmin12 + 1, lmax12
830  gg(1:nr, l) = grid_atom%rad(1:nr)*gg(1:nr, l - 1)
831  gg_lm1(1:nr, l) = gg(1:nr, l - 1)
832  END DO
833  !
834  DO l = lmin12, lmax12
835  dgg_1(1:nr, l) = 2.0_dp*(zet(ipgf1, iset1) - zet(ipgf2, iset2))&
836  & *gg(1:nr, l)*grid_atom%rad(1:nr)
837  END DO
838  ELSE
839  cycle
840  END IF ! lmin12
841  !
842  DO iat = bo(1), bo(2)
843  iatom = atom_list(iat)
844  !
845  DO ispin = 1, nspins
846  !------------------------------------------------------------------
847  ! P_\alpha\alpha'
848  cjc0_h_block = huge(1.0_dp)
849  cjc0_s_block = huge(1.0_dp)
850  !
851  ! Hard term
852  coeff => jrho1_atom_set(iatom)%cjc0_h(ispin)%r_coef
853  cjc0_h_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = &
854  coeff(iso1_first:iso1_last, iso2_first:iso2_last)
855  !
856  ! Soft term
857  coeff => jrho1_atom_set(iatom)%cjc0_s(ispin)%r_coef
858  cjc0_s_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = &
859  coeff(iso1_first:iso1_last, iso2_first:iso2_last)
860  !------------------------------------------------------------------
861  ! mQai_\alpha\alpha'
862  cjc_h_block = huge(1.0_dp)
863  cjc_s_block = huge(1.0_dp)
864  !
865  ! Hard term
866  coeff => jrho1_atom_set(iatom)%cjc_h(ispin)%r_coef
867  cjc_h_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = &
868  coeff(iso1_first:iso1_last, iso2_first:iso2_last)
869  !
870  ! Soft term
871  coeff => jrho1_atom_set(iatom)%cjc_s(ispin)%r_coef
872  cjc_s_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = &
873  coeff(iso1_first:iso1_last, iso2_first:iso2_last)
874  !------------------------------------------------------------------
875  ! Qci_\alpha\alpha'
876  cjc_ii_h_block = huge(1.0_dp)
877  cjc_ii_s_block = huge(1.0_dp)
878  !
879  ! Hard term
880  coeff => jrho1_atom_set(iatom)%cjc_ii_h(ispin)%r_coef
881  cjc_ii_h_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = &
882  coeff(iso1_first:iso1_last, iso2_first:iso2_last)
883  !
884  ! Soft term
885  coeff => jrho1_atom_set(iatom)%cjc_ii_s(ispin)%r_coef
886  cjc_ii_s_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = &
887  coeff(iso1_first:iso1_last, iso2_first:iso2_last)
888  !------------------------------------------------------------------
889  ! Qbi_\alpha\alpha'
890  cjc_iii_h_block = huge(1.0_dp)
891  cjc_iii_s_block = huge(1.0_dp)
892  !
893  !
894  ! Hard term
895  coeff => jrho1_atom_set(iatom)%cjc_iii_h(ispin)%r_coef
896  cjc_iii_h_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = &
897  coeff(iso1_first:iso1_last, iso2_first:iso2_last)
898  !
899  ! Soft term
900  coeff => jrho1_atom_set(iatom)%cjc_iii_s(ispin)%r_coef
901  cjc_iii_s_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = &
902  coeff(iso1_first:iso1_last, iso2_first:iso2_last)
903  !------------------------------------------------------------------
904  !
905  ! Allocation radial functions
906  jrho1_atom => jrho1_atom_set(iatom)
907  IF (.NOT. ASSOCIATED(jrho1_atom%jrho_a_h(ispin)%r_coef)) THEN
908  CALL allocate_jrho_atom_rad(jrho1_atom, ispin, nr, na, &
909  max_max_iso_not0)
910  is_set_to_0(iatom, ispin) = .true.
911  ELSE
912  IF (.NOT. is_set_to_0(iatom, ispin)) THEN
913  CALL set2zero_jrho_atom_rad(jrho1_atom, ispin)
914  is_set_to_0(iatom, ispin) = .true.
915  END IF
916  END IF
917  !------------------------------------------------------------------
918  !
919  fr_h => jrho1_atom%jrho_h(ispin)%r_coef
920  fr_s => jrho1_atom%jrho_s(ispin)%r_coef
921  !------------------------------------------------------------------
922  !
923  fr_a_h => jrho1_atom%jrho_a_h(ispin)%r_coef
924  fr_a_s => jrho1_atom%jrho_a_s(ispin)%r_coef
925  fr_b_h => jrho1_atom%jrho_b_h(ispin)%r_coef
926  fr_b_s => jrho1_atom%jrho_b_s(ispin)%r_coef
927  !------------------------------------------------------------------
928  !
929  fr_a_h_ii => jrho1_atom%jrho_a_h_ii(ispin)%r_coef
930  fr_a_s_ii => jrho1_atom%jrho_a_s_ii(ispin)%r_coef
931  fr_b_h_ii => jrho1_atom%jrho_b_h_ii(ispin)%r_coef
932  fr_b_s_ii => jrho1_atom%jrho_b_s_ii(ispin)%r_coef
933  !------------------------------------------------------------------
934  !
935  fr_a_h_iii => jrho1_atom%jrho_a_h_iii(ispin)%r_coef
936  fr_a_s_iii => jrho1_atom%jrho_a_s_iii(ispin)%r_coef
937  fr_b_h_iii => jrho1_atom%jrho_b_h_iii(ispin)%r_coef
938  fr_b_s_iii => jrho1_atom%jrho_b_s_iii(ispin)%r_coef
939  !------------------------------------------------------------------
940  !
941  DO iso = 1, max_iso_not0_local
942  l_iso = indso(1, iso) ! not needed
943  m_iso = indso(2, iso) ! not needed
944  !
945  DO icg = 1, cg_n_list(iso)
946  !
947  iso1 = cg_list(1, icg, iso)
948  iso2 = cg_list(2, icg, iso)
949  !
950  IF (.NOT. (iso2 > 0 .AND. iso1 > 0)) THEN
951  WRITE (*, *) 'iso1=', iso1, ' iso2=', iso2, ' iso=', iso, ' icg=', icg
952  WRITE (*, *) '.... will stop!'
953  END IF
954  cpassert(iso2 > 0 .AND. iso1 > 0)
955  !
956  l = indso(1, iso1) + indso(1, iso2)
957  IF (l .GT. lmax_expansion .OR. l .LT. .0) THEN
958  WRITE (*, *) 'calculate_jrho_atom_rad: 1 l', l
959  WRITE (*, *) 'calculate_jrho_atom_rad: 1 lmax_expansion', lmax_expansion
960  WRITE (*, *) '.... will stop!'
961  END IF
962  cpassert(l <= lmax_expansion)
963  !------------------------------------------------------------------
964  ! P0
965  !
966  IF (current_env%gauge .EQ. current_gauge_atom) THEN
967  ! Hard term
968  fr_h(1:nr, iso) = fr_h(1:nr, iso) + &
969  gg(1:nr, l)*cjc0_h_block(iso1, iso2)* &
970  my_cg(iso1, iso2, iso)
971  ! Soft term
972  fr_s(1:nr, iso) = fr_s(1:nr, iso) + &
973  gg(1:nr, l)*cjc0_s_block(iso1, iso2)* &
974  my_cg(iso1, iso2, iso)
975  ELSE
976  ! Hard term
977  fr_h(1:nr, iso) = fr_h(1:nr, iso) + &
978  gg(1:nr, l)*cjc0_h_block(iso1, iso2)* &
979  my_cg(iso1, iso2, iso)*(grid_atom%rad(1:nr) - gauge_h(1:nr))
980  ! Soft term
981  fr_s(1:nr, iso) = fr_s(1:nr, iso) + &
982  gg(1:nr, l)*cjc0_s_block(iso1, iso2)* &
983  my_cg(iso1, iso2, iso)*(grid_atom%rad(1:nr) - gauge_s(1:nr))
984  END IF
985  !------------------------------------------------------------------
986  ! Rai
987  !
988  ! Hard term
989  fr_a_h(1:nr, iso) = fr_a_h(1:nr, iso) + &
990  dgg_1(1:nr, l)*cjc_h_block(iso1, iso2)* &
991  my_cg(iso1, iso2, iso)
992  !
993  ! Soft term
994  fr_a_s(1:nr, iso) = fr_a_s(1:nr, iso) + &
995  dgg_1(1:nr, l)*cjc_s_block(iso1, iso2)* &
996  my_cg(iso1, iso2, iso)
997  !------------------------------------------------------------------
998  ! Rci
999  !
1000  IF (current_env%gauge .EQ. current_gauge_atom) THEN
1001  ! Hard term
1002  fr_a_h_ii(1:nr, iso) = fr_a_h_ii(1:nr, iso) + &
1003  dgg_1(1:nr, l)* &
1004  cjc_ii_h_block(iso1, iso2)* &
1005  my_cg(iso1, iso2, iso)
1006  ! Soft term
1007  fr_a_s_ii(1:nr, iso) = fr_a_s_ii(1:nr, iso) + &
1008  dgg_1(1:nr, l)* &
1009  cjc_ii_s_block(iso1, iso2)* &
1010  my_cg(iso1, iso2, iso)
1011  ELSE
1012  ! Hard term
1013  fr_a_h_ii(1:nr, iso) = fr_a_h_ii(1:nr, iso) + &
1014  dgg_1(1:nr, l)*gauge_h(1:nr)* &
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)*gauge_s(1:nr)* &
1020  cjc_ii_s_block(iso1, iso2)* &
1021  my_cg(iso1, iso2, iso)
1022  END IF
1023  !------------------------------------------------------------------
1024  ! Rbi
1025  !
1026  IF (current_env%gauge .EQ. current_gauge_atom) THEN
1027  ! Hard term
1028  fr_a_h_iii(1:nr, iso) = fr_a_h_iii(1:nr, iso) + &
1029  dgg_1(1:nr, l)* &
1030  cjc_iii_h_block(iso1, iso2)* &
1031  my_cg(iso1, iso2, iso)
1032  ! Soft term
1033  fr_a_s_iii(1:nr, iso) = fr_a_s_iii(1:nr, iso) + &
1034  dgg_1(1:nr, l)* &
1035  cjc_iii_s_block(iso1, iso2)* &
1036  my_cg(iso1, iso2, iso)
1037  ELSE
1038  ! Hard term
1039  fr_a_h_iii(1:nr, iso) = fr_a_h_iii(1:nr, iso) + &
1040  dgg_1(1:nr, l)*gauge_h(1:nr)* &
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)*gauge_s(1:nr)* &
1046  cjc_iii_s_block(iso1, iso2)* &
1047  my_cg(iso1, iso2, iso)
1048  END IF
1049  !------------------------------------------------------------------
1050  END DO !icg
1051  !
1052  END DO ! iso
1053  !
1054  !
1055  DO iso = 1, damax_iso_not0_local
1056  l_iso = indso(1, iso) ! not needed
1057  m_iso = indso(2, iso) ! not needed
1058  !
1059  DO icg = 1, dacg_n_list(iso)
1060  !
1061  iso1 = dacg_list(1, icg, iso)
1062  iso2 = dacg_list(2, icg, iso)
1063  !
1064  IF (.NOT. (iso2 > 0 .AND. iso1 > 0)) THEN
1065  WRITE (*, *) 'iso1=', iso1, ' iso2=', iso2, ' iso=', iso, ' icg=', icg
1066  WRITE (*, *) '.... will stop!'
1067  END IF
1068  cpassert(iso2 > 0 .AND. iso1 > 0)
1069  !
1070  l = indso(1, iso1) + indso(1, iso2)
1071  IF (l .GT. lmax_expansion) THEN
1072  WRITE (*, *) 'calculate_jrho_atom_rad: 1 l', l
1073  WRITE (*, *) 'calculate_jrho_atom_rad: 1 lmax_expansion', lmax_expansion
1074  WRITE (*, *) '.... will stop!'
1075  END IF
1076  cpassert(l <= lmax_expansion)
1077  !------------------------------------------------------------------
1078  ! Daij
1079  !
1080  ! Hard term
1081  fr_b_h(1:nr, iso) = fr_b_h(1:nr, iso) + &
1082  gg_lm1(1:nr, l)*cjc_h_block(iso1, iso2)* &
1083  my_cg_dxyz_asym(idir, iso1, iso2, iso)
1084  !
1085  ! Soft term
1086  fr_b_s(1:nr, iso) = fr_b_s(1:nr, iso) + &
1087  gg_lm1(1:nr, l)*cjc_s_block(iso1, iso2)* &
1088  my_cg_dxyz_asym(idir, iso1, iso2, iso)
1089  !
1090  !------------------------------------------------------------------
1091  ! Dcij
1092  !
1093  IF (current_env%gauge .EQ. current_gauge_atom) THEN
1094  ! Hard term
1095  fr_b_h_ii(1:nr, iso) = fr_b_h_ii(1:nr, iso) + &
1096  gg_lm1(1:nr, l)* &
1097  cjc_ii_h_block(iso1, iso2)* &
1098  my_cg_dxyz_asym(idir, iso1, iso2, iso)
1099  ! Soft term
1100  fr_b_s_ii(1:nr, iso) = fr_b_s_ii(1:nr, iso) + &
1101  gg_lm1(1:nr, l)* &
1102  cjc_ii_s_block(iso1, iso2)* &
1103  my_cg_dxyz_asym(idir, iso1, iso2, iso)
1104  ELSE
1105  ! Hard term
1106  fr_b_h_ii(1:nr, iso) = fr_b_h_ii(1:nr, iso) + &
1107  gg_lm1(1:nr, l)*gauge_h(1:nr)* &
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)*gauge_s(1:nr)* &
1113  cjc_ii_s_block(iso1, iso2)* &
1114  my_cg_dxyz_asym(idir, iso1, iso2, iso)
1115  END IF
1116  !------------------------------------------------------------------
1117  ! Dbij
1118  !
1119  IF (current_env%gauge .EQ. current_gauge_atom) THEN
1120  ! Hard term
1121  fr_b_h_iii(1:nr, iso) = fr_b_h_iii(1:nr, iso) + &
1122  gg_lm1(1:nr, l)* &
1123  cjc_iii_h_block(iso1, iso2)* &
1124  my_cg_dxyz_asym(idir, iso1, iso2, iso)
1125  ! Soft term
1126  fr_b_s_iii(1:nr, iso) = fr_b_s_iii(1:nr, iso) + &
1127  gg_lm1(1:nr, l)* &
1128  cjc_iii_s_block(iso1, iso2)* &
1129  my_cg_dxyz_asym(idir, iso1, iso2, iso)
1130  ELSE
1131  ! Hard term
1132  fr_b_h_iii(1:nr, iso) = fr_b_h_iii(1:nr, iso) + &
1133  gg_lm1(1:nr, l)*gauge_h(1:nr)* &
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)*gauge_s(1:nr)* &
1139  cjc_iii_s_block(iso1, iso2)* &
1140  my_cg_dxyz_asym(idir, iso1, iso2, iso)
1141  END IF
1142  !------------------------------------------------------------------
1143  END DO ! icg
1144  END DO ! iso
1145  !
1146  END DO ! ispin
1147  END DO ! iat
1148  !
1149  !------------------------------------------------------------------
1150  !
1151  END DO !ipgf2
1152  END DO ! ipgf1
1153  m2s = m2s + maxso
1154  END DO ! iset2
1155  m1s = m1s + maxso
1156  END DO ! iset1
1157  !
1158  DEALLOCATE (cjc0_h_block, cjc0_s_block, cjc_h_block, cjc_s_block, cjc_ii_h_block, cjc_ii_s_block, &
1159  cjc_iii_h_block, cjc_iii_s_block, g1, g2, gg, gg_lm1, dgg_1, gauge_h, gauge_s, &
1160  cg_list, cg_n_list, dacg_list, dacg_n_list)
1161  DEALLOCATE (o2nindex)
1162  END DO ! ikind
1163  !
1164  !
1165  DEALLOCATE (is_set_to_0)
1166  !
1167  CALL timestop(handle)
1168  !
1169  END SUBROUTINE calculate_jrho_atom_rad
1170 
1171 ! **************************************************************************************************
1172 !> \brief ...
1173 !> \param jrho1_atom ...
1174 !> \param jrho_h ...
1175 !> \param jrho_s ...
1176 !> \param grid_atom ...
1177 !> \param harmonics ...
1178 !> \param do_igaim ...
1179 !> \param ratom ...
1180 !> \param natm_gauge ...
1181 !> \param iB ...
1182 !> \param idir ...
1183 !> \param ispin ...
1184 ! **************************************************************************************************
1185  SUBROUTINE calculate_jrho_atom_ang(jrho1_atom, jrho_h, jrho_s, grid_atom, &
1186  harmonics, do_igaim, ratom, natm_gauge, &
1187  iB, idir, ispin)
1188  !
1189  TYPE(jrho_atom_type), POINTER :: jrho1_atom
1190  REAL(dp), DIMENSION(:, :), POINTER :: jrho_h, jrho_s
1191  TYPE(grid_atom_type), POINTER :: grid_atom
1192  TYPE(harmonics_atom_type), POINTER :: harmonics
1193  LOGICAL, INTENT(IN) :: do_igaim
1194  INTEGER, INTENT(IN) :: ib, idir, ispin, natm_gauge
1195  REAL(dp), INTENT(IN) :: ratom(:, :)
1196 
1197  INTEGER :: ia, idir2, iib, iiib, ir, &
1198  iso, max_max_iso_not0, na, nr
1199  REAL(dp) :: rad_part, scale
1200  REAL(dp), DIMENSION(:, :), POINTER :: a, fr_a_h, fr_a_h_ii, fr_a_h_iii, &
1201  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, &
1202  fr_b_s_ii, fr_b_s_iii, fr_h, fr_s, slm
1203  REAL(dp), DIMENSION(:), POINTER :: r
1204  REAL(dp), DIMENSION(:, :, :), ALLOCATABLE :: g
1205 !
1206 !
1207  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, &
1208  fr_b_h, fr_b_s, fr_b_h_ii, fr_b_s_ii, fr_b_h_iii, fr_b_s_iii, &
1209  a, slm)
1210  !
1211  cpassert(ASSOCIATED(jrho_h))
1212  cpassert(ASSOCIATED(jrho_s))
1213  cpassert(ASSOCIATED(jrho1_atom))
1214  ! just to be sure...
1215  cpassert(ASSOCIATED(jrho1_atom%jrho_a_h))
1216  cpassert(ASSOCIATED(jrho1_atom%jrho_a_s))
1217  cpassert(ASSOCIATED(jrho1_atom%jrho_b_h))
1218  cpassert(ASSOCIATED(jrho1_atom%jrho_b_s))
1219  cpassert(ASSOCIATED(jrho1_atom%jrho_a_h(ispin)%r_coef))
1220  cpassert(ASSOCIATED(jrho1_atom%jrho_a_s(ispin)%r_coef))
1221  cpassert(ASSOCIATED(jrho1_atom%jrho_b_h(ispin)%r_coef))
1222  cpassert(ASSOCIATED(jrho1_atom%jrho_b_s(ispin)%r_coef))
1223  cpassert(ASSOCIATED(jrho1_atom%jrho_a_h_ii))
1224  cpassert(ASSOCIATED(jrho1_atom%jrho_a_s_ii))
1225  cpassert(ASSOCIATED(jrho1_atom%jrho_b_h_ii))
1226  cpassert(ASSOCIATED(jrho1_atom%jrho_b_s_ii))
1227  cpassert(ASSOCIATED(jrho1_atom%jrho_a_h_ii(ispin)%r_coef))
1228  cpassert(ASSOCIATED(jrho1_atom%jrho_a_s_ii(ispin)%r_coef))
1229  cpassert(ASSOCIATED(jrho1_atom%jrho_b_h_ii(ispin)%r_coef))
1230  cpassert(ASSOCIATED(jrho1_atom%jrho_b_s_ii(ispin)%r_coef))
1231  cpassert(ASSOCIATED(jrho1_atom%jrho_a_h_iii))
1232  cpassert(ASSOCIATED(jrho1_atom%jrho_a_s_iii))
1233  cpassert(ASSOCIATED(jrho1_atom%jrho_b_h_iii))
1234  cpassert(ASSOCIATED(jrho1_atom%jrho_b_s_iii))
1235  cpassert(ASSOCIATED(jrho1_atom%jrho_a_h_iii(ispin)%r_coef))
1236  cpassert(ASSOCIATED(jrho1_atom%jrho_a_s_iii(ispin)%r_coef))
1237  cpassert(ASSOCIATED(jrho1_atom%jrho_b_h_iii(ispin)%r_coef))
1238  cpassert(ASSOCIATED(jrho1_atom%jrho_b_s_iii(ispin)%r_coef))
1239  !
1240  !
1241  nr = grid_atom%nr
1242  na = grid_atom%ng_sphere
1243  max_max_iso_not0 = max(harmonics%max_iso_not0, harmonics%damax_iso_not0)
1244  ALLOCATE (g(3, nr, na))
1245  !------------------------------------------------------------------
1246  !
1247  fr_h => jrho1_atom%jrho_h(ispin)%r_coef
1248  fr_s => jrho1_atom%jrho_s(ispin)%r_coef
1249  !------------------------------------------------------------------
1250  !
1251  fr_a_h => jrho1_atom%jrho_a_h(ispin)%r_coef !Rai
1252  fr_a_s => jrho1_atom%jrho_a_s(ispin)%r_coef
1253  fr_b_h => jrho1_atom%jrho_b_h(ispin)%r_coef !Daij
1254  fr_b_s => jrho1_atom%jrho_b_s(ispin)%r_coef
1255  !------------------------------------------------------------------
1256  !
1257  fr_a_h_ii => jrho1_atom%jrho_a_h_ii(ispin)%r_coef !Rci
1258  fr_a_s_ii => jrho1_atom%jrho_a_s_ii(ispin)%r_coef
1259  fr_b_h_ii => jrho1_atom%jrho_b_h_ii(ispin)%r_coef !Dcij
1260  fr_b_s_ii => jrho1_atom%jrho_b_s_ii(ispin)%r_coef
1261  !------------------------------------------------------------------
1262  !
1263  fr_a_h_iii => jrho1_atom%jrho_a_h_iii(ispin)%r_coef !Rbi
1264  fr_a_s_iii => jrho1_atom%jrho_a_s_iii(ispin)%r_coef
1265  fr_b_h_iii => jrho1_atom%jrho_b_h_iii(ispin)%r_coef !Dbij
1266  fr_b_s_iii => jrho1_atom%jrho_b_s_iii(ispin)%r_coef
1267  !------------------------------------------------------------------
1268  !
1269  a => harmonics%a
1270  slm => harmonics%slm
1271  r => grid_atom%rad
1272  !
1273  CALL set_vecp(ib, iib, iiib)
1274  !
1275  scale = 0.0_dp
1276  idir2 = 1
1277  IF (idir .NE. ib) THEN
1278  CALL set_vecp_rev(idir, ib, idir2)
1279  scale = fac_vecp(idir, ib, idir2)
1280  END IF
1281  !
1282  ! Set the gauge
1283  CALL get_gauge()
1284  !
1285  DO ir = 1, nr
1286  DO iso = 1, max_max_iso_not0
1287  DO ia = 1, na
1288  IF (do_igaim) THEN
1289  !------------------------------------------------------------------
1290  ! Hard current density response
1291  ! radial(ia,ir) = ( aj(ia) * Rai(ir,iso) + Daij
1292  ! - aii(ia) * ( aj(ia) * Rbi(ir,iso) + Dbij )
1293  ! + aiii(ia) * ( aj(ia) * Rci(ir,iso) + Dcij )
1294  ! ) * Ylm(ia)
1295  rad_part = a(idir, ia)*fr_a_h(ir, iso) + fr_b_h(ir, iso) &
1296  & - g(iib, ir, ia)*(a(idir, ia)*fr_a_h_iii(ir, iso) + fr_b_h_iii(ir, iso))&
1297  & + g(iiib, ir, ia)*(a(idir, ia)*fr_a_h_ii(ir, iso) + fr_b_h_ii(ir, iso))&
1298  & + scale*(a(idir2, ia)*r(ir) - g(idir2, ir, ia))*fr_h(ir, iso)
1299  !
1300  jrho_h(ir, ia) = jrho_h(ir, ia) + rad_part*slm(ia, iso)
1301  !------------------------------------------------------------------
1302  ! Soft current density response
1303  rad_part = a(idir, ia)*fr_a_s(ir, iso) + fr_b_s(ir, iso) &
1304  & - g(iib, ir, ia)*(a(idir, ia)*fr_a_s_iii(ir, iso) + fr_b_s_iii(ir, iso))&
1305  & + g(iiib, ir, ia)*(a(idir, ia)*fr_a_s_ii(ir, iso) + fr_b_s_ii(ir, iso))&
1306  & + scale*(a(idir2, ia)*r(ir) - g(idir2, ir, ia))*fr_s(ir, iso)
1307  !
1308  jrho_s(ir, ia) = jrho_s(ir, ia) + rad_part*slm(ia, iso)
1309  !------------------------------------------------------------------
1310  ELSE
1311  !------------------------------------------------------------------
1312  ! Hard current density response
1313  ! radial(ia,ir) = ( aj(ia) * Rai(ir,iso) + Daij
1314  ! - aii(ia) * ( aj(ia) * Rbi(ir,iso) + Dbij )
1315  ! + aiii(ia) * ( aj(ia) * Rci(ir,iso) + Dcij )
1316  ! ) * Ylm(ia)
1317  rad_part = a(idir, ia)*fr_a_h(ir, iso) + fr_b_h(ir, iso) &
1318  & - a(iib, ia)*(a(idir, ia)*fr_a_h_iii(ir, iso) + fr_b_h_iii(ir, iso))&
1319  & + a(iiib, ia)*(a(idir, ia)*fr_a_h_ii(ir, iso) + fr_b_h_ii(ir, iso))&
1320  & + scale*a(idir2, ia)*fr_h(ir, iso)
1321  !
1322  jrho_h(ir, ia) = jrho_h(ir, ia) + rad_part*slm(ia, iso)
1323  !------------------------------------------------------------------
1324  ! Soft current density response
1325  rad_part = a(idir, ia)*fr_a_s(ir, iso) + fr_b_s(ir, iso) &
1326  & - a(iib, ia)*(a(idir, ia)*fr_a_s_iii(ir, iso) + fr_b_s_iii(ir, iso))&
1327  & + a(iiib, ia)*(a(idir, ia)*fr_a_s_ii(ir, iso) + fr_b_s_ii(ir, iso))&
1328  & + scale*a(idir2, ia)*fr_s(ir, iso)
1329  !
1330  jrho_s(ir, ia) = jrho_s(ir, ia) + rad_part*slm(ia, iso)
1331  !------------------------------------------------------------------
1332  END IF
1333  END DO ! ia
1334  END DO ! iso
1335  END DO ! ir
1336  !
1337  DEALLOCATE (g)
1338  !
1339  CONTAINS
1340  !
1341 ! **************************************************************************************************
1342 !> \brief ...
1343 ! **************************************************************************************************
1344  SUBROUTINE get_gauge()
1345  INTEGER :: iatom, ixyz, jatom
1346  REAL(dp) :: ab, pa, pb, point(3), pra(3), prb(3), &
1347  res, tmp
1348  REAL(dp), ALLOCATABLE, DIMENSION(:) :: buf
1349 
1350  ALLOCATE (buf(natm_gauge))
1351  DO ir = 1, nr
1352  DO ia = 1, na
1353  DO ixyz = 1, 3
1354  g(ixyz, ir, ia) = 0.0_dp
1355  END DO
1356  point(1) = r(ir)*a(1, ia)
1357  point(2) = r(ir)*a(2, ia)
1358  point(3) = r(ir)*a(3, ia)
1359  DO iatom = 1, natm_gauge
1360  buf(iatom) = 1.0_dp
1361  pra = point - ratom(:, iatom)
1362  pa = sqrt(pra(1)**2 + pra(2)**2 + pra(3)**2)
1363  DO jatom = 1, natm_gauge
1364  IF (iatom .EQ. jatom) cycle
1365  prb = point - ratom(:, jatom)
1366  pb = sqrt(prb(1)**2 + prb(2)**2 + prb(3)**2)
1367  ab = sqrt((pra(1) - prb(1))**2 + (pra(2) - prb(2))**2 + (pra(3) - prb(3))**2)
1368  !
1369  tmp = (pa - pb)/ab
1370  tmp = 0.5_dp*(3.0_dp - tmp**2)*tmp
1371  tmp = 0.5_dp*(3.0_dp - tmp**2)*tmp
1372  tmp = 0.5_dp*(3.0_dp - tmp**2)*tmp
1373  buf(iatom) = buf(iatom)*0.5_dp*(1.0_dp - tmp)
1374  END DO
1375  END DO
1376  DO ixyz = 1, 3
1377  res = 0.0_dp
1378  DO iatom = 1, natm_gauge
1379  res = res + ratom(ixyz, iatom)*buf(iatom)
1380  END DO
1381  res = res/sum(buf(1:natm_gauge))
1382  !
1383  g(ixyz, ir, ia) = res
1384  END DO
1385  END DO
1386  END DO
1387  DEALLOCATE (buf)
1388  END SUBROUTINE get_gauge
1389  END SUBROUTINE calculate_jrho_atom_ang
1390 
1391 ! **************************************************************************************************
1392 !> \brief ...
1393 !> \param current_env ...
1394 !> \param qs_env ...
1395 !> \param iB ...
1396 !> \param idir ...
1397 ! **************************************************************************************************
1398  SUBROUTINE calculate_jrho_atom(current_env, qs_env, iB, idir)
1399  TYPE(current_env_type) :: current_env
1400  TYPE(qs_environment_type), POINTER :: qs_env
1401  INTEGER, INTENT(IN) :: ib, idir
1402 
1403  INTEGER :: iat, iatom, ikind, ispin, jatom, &
1404  natm_gauge, natm_tot, natom, nkind, &
1405  nspins
1406  INTEGER, DIMENSION(2) :: bo
1407  INTEGER, DIMENSION(:), POINTER :: atom_list
1408  LOGICAL :: do_igaim, gapw, paw_atom
1409  REAL(dp) :: hard_radius, r(3)
1410  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: ratom
1411  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1412  TYPE(cell_type), POINTER :: cell
1413  TYPE(mp_para_env_type), POINTER :: para_env
1414  TYPE(dft_control_type), POINTER :: dft_control
1415  TYPE(grid_atom_type), POINTER :: grid_atom
1416  TYPE(harmonics_atom_type), POINTER :: harmonics
1417  TYPE(jrho_atom_type), DIMENSION(:), POINTER :: jrho1_atom_set
1418  TYPE(jrho_atom_type), POINTER :: jrho1_atom
1419  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1420  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1421 
1422  NULLIFY (para_env, dft_control)
1423  NULLIFY (jrho1_atom_set, grid_atom, harmonics)
1424  NULLIFY (atomic_kind_set, qs_kind_set, atom_list)
1425 
1426  CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, &
1427  atomic_kind_set=atomic_kind_set, &
1428  qs_kind_set=qs_kind_set, &
1429  particle_set=particle_set, &
1430  cell=cell, &
1431  para_env=para_env)
1432 
1433  CALL get_current_env(current_env=current_env, &
1434  jrho1_atom_set=jrho1_atom_set)
1435 
1436  do_igaim = .false.
1437  IF (current_env%gauge .EQ. current_gauge_atom) do_igaim = .true.
1438 
1439  gapw = dft_control%qs_control%gapw
1440  nkind = SIZE(qs_kind_set, 1)
1441  nspins = dft_control%nspins
1442 
1443  natm_tot = SIZE(particle_set)
1444  ALLOCATE (ratom(3, natm_tot))
1445 
1446  IF (gapw) THEN
1447  DO ikind = 1, nkind
1448  NULLIFY (atom_list, grid_atom, harmonics)
1449  CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
1450  CALL get_qs_kind(qs_kind_set(ikind), &
1451  grid_atom=grid_atom, &
1452  harmonics=harmonics, &
1453  hard_radius=hard_radius, &
1454  paw_atom=paw_atom)
1455 
1456  IF (.NOT. paw_atom) cycle
1457 
1458  ! Distribute the atoms of this kind
1459 
1460  bo = get_limit(natom, para_env%num_pe, para_env%mepos)
1461 
1462  DO iat = bo(1), bo(2)
1463  iatom = atom_list(iat)
1464  NULLIFY (jrho1_atom)
1465  jrho1_atom => jrho1_atom_set(iatom)
1466 
1467  natm_gauge = 0
1468  DO jatom = 1, natm_tot
1469  r(:) = pbc(particle_set(jatom)%r(:) - particle_set(iatom)%r(:), cell)
1470  ! SQRT(SUM(r(:)**2)) .LE. 2.0_dp*hard_radius
1471  IF (sum(r(:)**2) .LE. (4.0_dp*hard_radius*hard_radius)) THEN
1472  natm_gauge = natm_gauge + 1
1473  ratom(:, natm_gauge) = r(:)
1474  END IF
1475  END DO
1476 
1477  DO ispin = 1, nspins
1478  jrho1_atom%jrho_vec_rad_h(idir, ispin)%r_coef = 0.0_dp
1479  jrho1_atom%jrho_vec_rad_s(idir, ispin)%r_coef = 0.0_dp
1480  CALL calculate_jrho_atom_ang(jrho1_atom, &
1481  jrho1_atom%jrho_vec_rad_h(idir, ispin)%r_coef, &
1482  jrho1_atom%jrho_vec_rad_s(idir, ispin)%r_coef, &
1483  grid_atom, harmonics, &
1484  do_igaim, &
1485  ratom, natm_gauge, ib, idir, ispin)
1486  END DO !ispin
1487  END DO !iat
1488  END DO !ikind
1489  END IF
1490 
1491  DEALLOCATE (ratom)
1492 
1493  END SUBROUTINE calculate_jrho_atom
1494 
1495 END MODULE qs_linres_atom_current
pbc
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition: atomic_kind_types.F:138
basis_set_types
Definition: basis_set_types.F:15
basis_set_types::get_gto_basis_set
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius)
...
Definition: basis_set_types.F:615
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition: cp_log_handling.F:515
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::current_gauge_atom
integer, parameter, public current_gauge_atom
Definition: input_constants.F:501
input_constants::current_gauge_r
integer, parameter, public current_gauge_r
Definition: input_constants.F:501
input_constants::current_gauge_r_and_step_func
integer, parameter, public current_gauge_r_and_step_func
Definition: input_constants.F:501
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition: orbital_pointers.F:15
orbital_pointers::nsoset
integer, dimension(:), allocatable, public nsoset
Definition: orbital_pointers.F:42
orbital_pointers::indso
integer, dimension(:, :), allocatable, public indso
Definition: orbital_pointers.F:43
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
paw_basis_types
Definition: paw_basis_types.F:13
paw_basis_types::get_paw_basis_info
subroutine, public get_paw_basis_info(basis_1c, o2nindex, n2oindex, nsatbas)
Return some info on the PAW basis derived from a GTO basis set.
Definition: paw_basis_types.F:40
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
qs_grid_atom
Definition: qs_grid_atom.F:8
qs_harmonics_atom
Definition: qs_harmonics_atom.F:8
qs_kind_types
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition: qs_kind_types.F:451
qs_kind_types::get_qs_kind_set
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr)
Get attributes of an atomic kind set.
Definition: qs_kind_types.F:864
qs_linres_atom_current
given the response wavefunctions obtained by the application of the (rxp), p, and ((dk-dl)xp) operato...
Definition: qs_linres_atom_current.F:17
qs_linres_atom_current::calculate_jrho_atom
subroutine, public calculate_jrho_atom(current_env, qs_env, iB, idir)
...
Definition: qs_linres_atom_current.F:1399
qs_linres_atom_current::calculate_jrho_atom_rad
subroutine, public calculate_jrho_atom_rad(qs_env, current_env, idir)
...
Definition: qs_linres_atom_current.F:626
qs_linres_atom_current::calculate_jrho_atom_coeff
subroutine, public calculate_jrho_atom_coeff(qs_env, current_env, mat_d0, mat_jp, mat_jp_rii, mat_jp_riii, iB, idir)
Calculate the expansion coefficients for the atomic terms of the current densitiy in GAPW.
Definition: qs_linres_atom_current.F:108
qs_linres_op
Calculate the operators p rxp and D needed in the optimization of the different contribution of the f...
Definition: qs_linres_op.F:18
qs_linres_op::set_vecp
subroutine, public set_vecp(i1, i2, i3)
...
Definition: qs_linres_op.F:1155
qs_linres_op::fac_vecp
real(dp) function, public fac_vecp(a, b, c)
...
Definition: qs_linres_op.F:1106
qs_linres_op::set_vecp_rev
subroutine, public set_vecp_rev(i1, i2, i3)
...
Definition: qs_linres_op.F:1179
qs_linres_types
Type definitiona for linear response calculations.
Definition: qs_linres_types.F:12
qs_linres_types::allocate_jrho_coeff
subroutine, public allocate_jrho_coeff(jrho1_atom_set, iatom, nsotot)
...
Definition: qs_linres_types.F:980
qs_linres_types::set2zero_jrho_atom_rad
subroutine, public set2zero_jrho_atom_rad(jrho1_atom, ispin)
...
Definition: qs_linres_types.F:943
qs_linres_types::get_current_env
subroutine, public get_current_env(current_env, simple_done, simple_converged, full_done, nao, nstates, gauge, list_cubes, statetrueindex, gauge_name, basisfun_center, nbr_center, center_list, centers_set, psi1_p, psi1_rxp, psi1_D, p_psi0, rxp_psi0, jrho1_atom_set, jrho1_set, chi_tensor, chi_tensor_loc, gauge_atom_radius, rs_gauge, use_old_gauge_atom, chi_pbc, psi0_order)
...
Definition: qs_linres_types.F:426
qs_linres_types::allocate_jrho_atom_rad
subroutine, public allocate_jrho_atom_rad(jrho1_atom, ispin, nr, na, max_iso_not0)
...
Definition: qs_linres_types.F:885
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition: qs_neighbor_list_types.F:17
qs_neighbor_list_types::neighbor_list_iterator_create
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
Definition: qs_neighbor_list_types.F:165
qs_neighbor_list_types::neighbor_list_iterator_release
subroutine, public neighbor_list_iterator_release(iterator_set)
...
Definition: qs_neighbor_list_types.F:251
qs_neighbor_list_types::neighbor_list_iterate
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
Definition: qs_neighbor_list_types.F:351
qs_neighbor_list_types::get_iterator_info
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
Definition: qs_neighbor_list_types.F:549
qs_oce_methods
Routines for the construction of the coefficients for the expansion of the atomic densities rho1_hard...
Definition: qs_oce_methods.F:16
qs_oce_methods::proj_blk
subroutine, public proj_blk(h_a, s_a, na, h_b, s_b, nb, blk, ldb, proj_h, proj_s, nso, len1, len2, fac, distab)
Project a matrix block onto the local atomic functions.
Definition: qs_oce_methods.F:784
qs_oce_types
Definition: qs_oce_types.F:8
qs_rho_atom_types
Definition: qs_rho_atom_types.F:8
sap_kind_types
General overlap type integrals containers.
Definition: sap_kind_types.F:12
sap_kind_types::alist_pre_align_blk
subroutine, public alist_pre_align_blk(blk_in, ldin, blk_out, ldout, ilist, in, jlist, jn)
...
Definition: sap_kind_types.F:172
sap_kind_types::get_alist
subroutine, public get_alist(sap_int, alist, atom)
...
Definition: sap_kind_types.F:140
util
All kind of helpful little routines.
Definition: util.F:14
util::get_limit
pure integer function, dimension(2), public get_limit(m, n, me)
divide m entries into n parts, return size of part me
Definition: util.F:333