Line data Source code
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 : ! **************************************************************************************************
17 : MODULE qs_linres_atom_current
18 : USE atomic_kind_types, ONLY: atomic_kind_type,&
19 : get_atomic_kind
20 : USE basis_set_types, ONLY: get_gto_basis_set,&
21 : gto_basis_set_p_type,&
22 : gto_basis_set_type
23 : USE cell_types, ONLY: cell_type,&
24 : pbc
25 : USE cp_control_types, ONLY: dft_control_type
26 : USE cp_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 :
87 : CONTAINS
88 :
89 : ! **************************************************************************************************
90 : !> \brief Calculate the expansion coefficients for the atomic terms
91 : !> of the current densitiy in GAPW
92 : !> \param qs_env ...
93 : !> \param current_env ...
94 : !> \param mat_d0 ...
95 : !> \param mat_jp ...
96 : !> \param mat_jp_rii ...
97 : !> \param mat_jp_riii ...
98 : !> \param iB ...
99 : !> \param idir ...
100 : !> \par History
101 : !> 07.2006 created [MI]
102 : !> 02.2009 using new setup of projector-basis overlap [jgh]
103 : !> 08.2016 add OpenMP [EPCC]
104 : !> 09.2016 use automatic arrays [M Tucker]
105 : ! **************************************************************************************************
106 864 : 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 864 : 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 864 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: proj_work1, proj_work2
127 864 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: a_matrix, b_matrix, c_matrix, d_matrix
128 864 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: C_coeff_hh_a, C_coeff_hh_b, &
129 864 : C_coeff_ss_a, C_coeff_ss_b, r_coef_h, &
130 864 : r_coef_s, tmp_coeff, zero_coeff
131 : TYPE(alist_type), POINTER :: alist_ac, alist_bc
132 864 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
133 864 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_a, mat_b, mat_c, mat_d
134 : TYPE(dft_control_type), POINTER :: dft_control
135 864 : 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 864 : 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 864 : DIMENSION(:), POINTER :: nl_iterator
141 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
142 864 : POINTER :: sab_all
143 : TYPE(oce_matrix_type), POINTER :: oce
144 864 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
145 864 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: a_block, b_block, c_block, d_block, &
146 864 : jp2_RARnu, jp_RARnu
147 :
148 864 : !$ INTEGER(kind=omp_lock_kind), ALLOCATABLE, DIMENSION(:) :: proj_blk_lock, alloc_lock
149 : !$ INTEGER :: lock
150 :
151 864 : CALL timeset(routineN, handle)
152 :
153 864 : NULLIFY (atomic_kind_set, qs_kind_set, dft_control, sab_all, jrho1_atom_set, oce, &
154 864 : 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 864 : para_env=para_env)
163 864 : CPASSERT(ASSOCIATED(oce))
164 :
165 864 : CALL get_current_env(current_env=current_env, jrho1_atom_set=jrho1_atom_set)
166 864 : CPASSERT(ASSOCIATED(jrho1_atom_set))
167 :
168 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
169 864 : maxsgf=max_nsgf, maxgtops=max_gau, basis_type="GAPW_1C")
170 :
171 864 : eps_cpc = dft_control%qs_control%gapw_control%eps_cpc
172 :
173 864 : idir2 = 1
174 864 : IF (idir /= iB) THEN
175 576 : CALL set_vecp_rev(idir, iB, idir2)
176 : END IF
177 864 : CALL set_vecp(iB, ii, iii)
178 :
179 : ! Set pointers for the different gauge
180 864 : mat_a => mat_d0
181 864 : mat_b => mat_jp
182 864 : mat_c => mat_jp_rii
183 864 : mat_d => mat_jp_riii
184 :
185 : ! Density-like matrices
186 864 : nkind = SIZE(qs_kind_set)
187 864 : natom = SIZE(jrho1_atom_set)
188 864 : nspins = dft_control%nspins
189 :
190 : ! Reset CJC coefficients and local density arrays
191 2466 : DO ikind = 1, nkind
192 1602 : NULLIFY (atom_list)
193 : CALL get_atomic_kind(atomic_kind_set(ikind), &
194 : atom_list=atom_list, &
195 1602 : natom=nat)
196 1602 : CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
197 :
198 : ! Quick cycle if needed.
199 1602 : IF (.NOT. paw_atom) CYCLE
200 :
201 : ! Initialize the density matrix-like arrays.
202 5400 : DO iat = 1, nat
203 1692 : iatom = atom_list(iat)
204 5940 : DO ispin = 1, nspins
205 4338 : IF (ASSOCIATED(jrho1_atom_set(iatom)%cjc0_h(1)%r_coef)) THEN
206 726688 : jrho1_atom_set(iatom)%cjc0_h(ispin)%r_coef = 0.0_dp
207 726688 : jrho1_atom_set(iatom)%cjc0_s(ispin)%r_coef = 0.0_dp
208 726688 : jrho1_atom_set(iatom)%cjc_h(ispin)%r_coef = 0.0_dp
209 726688 : jrho1_atom_set(iatom)%cjc_s(ispin)%r_coef = 0.0_dp
210 726688 : jrho1_atom_set(iatom)%cjc_ii_h(ispin)%r_coef = 0.0_dp
211 726688 : jrho1_atom_set(iatom)%cjc_ii_s(ispin)%r_coef = 0.0_dp
212 726688 : jrho1_atom_set(iatom)%cjc_iii_h(ispin)%r_coef = 0.0_dp
213 726688 : 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 4194 : ALLOCATE (basis_set_list(nkind))
221 2466 : DO ikind = 1, nkind
222 1602 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set_a)
223 2466 : IF (ASSOCIATED(basis_set_a)) THEN
224 1602 : basis_set_list(ikind)%gto_basis_set => basis_set_a
225 : ELSE
226 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
227 : END IF
228 : END DO
229 :
230 864 : len_PC1 = max_nsgf*max_gau
231 864 : len_CPC = max_gau*max_gau
232 :
233 : num_pe = 1
234 864 : !$ num_pe = omp_get_max_threads()
235 864 : 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 864 : !$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 :
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 : ALLOCATE (proj_work1(len_PC1), proj_work2(len_CPC))
274 :
275 : !$OMP SINGLE
276 : !$ ALLOCATE (alloc_lock(natom))
277 : !$ ALLOCATE (proj_blk_lock(nspins*natom))
278 : !$OMP END SINGLE
279 :
280 : !$OMP DO
281 : !$ DO lock = 1, natom
282 : !$ call omp_init_lock(alloc_lock(lock))
283 : !$ END DO
284 : !$OMP END DO
285 :
286 : !$OMP DO
287 : !$ DO lock = 1, nspins*natom
288 : !$ call omp_init_lock(proj_blk_lock(lock))
289 : !$ END DO
290 : !$OMP END DO
291 :
292 : mepos = 0
293 : !$ mepos = omp_get_thread_num()
294 :
295 : DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
296 :
297 : CALL get_iterator_info(nl_iterator, mepos=mepos, &
298 : ikind=ikind, jkind=jkind, &
299 : iatom=iatom, jatom=jatom, cell=cell_b, r=rab)
300 : basis_set_a => basis_set_list(ikind)%gto_basis_set
301 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
302 : basis_set_b => basis_set_list(jkind)%gto_basis_set
303 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
304 : nsgfa = basis_set_a%nsgf
305 : nsgfb = basis_set_b%nsgf
306 : DO ispin = 1, nspins
307 : NULLIFY (jp_RARnu(ispin)%r_coef, jp2_RARnu(ispin)%r_coef)
308 : ALLOCATE (jp_RARnu(ispin)%r_coef(nsgfa, nsgfb), &
309 : jp2_RARnu(ispin)%r_coef(nsgfa, nsgfb))
310 : END DO
311 :
312 : ! Take the block \mu\nu of jpab, jpab_ii and jpab_iii
313 : jmax = 0._dp
314 : DO ispin = 1, nspins
315 : NULLIFY (a_block(ispin)%r_coef)
316 : NULLIFY (b_block(ispin)%r_coef)
317 : NULLIFY (c_block(ispin)%r_coef)
318 : NULLIFY (d_block(ispin)%r_coef)
319 : CALL dbcsr_get_block_p(matrix=mat_a(ispin)%matrix, &
320 : row=iatom, col=jatom, block=a_block(ispin)%r_coef, &
321 : found=den_found)
322 : jmax = jmax + MAXVAL(ABS(a_block(ispin)%r_coef))
323 : CALL dbcsr_get_block_p(matrix=mat_b(ispin)%matrix, &
324 : row=iatom, col=jatom, block=b_block(ispin)%r_coef, &
325 : found=den_found)
326 : jmax = jmax + MAXVAL(ABS(b_block(ispin)%r_coef))
327 : CALL dbcsr_get_block_p(matrix=mat_c(ispin)%matrix, &
328 : row=iatom, col=jatom, block=c_block(ispin)%r_coef, &
329 : found=den_found)
330 : jmax = jmax + MAXVAL(ABS(c_block(ispin)%r_coef))
331 : CALL dbcsr_get_block_p(matrix=mat_d(ispin)%matrix, &
332 : row=iatom, col=jatom, block=d_block(ispin)%r_coef, &
333 : found=den_found)
334 : jmax = jmax + MAXVAL(ABS(d_block(ispin)%r_coef))
335 : END DO
336 :
337 : ! Loop over atoms
338 : DO kkind = 1, nkind
339 : CALL get_qs_kind(qs_kind_set(kkind), &
340 : basis_set=basis_1c_set, basis_type="GAPW_1C", &
341 : paw_atom=paw_atom)
342 :
343 : ! Quick cycle if needed.
344 : IF (.NOT. paw_atom) CYCLE
345 :
346 : CALL get_paw_basis_info(basis_1c_set, nsatbas=nsatbas)
347 : nsoctot = nsatbas
348 :
349 : iac = ikind + nkind*(kkind - 1)
350 : ibc = jkind + nkind*(kkind - 1)
351 : IF (.NOT. ASSOCIATED(oce%intac(iac)%alist)) CYCLE
352 : IF (.NOT. ASSOCIATED(oce%intac(ibc)%alist)) CYCLE
353 :
354 : CALL get_alist(oce%intac(iac), alist_ac, iatom)
355 : CALL get_alist(oce%intac(ibc), alist_bc, jatom)
356 : IF (.NOT. ASSOCIATED(alist_ac)) CYCLE
357 : IF (.NOT. ASSOCIATED(alist_bc)) CYCLE
358 :
359 : DO kac = 1, alist_ac%nclist
360 : DO kbc = 1, alist_bc%nclist
361 : IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) CYCLE
362 :
363 : IF (ALL(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
364 : IF (jmax*alist_bc%clist(kbc)%maxac*alist_ac%clist(kac)%maxac < eps_cpc) CYCLE
365 :
366 : n_cont_a = alist_ac%clist(kac)%nsgf_cnt
367 : n_cont_b = alist_bc%clist(kbc)%nsgf_cnt
368 : sgf_soft_only_a = alist_ac%clist(kac)%sgf_soft_only
369 : sgf_soft_only_b = alist_bc%clist(kbc)%sgf_soft_only
370 : IF (n_cont_a == 0 .OR. n_cont_b == 0) CYCLE
371 :
372 : ! thanks to the linearity of the response, we
373 : ! can avoid computing soft-soft interactions.
374 : ! those terms are already included in the
375 : ! regular grid.
376 : IF (sgf_soft_only_a .AND. sgf_soft_only_b) CYCLE
377 :
378 : list_a => alist_ac%clist(kac)%sgf_list
379 : list_b => alist_bc%clist(kbc)%sgf_list
380 :
381 : katom = alist_ac%clist(kac)%catom
382 : !$ CALL omp_set_lock(alloc_lock(katom))
383 : IF (.NOT. ASSOCIATED(jrho1_atom_set(katom)%cjc0_h(1)%r_coef)) THEN
384 : CALL allocate_jrho_coeff(jrho1_atom_set, katom, nsoctot)
385 : END IF
386 : !$ CALL omp_unset_lock(alloc_lock(katom))
387 :
388 : ! Compute the modified Qai matrix as
389 : ! mQai_\mu\nu = Qai_\mu\nu - Qbi_\mu\nu * (R_A-R_\nu)_ii
390 : ! + Qci_\mu\nu * (R_A-R_\nu)_iii
391 : rbc = alist_bc%clist(kbc)%rac
392 : DO ispin = 1, nspins
393 : CALL DCOPY(nsgfa*nsgfb, b_block(ispin)%r_coef(1, 1), 1, &
394 : jp_RARnu(ispin)%r_coef(1, 1), 1)
395 : CALL DAXPY(nsgfa*nsgfb, -rbc(ii), d_block(ispin)%r_coef(1, 1), 1, &
396 : jp_RARnu(ispin)%r_coef(1, 1), 1)
397 : CALL DAXPY(nsgfa*nsgfb, rbc(iii), c_block(ispin)%r_coef(1, 1), 1, &
398 : jp_RARnu(ispin)%r_coef(1, 1), 1)
399 : END DO
400 :
401 : ! Get the d_A's for the hard and soft densities.
402 : IF (iatom == katom .AND. ALL(alist_ac%clist(kac)%cell == 0)) THEN
403 : C_coeff_hh_a => alist_ac%clist(kac)%achint(:, :, 1)
404 : C_coeff_ss_a => alist_ac%clist(kac)%acint(:, :, 1)
405 : dista = .FALSE.
406 : ELSE
407 : C_coeff_hh_a => alist_ac%clist(kac)%acint(:, :, 1)
408 : C_coeff_ss_a => alist_ac%clist(kac)%acint(:, :, 1)
409 : dista = .TRUE.
410 : END IF
411 : ! Get the d_B's for the hard and soft densities.
412 : IF (jatom == katom .AND. ALL(alist_bc%clist(kbc)%cell == 0)) THEN
413 : C_coeff_hh_b => alist_bc%clist(kbc)%achint(:, :, 1)
414 : C_coeff_ss_b => alist_bc%clist(kbc)%acint(:, :, 1)
415 : distb = .FALSE.
416 : ELSE
417 : C_coeff_hh_b => alist_bc%clist(kbc)%acint(:, :, 1)
418 : C_coeff_ss_b => alist_bc%clist(kbc)%acint(:, :, 1)
419 : distb = .TRUE.
420 : END IF
421 :
422 : distab = dista .AND. distb
423 :
424 : nso = nsoctot
425 :
426 : DO ispin = 1, nspins
427 :
428 : ! align the blocks
429 : CALL alist_pre_align_blk(a_block(ispin)%r_coef, SIZE(a_block(ispin)%r_coef, 1), &
430 : a_matrix, SIZE(a_matrix, 1), &
431 : list_a, n_cont_a, list_b, n_cont_b)
432 :
433 : CALL alist_pre_align_blk(jp_RARnu(ispin)%r_coef, SIZE(jp_RARnu(ispin)%r_coef, 1), &
434 : b_matrix, SIZE(b_matrix, 1), &
435 : list_a, n_cont_a, list_b, n_cont_b)
436 :
437 : CALL alist_pre_align_blk(c_block(ispin)%r_coef, SIZE(c_block(ispin)%r_coef, 1), &
438 : c_matrix, SIZE(c_matrix, 1), &
439 : list_a, n_cont_a, list_b, n_cont_b)
440 : CALL alist_pre_align_blk(d_block(ispin)%r_coef, SIZE(d_block(ispin)%r_coef, 1), &
441 : d_matrix, SIZE(d_matrix, 1), &
442 : list_a, n_cont_a, list_b, n_cont_b)
443 :
444 : !$ CALL omp_set_lock(proj_blk_lock((katom - 1)*nspins + ispin))
445 : !------------------------------------------------------------------
446 : ! P_\alpha\alpha'
447 : r_coef_h => jrho1_atom_set(katom)%cjc0_h(ispin)%r_coef
448 : r_coef_s => jrho1_atom_set(katom)%cjc0_s(ispin)%r_coef
449 : CALL proj_blk(C_coeff_hh_a, C_coeff_ss_a, n_cont_a, &
450 : C_coeff_hh_b, C_coeff_ss_b, n_cont_b, &
451 : a_matrix, max_nsgf, r_coef_h, r_coef_s, nso, &
452 : len_PC1, len_CPC, 1.0_dp, distab, proj_work1, proj_work2)
453 : !------------------------------------------------------------------
454 : ! mQai_\alpha\alpha'
455 : r_coef_h => jrho1_atom_set(katom)%cjc_h(ispin)%r_coef
456 : r_coef_s => jrho1_atom_set(katom)%cjc_s(ispin)%r_coef
457 : CALL proj_blk(C_coeff_hh_a, C_coeff_ss_a, n_cont_a, &
458 : C_coeff_hh_b, C_coeff_ss_b, n_cont_b, &
459 : b_matrix, max_nsgf, r_coef_h, r_coef_s, nso, &
460 : len_PC1, len_CPC, 1.0_dp, distab, proj_work1, proj_work2)
461 : !------------------------------------------------------------------
462 : ! Qci_\alpha\alpha'
463 : r_coef_h => jrho1_atom_set(katom)%cjc_ii_h(ispin)%r_coef
464 : r_coef_s => jrho1_atom_set(katom)%cjc_ii_s(ispin)%r_coef
465 : CALL proj_blk(C_coeff_hh_a, C_coeff_ss_a, n_cont_a, &
466 : C_coeff_hh_b, C_coeff_ss_b, n_cont_b, &
467 : c_matrix, max_nsgf, r_coef_h, r_coef_s, nso, &
468 : len_PC1, len_CPC, 1.0_dp, distab, proj_work1, proj_work2)
469 : !------------------------------------------------------------------
470 : ! Qbi_\alpha\alpha'
471 : r_coef_h => jrho1_atom_set(katom)%cjc_iii_h(ispin)%r_coef
472 : r_coef_s => jrho1_atom_set(katom)%cjc_iii_s(ispin)%r_coef
473 : CALL proj_blk(C_coeff_hh_a, C_coeff_ss_a, n_cont_a, &
474 : C_coeff_hh_b, C_coeff_ss_b, n_cont_b, &
475 : d_matrix, max_nsgf, r_coef_h, r_coef_s, nso, &
476 : len_PC1, len_CPC, 1.0_dp, distab, proj_work1, proj_work2)
477 : !------------------------------------------------------------------
478 : !$ CALL omp_unset_lock(proj_blk_lock((katom - 1)*nspins + ispin))
479 :
480 : END DO ! ispin
481 :
482 : EXIT !search loop over jatom-katom list
483 : END IF
484 : END DO
485 : END DO
486 :
487 : END DO ! kkind
488 : DO ispin = 1, nspins
489 : DEALLOCATE (jp_RARnu(ispin)%r_coef, jp2_RARnu(ispin)%r_coef)
490 : END DO
491 : END DO
492 : ! Wait for all threads to finish the loop before locks can be freed
493 : !$OMP BARRIER
494 :
495 : !$OMP DO
496 : !$ DO lock = 1, natom
497 : !$ call omp_destroy_lock(alloc_lock(lock))
498 : !$ END DO
499 : !$OMP END DO
500 :
501 : !$OMP DO
502 : !$ DO lock = 1, nspins*natom
503 : !$ call omp_destroy_lock(proj_blk_lock(lock))
504 : !$ END DO
505 : !$OMP END DO
506 :
507 : !$OMP SINGLE
508 : !$ DEALLOCATE (alloc_lock)
509 : !$ DEALLOCATE (proj_blk_lock)
510 : !$OMP END SINGLE NOWAIT
511 :
512 : DEALLOCATE (a_matrix, b_matrix, c_matrix, d_matrix, proj_work1, proj_work2, &
513 : a_block, b_block, c_block, d_block, &
514 : jp_RARnu, jp2_RARnu &
515 : )
516 :
517 : !$OMP END PARALLEL
518 :
519 864 : CALL neighbor_list_iterator_release(nl_iterator)
520 864 : DEALLOCATE (basis_set_list)
521 :
522 : ! parallel sum up
523 864 : nbr_dbl = 0.0_dp
524 2466 : DO ikind = 1, nkind
525 : CALL get_atomic_kind(atomic_kind_set(ikind), &
526 : atom_list=atom_list, &
527 1602 : natom=nat)
528 : CALL get_qs_kind(qs_kind_set(ikind), &
529 : basis_set=basis_1c_set, basis_type="GAPW_1C", &
530 1602 : paw_atom=paw_atom)
531 :
532 1602 : IF (.NOT. paw_atom) CYCLE
533 :
534 1242 : CALL get_paw_basis_info(basis_1c_set, nsatbas=nsatbas)
535 1242 : nsoctot = nsatbas
536 :
537 1242 : num_pe = para_env%num_pe
538 1242 : mepos = para_env%mepos
539 1242 : bo = get_limit(nat, num_pe, mepos)
540 :
541 4968 : ALLOCATE (zero_coeff(nsoctot, nsoctot))
542 2934 : DO iat = 1, nat
543 1692 : iatom = atom_list(iat)
544 1692 : is_not_associated = .NOT. ASSOCIATED(jrho1_atom_set(iatom)%cjc0_h(1)%r_coef)
545 :
546 1692 : IF (iat >= bo(1) .AND. iat <= bo(2) .AND. is_not_associated) THEN
547 0 : CALL allocate_jrho_coeff(jrho1_atom_set, iatom, nsoctot)
548 : END IF
549 :
550 5580 : DO ispin = 1, nspins
551 :
552 2646 : tmp_coeff => jrho1_atom_set(iatom)%cjc0_h(ispin)%r_coef
553 2646 : IF (is_not_associated) THEN
554 1026 : zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
555 : END IF
556 1634454 : CALL para_env%sum(tmp_coeff)
557 :
558 2646 : tmp_coeff => jrho1_atom_set(iatom)%cjc0_s(ispin)%r_coef
559 2646 : IF (is_not_associated) THEN
560 1026 : zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
561 : END IF
562 1634454 : CALL para_env%sum(tmp_coeff)
563 :
564 2646 : tmp_coeff => jrho1_atom_set(iatom)%cjc_h(ispin)%r_coef
565 2646 : IF (is_not_associated) THEN
566 1026 : zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
567 : END IF
568 :
569 1634454 : CALL para_env%sum(tmp_coeff)
570 2646 : tmp_coeff => jrho1_atom_set(iatom)%cjc_s(ispin)%r_coef
571 2646 : IF (is_not_associated) THEN
572 1026 : zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
573 : END IF
574 1634454 : CALL para_env%sum(tmp_coeff)
575 :
576 2646 : tmp_coeff => jrho1_atom_set(iatom)%cjc_ii_h(ispin)%r_coef
577 2646 : IF (is_not_associated) THEN
578 1026 : zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
579 : END IF
580 1634454 : CALL para_env%sum(tmp_coeff)
581 :
582 2646 : tmp_coeff => jrho1_atom_set(iatom)%cjc_ii_s(ispin)%r_coef
583 2646 : IF (is_not_associated) THEN
584 1026 : zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
585 : END IF
586 1634454 : CALL para_env%sum(tmp_coeff)
587 :
588 2646 : tmp_coeff => jrho1_atom_set(iatom)%cjc_iii_h(ispin)%r_coef
589 2646 : IF (is_not_associated) THEN
590 1026 : zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
591 : END IF
592 1634454 : CALL para_env%sum(tmp_coeff)
593 :
594 2646 : tmp_coeff => jrho1_atom_set(iatom)%cjc_iii_s(ispin)%r_coef
595 2646 : IF (is_not_associated) THEN
596 1026 : zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
597 : END IF
598 1634454 : CALL para_env%sum(tmp_coeff)
599 :
600 4338 : IF (ASSOCIATED(jrho1_atom_set(iatom)%cjc0_h(ispin)%r_coef)) THEN
601 7884 : nbr_dbl = nbr_dbl + 8.0_dp*REAL(SIZE(jrho1_atom_set(iatom)%cjc0_h(ispin)%r_coef), dp)
602 : END IF
603 : END DO ! ispin
604 : END DO ! iat
605 :
606 5310 : DEALLOCATE (zero_coeff)
607 :
608 : END DO ! ikind
609 :
610 864 : output_unit = cp_logger_get_default_io_unit()
611 864 : IF (output_unit > 0) THEN
612 432 : WRITE (output_unit, '(A,E8.2)') 'calculate_jrho_atom_coeff: nbr_dbl=', nbr_dbl
613 : END IF
614 :
615 864 : CALL timestop(handle)
616 :
617 3456 : END SUBROUTINE calculate_jrho_atom_coeff
618 :
619 : ! **************************************************************************************************
620 : !> \brief ...
621 : !> \param qs_env ...
622 : !> \param current_env ...
623 : !> \param idir ...
624 : ! **************************************************************************************************
625 864 : 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 864 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list, dacg_n_list
640 864 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list, dacg_list
641 : INTEGER, DIMENSION(2) :: bo
642 864 : INTEGER, DIMENSION(:), POINTER :: atom_list, lmax, lmin, npgf, o2nindex
643 : LOGICAL :: paw_atom
644 864 : LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: is_set_to_0
645 : REAL(dp) :: hard_radius
646 864 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: g1, g2, gauge_h, gauge_s
647 864 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: cjc0_h_block, cjc0_s_block, cjc_h_block, &
648 864 : cjc_ii_h_block, cjc_ii_s_block, cjc_iii_h_block, cjc_iii_s_block, cjc_s_block, dgg_1, gg, &
649 864 : gg_lm1
650 864 : REAL(dp), DIMENSION(:, :), POINTER :: coeff, Fr_a_h, Fr_a_h_ii, Fr_a_h_iii, Fr_a_s, &
651 864 : 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 864 : Fr_h, Fr_s, zet
653 : REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG
654 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: my_CG_dxyz_asym
655 864 : 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 864 : 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 864 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
664 :
665 864 : CALL timeset(routineN, handle)
666 : !
667 864 : NULLIFY (atomic_kind_set, qs_kind_set, dft_control, para_env, &
668 864 : coeff, Fr_h, Fr_s, Fr_a_h, Fr_a_s, Fr_a_h_ii, Fr_a_s_ii, &
669 864 : Fr_a_h_iii, Fr_a_s_iii, Fr_b_h, Fr_b_s, Fr_b_h_ii, &
670 864 : Fr_b_s_ii, Fr_b_h_iii, Fr_b_s_iii, jrho1_atom_set, &
671 864 : 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 864 : para_env=para_env)
678 :
679 864 : 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 864 : jrho1_atom_set=jrho1_atom_set)
684 : !
685 :
686 864 : nkind = SIZE(qs_kind_set)
687 864 : nspins = dft_control%nspins
688 : !
689 864 : natom_tot = SIZE(jrho1_atom_set, 1)
690 3456 : ALLOCATE (is_set_to_0(natom_tot, nspins))
691 864 : is_set_to_0(:, :) = .FALSE.
692 :
693 : !
694 2466 : DO ikind = 1, nkind
695 1602 : NULLIFY (atom_list, grid_atom, harmonics, basis_1c_set, &
696 1602 : 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 1602 : 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 1602 : basis_type="GAPW_1C")
708 : !
709 : ! Quick cycle if needed.
710 1602 : 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 1242 : maxso=maxso)
717 1242 : CALL get_paw_basis_info(basis_1c_set, o2nindex=o2nindex)
718 : !
719 1242 : nr = grid_atom%nr
720 1242 : na = grid_atom%ng_sphere
721 1242 : max_iso_not0 = harmonics%max_iso_not0
722 1242 : damax_iso_not0 = harmonics%damax_iso_not0
723 1242 : max_max_iso_not0 = MAX(max_iso_not0, damax_iso_not0)
724 1242 : lmax_expansion = indso(1, max_max_iso_not0)
725 1242 : max_s_harm = harmonics%max_s_harm
726 1242 : llmax = harmonics%llmax
727 : !
728 : ! Distribute the atoms of this kind
729 1242 : num_pe = para_env%num_pe
730 1242 : mepos = para_env%mepos
731 1242 : bo = get_limit(natom, num_pe, mepos)
732 : !
733 1242 : my_CG => harmonics%my_CG
734 1242 : my_CG_dxyz_asym => harmonics%my_CG_dxyz_asym
735 : !
736 : ! Allocate some arrays.
737 1242 : max_nso = nsoset(maxl)
738 0 : ALLOCATE (g1(nr), g2(nr), gg(nr, 0:2*maxl), gg_lm1(nr, 0:2*maxl), dgg_1(nr, 0:2*maxl), &
739 0 : cjc0_h_block(max_nso, max_nso), cjc0_s_block(max_nso, max_nso), &
740 0 : cjc_h_block(max_nso, max_nso), cjc_s_block(max_nso, max_nso), &
741 0 : cjc_ii_h_block(max_nso, max_nso), cjc_ii_s_block(max_nso, max_nso), &
742 0 : cjc_iii_h_block(max_nso, max_nso), cjc_iii_s_block(max_nso, max_nso), &
743 0 : cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm), &
744 0 : dacg_list(2, nsoset(maxl)**2, max_s_harm), dacg_n_list(max_s_harm), &
745 47196 : gauge_h(nr), gauge_s(nr))
746 : !
747 : ! Compute the gauge
748 1296 : SELECT CASE (current_env%gauge)
749 : CASE (current_gauge_r)
750 : ! d(r)=r
751 2754 : gauge_h(1:nr) = grid_atom%rad(1:nr)
752 2754 : gauge_s(1:nr) = grid_atom%rad(1:nr)
753 : CASE (current_gauge_r_and_step_func)
754 : ! step function
755 43740 : gauge_h(1:nr) = 0e0_dp
756 43740 : DO ir = 1, nr
757 43740 : IF (grid_atom%rad(ir) <= hard_radius) THEN
758 24138 : gauge_s(ir) = grid_atom%rad(ir)
759 : ELSE
760 18702 : gauge_s(ir) = gauge_h(ir)
761 : END IF
762 : END DO
763 : CASE (current_gauge_atom)
764 : ! d(r)=A
765 15588 : gauge_h(1:nr) = HUGE(0e0_dp) !0e0_dp
766 15588 : gauge_s(1:nr) = HUGE(0e0_dp) !0e0_dp
767 : CASE DEFAULT
768 1242 : CPABORT("Unknown gauge, try again...")
769 : END SELECT
770 : !
771 : !
772 1242 : m1s = 0
773 4968 : DO iset1 = 1, nset
774 : m2s = 0
775 15300 : DO iset2 = 1, nset
776 : CALL get_none0_cg_list(my_CG, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
777 11574 : max_s_harm, lmax_expansion, cg_list, cg_n_list, max_iso_not0_local)
778 11574 : CPASSERT(max_iso_not0_local <= max_iso_not0)
779 : CALL get_none0_cg_list(my_CG_dxyz_asym, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
780 11574 : max_s_harm, lmax_expansion, dacg_list, dacg_n_list, damax_iso_not0_local)
781 11574 : CPASSERT(damax_iso_not0_local <= damax_iso_not0)
782 :
783 11574 : n1s = nsoset(lmax(iset1))
784 36504 : DO ipgf1 = 1, npgf(iset1)
785 24930 : iso1_first = nsoset(lmin(iset1) - 1) + 1 + n1s*(ipgf1 - 1) + m1s
786 24930 : iso1_last = nsoset(lmax(iset1)) + n1s*(ipgf1 - 1) + m1s
787 24930 : size1 = iso1_last - iso1_first + 1
788 24930 : iso1_first = o2nindex(iso1_first)
789 24930 : iso1_last = o2nindex(iso1_last)
790 24930 : i1 = iso1_last - iso1_first + 1
791 24930 : CPASSERT(size1 == i1)
792 24930 : i1 = nsoset(lmin(iset1) - 1) + 1
793 : !
794 1244430 : g1(1:nr) = EXP(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))
795 : !
796 24930 : n2s = nsoset(lmax(iset2))
797 91674 : DO ipgf2 = 1, npgf(iset2)
798 55170 : iso2_first = nsoset(lmin(iset2) - 1) + 1 + n2s*(ipgf2 - 1) + m2s
799 55170 : iso2_last = nsoset(lmax(iset2)) + n2s*(ipgf2 - 1) + m2s
800 55170 : size2 = iso2_last - iso2_first + 1
801 55170 : iso2_first = o2nindex(iso2_first)
802 55170 : iso2_last = o2nindex(iso2_last)
803 55170 : i2 = iso2_last - iso2_first + 1
804 55170 : CPASSERT(size2 == i2)
805 55170 : i2 = nsoset(lmin(iset2) - 1) + 1
806 : !
807 2730510 : g2(1:nr) = EXP(-zet(ipgf2, iset2)*grid_atom%rad2(1:nr))
808 : !
809 55170 : lmin12 = lmin(iset1) + lmin(iset2)
810 55170 : lmax12 = lmax(iset1) + lmax(iset2)
811 : !
812 55170 : gg = 0.0_dp
813 55170 : gg_lm1 = 0.0_dp
814 55170 : dgg_1 = 0.0_dp
815 : !
816 : ! Take only the terms of expansion for L < lmax_expansion
817 55170 : IF (lmin12 <= lmax_expansion) THEN
818 : !
819 55170 : IF (lmax12 > lmax_expansion) lmax12 = lmax_expansion
820 : !
821 55170 : IF (lmin12 == 0) THEN
822 2510514 : gg(1:nr, lmin12) = g1(1:nr)*g2(1:nr)
823 2510514 : gg_lm1(1:nr, lmin12) = 0.0_dp
824 : ELSE
825 219996 : gg(1:nr, lmin12) = grid_atom%rad2l(1:nr, lmin12)*g1(1:nr)*g2(1:nr)
826 219996 : gg_lm1(1:nr, lmin12) = grid_atom%rad2l(1:nr, lmin12 - 1)*g1(1:nr)*g2(1:nr)
827 : END IF
828 : !
829 103122 : DO l = lmin12 + 1, lmax12
830 2363472 : gg(1:nr, l) = grid_atom%rad(1:nr)*gg(1:nr, l - 1)
831 2418642 : gg_lm1(1:nr, l) = gg(1:nr, l - 1)
832 : END DO
833 : !
834 158292 : DO l = lmin12, lmax12
835 : dgg_1(1:nr, l) = 2.0_dp*(zet(ipgf1, iset1) - zet(ipgf2, iset2))&
836 5149152 : & *gg(1:nr, l)*grid_atom%rad(1:nr)
837 : END DO
838 : ELSE
839 : CYCLE
840 : END IF ! lmin12
841 : !
842 118251 : DO iat = bo(1), bo(2)
843 38151 : iatom = atom_list(iat)
844 : !
845 155052 : DO ispin = 1, nspins
846 : !------------------------------------------------------------------
847 : ! P_\alpha\alpha'
848 3125673 : cjc0_h_block = HUGE(1.0_dp)
849 3125673 : cjc0_s_block = HUGE(1.0_dp)
850 : !
851 : ! Hard term
852 61731 : coeff => jrho1_atom_set(iatom)%cjc0_h(ispin)%r_coef
853 : cjc0_h_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = &
854 601002 : coeff(iso1_first:iso1_last, iso2_first:iso2_last)
855 : !
856 : ! Soft term
857 61731 : coeff => jrho1_atom_set(iatom)%cjc0_s(ispin)%r_coef
858 : cjc0_s_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = &
859 601002 : coeff(iso1_first:iso1_last, iso2_first:iso2_last)
860 : !------------------------------------------------------------------
861 : ! mQai_\alpha\alpha'
862 3125673 : cjc_h_block = HUGE(1.0_dp)
863 3125673 : cjc_s_block = HUGE(1.0_dp)
864 : !
865 : ! Hard term
866 61731 : coeff => jrho1_atom_set(iatom)%cjc_h(ispin)%r_coef
867 : cjc_h_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = &
868 601002 : coeff(iso1_first:iso1_last, iso2_first:iso2_last)
869 : !
870 : ! Soft term
871 61731 : coeff => jrho1_atom_set(iatom)%cjc_s(ispin)%r_coef
872 : cjc_s_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = &
873 601002 : coeff(iso1_first:iso1_last, iso2_first:iso2_last)
874 : !------------------------------------------------------------------
875 : ! Qci_\alpha\alpha'
876 3125673 : cjc_ii_h_block = HUGE(1.0_dp)
877 3125673 : cjc_ii_s_block = HUGE(1.0_dp)
878 : !
879 : ! Hard term
880 61731 : 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 601002 : coeff(iso1_first:iso1_last, iso2_first:iso2_last)
883 : !
884 : ! Soft term
885 61731 : 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 601002 : coeff(iso1_first:iso1_last, iso2_first:iso2_last)
888 : !------------------------------------------------------------------
889 : ! Qbi_\alpha\alpha'
890 3125673 : cjc_iii_h_block = HUGE(1.0_dp)
891 3125673 : cjc_iii_s_block = HUGE(1.0_dp)
892 : !
893 : !
894 : ! Hard term
895 61731 : 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 601002 : coeff(iso1_first:iso1_last, iso2_first:iso2_last)
898 : !
899 : ! Soft term
900 61731 : 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 601002 : coeff(iso1_first:iso1_last, iso2_first:iso2_last)
903 : !------------------------------------------------------------------
904 : !
905 : ! Allocation radial functions
906 61731 : jrho1_atom => jrho1_atom_set(iatom)
907 61731 : IF (.NOT. ASSOCIATED(jrho1_atom%jrho_a_h(ispin)%r_coef)) THEN
908 : CALL allocate_jrho_atom_rad(jrho1_atom, ispin, nr, na, &
909 147 : max_max_iso_not0)
910 147 : is_set_to_0(iatom, ispin) = .TRUE.
911 : ELSE
912 61584 : IF (.NOT. is_set_to_0(iatom, ispin)) THEN
913 1176 : CALL set2zero_jrho_atom_rad(jrho1_atom, ispin)
914 1176 : is_set_to_0(iatom, ispin) = .TRUE.
915 : END IF
916 : END IF
917 : !------------------------------------------------------------------
918 : !
919 61731 : Fr_h => jrho1_atom%jrho_h(ispin)%r_coef
920 61731 : Fr_s => jrho1_atom%jrho_s(ispin)%r_coef
921 : !------------------------------------------------------------------
922 : !
923 61731 : Fr_a_h => jrho1_atom%jrho_a_h(ispin)%r_coef
924 61731 : Fr_a_s => jrho1_atom%jrho_a_s(ispin)%r_coef
925 61731 : Fr_b_h => jrho1_atom%jrho_b_h(ispin)%r_coef
926 61731 : Fr_b_s => jrho1_atom%jrho_b_s(ispin)%r_coef
927 : !------------------------------------------------------------------
928 : !
929 61731 : Fr_a_h_ii => jrho1_atom%jrho_a_h_ii(ispin)%r_coef
930 61731 : Fr_a_s_ii => jrho1_atom%jrho_a_s_ii(ispin)%r_coef
931 61731 : Fr_b_h_ii => jrho1_atom%jrho_b_h_ii(ispin)%r_coef
932 61731 : Fr_b_s_ii => jrho1_atom%jrho_b_s_ii(ispin)%r_coef
933 : !------------------------------------------------------------------
934 : !
935 61731 : Fr_a_h_iii => jrho1_atom%jrho_a_h_iii(ispin)%r_coef
936 61731 : Fr_a_s_iii => jrho1_atom%jrho_a_s_iii(ispin)%r_coef
937 61731 : Fr_b_h_iii => jrho1_atom%jrho_b_h_iii(ispin)%r_coef
938 61731 : Fr_b_s_iii => jrho1_atom%jrho_b_s_iii(ispin)%r_coef
939 : !------------------------------------------------------------------
940 : !
941 362160 : DO iso = 1, max_iso_not0_local
942 300429 : l_iso = indso(1, iso) ! not needed
943 300429 : m_iso = indso(2, iso) ! not needed
944 : !
945 858771 : DO icg = 1, cg_n_list(iso)
946 : !
947 496611 : iso1 = cg_list(1, icg, iso)
948 496611 : iso2 = cg_list(2, icg, iso)
949 : !
950 496611 : IF (.NOT. (iso2 > 0 .AND. iso1 > 0)) THEN
951 0 : WRITE (*, *) 'iso1=', iso1, ' iso2=', iso2, ' iso=', iso, ' icg=', icg
952 0 : WRITE (*, *) '.... will stop!'
953 : END IF
954 496611 : CPASSERT(iso2 > 0 .AND. iso1 > 0)
955 : !
956 496611 : l = indso(1, iso1) + indso(1, iso2)
957 496611 : IF (l > lmax_expansion .OR. l < .0) THEN
958 0 : WRITE (*, *) 'calculate_jrho_atom_rad: 1 l', l
959 0 : WRITE (*, *) 'calculate_jrho_atom_rad: 1 lmax_expansion', lmax_expansion
960 0 : WRITE (*, *) '.... will stop!'
961 : END IF
962 496611 : CPASSERT(l <= lmax_expansion)
963 : !------------------------------------------------------------------
964 : ! P0
965 : !
966 496611 : IF (current_env%gauge == 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 4477842 : 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 4477842 : 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 39661029 : 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 39661029 : 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 24512841 : 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 24512841 : my_CG(iso1, iso2, iso)
997 : !------------------------------------------------------------------
998 : ! Rci
999 : !
1000 496611 : IF (current_env%gauge == 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 4477842 : 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 4477842 : 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 20034999 : 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 20034999 : my_CG(iso1, iso2, iso)
1022 : END IF
1023 : !------------------------------------------------------------------
1024 : ! Rbi
1025 : !
1026 797040 : IF (current_env%gauge == 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 4477842 : 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 4477842 : 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 20034999 : 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 20034999 : my_CG(iso1, iso2, iso)
1048 : END IF
1049 : !------------------------------------------------------------------
1050 : END DO !icg
1051 : !
1052 : END DO ! iso
1053 : !
1054 : !
1055 212436 : DO iso = 1, damax_iso_not0_local
1056 112554 : l_iso = indso(1, iso) ! not needed
1057 112554 : m_iso = indso(2, iso) ! not needed
1058 : !
1059 669321 : DO icg = 1, dacg_n_list(iso)
1060 : !
1061 495036 : iso1 = dacg_list(1, icg, iso)
1062 495036 : iso2 = dacg_list(2, icg, iso)
1063 : !
1064 495036 : IF (.NOT. (iso2 > 0 .AND. iso1 > 0)) THEN
1065 0 : WRITE (*, *) 'iso1=', iso1, ' iso2=', iso2, ' iso=', iso, ' icg=', icg
1066 0 : WRITE (*, *) '.... will stop!'
1067 : END IF
1068 495036 : CPASSERT(iso2 > 0 .AND. iso1 > 0)
1069 : !
1070 495036 : l = indso(1, iso1) + indso(1, iso2)
1071 495036 : IF (l > lmax_expansion) THEN
1072 0 : WRITE (*, *) 'calculate_jrho_atom_rad: 1 l', l
1073 0 : WRITE (*, *) 'calculate_jrho_atom_rad: 1 lmax_expansion', lmax_expansion
1074 0 : WRITE (*, *) '.... will stop!'
1075 : END IF
1076 495036 : 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 24139836 : 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 24139836 : my_CG_dxyz_asym(idir, iso1, iso2, iso)
1089 : !
1090 : !------------------------------------------------------------------
1091 : ! Dcij
1092 : !
1093 495036 : IF (current_env%gauge == 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 3569184 : 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 3569184 : 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 20570652 : 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 20570652 : my_CG_dxyz_asym(idir, iso1, iso2, iso)
1115 : END IF
1116 : !------------------------------------------------------------------
1117 : ! Dbij
1118 : !
1119 607590 : IF (current_env%gauge == 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 3569184 : 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 3569184 : 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 20570652 : 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 20570652 : 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 38448 : m2s = m2s + maxso
1154 : END DO ! iset2
1155 4968 : m1s = m1s + maxso
1156 : END DO ! iset1
1157 : !
1158 0 : DEALLOCATE (cjc0_h_block, cjc0_s_block, cjc_h_block, cjc_s_block, cjc_ii_h_block, cjc_ii_s_block, &
1159 0 : cjc_iii_h_block, cjc_iii_s_block, g1, g2, gg, gg_lm1, dgg_1, gauge_h, gauge_s, &
1160 1242 : cg_list, cg_n_list, dacg_list, dacg_n_list)
1161 5310 : DEALLOCATE (o2nindex)
1162 : END DO ! ikind
1163 : !
1164 : !
1165 864 : DEALLOCATE (is_set_to_0)
1166 : !
1167 864 : CALL timestop(handle)
1168 : !
1169 2592 : 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 1323 : SUBROUTINE calculate_jrho_atom_ang(jrho1_atom, jrho_h, jrho_s, grid_atom, &
1186 1323 : 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 1323 : REAL(dp), DIMENSION(:, :), POINTER :: a, Fr_a_h, Fr_a_h_ii, Fr_a_h_iii, &
1201 1323 : 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 1323 : Fr_b_s_ii, Fr_b_s_iii, Fr_h, Fr_s, slm
1203 1323 : REAL(dp), DIMENSION(:), POINTER :: r
1204 : REAL(dp), DIMENSION(:, :, :), ALLOCATABLE :: g
1205 : !
1206 : !
1207 1323 : 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 1323 : Fr_b_h, Fr_b_s, Fr_b_h_ii, Fr_b_s_ii, Fr_b_h_iii, Fr_b_s_iii, &
1209 1323 : a, slm)
1210 : !
1211 0 : CPASSERT(ASSOCIATED(jrho_h))
1212 1323 : CPASSERT(ASSOCIATED(jrho_s))
1213 1323 : CPASSERT(ASSOCIATED(jrho1_atom))
1214 : ! just to be sure...
1215 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_h))
1216 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_s))
1217 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_h))
1218 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_s))
1219 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_h(ispin)%r_coef))
1220 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_s(ispin)%r_coef))
1221 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_h(ispin)%r_coef))
1222 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_s(ispin)%r_coef))
1223 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_h_ii))
1224 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_s_ii))
1225 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_h_ii))
1226 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_s_ii))
1227 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_h_ii(ispin)%r_coef))
1228 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_s_ii(ispin)%r_coef))
1229 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_h_ii(ispin)%r_coef))
1230 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_s_ii(ispin)%r_coef))
1231 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_h_iii))
1232 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_s_iii))
1233 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_h_iii))
1234 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_s_iii))
1235 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_h_iii(ispin)%r_coef))
1236 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_s_iii(ispin)%r_coef))
1237 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_h_iii(ispin)%r_coef))
1238 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_s_iii(ispin)%r_coef))
1239 : !
1240 : !
1241 1323 : nr = grid_atom%nr
1242 1323 : na = grid_atom%ng_sphere
1243 1323 : max_max_iso_not0 = MAX(harmonics%max_iso_not0, harmonics%damax_iso_not0)
1244 5292 : ALLOCATE (g(3, nr, na))
1245 : !------------------------------------------------------------------
1246 : !
1247 1323 : Fr_h => jrho1_atom%jrho_h(ispin)%r_coef
1248 1323 : Fr_s => jrho1_atom%jrho_s(ispin)%r_coef
1249 : !------------------------------------------------------------------
1250 : !
1251 1323 : Fr_a_h => jrho1_atom%jrho_a_h(ispin)%r_coef !Rai
1252 1323 : Fr_a_s => jrho1_atom%jrho_a_s(ispin)%r_coef
1253 1323 : Fr_b_h => jrho1_atom%jrho_b_h(ispin)%r_coef !Daij
1254 1323 : Fr_b_s => jrho1_atom%jrho_b_s(ispin)%r_coef
1255 : !------------------------------------------------------------------
1256 : !
1257 1323 : Fr_a_h_ii => jrho1_atom%jrho_a_h_ii(ispin)%r_coef !Rci
1258 1323 : Fr_a_s_ii => jrho1_atom%jrho_a_s_ii(ispin)%r_coef
1259 1323 : Fr_b_h_ii => jrho1_atom%jrho_b_h_ii(ispin)%r_coef !Dcij
1260 1323 : Fr_b_s_ii => jrho1_atom%jrho_b_s_ii(ispin)%r_coef
1261 : !------------------------------------------------------------------
1262 : !
1263 1323 : Fr_a_h_iii => jrho1_atom%jrho_a_h_iii(ispin)%r_coef !Rbi
1264 1323 : Fr_a_s_iii => jrho1_atom%jrho_a_s_iii(ispin)%r_coef
1265 1323 : Fr_b_h_iii => jrho1_atom%jrho_b_h_iii(ispin)%r_coef !Dbij
1266 1323 : Fr_b_s_iii => jrho1_atom%jrho_b_s_iii(ispin)%r_coef
1267 : !------------------------------------------------------------------
1268 : !
1269 1323 : a => harmonics%a
1270 1323 : slm => harmonics%slm
1271 1323 : r => grid_atom%rad
1272 : !
1273 1323 : CALL set_vecp(iB, iiB, iiiB)
1274 : !
1275 1323 : scale = 0.0_dp
1276 1323 : idir2 = 1
1277 1323 : IF (idir /= iB) THEN
1278 882 : CALL set_vecp_rev(idir, iB, idir2)
1279 882 : scale = fac_vecp(idir, iB, idir2)
1280 : END IF
1281 : !
1282 : ! Set the gauge
1283 1323 : CALL get_gauge()
1284 : !
1285 65943 : DO ir = 1, nr
1286 820323 : DO iso = 1, max_max_iso_not0
1287 38013120 : DO ia = 1, na
1288 37948500 : 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 7380000 : & + scale*(a(idir2, ia)*r(ir) - g(idir2, ir, ia))*Fr_h(ir, iso)
1299 : !
1300 7380000 : 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 7380000 : & + scale*(a(idir2, ia)*r(ir) - g(idir2, ir, ia))*Fr_s(ir, iso)
1307 : !
1308 7380000 : 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 29814120 : & + scale*a(idir2, ia)*Fr_h(ir, iso)
1321 : !
1322 29814120 : 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 29814120 : & + scale*a(idir2, ia)*Fr_s(ir, iso)
1329 : !
1330 29814120 : 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 3969 : DEALLOCATE (g)
1338 : !
1339 : CONTAINS
1340 : !
1341 : ! **************************************************************************************************
1342 : !> \brief ...
1343 : ! **************************************************************************************************
1344 1323 : SUBROUTINE get_gauge()
1345 : INTEGER :: iatom, ixyz, jatom
1346 : REAL(dp) :: ab, pa, pb, point(3), pra(3), prb(3), &
1347 : res, tmp
1348 1323 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: buf
1349 :
1350 3969 : ALLOCATE (buf(natm_gauge))
1351 65943 : DO ir = 1, nr
1352 3238623 : DO ia = 1, na
1353 12690720 : DO ixyz = 1, 3
1354 12690720 : g(ixyz, ir, ia) = 0.0_dp
1355 : END DO
1356 3172680 : point(1) = r(ir)*a(1, ia)
1357 3172680 : point(2) = r(ir)*a(2, ia)
1358 3172680 : point(3) = r(ir)*a(3, ia)
1359 10785600 : DO iatom = 1, natm_gauge
1360 7612920 : buf(iatom) = 1.0_dp
1361 30451680 : pra = point - ratom(:, iatom)
1362 7612920 : pa = SQRT(pra(1)**2 + pra(2)**2 + pra(3)**2)
1363 31404240 : DO jatom = 1, natm_gauge
1364 20618640 : IF (iatom == jatom) CYCLE
1365 52022880 : prb = point - ratom(:, jatom)
1366 13005720 : pb = SQRT(prb(1)**2 + prb(2)**2 + prb(3)**2)
1367 13005720 : ab = SQRT((pra(1) - prb(1))**2 + (pra(2) - prb(2))**2 + (pra(3) - prb(3))**2)
1368 : !
1369 13005720 : tmp = (pa - pb)/ab
1370 13005720 : tmp = 0.5_dp*(3.0_dp - tmp**2)*tmp
1371 13005720 : tmp = 0.5_dp*(3.0_dp - tmp**2)*tmp
1372 13005720 : tmp = 0.5_dp*(3.0_dp - tmp**2)*tmp
1373 28231560 : buf(iatom) = buf(iatom)*0.5_dp*(1.0_dp - tmp)
1374 : END DO
1375 : END DO
1376 12755340 : DO ixyz = 1, 3
1377 9518040 : res = 0.0_dp
1378 32356800 : DO iatom = 1, natm_gauge
1379 32356800 : res = res + ratom(ixyz, iatom)*buf(iatom)
1380 : END DO
1381 32356800 : res = res/SUM(buf(1:natm_gauge))
1382 : !
1383 12690720 : g(ixyz, ir, ia) = res
1384 : END DO
1385 : END DO
1386 : END DO
1387 1323 : DEALLOCATE (buf)
1388 1323 : 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 864 : 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 864 : 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 864 : 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 864 : TYPE(jrho_atom_type), DIMENSION(:), POINTER :: jrho1_atom_set
1418 : TYPE(jrho_atom_type), POINTER :: jrho1_atom
1419 864 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1420 864 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1421 :
1422 864 : NULLIFY (para_env, dft_control)
1423 864 : NULLIFY (jrho1_atom_set, grid_atom, harmonics)
1424 864 : 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 864 : para_env=para_env)
1432 :
1433 : CALL get_current_env(current_env=current_env, &
1434 864 : jrho1_atom_set=jrho1_atom_set)
1435 :
1436 864 : do_igaim = .FALSE.
1437 864 : IF (current_env%gauge == current_gauge_atom) do_igaim = .TRUE.
1438 :
1439 864 : gapw = dft_control%qs_control%gapw
1440 864 : nkind = SIZE(qs_kind_set, 1)
1441 864 : nspins = dft_control%nspins
1442 :
1443 864 : natm_tot = SIZE(particle_set)
1444 2592 : ALLOCATE (ratom(3, natm_tot))
1445 :
1446 864 : IF (gapw) THEN
1447 2466 : DO ikind = 1, nkind
1448 1602 : NULLIFY (atom_list, grid_atom, harmonics)
1449 1602 : 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 1602 : paw_atom=paw_atom)
1455 :
1456 1602 : IF (.NOT. paw_atom) CYCLE
1457 :
1458 : ! Distribute the atoms of this kind
1459 :
1460 1242 : bo = get_limit(natom, para_env%num_pe, para_env%mepos)
1461 :
1462 4554 : DO iat = bo(1), bo(2)
1463 846 : iatom = atom_list(iat)
1464 : NULLIFY (jrho1_atom)
1465 846 : jrho1_atom => jrho1_atom_set(iatom)
1466 :
1467 846 : natm_gauge = 0
1468 3690 : DO jatom = 1, natm_tot
1469 11376 : r(:) = pbc(particle_set(jatom)%r(:) - particle_set(iatom)%r(:), cell)
1470 : ! SQRT(SUM(r(:)**2)) <= 2.0_dp*hard_radius
1471 12222 : IF (SUM(r(:)**2) <= (4.0_dp*hard_radius*hard_radius)) THEN
1472 2160 : natm_gauge = natm_gauge + 1
1473 8640 : ratom(:, natm_gauge) = r(:)
1474 : END IF
1475 : END DO
1476 :
1477 3771 : DO ispin = 1, nspins
1478 3237237 : jrho1_atom%jrho_vec_rad_h(idir, ispin)%r_coef = 0.0_dp
1479 3237237 : 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 2169 : ratom, natm_gauge, iB, idir, ispin)
1486 : END DO !ispin
1487 : END DO !iat
1488 : END DO !ikind
1489 : END IF
1490 :
1491 864 : DEALLOCATE (ratom)
1492 :
1493 864 : END SUBROUTINE calculate_jrho_atom
1494 :
1495 : END MODULE qs_linres_atom_current
|
