(git:ed6f26b)
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-2025 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(:, :) :: 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
1495END MODULE qs_linres_atom_current
cell_types::pbc
Definition cell_types.F:92
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:624
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:679
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:502
input_constants::current_gauge_r
integer, parameter, public current_gauge_r
Definition input_constants.F:502
input_constants::current_gauge_r_and_step_func
integer, parameter, public current_gauge_r_and_step_func
Definition input_constants.F:502
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, 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, 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, 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)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:514
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, 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, 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:454
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)
Get attributes of an atomic kind set.
Definition qs_kind_types.F:880
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_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:626
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)
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
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:55
cp_control_types::dft_control_type
Definition cp_control_types.F:570
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:170
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:217
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:171
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