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 : STAT=istat)
271 : CPASSERT(istat == 0)
272 :
273 : ALLOCATE (a_matrix(max_nsgf, max_nsgf), b_matrix(max_nsgf, max_nsgf), &
274 : c_matrix(max_nsgf, max_nsgf), d_matrix(max_nsgf, max_nsgf), &
275 : STAT=istat)
276 : CPASSERT(istat == 0)
277 : ALLOCATE (proj_work1(len_PC1), proj_work2(len_CPC), STAT=istat)
278 : CPASSERT(istat == 0)
279 :
280 : !$OMP SINGLE
281 : !$ ALLOCATE (alloc_lock(natom))
282 : !$ ALLOCATE (proj_blk_lock(nspins*natom))
283 : !$OMP END SINGLE
284 :
285 : !$OMP DO
286 : !$ DO lock = 1, natom
287 : !$ call omp_init_lock(alloc_lock(lock))
288 : !$ END DO
289 : !$OMP END DO
290 :
291 : !$OMP DO
292 : !$ DO lock = 1, nspins*natom
293 : !$ call omp_init_lock(proj_blk_lock(lock))
294 : !$ END DO
295 : !$OMP END DO
296 :
297 : mepos = 0
298 : !$ mepos = omp_get_thread_num()
299 :
300 : DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
301 :
302 : CALL get_iterator_info(nl_iterator, mepos=mepos, &
303 : ikind=ikind, jkind=jkind, &
304 : iatom=iatom, jatom=jatom, cell=cell_b, r=rab)
305 : basis_set_a => basis_set_list(ikind)%gto_basis_set
306 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
307 : basis_set_b => basis_set_list(jkind)%gto_basis_set
308 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
309 : nsgfa = basis_set_a%nsgf
310 : nsgfb = basis_set_b%nsgf
311 : DO ispin = 1, nspins
312 : NULLIFY (jp_RARnu(ispin)%r_coef, jp2_RARnu(ispin)%r_coef)
313 : ALLOCATE (jp_RARnu(ispin)%r_coef(nsgfa, nsgfb), &
314 : jp2_RARnu(ispin)%r_coef(nsgfa, nsgfb))
315 : END DO
316 :
317 : ! Take the block \mu\nu of jpab, jpab_ii and jpab_iii
318 : jmax = 0._dp
319 : DO ispin = 1, nspins
320 : NULLIFY (a_block(ispin)%r_coef)
321 : NULLIFY (b_block(ispin)%r_coef)
322 : NULLIFY (c_block(ispin)%r_coef)
323 : NULLIFY (d_block(ispin)%r_coef)
324 : CALL dbcsr_get_block_p(matrix=mat_a(ispin)%matrix, &
325 : row=iatom, col=jatom, block=a_block(ispin)%r_coef, &
326 : found=den_found)
327 : jmax = jmax + MAXVAL(ABS(a_block(ispin)%r_coef))
328 : CALL dbcsr_get_block_p(matrix=mat_b(ispin)%matrix, &
329 : row=iatom, col=jatom, block=b_block(ispin)%r_coef, &
330 : found=den_found)
331 : jmax = jmax + MAXVAL(ABS(b_block(ispin)%r_coef))
332 : CALL dbcsr_get_block_p(matrix=mat_c(ispin)%matrix, &
333 : row=iatom, col=jatom, block=c_block(ispin)%r_coef, &
334 : found=den_found)
335 : jmax = jmax + MAXVAL(ABS(c_block(ispin)%r_coef))
336 : CALL dbcsr_get_block_p(matrix=mat_d(ispin)%matrix, &
337 : row=iatom, col=jatom, block=d_block(ispin)%r_coef, &
338 : found=den_found)
339 : jmax = jmax + MAXVAL(ABS(d_block(ispin)%r_coef))
340 : END DO
341 :
342 : ! Loop over atoms
343 : DO kkind = 1, nkind
344 : CALL get_qs_kind(qs_kind_set(kkind), &
345 : basis_set=basis_1c_set, basis_type="GAPW_1C", &
346 : paw_atom=paw_atom)
347 :
348 : ! Quick cycle if needed.
349 : IF (.NOT. paw_atom) CYCLE
350 :
351 : CALL get_paw_basis_info(basis_1c_set, nsatbas=nsatbas)
352 : nsoctot = nsatbas
353 :
354 : iac = ikind + nkind*(kkind - 1)
355 : ibc = jkind + nkind*(kkind - 1)
356 : IF (.NOT. ASSOCIATED(oce%intac(iac)%alist)) CYCLE
357 : IF (.NOT. ASSOCIATED(oce%intac(ibc)%alist)) CYCLE
358 :
359 : CALL get_alist(oce%intac(iac), alist_ac, iatom)
360 : CALL get_alist(oce%intac(ibc), alist_bc, jatom)
361 : IF (.NOT. ASSOCIATED(alist_ac)) CYCLE
362 : IF (.NOT. ASSOCIATED(alist_bc)) CYCLE
363 :
364 : DO kac = 1, alist_ac%nclist
365 : DO kbc = 1, alist_bc%nclist
366 : IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) CYCLE
367 :
368 : IF (ALL(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
369 : IF (jmax*alist_bc%clist(kbc)%maxac*alist_ac%clist(kac)%maxac < eps_cpc) CYCLE
370 :
371 : n_cont_a = alist_ac%clist(kac)%nsgf_cnt
372 : n_cont_b = alist_bc%clist(kbc)%nsgf_cnt
373 : sgf_soft_only_a = alist_ac%clist(kac)%sgf_soft_only
374 : sgf_soft_only_b = alist_bc%clist(kbc)%sgf_soft_only
375 : IF (n_cont_a == 0 .OR. n_cont_b == 0) CYCLE
376 :
377 : ! thanks to the linearity of the response, we
378 : ! can avoid computing soft-soft interactions.
379 : ! those terms are already included in the
380 : ! regular grid.
381 : IF (sgf_soft_only_a .AND. sgf_soft_only_b) CYCLE
382 :
383 : list_a => alist_ac%clist(kac)%sgf_list
384 : list_b => alist_bc%clist(kbc)%sgf_list
385 :
386 : katom = alist_ac%clist(kac)%catom
387 : !$ CALL omp_set_lock(alloc_lock(katom))
388 : IF (.NOT. ASSOCIATED(jrho1_atom_set(katom)%cjc0_h(1)%r_coef)) THEN
389 : CALL allocate_jrho_coeff(jrho1_atom_set, katom, nsoctot)
390 : END IF
391 : !$ CALL omp_unset_lock(alloc_lock(katom))
392 :
393 : ! Compute the modified Qai matrix as
394 : ! mQai_\mu\nu = Qai_\mu\nu - Qbi_\mu\nu * (R_A-R_\nu)_ii
395 : ! + Qci_\mu\nu * (R_A-R_\nu)_iii
396 : rbc = alist_bc%clist(kbc)%rac
397 : DO ispin = 1, nspins
398 : CALL DCOPY(nsgfa*nsgfb, b_block(ispin)%r_coef(1, 1), 1, &
399 : jp_RARnu(ispin)%r_coef(1, 1), 1)
400 : CALL DAXPY(nsgfa*nsgfb, -rbc(ii), d_block(ispin)%r_coef(1, 1), 1, &
401 : jp_RARnu(ispin)%r_coef(1, 1), 1)
402 : CALL DAXPY(nsgfa*nsgfb, rbc(iii), c_block(ispin)%r_coef(1, 1), 1, &
403 : jp_RARnu(ispin)%r_coef(1, 1), 1)
404 : END DO
405 :
406 : ! Get the d_A's for the hard and soft densities.
407 : IF (iatom == katom .AND. ALL(alist_ac%clist(kac)%cell == 0)) THEN
408 : C_coeff_hh_a => alist_ac%clist(kac)%achint(:, :, 1)
409 : C_coeff_ss_a => alist_ac%clist(kac)%acint(:, :, 1)
410 : dista = .FALSE.
411 : ELSE
412 : C_coeff_hh_a => alist_ac%clist(kac)%acint(:, :, 1)
413 : C_coeff_ss_a => alist_ac%clist(kac)%acint(:, :, 1)
414 : dista = .TRUE.
415 : END IF
416 : ! Get the d_B's for the hard and soft densities.
417 : IF (jatom == katom .AND. ALL(alist_bc%clist(kbc)%cell == 0)) THEN
418 : C_coeff_hh_b => alist_bc%clist(kbc)%achint(:, :, 1)
419 : C_coeff_ss_b => alist_bc%clist(kbc)%acint(:, :, 1)
420 : distb = .FALSE.
421 : ELSE
422 : C_coeff_hh_b => alist_bc%clist(kbc)%acint(:, :, 1)
423 : C_coeff_ss_b => alist_bc%clist(kbc)%acint(:, :, 1)
424 : distb = .TRUE.
425 : END IF
426 :
427 : distab = dista .AND. distb
428 :
429 : nso = nsoctot
430 :
431 : DO ispin = 1, nspins
432 :
433 : ! align the blocks
434 : CALL alist_pre_align_blk(a_block(ispin)%r_coef, SIZE(a_block(ispin)%r_coef, 1), &
435 : a_matrix, SIZE(a_matrix, 1), &
436 : list_a, n_cont_a, list_b, n_cont_b)
437 :
438 : CALL alist_pre_align_blk(jp_RARnu(ispin)%r_coef, SIZE(jp_RARnu(ispin)%r_coef, 1), &
439 : b_matrix, SIZE(b_matrix, 1), &
440 : list_a, n_cont_a, list_b, n_cont_b)
441 :
442 : CALL alist_pre_align_blk(c_block(ispin)%r_coef, SIZE(c_block(ispin)%r_coef, 1), &
443 : c_matrix, SIZE(c_matrix, 1), &
444 : list_a, n_cont_a, list_b, n_cont_b)
445 : CALL alist_pre_align_blk(d_block(ispin)%r_coef, SIZE(d_block(ispin)%r_coef, 1), &
446 : d_matrix, SIZE(d_matrix, 1), &
447 : list_a, n_cont_a, list_b, n_cont_b)
448 :
449 : !$ CALL omp_set_lock(proj_blk_lock((katom - 1)*nspins + ispin))
450 : !------------------------------------------------------------------
451 : ! P_\alpha\alpha'
452 : r_coef_h => jrho1_atom_set(katom)%cjc0_h(ispin)%r_coef
453 : r_coef_s => jrho1_atom_set(katom)%cjc0_s(ispin)%r_coef
454 : CALL proj_blk(C_coeff_hh_a, C_coeff_ss_a, n_cont_a, &
455 : C_coeff_hh_b, C_coeff_ss_b, n_cont_b, &
456 : a_matrix, max_nsgf, r_coef_h, r_coef_s, nso, &
457 : len_PC1, len_CPC, 1.0_dp, distab, proj_work1, proj_work2)
458 : !------------------------------------------------------------------
459 : ! mQai_\alpha\alpha'
460 : r_coef_h => jrho1_atom_set(katom)%cjc_h(ispin)%r_coef
461 : r_coef_s => jrho1_atom_set(katom)%cjc_s(ispin)%r_coef
462 : CALL proj_blk(C_coeff_hh_a, C_coeff_ss_a, n_cont_a, &
463 : C_coeff_hh_b, C_coeff_ss_b, n_cont_b, &
464 : b_matrix, max_nsgf, r_coef_h, r_coef_s, nso, &
465 : len_PC1, len_CPC, 1.0_dp, distab, proj_work1, proj_work2)
466 : !------------------------------------------------------------------
467 : ! Qci_\alpha\alpha'
468 : r_coef_h => jrho1_atom_set(katom)%cjc_ii_h(ispin)%r_coef
469 : r_coef_s => jrho1_atom_set(katom)%cjc_ii_s(ispin)%r_coef
470 : CALL proj_blk(C_coeff_hh_a, C_coeff_ss_a, n_cont_a, &
471 : C_coeff_hh_b, C_coeff_ss_b, n_cont_b, &
472 : c_matrix, max_nsgf, r_coef_h, r_coef_s, nso, &
473 : len_PC1, len_CPC, 1.0_dp, distab, proj_work1, proj_work2)
474 : !------------------------------------------------------------------
475 : ! Qbi_\alpha\alpha'
476 : r_coef_h => jrho1_atom_set(katom)%cjc_iii_h(ispin)%r_coef
477 : r_coef_s => jrho1_atom_set(katom)%cjc_iii_s(ispin)%r_coef
478 : CALL proj_blk(C_coeff_hh_a, C_coeff_ss_a, n_cont_a, &
479 : C_coeff_hh_b, C_coeff_ss_b, n_cont_b, &
480 : d_matrix, max_nsgf, r_coef_h, r_coef_s, nso, &
481 : len_PC1, len_CPC, 1.0_dp, distab, proj_work1, proj_work2)
482 : !------------------------------------------------------------------
483 : !$ CALL omp_unset_lock(proj_blk_lock((katom - 1)*nspins + ispin))
484 :
485 : END DO ! ispin
486 :
487 : EXIT !search loop over jatom-katom list
488 : END IF
489 : END DO
490 : END DO
491 :
492 : END DO ! kkind
493 : DO ispin = 1, nspins
494 : DEALLOCATE (jp_RARnu(ispin)%r_coef, jp2_RARnu(ispin)%r_coef)
495 : END DO
496 : END DO
497 : ! Wait for all threads to finish the loop before locks can be freed
498 : !$OMP BARRIER
499 :
500 : !$OMP DO
501 : !$ DO lock = 1, natom
502 : !$ call omp_destroy_lock(alloc_lock(lock))
503 : !$ END DO
504 : !$OMP END DO
505 :
506 : !$OMP DO
507 : !$ DO lock = 1, nspins*natom
508 : !$ call omp_destroy_lock(proj_blk_lock(lock))
509 : !$ END DO
510 : !$OMP END DO
511 :
512 : !$OMP SINGLE
513 : !$ DEALLOCATE (alloc_lock)
514 : !$ DEALLOCATE (proj_blk_lock)
515 : !$OMP END SINGLE NOWAIT
516 :
517 : DEALLOCATE (a_matrix, b_matrix, c_matrix, d_matrix, proj_work1, proj_work2, &
518 : a_block, b_block, c_block, d_block, &
519 : jp_RARnu, jp2_RARnu &
520 : )
521 :
522 : !$OMP END PARALLEL
523 :
524 864 : CALL neighbor_list_iterator_release(nl_iterator)
525 864 : DEALLOCATE (basis_set_list)
526 :
527 : ! parallel sum up
528 864 : nbr_dbl = 0.0_dp
529 2466 : DO ikind = 1, nkind
530 : CALL get_atomic_kind(atomic_kind_set(ikind), &
531 : atom_list=atom_list, &
532 1602 : natom=nat)
533 : CALL get_qs_kind(qs_kind_set(ikind), &
534 : basis_set=basis_1c_set, basis_type="GAPW_1C", &
535 1602 : paw_atom=paw_atom)
536 :
537 1602 : IF (.NOT. paw_atom) CYCLE
538 :
539 1242 : CALL get_paw_basis_info(basis_1c_set, nsatbas=nsatbas)
540 1242 : nsoctot = nsatbas
541 :
542 1242 : num_pe = para_env%num_pe
543 1242 : mepos = para_env%mepos
544 1242 : bo = get_limit(nat, num_pe, mepos)
545 :
546 4968 : ALLOCATE (zero_coeff(nsoctot, nsoctot))
547 2934 : DO iat = 1, nat
548 1692 : iatom = atom_list(iat)
549 1692 : is_not_associated = .NOT. ASSOCIATED(jrho1_atom_set(iatom)%cjc0_h(1)%r_coef)
550 :
551 1692 : IF (iat >= bo(1) .AND. iat <= bo(2) .AND. is_not_associated) THEN
552 0 : CALL allocate_jrho_coeff(jrho1_atom_set, iatom, nsoctot)
553 : END IF
554 :
555 5580 : DO ispin = 1, nspins
556 :
557 2646 : tmp_coeff => jrho1_atom_set(iatom)%cjc0_h(ispin)%r_coef
558 2646 : IF (is_not_associated) THEN
559 1026 : zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
560 : END IF
561 1634454 : CALL para_env%sum(tmp_coeff)
562 :
563 2646 : tmp_coeff => jrho1_atom_set(iatom)%cjc0_s(ispin)%r_coef
564 2646 : IF (is_not_associated) THEN
565 1026 : zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
566 : END IF
567 1634454 : CALL para_env%sum(tmp_coeff)
568 :
569 2646 : tmp_coeff => jrho1_atom_set(iatom)%cjc_h(ispin)%r_coef
570 2646 : IF (is_not_associated) THEN
571 1026 : zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
572 : END IF
573 :
574 1634454 : CALL para_env%sum(tmp_coeff)
575 2646 : tmp_coeff => jrho1_atom_set(iatom)%cjc_s(ispin)%r_coef
576 2646 : IF (is_not_associated) THEN
577 1026 : zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
578 : END IF
579 1634454 : CALL para_env%sum(tmp_coeff)
580 :
581 2646 : tmp_coeff => jrho1_atom_set(iatom)%cjc_ii_h(ispin)%r_coef
582 2646 : IF (is_not_associated) THEN
583 1026 : zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
584 : END IF
585 1634454 : CALL para_env%sum(tmp_coeff)
586 :
587 2646 : tmp_coeff => jrho1_atom_set(iatom)%cjc_ii_s(ispin)%r_coef
588 2646 : IF (is_not_associated) THEN
589 1026 : zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
590 : END IF
591 1634454 : CALL para_env%sum(tmp_coeff)
592 :
593 2646 : tmp_coeff => jrho1_atom_set(iatom)%cjc_iii_h(ispin)%r_coef
594 2646 : IF (is_not_associated) THEN
595 1026 : zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
596 : END IF
597 1634454 : CALL para_env%sum(tmp_coeff)
598 :
599 2646 : tmp_coeff => jrho1_atom_set(iatom)%cjc_iii_s(ispin)%r_coef
600 2646 : IF (is_not_associated) THEN
601 1026 : zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
602 : END IF
603 1634454 : CALL para_env%sum(tmp_coeff)
604 :
605 2646 : IF (ASSOCIATED(jrho1_atom_set(iatom)%cjc0_h(ispin)%r_coef)) &
606 9576 : nbr_dbl = nbr_dbl + 8.0_dp*REAL(SIZE(jrho1_atom_set(iatom)%cjc0_h(ispin)%r_coef), dp)
607 : END DO ! ispin
608 : END DO ! iat
609 :
610 5310 : DEALLOCATE (zero_coeff)
611 :
612 : END DO ! ikind
613 :
614 864 : output_unit = cp_logger_get_default_io_unit()
615 864 : IF (output_unit > 0) THEN
616 432 : WRITE (output_unit, '(A,E8.2)') 'calculate_jrho_atom_coeff: nbr_dbl=', nbr_dbl
617 : END IF
618 :
619 864 : CALL timestop(handle)
620 :
621 3456 : END SUBROUTINE calculate_jrho_atom_coeff
622 :
623 : ! **************************************************************************************************
624 : !> \brief ...
625 : !> \param qs_env ...
626 : !> \param current_env ...
627 : !> \param idir ...
628 : ! **************************************************************************************************
629 864 : SUBROUTINE calculate_jrho_atom_rad(qs_env, current_env, idir)
630 : !
631 : TYPE(qs_environment_type), POINTER :: qs_env
632 : TYPE(current_env_type) :: current_env
633 : INTEGER, INTENT(IN) :: idir
634 :
635 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_jrho_atom_rad'
636 :
637 : INTEGER :: damax_iso_not0, damax_iso_not0_local, handle, i1, i2, iat, iatom, icg, ikind, &
638 : ipgf1, ipgf2, ir, iset1, iset2, iso, iso1, iso1_first, iso1_last, iso2, iso2_first, &
639 : iso2_last, ispin, l, l_iso, llmax, lmax12, lmax_expansion, lmin12, m1s, m2s, m_iso, &
640 : max_iso_not0, max_iso_not0_local, max_max_iso_not0, max_nso, max_s_harm, maxl, maxlgto, &
641 : maxso, mepos, n1s, n2s, na, natom, natom_tot, nkind, nr, nset, nspins, num_pe, size1, &
642 : size2
643 864 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list, dacg_n_list
644 864 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list, dacg_list
645 : INTEGER, DIMENSION(2) :: bo
646 864 : INTEGER, DIMENSION(:), POINTER :: atom_list, lmax, lmin, npgf, o2nindex
647 : LOGICAL :: paw_atom
648 864 : LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: is_set_to_0
649 : REAL(dp) :: hard_radius
650 864 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: g1, g2, gauge_h, gauge_s
651 864 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: cjc0_h_block, cjc0_s_block, cjc_h_block, &
652 864 : cjc_ii_h_block, cjc_ii_s_block, cjc_iii_h_block, cjc_iii_s_block, cjc_s_block, dgg_1, gg, &
653 864 : gg_lm1
654 864 : REAL(dp), DIMENSION(:, :), POINTER :: coeff, Fr_a_h, Fr_a_h_ii, Fr_a_h_iii, Fr_a_s, &
655 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, &
656 864 : Fr_h, Fr_s, zet
657 864 : REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG
658 864 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: my_CG_dxyz_asym
659 864 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
660 : TYPE(dft_control_type), POINTER :: dft_control
661 : TYPE(grid_atom_type), POINTER :: grid_atom
662 : TYPE(gto_basis_set_type), POINTER :: basis_1c_set
663 : TYPE(harmonics_atom_type), POINTER :: harmonics
664 864 : TYPE(jrho_atom_type), DIMENSION(:), POINTER :: jrho1_atom_set
665 : TYPE(jrho_atom_type), POINTER :: jrho1_atom
666 : TYPE(mp_para_env_type), POINTER :: para_env
667 864 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
668 :
669 864 : CALL timeset(routineN, handle)
670 : !
671 864 : NULLIFY (atomic_kind_set, qs_kind_set, dft_control, para_env, &
672 864 : coeff, Fr_h, Fr_s, Fr_a_h, Fr_a_s, Fr_a_h_ii, Fr_a_s_ii, &
673 864 : Fr_a_h_iii, Fr_a_s_iii, Fr_b_h, Fr_b_s, Fr_b_h_ii, &
674 864 : Fr_b_s_ii, Fr_b_h_iii, Fr_b_s_iii, jrho1_atom_set, &
675 864 : jrho1_atom)
676 : !
677 : CALL get_qs_env(qs_env=qs_env, &
678 : atomic_kind_set=atomic_kind_set, &
679 : qs_kind_set=qs_kind_set, &
680 : dft_control=dft_control, &
681 864 : para_env=para_env)
682 :
683 864 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxlgto=maxlgto)
684 :
685 : !
686 : CALL get_current_env(current_env=current_env, &
687 864 : jrho1_atom_set=jrho1_atom_set)
688 : !
689 :
690 864 : nkind = SIZE(qs_kind_set)
691 864 : nspins = dft_control%nspins
692 : !
693 864 : natom_tot = SIZE(jrho1_atom_set, 1)
694 3456 : ALLOCATE (is_set_to_0(natom_tot, nspins))
695 5994 : is_set_to_0(:, :) = .FALSE.
696 :
697 : !
698 2466 : DO ikind = 1, nkind
699 1602 : NULLIFY (atom_list, grid_atom, harmonics, basis_1c_set, &
700 1602 : lmax, lmin, npgf, zet, grid_atom, harmonics, my_CG, my_CG_dxyz_asym)
701 : !
702 : CALL get_atomic_kind(atomic_kind_set(ikind), &
703 : atom_list=atom_list, &
704 1602 : natom=natom)
705 : CALL get_qs_kind(qs_kind_set(ikind), &
706 : grid_atom=grid_atom, &
707 : paw_atom=paw_atom, &
708 : harmonics=harmonics, &
709 : hard_radius=hard_radius, &
710 : basis_set=basis_1c_set, &
711 1602 : basis_type="GAPW_1C")
712 : !
713 : ! Quick cycle if needed.
714 1602 : IF (.NOT. paw_atom) CYCLE
715 : !
716 : CALL get_gto_basis_set(gto_basis_set=basis_1c_set, &
717 : lmax=lmax, lmin=lmin, &
718 : maxl=maxl, npgf=npgf, &
719 : nset=nset, zet=zet, &
720 1242 : maxso=maxso)
721 1242 : CALL get_paw_basis_info(basis_1c_set, o2nindex=o2nindex)
722 : !
723 1242 : nr = grid_atom%nr
724 1242 : na = grid_atom%ng_sphere
725 1242 : max_iso_not0 = harmonics%max_iso_not0
726 1242 : damax_iso_not0 = harmonics%damax_iso_not0
727 1242 : max_max_iso_not0 = MAX(max_iso_not0, damax_iso_not0)
728 1242 : lmax_expansion = indso(1, max_max_iso_not0)
729 1242 : max_s_harm = harmonics%max_s_harm
730 1242 : llmax = harmonics%llmax
731 : !
732 : ! Distribute the atoms of this kind
733 1242 : num_pe = para_env%num_pe
734 1242 : mepos = para_env%mepos
735 1242 : bo = get_limit(natom, num_pe, mepos)
736 : !
737 1242 : my_CG => harmonics%my_CG
738 1242 : my_CG_dxyz_asym => harmonics%my_CG_dxyz_asym
739 : !
740 : ! Allocate some arrays.
741 1242 : max_nso = nsoset(maxl)
742 0 : ALLOCATE (g1(nr), g2(nr), gg(nr, 0:2*maxl), gg_lm1(nr, 0:2*maxl), dgg_1(nr, 0:2*maxl), &
743 0 : cjc0_h_block(max_nso, max_nso), cjc0_s_block(max_nso, max_nso), &
744 0 : cjc_h_block(max_nso, max_nso), cjc_s_block(max_nso, max_nso), &
745 0 : cjc_ii_h_block(max_nso, max_nso), cjc_ii_s_block(max_nso, max_nso), &
746 0 : cjc_iii_h_block(max_nso, max_nso), cjc_iii_s_block(max_nso, max_nso), &
747 0 : cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm), &
748 0 : dacg_list(2, nsoset(maxl)**2, max_s_harm), dacg_n_list(max_s_harm), &
749 47196 : gauge_h(nr), gauge_s(nr))
750 : !
751 : ! Compute the gauge
752 1296 : SELECT CASE (current_env%gauge)
753 : CASE (current_gauge_r)
754 : ! d(r)=r
755 2754 : gauge_h(1:nr) = grid_atom%rad(1:nr)
756 2754 : gauge_s(1:nr) = grid_atom%rad(1:nr)
757 : CASE (current_gauge_r_and_step_func)
758 : ! step function
759 43740 : gauge_h(1:nr) = 0e0_dp
760 43740 : DO ir = 1, nr
761 43740 : IF (grid_atom%rad(ir) <= hard_radius) THEN
762 24138 : gauge_s(ir) = grid_atom%rad(ir)
763 : ELSE
764 18702 : gauge_s(ir) = gauge_h(ir)
765 : END IF
766 : END DO
767 : CASE (current_gauge_atom)
768 : ! d(r)=A
769 15588 : gauge_h(1:nr) = HUGE(0e0_dp) !0e0_dp
770 15588 : gauge_s(1:nr) = HUGE(0e0_dp) !0e0_dp
771 : CASE DEFAULT
772 1242 : CPABORT("Unknown gauge, try again...")
773 : END SELECT
774 : !
775 : !
776 1242 : m1s = 0
777 4968 : DO iset1 = 1, nset
778 : m2s = 0
779 15300 : DO iset2 = 1, nset
780 : CALL get_none0_cg_list(my_CG, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
781 11574 : max_s_harm, lmax_expansion, cg_list, cg_n_list, max_iso_not0_local)
782 11574 : CPASSERT(max_iso_not0_local <= max_iso_not0)
783 : CALL get_none0_cg_list(my_CG_dxyz_asym, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
784 11574 : max_s_harm, lmax_expansion, dacg_list, dacg_n_list, damax_iso_not0_local)
785 11574 : CPASSERT(damax_iso_not0_local <= damax_iso_not0)
786 :
787 11574 : n1s = nsoset(lmax(iset1))
788 36504 : DO ipgf1 = 1, npgf(iset1)
789 24930 : iso1_first = nsoset(lmin(iset1) - 1) + 1 + n1s*(ipgf1 - 1) + m1s
790 24930 : iso1_last = nsoset(lmax(iset1)) + n1s*(ipgf1 - 1) + m1s
791 24930 : size1 = iso1_last - iso1_first + 1
792 24930 : iso1_first = o2nindex(iso1_first)
793 24930 : iso1_last = o2nindex(iso1_last)
794 24930 : i1 = iso1_last - iso1_first + 1
795 24930 : CPASSERT(size1 == i1)
796 24930 : i1 = nsoset(lmin(iset1) - 1) + 1
797 : !
798 1244430 : g1(1:nr) = EXP(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))
799 : !
800 24930 : n2s = nsoset(lmax(iset2))
801 91674 : DO ipgf2 = 1, npgf(iset2)
802 55170 : iso2_first = nsoset(lmin(iset2) - 1) + 1 + n2s*(ipgf2 - 1) + m2s
803 55170 : iso2_last = nsoset(lmax(iset2)) + n2s*(ipgf2 - 1) + m2s
804 55170 : size2 = iso2_last - iso2_first + 1
805 55170 : iso2_first = o2nindex(iso2_first)
806 55170 : iso2_last = o2nindex(iso2_last)
807 55170 : i2 = iso2_last - iso2_first + 1
808 55170 : CPASSERT(size2 == i2)
809 55170 : i2 = nsoset(lmin(iset2) - 1) + 1
810 : !
811 2730510 : g2(1:nr) = EXP(-zet(ipgf2, iset2)*grid_atom%rad2(1:nr))
812 : !
813 55170 : lmin12 = lmin(iset1) + lmin(iset2)
814 55170 : lmax12 = lmax(iset1) + lmax(iset2)
815 : !
816 10133028 : gg = 0.0_dp
817 10133028 : gg_lm1 = 0.0_dp
818 10133028 : dgg_1 = 0.0_dp
819 : !
820 : ! Take only the terms of expansion for L < lmax_expansion
821 55170 : IF (lmin12 <= lmax_expansion) THEN
822 : !
823 55170 : IF (lmax12 > lmax_expansion) lmax12 = lmax_expansion
824 : !
825 55170 : IF (lmin12 == 0) THEN
826 2510514 : gg(1:nr, lmin12) = g1(1:nr)*g2(1:nr)
827 2510514 : gg_lm1(1:nr, lmin12) = 0.0_dp
828 : ELSE
829 219996 : gg(1:nr, lmin12) = grid_atom%rad2l(1:nr, lmin12)*g1(1:nr)*g2(1:nr)
830 219996 : gg_lm1(1:nr, lmin12) = grid_atom%rad2l(1:nr, lmin12 - 1)*g1(1:nr)*g2(1:nr)
831 : END IF
832 : !
833 103122 : DO l = lmin12 + 1, lmax12
834 2363472 : gg(1:nr, l) = grid_atom%rad(1:nr)*gg(1:nr, l - 1)
835 2418642 : gg_lm1(1:nr, l) = gg(1:nr, l - 1)
836 : END DO
837 : !
838 158292 : DO l = lmin12, lmax12
839 : dgg_1(1:nr, l) = 2.0_dp*(zet(ipgf1, iset1) - zet(ipgf2, iset2))&
840 5149152 : & *gg(1:nr, l)*grid_atom%rad(1:nr)
841 : END DO
842 : ELSE
843 : CYCLE
844 : END IF ! lmin12
845 : !
846 118251 : DO iat = bo(1), bo(2)
847 38151 : iatom = atom_list(iat)
848 : !
849 155052 : DO ispin = 1, nspins
850 : !------------------------------------------------------------------
851 : ! P_\alpha\alpha'
852 3125673 : cjc0_h_block = HUGE(1.0_dp)
853 3125673 : cjc0_s_block = HUGE(1.0_dp)
854 : !
855 : ! Hard term
856 61731 : coeff => jrho1_atom_set(iatom)%cjc0_h(ispin)%r_coef
857 : cjc0_h_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = &
858 601002 : coeff(iso1_first:iso1_last, iso2_first:iso2_last)
859 : !
860 : ! Soft term
861 61731 : coeff => jrho1_atom_set(iatom)%cjc0_s(ispin)%r_coef
862 : cjc0_s_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = &
863 601002 : coeff(iso1_first:iso1_last, iso2_first:iso2_last)
864 : !------------------------------------------------------------------
865 : ! mQai_\alpha\alpha'
866 3125673 : cjc_h_block = HUGE(1.0_dp)
867 3125673 : cjc_s_block = HUGE(1.0_dp)
868 : !
869 : ! Hard term
870 61731 : coeff => jrho1_atom_set(iatom)%cjc_h(ispin)%r_coef
871 : cjc_h_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = &
872 601002 : coeff(iso1_first:iso1_last, iso2_first:iso2_last)
873 : !
874 : ! Soft term
875 61731 : coeff => jrho1_atom_set(iatom)%cjc_s(ispin)%r_coef
876 : cjc_s_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = &
877 601002 : coeff(iso1_first:iso1_last, iso2_first:iso2_last)
878 : !------------------------------------------------------------------
879 : ! Qci_\alpha\alpha'
880 3125673 : cjc_ii_h_block = HUGE(1.0_dp)
881 3125673 : cjc_ii_s_block = HUGE(1.0_dp)
882 : !
883 : ! Hard term
884 61731 : coeff => jrho1_atom_set(iatom)%cjc_ii_h(ispin)%r_coef
885 : cjc_ii_h_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = &
886 601002 : coeff(iso1_first:iso1_last, iso2_first:iso2_last)
887 : !
888 : ! Soft term
889 61731 : coeff => jrho1_atom_set(iatom)%cjc_ii_s(ispin)%r_coef
890 : cjc_ii_s_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = &
891 601002 : coeff(iso1_first:iso1_last, iso2_first:iso2_last)
892 : !------------------------------------------------------------------
893 : ! Qbi_\alpha\alpha'
894 3125673 : cjc_iii_h_block = HUGE(1.0_dp)
895 3125673 : cjc_iii_s_block = HUGE(1.0_dp)
896 : !
897 : !
898 : ! Hard term
899 61731 : coeff => jrho1_atom_set(iatom)%cjc_iii_h(ispin)%r_coef
900 : cjc_iii_h_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = &
901 601002 : coeff(iso1_first:iso1_last, iso2_first:iso2_last)
902 : !
903 : ! Soft term
904 61731 : coeff => jrho1_atom_set(iatom)%cjc_iii_s(ispin)%r_coef
905 : cjc_iii_s_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = &
906 601002 : coeff(iso1_first:iso1_last, iso2_first:iso2_last)
907 : !------------------------------------------------------------------
908 : !
909 : ! Allocation radial functions
910 61731 : jrho1_atom => jrho1_atom_set(iatom)
911 61731 : IF (.NOT. ASSOCIATED(jrho1_atom%jrho_a_h(ispin)%r_coef)) THEN
912 : CALL allocate_jrho_atom_rad(jrho1_atom, ispin, nr, na, &
913 147 : max_max_iso_not0)
914 147 : is_set_to_0(iatom, ispin) = .TRUE.
915 : ELSE
916 61584 : IF (.NOT. is_set_to_0(iatom, ispin)) THEN
917 1176 : CALL set2zero_jrho_atom_rad(jrho1_atom, ispin)
918 1176 : is_set_to_0(iatom, ispin) = .TRUE.
919 : END IF
920 : END IF
921 : !------------------------------------------------------------------
922 : !
923 61731 : Fr_h => jrho1_atom%jrho_h(ispin)%r_coef
924 61731 : Fr_s => jrho1_atom%jrho_s(ispin)%r_coef
925 : !------------------------------------------------------------------
926 : !
927 61731 : Fr_a_h => jrho1_atom%jrho_a_h(ispin)%r_coef
928 61731 : Fr_a_s => jrho1_atom%jrho_a_s(ispin)%r_coef
929 61731 : Fr_b_h => jrho1_atom%jrho_b_h(ispin)%r_coef
930 61731 : Fr_b_s => jrho1_atom%jrho_b_s(ispin)%r_coef
931 : !------------------------------------------------------------------
932 : !
933 61731 : Fr_a_h_ii => jrho1_atom%jrho_a_h_ii(ispin)%r_coef
934 61731 : Fr_a_s_ii => jrho1_atom%jrho_a_s_ii(ispin)%r_coef
935 61731 : Fr_b_h_ii => jrho1_atom%jrho_b_h_ii(ispin)%r_coef
936 61731 : Fr_b_s_ii => jrho1_atom%jrho_b_s_ii(ispin)%r_coef
937 : !------------------------------------------------------------------
938 : !
939 61731 : Fr_a_h_iii => jrho1_atom%jrho_a_h_iii(ispin)%r_coef
940 61731 : Fr_a_s_iii => jrho1_atom%jrho_a_s_iii(ispin)%r_coef
941 61731 : Fr_b_h_iii => jrho1_atom%jrho_b_h_iii(ispin)%r_coef
942 61731 : Fr_b_s_iii => jrho1_atom%jrho_b_s_iii(ispin)%r_coef
943 : !------------------------------------------------------------------
944 : !
945 362160 : DO iso = 1, max_iso_not0_local
946 300429 : l_iso = indso(1, iso) ! not needed
947 300429 : m_iso = indso(2, iso) ! not needed
948 : !
949 858771 : DO icg = 1, cg_n_list(iso)
950 : !
951 496611 : iso1 = cg_list(1, icg, iso)
952 496611 : iso2 = cg_list(2, icg, iso)
953 : !
954 496611 : IF (.NOT. (iso2 > 0 .AND. iso1 > 0)) THEN
955 0 : WRITE (*, *) 'iso1=', iso1, ' iso2=', iso2, ' iso=', iso, ' icg=', icg
956 0 : WRITE (*, *) '.... will stop!'
957 : END IF
958 496611 : CPASSERT(iso2 > 0 .AND. iso1 > 0)
959 : !
960 496611 : l = indso(1, iso1) + indso(1, iso2)
961 496611 : IF (l > lmax_expansion .OR. l < .0) THEN
962 0 : WRITE (*, *) 'calculate_jrho_atom_rad: 1 l', l
963 0 : WRITE (*, *) 'calculate_jrho_atom_rad: 1 lmax_expansion', lmax_expansion
964 0 : WRITE (*, *) '.... will stop!'
965 : END IF
966 496611 : CPASSERT(l <= lmax_expansion)
967 : !------------------------------------------------------------------
968 : ! P0
969 : !
970 496611 : IF (current_env%gauge == current_gauge_atom) THEN
971 : ! Hard term
972 : Fr_h(1:nr, iso) = Fr_h(1:nr, iso) + &
973 : gg(1:nr, l)*cjc0_h_block(iso1, iso2)* &
974 4477842 : my_CG(iso1, iso2, iso)
975 : ! Soft term
976 : Fr_s(1:nr, iso) = Fr_s(1:nr, iso) + &
977 : gg(1:nr, l)*cjc0_s_block(iso1, iso2)* &
978 4477842 : my_CG(iso1, iso2, iso)
979 : ELSE
980 : ! Hard term
981 : Fr_h(1:nr, iso) = Fr_h(1:nr, iso) + &
982 : gg(1:nr, l)*cjc0_h_block(iso1, iso2)* &
983 39661029 : my_CG(iso1, iso2, iso)*(grid_atom%rad(1:nr) - gauge_h(1:nr))
984 : ! Soft term
985 : Fr_s(1:nr, iso) = Fr_s(1:nr, iso) + &
986 : gg(1:nr, l)*cjc0_s_block(iso1, iso2)* &
987 39661029 : my_CG(iso1, iso2, iso)*(grid_atom%rad(1:nr) - gauge_s(1:nr))
988 : END IF
989 : !------------------------------------------------------------------
990 : ! Rai
991 : !
992 : ! Hard term
993 : Fr_a_h(1:nr, iso) = Fr_a_h(1:nr, iso) + &
994 : dgg_1(1:nr, l)*cjc_h_block(iso1, iso2)* &
995 24512841 : my_CG(iso1, iso2, iso)
996 : !
997 : ! Soft term
998 : Fr_a_s(1:nr, iso) = Fr_a_s(1:nr, iso) + &
999 : dgg_1(1:nr, l)*cjc_s_block(iso1, iso2)* &
1000 24512841 : my_CG(iso1, iso2, iso)
1001 : !------------------------------------------------------------------
1002 : ! Rci
1003 : !
1004 496611 : IF (current_env%gauge == current_gauge_atom) THEN
1005 : ! Hard term
1006 : Fr_a_h_ii(1:nr, iso) = Fr_a_h_ii(1:nr, iso) + &
1007 : dgg_1(1:nr, l)* &
1008 : cjc_ii_h_block(iso1, iso2)* &
1009 4477842 : my_CG(iso1, iso2, iso)
1010 : ! Soft term
1011 : Fr_a_s_ii(1:nr, iso) = Fr_a_s_ii(1:nr, iso) + &
1012 : dgg_1(1:nr, l)* &
1013 : cjc_ii_s_block(iso1, iso2)* &
1014 4477842 : my_CG(iso1, iso2, iso)
1015 : ELSE
1016 : ! Hard term
1017 : Fr_a_h_ii(1:nr, iso) = Fr_a_h_ii(1:nr, iso) + &
1018 : dgg_1(1:nr, l)*gauge_h(1:nr)* &
1019 : cjc_ii_h_block(iso1, iso2)* &
1020 20034999 : my_CG(iso1, iso2, iso)
1021 : ! Soft term
1022 : Fr_a_s_ii(1:nr, iso) = Fr_a_s_ii(1:nr, iso) + &
1023 : dgg_1(1:nr, l)*gauge_s(1:nr)* &
1024 : cjc_ii_s_block(iso1, iso2)* &
1025 20034999 : my_CG(iso1, iso2, iso)
1026 : END IF
1027 : !------------------------------------------------------------------
1028 : ! Rbi
1029 : !
1030 797040 : IF (current_env%gauge == current_gauge_atom) THEN
1031 : ! Hard term
1032 : Fr_a_h_iii(1:nr, iso) = Fr_a_h_iii(1:nr, iso) + &
1033 : dgg_1(1:nr, l)* &
1034 : cjc_iii_h_block(iso1, iso2)* &
1035 4477842 : my_CG(iso1, iso2, iso)
1036 : ! Soft term
1037 : Fr_a_s_iii(1:nr, iso) = Fr_a_s_iii(1:nr, iso) + &
1038 : dgg_1(1:nr, l)* &
1039 : cjc_iii_s_block(iso1, iso2)* &
1040 4477842 : my_CG(iso1, iso2, iso)
1041 : ELSE
1042 : ! Hard term
1043 : Fr_a_h_iii(1:nr, iso) = Fr_a_h_iii(1:nr, iso) + &
1044 : dgg_1(1:nr, l)*gauge_h(1:nr)* &
1045 : cjc_iii_h_block(iso1, iso2)* &
1046 20034999 : my_CG(iso1, iso2, iso)
1047 : ! Soft term
1048 : Fr_a_s_iii(1:nr, iso) = Fr_a_s_iii(1:nr, iso) + &
1049 : dgg_1(1:nr, l)*gauge_s(1:nr)* &
1050 : cjc_iii_s_block(iso1, iso2)* &
1051 20034999 : my_CG(iso1, iso2, iso)
1052 : END IF
1053 : !------------------------------------------------------------------
1054 : END DO !icg
1055 : !
1056 : END DO ! iso
1057 : !
1058 : !
1059 212436 : DO iso = 1, damax_iso_not0_local
1060 112554 : l_iso = indso(1, iso) ! not needed
1061 112554 : m_iso = indso(2, iso) ! not needed
1062 : !
1063 669321 : DO icg = 1, dacg_n_list(iso)
1064 : !
1065 495036 : iso1 = dacg_list(1, icg, iso)
1066 495036 : iso2 = dacg_list(2, icg, iso)
1067 : !
1068 495036 : IF (.NOT. (iso2 > 0 .AND. iso1 > 0)) THEN
1069 0 : WRITE (*, *) 'iso1=', iso1, ' iso2=', iso2, ' iso=', iso, ' icg=', icg
1070 0 : WRITE (*, *) '.... will stop!'
1071 : END IF
1072 495036 : CPASSERT(iso2 > 0 .AND. iso1 > 0)
1073 : !
1074 495036 : l = indso(1, iso1) + indso(1, iso2)
1075 495036 : IF (l > lmax_expansion) THEN
1076 0 : WRITE (*, *) 'calculate_jrho_atom_rad: 1 l', l
1077 0 : WRITE (*, *) 'calculate_jrho_atom_rad: 1 lmax_expansion', lmax_expansion
1078 0 : WRITE (*, *) '.... will stop!'
1079 : END IF
1080 495036 : CPASSERT(l <= lmax_expansion)
1081 : !------------------------------------------------------------------
1082 : ! Daij
1083 : !
1084 : ! Hard term
1085 : Fr_b_h(1:nr, iso) = Fr_b_h(1:nr, iso) + &
1086 : gg_lm1(1:nr, l)*cjc_h_block(iso1, iso2)* &
1087 24139836 : my_CG_dxyz_asym(idir, iso1, iso2, iso)
1088 : !
1089 : ! Soft term
1090 : Fr_b_s(1:nr, iso) = Fr_b_s(1:nr, iso) + &
1091 : gg_lm1(1:nr, l)*cjc_s_block(iso1, iso2)* &
1092 24139836 : my_CG_dxyz_asym(idir, iso1, iso2, iso)
1093 : !
1094 : !------------------------------------------------------------------
1095 : ! Dcij
1096 : !
1097 495036 : IF (current_env%gauge == current_gauge_atom) THEN
1098 : ! Hard term
1099 : Fr_b_h_ii(1:nr, iso) = Fr_b_h_ii(1:nr, iso) + &
1100 : gg_lm1(1:nr, l)* &
1101 : cjc_ii_h_block(iso1, iso2)* &
1102 3569184 : my_CG_dxyz_asym(idir, iso1, iso2, iso)
1103 : ! Soft term
1104 : Fr_b_s_ii(1:nr, iso) = Fr_b_s_ii(1:nr, iso) + &
1105 : gg_lm1(1:nr, l)* &
1106 : cjc_ii_s_block(iso1, iso2)* &
1107 3569184 : my_CG_dxyz_asym(idir, iso1, iso2, iso)
1108 : ELSE
1109 : ! Hard term
1110 : Fr_b_h_ii(1:nr, iso) = Fr_b_h_ii(1:nr, iso) + &
1111 : gg_lm1(1:nr, l)*gauge_h(1:nr)* &
1112 : cjc_ii_h_block(iso1, iso2)* &
1113 20570652 : my_CG_dxyz_asym(idir, iso1, iso2, iso)
1114 : ! Soft term
1115 : Fr_b_s_ii(1:nr, iso) = Fr_b_s_ii(1:nr, iso) + &
1116 : gg_lm1(1:nr, l)*gauge_s(1:nr)* &
1117 : cjc_ii_s_block(iso1, iso2)* &
1118 20570652 : my_CG_dxyz_asym(idir, iso1, iso2, iso)
1119 : END IF
1120 : !------------------------------------------------------------------
1121 : ! Dbij
1122 : !
1123 607590 : IF (current_env%gauge == current_gauge_atom) THEN
1124 : ! Hard term
1125 : Fr_b_h_iii(1:nr, iso) = Fr_b_h_iii(1:nr, iso) + &
1126 : gg_lm1(1:nr, l)* &
1127 : cjc_iii_h_block(iso1, iso2)* &
1128 3569184 : my_CG_dxyz_asym(idir, iso1, iso2, iso)
1129 : ! Soft term
1130 : Fr_b_s_iii(1:nr, iso) = Fr_b_s_iii(1:nr, iso) + &
1131 : gg_lm1(1:nr, l)* &
1132 : cjc_iii_s_block(iso1, iso2)* &
1133 3569184 : my_CG_dxyz_asym(idir, iso1, iso2, iso)
1134 : ELSE
1135 : ! Hard term
1136 : Fr_b_h_iii(1:nr, iso) = Fr_b_h_iii(1:nr, iso) + &
1137 : gg_lm1(1:nr, l)*gauge_h(1:nr)* &
1138 : cjc_iii_h_block(iso1, iso2)* &
1139 20570652 : my_CG_dxyz_asym(idir, iso1, iso2, iso)
1140 : ! Soft term
1141 : Fr_b_s_iii(1:nr, iso) = Fr_b_s_iii(1:nr, iso) + &
1142 : gg_lm1(1:nr, l)*gauge_s(1:nr)* &
1143 : cjc_iii_s_block(iso1, iso2)* &
1144 20570652 : my_CG_dxyz_asym(idir, iso1, iso2, iso)
1145 : END IF
1146 : !------------------------------------------------------------------
1147 : END DO ! icg
1148 : END DO ! iso
1149 : !
1150 : END DO ! ispin
1151 : END DO ! iat
1152 : !
1153 : !------------------------------------------------------------------
1154 : !
1155 : END DO !ipgf2
1156 : END DO ! ipgf1
1157 38448 : m2s = m2s + maxso
1158 : END DO ! iset2
1159 4968 : m1s = m1s + maxso
1160 : END DO ! iset1
1161 : !
1162 0 : DEALLOCATE (cjc0_h_block, cjc0_s_block, cjc_h_block, cjc_s_block, cjc_ii_h_block, cjc_ii_s_block, &
1163 0 : cjc_iii_h_block, cjc_iii_s_block, g1, g2, gg, gg_lm1, dgg_1, gauge_h, gauge_s, &
1164 1242 : cg_list, cg_n_list, dacg_list, dacg_n_list)
1165 5310 : DEALLOCATE (o2nindex)
1166 : END DO ! ikind
1167 : !
1168 : !
1169 864 : DEALLOCATE (is_set_to_0)
1170 : !
1171 864 : CALL timestop(handle)
1172 : !
1173 2592 : END SUBROUTINE calculate_jrho_atom_rad
1174 :
1175 : ! **************************************************************************************************
1176 : !> \brief ...
1177 : !> \param jrho1_atom ...
1178 : !> \param jrho_h ...
1179 : !> \param jrho_s ...
1180 : !> \param grid_atom ...
1181 : !> \param harmonics ...
1182 : !> \param do_igaim ...
1183 : !> \param ratom ...
1184 : !> \param natm_gauge ...
1185 : !> \param iB ...
1186 : !> \param idir ...
1187 : !> \param ispin ...
1188 : ! **************************************************************************************************
1189 1323 : SUBROUTINE calculate_jrho_atom_ang(jrho1_atom, jrho_h, jrho_s, grid_atom, &
1190 1323 : harmonics, do_igaim, ratom, natm_gauge, &
1191 : iB, idir, ispin)
1192 : !
1193 : TYPE(jrho_atom_type), POINTER :: jrho1_atom
1194 : REAL(dp), DIMENSION(:, :), POINTER :: jrho_h, jrho_s
1195 : TYPE(grid_atom_type), POINTER :: grid_atom
1196 : TYPE(harmonics_atom_type), POINTER :: harmonics
1197 : LOGICAL, INTENT(IN) :: do_igaim
1198 : INTEGER, INTENT(IN) :: iB, idir, ispin, natm_gauge
1199 : REAL(dp), INTENT(IN) :: ratom(:, :)
1200 :
1201 : INTEGER :: ia, idir2, iiB, iiiB, ir, &
1202 : iso, max_max_iso_not0, na, nr
1203 : REAL(dp) :: rad_part, scale
1204 1323 : REAL(dp), DIMENSION(:, :), POINTER :: a, Fr_a_h, Fr_a_h_ii, Fr_a_h_iii, &
1205 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, &
1206 1323 : Fr_b_s_ii, Fr_b_s_iii, Fr_h, Fr_s, slm
1207 1323 : REAL(dp), DIMENSION(:), POINTER :: r
1208 : REAL(dp), DIMENSION(:, :, :), ALLOCATABLE :: g
1209 : !
1210 : !
1211 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, &
1212 1323 : Fr_b_h, Fr_b_s, Fr_b_h_ii, Fr_b_s_ii, Fr_b_h_iii, Fr_b_s_iii, &
1213 1323 : a, slm)
1214 : !
1215 0 : CPASSERT(ASSOCIATED(jrho_h))
1216 1323 : CPASSERT(ASSOCIATED(jrho_s))
1217 1323 : CPASSERT(ASSOCIATED(jrho1_atom))
1218 : ! just to be sure...
1219 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_h))
1220 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_s))
1221 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_h))
1222 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_s))
1223 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_h(ispin)%r_coef))
1224 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_s(ispin)%r_coef))
1225 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_h(ispin)%r_coef))
1226 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_s(ispin)%r_coef))
1227 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_h_ii))
1228 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_s_ii))
1229 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_h_ii))
1230 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_s_ii))
1231 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_h_ii(ispin)%r_coef))
1232 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_s_ii(ispin)%r_coef))
1233 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_h_ii(ispin)%r_coef))
1234 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_s_ii(ispin)%r_coef))
1235 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_h_iii))
1236 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_s_iii))
1237 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_h_iii))
1238 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_s_iii))
1239 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_h_iii(ispin)%r_coef))
1240 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_s_iii(ispin)%r_coef))
1241 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_h_iii(ispin)%r_coef))
1242 1323 : CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_s_iii(ispin)%r_coef))
1243 : !
1244 : !
1245 1323 : nr = grid_atom%nr
1246 1323 : na = grid_atom%ng_sphere
1247 1323 : max_max_iso_not0 = MAX(harmonics%max_iso_not0, harmonics%damax_iso_not0)
1248 5292 : ALLOCATE (g(3, nr, na))
1249 : !------------------------------------------------------------------
1250 : !
1251 1323 : Fr_h => jrho1_atom%jrho_h(ispin)%r_coef
1252 1323 : Fr_s => jrho1_atom%jrho_s(ispin)%r_coef
1253 : !------------------------------------------------------------------
1254 : !
1255 1323 : Fr_a_h => jrho1_atom%jrho_a_h(ispin)%r_coef !Rai
1256 1323 : Fr_a_s => jrho1_atom%jrho_a_s(ispin)%r_coef
1257 1323 : Fr_b_h => jrho1_atom%jrho_b_h(ispin)%r_coef !Daij
1258 1323 : Fr_b_s => jrho1_atom%jrho_b_s(ispin)%r_coef
1259 : !------------------------------------------------------------------
1260 : !
1261 1323 : Fr_a_h_ii => jrho1_atom%jrho_a_h_ii(ispin)%r_coef !Rci
1262 1323 : Fr_a_s_ii => jrho1_atom%jrho_a_s_ii(ispin)%r_coef
1263 1323 : Fr_b_h_ii => jrho1_atom%jrho_b_h_ii(ispin)%r_coef !Dcij
1264 1323 : Fr_b_s_ii => jrho1_atom%jrho_b_s_ii(ispin)%r_coef
1265 : !------------------------------------------------------------------
1266 : !
1267 1323 : Fr_a_h_iii => jrho1_atom%jrho_a_h_iii(ispin)%r_coef !Rbi
1268 1323 : Fr_a_s_iii => jrho1_atom%jrho_a_s_iii(ispin)%r_coef
1269 1323 : Fr_b_h_iii => jrho1_atom%jrho_b_h_iii(ispin)%r_coef !Dbij
1270 1323 : Fr_b_s_iii => jrho1_atom%jrho_b_s_iii(ispin)%r_coef
1271 : !------------------------------------------------------------------
1272 : !
1273 1323 : a => harmonics%a
1274 1323 : slm => harmonics%slm
1275 1323 : r => grid_atom%rad
1276 : !
1277 1323 : CALL set_vecp(iB, iiB, iiiB)
1278 : !
1279 1323 : scale = 0.0_dp
1280 1323 : idir2 = 1
1281 1323 : IF (idir /= iB) THEN
1282 882 : CALL set_vecp_rev(idir, iB, idir2)
1283 882 : scale = fac_vecp(idir, iB, idir2)
1284 : END IF
1285 : !
1286 : ! Set the gauge
1287 1323 : CALL get_gauge()
1288 : !
1289 65943 : DO ir = 1, nr
1290 820323 : DO iso = 1, max_max_iso_not0
1291 38013120 : DO ia = 1, na
1292 37948500 : IF (do_igaim) THEN
1293 : !------------------------------------------------------------------
1294 : ! Hard current density response
1295 : ! radial(ia,ir) = ( aj(ia) * Rai(ir,iso) + Daij
1296 : ! - aii(ia) * ( aj(ia) * Rbi(ir,iso) + Dbij )
1297 : ! + aiii(ia) * ( aj(ia) * Rci(ir,iso) + Dcij )
1298 : ! ) * Ylm(ia)
1299 : rad_part = a(idir, ia)*Fr_a_h(ir, iso) + Fr_b_h(ir, iso) &
1300 : & - g(iiB, ir, ia)*(a(idir, ia)*Fr_a_h_iii(ir, iso) + Fr_b_h_iii(ir, iso))&
1301 : & + g(iiiB, ir, ia)*(a(idir, ia)*Fr_a_h_ii(ir, iso) + Fr_b_h_ii(ir, iso))&
1302 7380000 : & + scale*(a(idir2, ia)*r(ir) - g(idir2, ir, ia))*Fr_h(ir, iso)
1303 : !
1304 7380000 : jrho_h(ir, ia) = jrho_h(ir, ia) + rad_part*slm(ia, iso)
1305 : !------------------------------------------------------------------
1306 : ! Soft current density response
1307 : rad_part = a(idir, ia)*Fr_a_s(ir, iso) + Fr_b_s(ir, iso) &
1308 : & - g(iiB, ir, ia)*(a(idir, ia)*Fr_a_s_iii(ir, iso) + Fr_b_s_iii(ir, iso))&
1309 : & + g(iiiB, ir, ia)*(a(idir, ia)*Fr_a_s_ii(ir, iso) + Fr_b_s_ii(ir, iso))&
1310 7380000 : & + scale*(a(idir2, ia)*r(ir) - g(idir2, ir, ia))*Fr_s(ir, iso)
1311 : !
1312 7380000 : jrho_s(ir, ia) = jrho_s(ir, ia) + rad_part*slm(ia, iso)
1313 : !------------------------------------------------------------------
1314 : ELSE
1315 : !------------------------------------------------------------------
1316 : ! Hard current density response
1317 : ! radial(ia,ir) = ( aj(ia) * Rai(ir,iso) + Daij
1318 : ! - aii(ia) * ( aj(ia) * Rbi(ir,iso) + Dbij )
1319 : ! + aiii(ia) * ( aj(ia) * Rci(ir,iso) + Dcij )
1320 : ! ) * Ylm(ia)
1321 : rad_part = a(idir, ia)*Fr_a_h(ir, iso) + Fr_b_h(ir, iso) &
1322 : & - a(iiB, ia)*(a(idir, ia)*Fr_a_h_iii(ir, iso) + Fr_b_h_iii(ir, iso))&
1323 : & + a(iiiB, ia)*(a(idir, ia)*Fr_a_h_ii(ir, iso) + Fr_b_h_ii(ir, iso))&
1324 29814120 : & + scale*a(idir2, ia)*Fr_h(ir, iso)
1325 : !
1326 29814120 : jrho_h(ir, ia) = jrho_h(ir, ia) + rad_part*slm(ia, iso)
1327 : !------------------------------------------------------------------
1328 : ! Soft current density response
1329 : rad_part = a(idir, ia)*Fr_a_s(ir, iso) + Fr_b_s(ir, iso) &
1330 : & - a(iiB, ia)*(a(idir, ia)*Fr_a_s_iii(ir, iso) + Fr_b_s_iii(ir, iso))&
1331 : & + a(iiiB, ia)*(a(idir, ia)*Fr_a_s_ii(ir, iso) + Fr_b_s_ii(ir, iso))&
1332 29814120 : & + scale*a(idir2, ia)*Fr_s(ir, iso)
1333 : !
1334 29814120 : jrho_s(ir, ia) = jrho_s(ir, ia) + rad_part*slm(ia, iso)
1335 : !------------------------------------------------------------------
1336 : END IF
1337 : END DO ! ia
1338 : END DO ! iso
1339 : END DO ! ir
1340 : !
1341 3969 : DEALLOCATE (g)
1342 : !
1343 : CONTAINS
1344 : !
1345 : ! **************************************************************************************************
1346 : !> \brief ...
1347 : ! **************************************************************************************************
1348 1323 : SUBROUTINE get_gauge()
1349 : INTEGER :: iatom, ixyz, jatom
1350 : REAL(dp) :: ab, pa, pb, point(3), pra(3), prb(3), &
1351 : res, tmp
1352 1323 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: buf
1353 :
1354 3969 : ALLOCATE (buf(natm_gauge))
1355 65943 : DO ir = 1, nr
1356 3238623 : DO ia = 1, na
1357 12690720 : DO ixyz = 1, 3
1358 12690720 : g(ixyz, ir, ia) = 0.0_dp
1359 : END DO
1360 3172680 : point(1) = r(ir)*a(1, ia)
1361 3172680 : point(2) = r(ir)*a(2, ia)
1362 3172680 : point(3) = r(ir)*a(3, ia)
1363 10785600 : DO iatom = 1, natm_gauge
1364 7612920 : buf(iatom) = 1.0_dp
1365 30451680 : pra = point - ratom(:, iatom)
1366 7612920 : pa = SQRT(pra(1)**2 + pra(2)**2 + pra(3)**2)
1367 31404240 : DO jatom = 1, natm_gauge
1368 20618640 : IF (iatom == jatom) CYCLE
1369 52022880 : prb = point - ratom(:, jatom)
1370 13005720 : pb = SQRT(prb(1)**2 + prb(2)**2 + prb(3)**2)
1371 13005720 : ab = SQRT((pra(1) - prb(1))**2 + (pra(2) - prb(2))**2 + (pra(3) - prb(3))**2)
1372 : !
1373 13005720 : tmp = (pa - pb)/ab
1374 13005720 : tmp = 0.5_dp*(3.0_dp - tmp**2)*tmp
1375 13005720 : tmp = 0.5_dp*(3.0_dp - tmp**2)*tmp
1376 13005720 : tmp = 0.5_dp*(3.0_dp - tmp**2)*tmp
1377 28231560 : buf(iatom) = buf(iatom)*0.5_dp*(1.0_dp - tmp)
1378 : END DO
1379 : END DO
1380 12755340 : DO ixyz = 1, 3
1381 9518040 : res = 0.0_dp
1382 32356800 : DO iatom = 1, natm_gauge
1383 32356800 : res = res + ratom(ixyz, iatom)*buf(iatom)
1384 : END DO
1385 32356800 : res = res/SUM(buf(1:natm_gauge))
1386 : !
1387 12690720 : g(ixyz, ir, ia) = res
1388 : END DO
1389 : END DO
1390 : END DO
1391 1323 : DEALLOCATE (buf)
1392 1323 : END SUBROUTINE get_gauge
1393 : END SUBROUTINE calculate_jrho_atom_ang
1394 :
1395 : ! **************************************************************************************************
1396 : !> \brief ...
1397 : !> \param current_env ...
1398 : !> \param qs_env ...
1399 : !> \param iB ...
1400 : !> \param idir ...
1401 : ! **************************************************************************************************
1402 864 : SUBROUTINE calculate_jrho_atom(current_env, qs_env, iB, idir)
1403 : TYPE(current_env_type) :: current_env
1404 : TYPE(qs_environment_type), POINTER :: qs_env
1405 : INTEGER, INTENT(IN) :: iB, idir
1406 :
1407 : INTEGER :: iat, iatom, ikind, ispin, jatom, &
1408 : natm_gauge, natm_tot, natom, nkind, &
1409 : nspins
1410 : INTEGER, DIMENSION(2) :: bo
1411 864 : INTEGER, DIMENSION(:), POINTER :: atom_list
1412 : LOGICAL :: do_igaim, gapw, paw_atom
1413 : REAL(dp) :: hard_radius, r(3)
1414 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: ratom
1415 864 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1416 : TYPE(cell_type), POINTER :: cell
1417 : TYPE(mp_para_env_type), POINTER :: para_env
1418 : TYPE(dft_control_type), POINTER :: dft_control
1419 : TYPE(grid_atom_type), POINTER :: grid_atom
1420 : TYPE(harmonics_atom_type), POINTER :: harmonics
1421 864 : TYPE(jrho_atom_type), DIMENSION(:), POINTER :: jrho1_atom_set
1422 : TYPE(jrho_atom_type), POINTER :: jrho1_atom
1423 864 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1424 864 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1425 :
1426 864 : NULLIFY (para_env, dft_control)
1427 864 : NULLIFY (jrho1_atom_set, grid_atom, harmonics)
1428 864 : NULLIFY (atomic_kind_set, qs_kind_set, atom_list)
1429 :
1430 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, &
1431 : atomic_kind_set=atomic_kind_set, &
1432 : qs_kind_set=qs_kind_set, &
1433 : particle_set=particle_set, &
1434 : cell=cell, &
1435 864 : para_env=para_env)
1436 :
1437 : CALL get_current_env(current_env=current_env, &
1438 864 : jrho1_atom_set=jrho1_atom_set)
1439 :
1440 864 : do_igaim = .FALSE.
1441 864 : IF (current_env%gauge == current_gauge_atom) do_igaim = .TRUE.
1442 :
1443 864 : gapw = dft_control%qs_control%gapw
1444 864 : nkind = SIZE(qs_kind_set, 1)
1445 864 : nspins = dft_control%nspins
1446 :
1447 864 : natm_tot = SIZE(particle_set)
1448 2592 : ALLOCATE (ratom(3, natm_tot))
1449 :
1450 864 : IF (gapw) THEN
1451 2466 : DO ikind = 1, nkind
1452 1602 : NULLIFY (atom_list, grid_atom, harmonics)
1453 1602 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
1454 : CALL get_qs_kind(qs_kind_set(ikind), &
1455 : grid_atom=grid_atom, &
1456 : harmonics=harmonics, &
1457 : hard_radius=hard_radius, &
1458 1602 : paw_atom=paw_atom)
1459 :
1460 1602 : IF (.NOT. paw_atom) CYCLE
1461 :
1462 : ! Distribute the atoms of this kind
1463 :
1464 1242 : bo = get_limit(natom, para_env%num_pe, para_env%mepos)
1465 :
1466 4554 : DO iat = bo(1), bo(2)
1467 846 : iatom = atom_list(iat)
1468 : NULLIFY (jrho1_atom)
1469 846 : jrho1_atom => jrho1_atom_set(iatom)
1470 :
1471 846 : natm_gauge = 0
1472 3690 : DO jatom = 1, natm_tot
1473 11376 : r(:) = pbc(particle_set(jatom)%r(:) - particle_set(iatom)%r(:), cell)
1474 : ! SQRT(SUM(r(:)**2)) <= 2.0_dp*hard_radius
1475 12222 : IF (SUM(r(:)**2) <= (4.0_dp*hard_radius*hard_radius)) THEN
1476 2160 : natm_gauge = natm_gauge + 1
1477 8640 : ratom(:, natm_gauge) = r(:)
1478 : END IF
1479 : END DO
1480 :
1481 3771 : DO ispin = 1, nspins
1482 3237237 : jrho1_atom%jrho_vec_rad_h(idir, ispin)%r_coef = 0.0_dp
1483 3237237 : jrho1_atom%jrho_vec_rad_s(idir, ispin)%r_coef = 0.0_dp
1484 : CALL calculate_jrho_atom_ang(jrho1_atom, &
1485 : jrho1_atom%jrho_vec_rad_h(idir, ispin)%r_coef, &
1486 : jrho1_atom%jrho_vec_rad_s(idir, ispin)%r_coef, &
1487 : grid_atom, harmonics, &
1488 : do_igaim, &
1489 2169 : ratom, natm_gauge, iB, idir, ispin)
1490 : END DO !ispin
1491 : END DO !iat
1492 : END DO !ikind
1493 : END IF
1494 :
1495 864 : DEALLOCATE (ratom)
1496 :
1497 864 : END SUBROUTINE calculate_jrho_atom
1498 :
1499 : END MODULE qs_linres_atom_current
|
