Line data Source code
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 : ! **************************************************************************************************
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(:, :) :: a_matrix, b_matrix, c_matrix, d_matrix
127 864 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: C_coeff_hh_a, C_coeff_hh_b, &
128 864 : C_coeff_ss_a, C_coeff_ss_b, r_coef_h, &
129 864 : r_coef_s, tmp_coeff, zero_coeff
130 : TYPE(alist_type), POINTER :: alist_ac, alist_bc
131 864 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
132 864 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_a, mat_b, mat_c, mat_d
133 : TYPE(dft_control_type), POINTER :: dft_control
134 864 : 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 864 : 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 864 : DIMENSION(:), POINTER :: nl_iterator
140 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
141 864 : POINTER :: sab_all
142 : TYPE(oce_matrix_type), POINTER :: oce
143 864 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
144 864 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: a_block, b_block, c_block, d_block, &
145 864 : jp2_RARnu, jp_RARnu
146 :
147 864 : !$ INTEGER(kind=omp_lock_kind), ALLOCATABLE, DIMENSION(:) :: proj_blk_lock, alloc_lock
148 : !$ INTEGER :: lock
149 :
150 864 : CALL timeset(routineN, handle)
151 :
152 864 : NULLIFY (atomic_kind_set, qs_kind_set, dft_control, sab_all, jrho1_atom_set, oce, &
153 864 : 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 864 : para_env=para_env)
162 864 : CPASSERT(ASSOCIATED(oce))
163 :
164 864 : CALL get_current_env(current_env=current_env, jrho1_atom_set=jrho1_atom_set)
165 864 : CPASSERT(ASSOCIATED(jrho1_atom_set))
166 :
167 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
168 864 : maxsgf=max_nsgf, maxgtops=max_gau, basis_type="GAPW_1C")
169 :
170 864 : eps_cpc = dft_control%qs_control%gapw_control%eps_cpc
171 :
172 864 : idir2 = 1
173 864 : IF (idir .NE. iB) THEN
174 576 : CALL set_vecp_rev(idir, iB, idir2)
175 : END IF
176 864 : CALL set_vecp(iB, ii, iii)
177 :
178 : ! Set pointers for the different gauge
179 864 : mat_a => mat_d0
180 864 : mat_b => mat_jp
181 864 : mat_c => mat_jp_rii
182 864 : mat_d => mat_jp_riii
183 :
184 : ! Density-like matrices
185 864 : nkind = SIZE(qs_kind_set)
186 864 : natom = SIZE(jrho1_atom_set)
187 864 : nspins = dft_control%nspins
188 :
189 : ! Reset CJC coefficients and local density arrays
190 2466 : DO ikind = 1, nkind
191 1602 : NULLIFY (atom_list)
192 : CALL get_atomic_kind(atomic_kind_set(ikind), &
193 : atom_list=atom_list, &
194 1602 : natom=nat)
195 1602 : CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
196 :
197 : ! Quick cycle if needed.
198 1602 : IF (.NOT. paw_atom) CYCLE
199 :
200 : ! Initialize the density matrix-like arrays.
201 5400 : DO iat = 1, nat
202 1692 : iatom = atom_list(iat)
203 5940 : DO ispin = 1, nspins
204 4338 : IF (ASSOCIATED(jrho1_atom_set(iatom)%cjc0_h(1)%r_coef)) THEN
205 726688 : jrho1_atom_set(iatom)%cjc0_h(ispin)%r_coef = 0.0_dp
206 726688 : jrho1_atom_set(iatom)%cjc0_s(ispin)%r_coef = 0.0_dp
207 726688 : jrho1_atom_set(iatom)%cjc_h(ispin)%r_coef = 0.0_dp
208 726688 : jrho1_atom_set(iatom)%cjc_s(ispin)%r_coef = 0.0_dp
209 726688 : jrho1_atom_set(iatom)%cjc_ii_h(ispin)%r_coef = 0.0_dp
210 726688 : jrho1_atom_set(iatom)%cjc_ii_s(ispin)%r_coef = 0.0_dp
211 726688 : jrho1_atom_set(iatom)%cjc_iii_h(ispin)%r_coef = 0.0_dp
212 726688 : 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 4194 : ALLOCATE (basis_set_list(nkind))
220 2466 : DO ikind = 1, nkind
221 1602 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set_a)
222 2466 : IF (ASSOCIATED(basis_set_a)) THEN
223 1602 : basis_set_list(ikind)%gto_basis_set => basis_set_a
224 : ELSE
225 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
226 : END IF
227 : END DO
228 :
229 864 : len_PC1 = max_nsgf*max_gau
230 864 : len_CPC = max_gau*max_gau
231 :
232 : num_pe = 1
233 864 : !$ num_pe = omp_get_max_threads()
234 864 : 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 864 : !$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 864 : CALL neighbor_list_iterator_release(nl_iterator)
521 864 : DEALLOCATE (basis_set_list)
522 :
523 : ! parallel sum up
524 864 : nbr_dbl = 0.0_dp
525 2466 : DO ikind = 1, nkind
526 : CALL get_atomic_kind(atomic_kind_set(ikind), &
527 : atom_list=atom_list, &
528 1602 : natom=nat)
529 : CALL get_qs_kind(qs_kind_set(ikind), &
530 : basis_set=basis_1c_set, basis_type="GAPW_1C", &
531 1602 : paw_atom=paw_atom)
532 :
533 1602 : IF (.NOT. paw_atom) CYCLE
534 :
535 1242 : CALL get_paw_basis_info(basis_1c_set, nsatbas=nsatbas)
536 1242 : nsoctot = nsatbas
537 :
538 1242 : num_pe = para_env%num_pe
539 1242 : mepos = para_env%mepos
540 1242 : bo = get_limit(nat, num_pe, mepos)
541 :
542 4968 : ALLOCATE (zero_coeff(nsoctot, nsoctot))
543 2934 : DO iat = 1, nat
544 1692 : iatom = atom_list(iat)
545 1692 : is_not_associated = .NOT. ASSOCIATED(jrho1_atom_set(iatom)%cjc0_h(1)%r_coef)
546 :
547 1692 : IF (iat .GE. bo(1) .AND. iat .LE. bo(2) .AND. is_not_associated) THEN
548 0 : CALL allocate_jrho_coeff(jrho1_atom_set, iatom, nsoctot)
549 : END IF
550 :
551 5580 : DO ispin = 1, nspins
552 :
553 2646 : tmp_coeff => jrho1_atom_set(iatom)%cjc0_h(ispin)%r_coef
554 2646 : IF (is_not_associated) THEN
555 1026 : zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
556 : END IF
557 1634454 : CALL para_env%sum(tmp_coeff)
558 :
559 2646 : tmp_coeff => jrho1_atom_set(iatom)%cjc0_s(ispin)%r_coef
560 2646 : IF (is_not_associated) THEN
561 1026 : zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
562 : END IF
563 1634454 : CALL para_env%sum(tmp_coeff)
564 :
565 2646 : tmp_coeff => jrho1_atom_set(iatom)%cjc_h(ispin)%r_coef
566 2646 : IF (is_not_associated) THEN
567 1026 : zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
568 : END IF
569 :
570 1634454 : CALL para_env%sum(tmp_coeff)
571 2646 : tmp_coeff => jrho1_atom_set(iatom)%cjc_s(ispin)%r_coef
572 2646 : IF (is_not_associated) THEN
573 1026 : zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
574 : END IF
575 1634454 : CALL para_env%sum(tmp_coeff)
576 :
577 2646 : tmp_coeff => jrho1_atom_set(iatom)%cjc_ii_h(ispin)%r_coef
578 2646 : IF (is_not_associated) THEN
579 1026 : zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
580 : END IF
581 1634454 : CALL para_env%sum(tmp_coeff)
582 :
583 2646 : tmp_coeff => jrho1_atom_set(iatom)%cjc_ii_s(ispin)%r_coef
584 2646 : IF (is_not_associated) THEN
585 1026 : zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
586 : END IF
587 1634454 : CALL para_env%sum(tmp_coeff)
588 :
589 2646 : tmp_coeff => jrho1_atom_set(iatom)%cjc_iii_h(ispin)%r_coef
590 2646 : IF (is_not_associated) THEN
591 1026 : zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
592 : END IF
593 1634454 : CALL para_env%sum(tmp_coeff)
594 :
595 2646 : tmp_coeff => jrho1_atom_set(iatom)%cjc_iii_s(ispin)%r_coef
596 2646 : IF (is_not_associated) THEN
597 1026 : zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
598 : END IF
599 1634454 : CALL para_env%sum(tmp_coeff)
600 :
601 2646 : IF (ASSOCIATED(jrho1_atom_set(iatom)%cjc0_h(ispin)%r_coef)) &
602 9576 : 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 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 864 : REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG
654 864 : 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 5994 : 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) .LE. 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 .LE. 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 .LE. 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 10133028 : gg = 0.0_dp
813 10133028 : gg_lm1 = 0.0_dp
814 10133028 : dgg_1 = 0.0_dp
815 : !
816 : ! Take only the terms of expansion for L < lmax_expansion
817 55170 : IF (lmin12 .LE. lmax_expansion) THEN
818 : !
819 55170 : IF (lmax12 .GT. 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 .GT. lmax_expansion .OR. l .LT. .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 .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 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 .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 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 .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 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 .GT. 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 .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 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 .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 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 .NE. 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 .EQ. 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 .EQ. 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)) .LE. 2.0_dp*hard_radius
1471 12222 : IF (SUM(r(:)**2) .LE. (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
|
