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