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 Support routines for integrals of the Vxc/Fxc/Gxc potentials calculated
10 : !> for the atomic density in the basis set of spherical primitives
11 : ! **************************************************************************************************
12 : MODULE qs_vxc_atom_utils
13 : USE basis_set_types, ONLY: get_gto_basis_set,&
14 : gto_basis_set_type
15 : USE external_potential_types, ONLY: get_potential,&
16 : gth_potential_type,&
17 : sgp_potential_type
18 : USE kinds, ONLY: dp
19 : USE mathconstants, ONLY: fac
20 : USE orbital_pointers, ONLY: indso,&
21 : nsoset
22 : USE paw_basis_types, ONLY: get_paw_basis_info
23 : USE qs_grid_atom, ONLY: grid_atom_type
24 : USE qs_harmonics_atom, ONLY: get_none0_cg_list,&
25 : harmonics_atom_type
26 : USE qs_kind_types, ONLY: has_nlcc
27 : USE qs_rho_atom_types, ONLY: rho_atom_coeff,&
28 : rho_atom_type
29 : #include "./base/base_uses.f90"
30 :
31 : IMPLICIT NONE
32 :
33 : PRIVATE
34 :
35 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_vxc_atom_utils'
36 :
37 : TYPE tau_basis_cache_type
38 : INTEGER :: maxso = 0, na = 0, nr = 0, nsatbas = 0, &
39 : nset = 0
40 : INTEGER, DIMENSION(:), POINTER :: lmax => NULL(), lmin => NULL(), &
41 : n2oindex => NULL(), npgf => NULL(), &
42 : o2nindex => NULL()
43 : REAL(dp), DIMENSION(:, :), POINTER :: zet => NULL()
44 : REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: grad
45 : END TYPE tau_basis_cache_type
46 :
47 : PUBLIC :: evaluate_nlcc_primitive_fields, gapw_atom_grid_support_radius, &
48 : create_tau_basis_cache, release_tau_basis_cache, &
49 : calc_rho_angular, calc_tau_atom, calc_rho_nlcc, &
50 : gaVxcgb_noGC, gaVxcgb_GC, dgaVtaudgb
51 : PUBLIC :: tau_basis_cache_type
52 : PUBLIC :: calc_weight_function
53 :
54 : ! **************************************************************************************************
55 :
56 : CONTAINS
57 :
58 : ! **************************************************************************************************
59 : !> \brief Evaluate an NLCC density and its first two Cartesian derivatives.
60 : !> \param point evaluation point
61 : !> \param center pseudopotential center
62 : !> \param gth_potential optional GTH potential
63 : !> \param sgp_potential optional separable Gaussian potential
64 : !> \param rho core density
65 : !> \param gradient Cartesian density gradient
66 : !> \param hessian Cartesian density Hessian
67 : ! **************************************************************************************************
68 112500 : SUBROUTINE evaluate_nlcc_primitive_fields(point, center, gth_potential, sgp_potential, &
69 : rho, gradient, hessian)
70 : REAL(dp), DIMENSION(3), INTENT(IN) :: point, center
71 : TYPE(gth_potential_type), INTENT(IN), POINTER :: gth_potential
72 : TYPE(sgp_potential_type), INTENT(IN), POINTER :: sgp_potential
73 : REAL(dp), INTENT(OUT) :: rho
74 : REAL(dp), DIMENSION(3), INTENT(OUT) :: gradient
75 : REAL(dp), DIMENSION(3, 3), INTENT(OUT) :: hessian
76 :
77 : INTEGER :: ic, idir, iexp, jdir, n_nlcc, nexp_nlcc, &
78 : power
79 112500 : INTEGER, DIMENSION(:), POINTER :: nct_nlcc
80 : LOGICAL :: has_sgp_nlcc, nlcc_present
81 : REAL(dp) :: alpha, beta, d2poly, dpoly, exponential, &
82 : poly, r2, rho_x, rho_xx, scaled_r2
83 : REAL(dp), DIMENSION(3) :: displacement
84 112500 : REAL(dp), DIMENSION(:), POINTER :: a_nlcc, alpha_nlcc, c_nlcc
85 112500 : REAL(dp), DIMENSION(:, :), POINTER :: cval_nlcc
86 :
87 112500 : NULLIFY (a_nlcc, alpha_nlcc, c_nlcc, cval_nlcc, nct_nlcc)
88 112500 : rho = 0.0_dp
89 112500 : rho_x = 0.0_dp
90 112500 : rho_xx = 0.0_dp
91 450000 : displacement = point - center
92 450000 : r2 = DOT_PRODUCT(displacement, displacement)
93 :
94 112500 : IF (ASSOCIATED(gth_potential)) THEN
95 : CALL get_potential(gth_potential, nlcc_present=nlcc_present, &
96 : nexp_nlcc=nexp_nlcc, alpha_nlcc=alpha_nlcc, &
97 112500 : nct_nlcc=nct_nlcc, cval_nlcc=cval_nlcc)
98 112500 : IF (nlcc_present) THEN
99 105000 : DO iexp = 1, nexp_nlcc
100 52500 : alpha = alpha_nlcc(iexp)
101 52500 : beta = 0.5_dp/(alpha*alpha)
102 52500 : scaled_r2 = r2/(alpha*alpha)
103 52500 : exponential = EXP(-0.5_dp*scaled_r2)
104 157500 : DO ic = 1, nct_nlcc(iexp)
105 52500 : power = ic - 1
106 52500 : poly = cval_nlcc(ic, iexp)*scaled_r2**power
107 52500 : dpoly = 0.0_dp
108 52500 : IF (power > 0) THEN
109 : dpoly = cval_nlcc(ic, iexp)*REAL(power, dp)* &
110 0 : scaled_r2**(power - 1)/(alpha*alpha)
111 : END IF
112 0 : d2poly = 0.0_dp
113 0 : IF (power > 1) THEN
114 : d2poly = cval_nlcc(ic, iexp)*REAL(power*(power - 1), dp)* &
115 0 : scaled_r2**(power - 2)/(alpha**4)
116 : END IF
117 52500 : rho = rho + exponential*poly
118 52500 : rho_x = rho_x + exponential*(dpoly - beta*poly)
119 105000 : rho_xx = rho_xx + exponential*(d2poly - 2.0_dp*beta*dpoly + beta*beta*poly)
120 : END DO
121 : END DO
122 : END IF
123 0 : ELSE IF (ASSOCIATED(sgp_potential)) THEN
124 : CALL get_potential(sgp_potential, has_nlcc=has_sgp_nlcc, n_nlcc=n_nlcc, &
125 0 : a_nlcc=a_nlcc, c_nlcc=c_nlcc)
126 0 : IF (has_sgp_nlcc) THEN
127 0 : DO iexp = 1, n_nlcc
128 0 : exponential = EXP(-a_nlcc(iexp)*r2)
129 0 : rho = rho + c_nlcc(iexp)*exponential
130 0 : rho_x = rho_x - a_nlcc(iexp)*c_nlcc(iexp)*exponential
131 0 : rho_xx = rho_xx + a_nlcc(iexp)**2*c_nlcc(iexp)*exponential
132 : END DO
133 : END IF
134 : END IF
135 :
136 450000 : gradient = 2.0_dp*rho_x*displacement
137 450000 : DO idir = 1, 3
138 1350000 : DO jdir = 1, 3
139 1350000 : hessian(idir, jdir) = 4.0_dp*rho_xx*displacement(idir)*displacement(jdir)
140 : END DO
141 450000 : hessian(idir, idir) = hessian(idir, idir) + 2.0_dp*rho_x
142 : END DO
143 112500 : END SUBROUTINE evaluate_nlcc_primitive_fields
144 :
145 : ! **************************************************************************************************
146 : !> \brief Compact support radius of one hard-minus-soft atom-grid field.
147 : !> \param grid_atom radial and angular source grid
148 : !> \param rho_h hard one-center density values
149 : !> \param rho_s soft one-center density values
150 : !> \param drho_h hard one-center density-gradient values
151 : !> \param drho_s soft one-center density-gradient values
152 : !> \param tau_h hard one-center kinetic-energy-density values
153 : !> \param tau_s soft one-center kinetic-energy-density values
154 : !> \return outermost radius required by any primitive hard-minus-soft field
155 : ! **************************************************************************************************
156 574 : FUNCTION gapw_atom_grid_support_radius( &
157 574 : grid_atom, rho_h, rho_s, drho_h, drho_s, tau_h, tau_s) RESULT(cutoff)
158 : TYPE(grid_atom_type), POINTER :: grid_atom
159 : REAL(dp), DIMENSION(:, :, :), INTENT(IN) :: rho_h, rho_s
160 : REAL(dp), DIMENSION(:, :, :, :), INTENT(IN) :: drho_h, drho_s
161 : REAL(dp), DIMENSION(:, :, :), INTENT(IN) :: tau_h, tau_s
162 : REAL(dp) :: cutoff
163 :
164 : INTEGER :: ir, support_ir
165 : LOGICAL :: descending
166 : REAL(dp) :: amplitude
167 :
168 574 : cutoff = 0.0_dp
169 574 : descending = grid_atom%rad(1) > grid_atom%rad(grid_atom%nr)
170 574 : IF (descending) THEN
171 8232 : DO ir = 1, grid_atom%nr
172 : amplitude = MAX( &
173 : MAXVAL(ABS(rho_h(:, ir, :) - rho_s(:, ir, :))), &
174 : MAXVAL(ABS(drho_h(1:3, :, ir, :) - drho_s(1:3, :, ir, :))), &
175 1562832 : MAXVAL(ABS(tau_h(:, ir, :) - tau_s(:, ir, :))))
176 8232 : IF (amplitude > 1.0E-12_dp) THEN
177 514 : support_ir = MAX(1, ir - 2)
178 514 : cutoff = grid_atom%rad(support_ir)
179 514 : RETURN
180 : END IF
181 : END DO
182 : ELSE
183 0 : DO ir = grid_atom%nr, 1, -1
184 : amplitude = MAX( &
185 : MAXVAL(ABS(rho_h(:, ir, :) - rho_s(:, ir, :))), &
186 : MAXVAL(ABS(drho_h(1:3, :, ir, :) - drho_s(1:3, :, ir, :))), &
187 0 : MAXVAL(ABS(tau_h(:, ir, :) - tau_s(:, ir, :))))
188 0 : IF (amplitude > 1.0E-12_dp) THEN
189 0 : support_ir = MIN(grid_atom%nr, ir + 2)
190 0 : cutoff = grid_atom%rad(support_ir)
191 0 : RETURN
192 : END IF
193 : END DO
194 : END IF
195 :
196 : END FUNCTION gapw_atom_grid_support_radius
197 :
198 : ! **************************************************************************************************
199 : !> \brief ...
200 : !> \param grid_atom ...
201 : !> \param harmonics ...
202 : !> \param nspins ...
203 : !> \param grad_func ...
204 : !> \param ir ...
205 : !> \param r_h ...
206 : !> \param r_s ...
207 : !> \param rho_h ...
208 : !> \param rho_s ...
209 : !> \param dr_h ...
210 : !> \param dr_s ...
211 : !> \param r_h_d ...
212 : !> \param r_s_d ...
213 : !> \param drho_h ...
214 : !> \param drho_s ...
215 : ! **************************************************************************************************
216 3079270 : SUBROUTINE calc_rho_angular(grid_atom, harmonics, nspins, grad_func, &
217 : ir, r_h, r_s, rho_h, rho_s, &
218 : dr_h, dr_s, r_h_d, r_s_d, drho_h, drho_s)
219 :
220 : TYPE(grid_atom_type), POINTER :: grid_atom
221 : TYPE(harmonics_atom_type), POINTER :: harmonics
222 : INTEGER, INTENT(IN) :: nspins
223 : LOGICAL, INTENT(IN) :: grad_func
224 : INTEGER, INTENT(IN) :: ir
225 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: r_h, r_s
226 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: rho_h, rho_s
227 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: dr_h, dr_s
228 : TYPE(rho_atom_coeff), DIMENSION(:, :), POINTER :: r_h_d, r_s_d
229 : REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: drho_h, drho_s
230 :
231 : INTEGER :: ia, iso, ispin, na
232 : REAL(KIND=dp) :: rad, urad
233 :
234 3079270 : CPASSERT(ASSOCIATED(r_h))
235 3079270 : CPASSERT(ASSOCIATED(r_s))
236 3079270 : CPASSERT(ASSOCIATED(rho_h))
237 3079270 : CPASSERT(ASSOCIATED(rho_s))
238 3079270 : IF (grad_func) THEN
239 1962520 : CPASSERT(ASSOCIATED(dr_h))
240 1962520 : CPASSERT(ASSOCIATED(dr_s))
241 1962520 : CPASSERT(ASSOCIATED(r_h_d))
242 1962520 : CPASSERT(ASSOCIATED(r_s_d))
243 1962520 : CPASSERT(ASSOCIATED(drho_h))
244 1962520 : CPASSERT(ASSOCIATED(drho_s))
245 : END IF
246 :
247 3079270 : na = grid_atom%ng_sphere
248 3079270 : rad = grid_atom%rad(ir)
249 3079270 : urad = grid_atom%oorad2l(ir, 1)
250 6536680 : DO ispin = 1, nspins
251 51864490 : DO iso = 1, harmonics%max_iso_not0
252 2315124600 : DO ia = 1, na
253 : rho_h(ia, ir, ispin) = rho_h(ia, ir, ispin) + &
254 2266339380 : r_h(ispin)%r_coef(ir, iso)*harmonics%slm(ia, iso)
255 : rho_s(ia, ir, ispin) = rho_s(ia, ir, ispin) + &
256 2311667190 : r_s(ispin)%r_coef(ir, iso)*harmonics%slm(ia, iso)
257 : END DO ! ia
258 : END DO ! iso
259 : END DO ! ispin
260 :
261 3079270 : IF (grad_func) THEN
262 4141630 : DO ispin = 1, nspins
263 31537820 : DO iso = 1, harmonics%max_iso_not0
264 1499422200 : DO ia = 1, na
265 :
266 : ! components of the gradient of rho1 hard
267 : drho_h(1, ia, ir, ispin) = drho_h(1, ia, ir, ispin) + &
268 : dr_h(ispin)%r_coef(ir, iso)* &
269 : harmonics%a(1, ia)*harmonics%slm(ia, iso) + &
270 : r_h_d(1, ispin)%r_coef(ir, iso)* &
271 1467884380 : harmonics%slm(ia, iso)
272 :
273 : drho_h(2, ia, ir, ispin) = drho_h(2, ia, ir, ispin) + &
274 : dr_h(ispin)%r_coef(ir, iso)* &
275 : harmonics%a(2, ia)*harmonics%slm(ia, iso) + &
276 : r_h_d(2, ispin)%r_coef(ir, iso)* &
277 1467884380 : harmonics%slm(ia, iso)
278 :
279 : drho_h(3, ia, ir, ispin) = drho_h(3, ia, ir, ispin) + &
280 : dr_h(ispin)%r_coef(ir, iso)* &
281 : harmonics%a(3, ia)*harmonics%slm(ia, iso) + &
282 : r_h_d(3, ispin)%r_coef(ir, iso)* &
283 1467884380 : harmonics%slm(ia, iso)
284 :
285 : ! components of the gradient of rho1 soft
286 : drho_s(1, ia, ir, ispin) = drho_s(1, ia, ir, ispin) + &
287 : dr_s(ispin)%r_coef(ir, iso)* &
288 : harmonics%a(1, ia)*harmonics%slm(ia, iso) + &
289 : r_s_d(1, ispin)%r_coef(ir, iso)* &
290 1467884380 : harmonics%slm(ia, iso)
291 :
292 : drho_s(2, ia, ir, ispin) = drho_s(2, ia, ir, ispin) + &
293 : dr_s(ispin)%r_coef(ir, iso)* &
294 : harmonics%a(2, ia)*harmonics%slm(ia, iso) + &
295 : r_s_d(2, ispin)%r_coef(ir, iso)* &
296 1467884380 : harmonics%slm(ia, iso)
297 :
298 : drho_s(3, ia, ir, ispin) = drho_s(3, ia, ir, ispin) + &
299 : dr_s(ispin)%r_coef(ir, iso)* &
300 : harmonics%a(3, ia)*harmonics%slm(ia, iso) + &
301 : r_s_d(3, ispin)%r_coef(ir, iso)* &
302 1497243090 : harmonics%slm(ia, iso)
303 :
304 : END DO ! ia
305 : END DO ! iso
306 112827130 : DO ia = 1, na
307 : drho_h(4, ia, ir, ispin) = SQRT( &
308 : drho_h(1, ia, ir, ispin)*drho_h(1, ia, ir, ispin) + &
309 : drho_h(2, ia, ir, ispin)*drho_h(2, ia, ir, ispin) + &
310 108685500 : drho_h(3, ia, ir, ispin)*drho_h(3, ia, ir, ispin))
311 :
312 : drho_s(4, ia, ir, ispin) = SQRT( &
313 : drho_s(1, ia, ir, ispin)*drho_s(1, ia, ir, ispin) + &
314 : drho_s(2, ia, ir, ispin)*drho_s(2, ia, ir, ispin) + &
315 110864610 : drho_s(3, ia, ir, ispin)*drho_s(3, ia, ir, ispin))
316 : END DO ! ia
317 : END DO ! ispin
318 : END IF
319 :
320 3079270 : END SUBROUTINE calc_rho_angular
321 :
322 : ! **************************************************************************************************
323 : !> \brief Precompute radial and angular factors for GAPW meta-GGA tau contractions
324 : !> \param tau_cache precomputed compact one-center gradient basis
325 : !> \param grid_atom atom-centered integration grid
326 : !> \param basis_1c GAPW one-center basis
327 : !> \param harmonics spherical harmonics on the atom-centered grid
328 : ! **************************************************************************************************
329 1700 : SUBROUTINE create_tau_basis_cache(tau_cache, grid_atom, basis_1c, harmonics)
330 :
331 : TYPE(tau_basis_cache_type), INTENT(INOUT) :: tau_cache
332 : TYPE(grid_atom_type), POINTER :: grid_atom
333 : TYPE(gto_basis_set_type), POINTER :: basis_1c
334 : TYPE(harmonics_atom_type), POINTER :: harmonics
335 :
336 : INTEGER :: dir, ia, igrid, ip, ipgf, ir, iset, iso, &
337 : l, starti
338 1700 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: a1, a2, gexp, r1, r2
339 1700 : REAL(dp), DIMENSION(:, :), POINTER :: slm
340 1700 : REAL(dp), DIMENSION(:, :, :), POINTER :: dslm_dxyz
341 :
342 1700 : NULLIFY (slm, dslm_dxyz)
343 :
344 1700 : CALL release_tau_basis_cache(tau_cache)
345 :
346 : CALL get_gto_basis_set(gto_basis_set=basis_1c, lmax=tau_cache%lmax, &
347 : lmin=tau_cache%lmin, maxso=tau_cache%maxso, &
348 : npgf=tau_cache%npgf, nset=tau_cache%nset, &
349 1700 : zet=tau_cache%zet)
350 : CALL get_paw_basis_info(basis_1c, o2nindex=tau_cache%o2nindex, &
351 : n2oindex=tau_cache%n2oindex, &
352 1700 : nsatbas=tau_cache%nsatbas)
353 :
354 1700 : tau_cache%nr = grid_atom%nr
355 1700 : tau_cache%na = grid_atom%ng_sphere
356 1700 : slm => harmonics%slm
357 1700 : dslm_dxyz => harmonics%dslm_dxyz
358 :
359 8500 : ALLOCATE (tau_cache%grad(tau_cache%na*tau_cache%nr, tau_cache%nsatbas, 3))
360 : ALLOCATE (a1(tau_cache%na), a2(tau_cache%na), gexp(tau_cache%nr), &
361 13600 : r1(tau_cache%nr), r2(tau_cache%nr))
362 283398830 : tau_cache%grad = 0.0_dp
363 :
364 6142 : DO iset = 1, tau_cache%nset
365 19874 : DO ipgf = 1, tau_cache%npgf(iset)
366 : starti = (iset - 1)*tau_cache%maxso + &
367 13732 : (ipgf - 1)*nsoset(tau_cache%lmax(iset))
368 : gexp(1:tau_cache%nr) = EXP(-tau_cache%zet(ipgf, iset)* &
369 663772 : grid_atom%rad2(1:tau_cache%nr))
370 56704 : DO iso = nsoset(tau_cache%lmin(iset) - 1) + 1, nsoset(tau_cache%lmax(iset))
371 38530 : ip = tau_cache%o2nindex(starti + iso)
372 38530 : IF (ip == 0) CYCLE
373 38530 : l = indso(1, iso)
374 :
375 1915750 : r1(1:tau_cache%nr) = grid_atom%rad(1:tau_cache%nr)**(l - 1)*gexp(1:tau_cache%nr)
376 : r2(1:tau_cache%nr) = -2.0_dp*tau_cache%zet(ipgf, iset)* &
377 1915750 : grid_atom%rad2(1:tau_cache%nr)*r1(1:tau_cache%nr)
378 :
379 167852 : DO dir = 1, 3
380 5665266 : a1(1:tau_cache%na) = dslm_dxyz(dir, 1:tau_cache%na, iso)
381 5665266 : a2(1:tau_cache%na) = harmonics%a(dir, 1:tau_cache%na)*slm(1:tau_cache%na, iso)
382 5785780 : DO ir = 1, tau_cache%nr
383 289023690 : DO ia = 1, tau_cache%na
384 283276440 : igrid = ia + (ir - 1)*tau_cache%na
385 288908100 : tau_cache%grad(igrid, ip, dir) = r1(ir)*a1(ia) + r2(ir)*a2(ia)
386 : END DO
387 : END DO
388 : END DO
389 : END DO
390 : END DO
391 : END DO
392 :
393 1700 : DEALLOCATE (a1, a2, gexp, r1, r2)
394 :
395 1700 : END SUBROUTINE create_tau_basis_cache
396 :
397 : ! **************************************************************************************************
398 : !> \brief Release precomputed GAPW meta-GGA tau factors
399 : !> \param tau_cache precomputed compact one-center gradient basis
400 : ! **************************************************************************************************
401 3400 : SUBROUTINE release_tau_basis_cache(tau_cache)
402 :
403 : TYPE(tau_basis_cache_type), INTENT(INOUT) :: tau_cache
404 :
405 3400 : IF (ALLOCATED(tau_cache%grad)) DEALLOCATE (tau_cache%grad)
406 3400 : IF (ASSOCIATED(tau_cache%n2oindex)) DEALLOCATE (tau_cache%n2oindex)
407 3400 : IF (ASSOCIATED(tau_cache%o2nindex)) DEALLOCATE (tau_cache%o2nindex)
408 3400 : NULLIFY (tau_cache%lmax, tau_cache%lmin, tau_cache%n2oindex, tau_cache%npgf, &
409 3400 : tau_cache%zet, tau_cache%o2nindex)
410 3400 : tau_cache%maxso = 0
411 3400 : tau_cache%na = 0
412 3400 : tau_cache%nr = 0
413 3400 : tau_cache%nsatbas = 0
414 3400 : tau_cache%nset = 0
415 :
416 3400 : END SUBROUTINE release_tau_basis_cache
417 :
418 : ! **************************************************************************************************
419 : !> \brief Computes tau hard and soft on the atomic grids for meta-GGA calculations
420 : !> \param tau_h the hard part of tau
421 : !> \param tau_s the soft part of tau
422 : !> \param rho_atom atom-centered density matrices
423 : !> \param tau_cache precomputed compact one-center gradient basis
424 : !> \param nspins number of spin channels
425 : !> \note This is a rewrite to correct a meta-GGA GAPW bug. This is more brute force than the original,
426 : !> which was done along in qs_rho_atom_methods.F, but makes sure that no corner is cut in
427 : !> terms of accuracy (A. Bussy)
428 : ! **************************************************************************************************
429 1512 : SUBROUTINE calc_tau_atom(tau_h, tau_s, rho_atom, tau_cache, nspins)
430 :
431 : REAL(dp), DIMENSION(:, :, :), INTENT(INOUT) :: tau_h, tau_s
432 : TYPE(rho_atom_type), POINTER :: rho_atom
433 : TYPE(tau_basis_cache_type), INTENT(IN) :: tau_cache
434 : INTEGER, INTENT(IN) :: nspins
435 :
436 : CHARACTER(len=*), PARAMETER :: routineN = 'calc_tau_atom'
437 :
438 : INTEGER :: dir, handle, ia, ibas, igrid, ir, ispin, &
439 : na, nbas, ngrid, nr
440 : REAL(dp) :: tau_sum
441 1512 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: work
442 :
443 1512 : CALL timeset(routineN, handle)
444 :
445 1512 : CPASSERT(ALLOCATED(tau_cache%grad))
446 :
447 : !zeroing tau, assuming it is already allocated
448 2949932 : tau_h = 0.0_dp
449 2949932 : tau_s = 0.0_dp
450 :
451 1512 : nr = tau_cache%nr
452 1512 : na = tau_cache%na
453 1512 : nbas = tau_cache%nsatbas
454 1512 : ngrid = na*nr
455 6048 : ALLOCATE (work(ngrid, nbas))
456 :
457 3032 : DO ispin = 1, nspins
458 7592 : DO dir = 1, 3
459 : CALL dgemm('N', 'T', ngrid, nbas, nbas, 0.5_dp, tau_cache%grad(:, :, dir), &
460 4560 : ngrid, rho_atom%cpc_h(ispin)%r_coef, nbas, 0.0_dp, work, ngrid)
461 : !$OMP PARALLEL DO DEFAULT(NONE) COLLAPSE(2) SCHEDULE(STATIC) &
462 : !$OMP SHARED(dir, ispin, na, nbas, nr, tau_cache, tau_h, work) &
463 4560 : !$OMP PRIVATE(ia, ibas, igrid, ir, tau_sum)
464 : DO ir = 1, nr
465 : DO ia = 1, na
466 : igrid = ia + (ir - 1)*na
467 : tau_sum = 0.0_dp
468 : DO ibas = 1, nbas
469 : tau_sum = tau_sum + tau_cache%grad(igrid, ibas, dir)*work(igrid, ibas)
470 : END DO
471 : tau_h(ia, ir, ispin) = tau_h(ia, ir, ispin) + tau_sum
472 : END DO
473 : END DO
474 : !$OMP END PARALLEL DO
475 :
476 : CALL dgemm('N', 'T', ngrid, nbas, nbas, 0.5_dp, tau_cache%grad(:, :, dir), &
477 4560 : ngrid, rho_atom%cpc_s(ispin)%r_coef, nbas, 0.0_dp, work, ngrid)
478 : !$OMP PARALLEL DO DEFAULT(NONE) COLLAPSE(2) SCHEDULE(STATIC) &
479 : !$OMP SHARED(dir, ispin, na, nbas, nr, tau_cache, tau_s, work) &
480 6080 : !$OMP PRIVATE(ia, ibas, igrid, ir, tau_sum)
481 : DO ir = 1, nr
482 : DO ia = 1, na
483 : igrid = ia + (ir - 1)*na
484 : tau_sum = 0.0_dp
485 : DO ibas = 1, nbas
486 : tau_sum = tau_sum + tau_cache%grad(igrid, ibas, dir)*work(igrid, ibas)
487 : END DO
488 : tau_s(ia, ir, ispin) = tau_s(ia, ir, ispin) + tau_sum
489 : END DO
490 : END DO
491 : !$OMP END PARALLEL DO
492 : END DO
493 : END DO
494 :
495 1512 : DEALLOCATE (work)
496 :
497 1512 : CALL timestop(handle)
498 :
499 3024 : END SUBROUTINE calc_tau_atom
500 :
501 : ! **************************************************************************************************
502 : !> \brief ...
503 : !> \param grid_atom ...
504 : !> \param nspins ...
505 : !> \param grad_func ...
506 : !> \param ir ...
507 : !> \param rho_nlcc ...
508 : !> \param rho_h ...
509 : !> \param rho_s ...
510 : !> \param drho_nlcc ...
511 : !> \param drho_h ...
512 : !> \param drho_s ...
513 : ! **************************************************************************************************
514 19750 : SUBROUTINE calc_rho_nlcc(grid_atom, nspins, grad_func, &
515 19750 : ir, rho_nlcc, rho_h, rho_s, drho_nlcc, drho_h, drho_s)
516 :
517 : TYPE(grid_atom_type), POINTER :: grid_atom
518 : INTEGER, INTENT(IN) :: nspins
519 : LOGICAL, INTENT(IN) :: grad_func
520 : INTEGER, INTENT(IN) :: ir
521 : REAL(KIND=dp), DIMENSION(:) :: rho_nlcc
522 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: rho_h, rho_s
523 : REAL(KIND=dp), DIMENSION(:) :: drho_nlcc
524 : REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: drho_h, drho_s
525 :
526 : INTEGER :: ia, ispin, na
527 : REAL(KIND=dp) :: drho, dx, dy, dz, rad, rho, urad, xsp
528 :
529 19750 : CPASSERT(ASSOCIATED(rho_h))
530 19750 : CPASSERT(ASSOCIATED(rho_s))
531 19750 : IF (grad_func) THEN
532 10450 : CPASSERT(ASSOCIATED(drho_h))
533 10450 : CPASSERT(ASSOCIATED(drho_s))
534 : END IF
535 :
536 19750 : na = grid_atom%ng_sphere
537 19750 : rad = grid_atom%rad(ir)
538 19750 : urad = grid_atom%oorad2l(ir, 1)
539 :
540 19750 : xsp = REAL(nspins, KIND=dp)
541 19750 : rho = rho_nlcc(ir)/xsp
542 39500 : DO ispin = 1, nspins
543 1007250 : rho_h(1:na, ir, ispin) = rho_h(1:na, ir, ispin) + rho
544 1027000 : rho_s(1:na, ir, ispin) = rho_s(1:na, ir, ispin) + rho
545 : END DO ! ispin
546 :
547 19750 : IF (grad_func) THEN
548 10450 : drho = drho_nlcc(ir)/xsp
549 20900 : DO ispin = 1, nspins
550 543400 : DO ia = 1, na
551 522500 : IF (grid_atom%azi(ia) == 0.0_dp) THEN
552 : dx = 0.0_dp
553 : dy = 0.0_dp
554 : ELSE
555 470250 : dx = grid_atom%sin_pol(ia)*grid_atom%sin_azi(ia)
556 470250 : dy = grid_atom%sin_pol(ia)*grid_atom%cos_azi(ia)
557 : END IF
558 522500 : dz = grid_atom%cos_pol(ia)
559 : ! components of the gradient of rho1 hard
560 522500 : drho_h(1, ia, ir, ispin) = drho_h(1, ia, ir, ispin) + drho*dx
561 522500 : drho_h(2, ia, ir, ispin) = drho_h(2, ia, ir, ispin) + drho*dy
562 522500 : drho_h(3, ia, ir, ispin) = drho_h(3, ia, ir, ispin) + drho*dz
563 : ! components of the gradient of rho1 soft
564 522500 : drho_s(1, ia, ir, ispin) = drho_s(1, ia, ir, ispin) + drho*dx
565 522500 : drho_s(2, ia, ir, ispin) = drho_s(2, ia, ir, ispin) + drho*dy
566 522500 : drho_s(3, ia, ir, ispin) = drho_s(3, ia, ir, ispin) + drho*dz
567 : ! norm of gradient
568 : drho_h(4, ia, ir, ispin) = SQRT( &
569 : drho_h(1, ia, ir, ispin)*drho_h(1, ia, ir, ispin) + &
570 : drho_h(2, ia, ir, ispin)*drho_h(2, ia, ir, ispin) + &
571 522500 : drho_h(3, ia, ir, ispin)*drho_h(3, ia, ir, ispin))
572 :
573 : drho_s(4, ia, ir, ispin) = SQRT( &
574 : drho_s(1, ia, ir, ispin)*drho_s(1, ia, ir, ispin) + &
575 : drho_s(2, ia, ir, ispin)*drho_s(2, ia, ir, ispin) + &
576 532950 : drho_s(3, ia, ir, ispin)*drho_s(3, ia, ir, ispin))
577 : END DO ! ia
578 : END DO ! ispin
579 : END IF
580 :
581 19750 : END SUBROUTINE calc_rho_nlcc
582 :
583 : ! **************************************************************************************************
584 : !> \brief ...
585 : !> \param vxc_h ...
586 : !> \param vxc_s ...
587 : !> \param int_hh ...
588 : !> \param int_ss ...
589 : !> \param grid_atom ...
590 : !> \param basis_1c ...
591 : !> \param harmonics ...
592 : !> \param nspins ...
593 : ! **************************************************************************************************
594 16817 : SUBROUTINE gaVxcgb_noGC(vxc_h, vxc_s, int_hh, int_ss, grid_atom, basis_1c, harmonics, nspins)
595 :
596 : REAL(dp), DIMENSION(:, :, :), POINTER :: vxc_h, vxc_s
597 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: int_hh, int_ss
598 : TYPE(grid_atom_type), POINTER :: grid_atom
599 : TYPE(gto_basis_set_type), POINTER :: basis_1c
600 : TYPE(harmonics_atom_type), POINTER :: harmonics
601 : INTEGER, INTENT(IN) :: nspins
602 :
603 : CHARACTER(len=*), PARAMETER :: routineN = 'gaVxcgb_noGC'
604 :
605 : INTEGER :: handle, ia, ic, icg, ipgf1, ipgf2, ir, iset1, iset2, iso, iso1, iso2, ispin, l, &
606 : ld, lmax12, lmax_expansion, lmin12, m1, m2, max_iso_not0, max_iso_not0_local, max_s_harm, &
607 : maxl, maxso, n1, n2, na, ngau1, ngau2, nngau1, nr, nset, size1
608 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list
609 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list
610 16817 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf
611 16817 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: g1, g2
612 16817 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: gg, gVg_h, gVg_s, matso_h, matso_s, vx
613 16817 : REAL(dp), DIMENSION(:, :), POINTER :: zet
614 : REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG
615 :
616 16817 : CALL timeset(routineN, handle)
617 :
618 16817 : NULLIFY (lmin, lmax, npgf, zet, my_CG)
619 :
620 : CALL get_gto_basis_set(gto_basis_set=basis_1c, lmax=lmax, lmin=lmin, &
621 : maxso=maxso, maxl=maxl, npgf=npgf, &
622 16817 : nset=nset, zet=zet)
623 :
624 16817 : nr = grid_atom%nr
625 16817 : na = grid_atom%ng_sphere
626 16817 : my_CG => harmonics%my_CG
627 16817 : max_iso_not0 = harmonics%max_iso_not0
628 16817 : lmax_expansion = indso(1, max_iso_not0)
629 16817 : max_s_harm = harmonics%max_s_harm
630 :
631 117719 : ALLOCATE (g1(nr), g2(nr), gg(nr, 0:2*maxl))
632 100902 : ALLOCATE (gVg_h(na, 0:2*maxl), gVg_s(na, 0:2*maxl))
633 : ALLOCATE (matso_h(nsoset(maxl), nsoset(maxl)), &
634 100902 : matso_s(nsoset(maxl), nsoset(maxl)))
635 67268 : ALLOCATE (vx(na, nr))
636 100902 : ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm))
637 :
638 16817 : g1 = 0.0_dp
639 16817 : g2 = 0.0_dp
640 16817 : m1 = 0
641 56701 : DO iset1 = 1, nset
642 39884 : n1 = nsoset(lmax(iset1))
643 39884 : m2 = 0
644 163658 : DO iset2 = 1, nset
645 : CALL get_none0_cg_list(my_CG, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
646 123774 : max_s_harm, lmax_expansion, cg_list, cg_n_list, max_iso_not0_local)
647 123774 : CPASSERT(max_iso_not0_local <= max_iso_not0)
648 :
649 123774 : n2 = nsoset(lmax(iset2))
650 450171 : DO ipgf1 = 1, npgf(iset1)
651 326397 : ngau1 = n1*(ipgf1 - 1) + m1
652 326397 : size1 = nsoset(lmax(iset1)) - nsoset(lmin(iset1) - 1)
653 326397 : nngau1 = nsoset(lmin(iset1) - 1) + ngau1
654 :
655 18755607 : g1(1:nr) = EXP(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))
656 1500527 : DO ipgf2 = 1, npgf(iset2)
657 1050356 : ngau2 = n2*(ipgf2 - 1) + m2
658 :
659 60010016 : g2(1:nr) = EXP(-zet(ipgf2, iset2)*grid_atom%rad2(1:nr))
660 1050356 : lmin12 = lmin(iset1) + lmin(iset2)
661 1050356 : lmax12 = lmax(iset1) + lmax(iset2)
662 :
663 : ! reduce expansion local densities
664 1376753 : IF (lmin12 <= lmax_expansion) THEN
665 :
666 1049411 : gg = 0.0_dp
667 1049411 : IF (lmin12 == 0) THEN
668 33829782 : gg(1:nr, lmin12) = g1(1:nr)*g2(1:nr)
669 : ELSE
670 26132039 : gg(1:nr, lmin12) = grid_atom%rad2l(1:nr, lmin12)*g1(1:nr)*g2(1:nr)
671 : END IF
672 :
673 : ! limit the expansion of the local densities to a max L
674 1049411 : IF (lmax12 > lmax_expansion) lmax12 = lmax_expansion
675 :
676 1522423 : DO l = lmin12 + 1, lmax12
677 29661823 : gg(1:nr, l) = grid_atom%rad(1:nr)*gg(:, l - 1)
678 : END DO
679 :
680 2275012 : DO ispin = 1, nspins
681 1225601 : ld = lmax12 + 1
682 73959811 : DO ir = 1, nr
683 3709498151 : vx(1:na, ir) = vxc_h(1:na, ir, ispin)
684 : END DO
685 : CALL dgemm('N', 'N', na, ld, nr, 1.0_dp, vx(1:na, 1:nr), na, &
686 1225601 : gg(1:nr, 0:lmax12), nr, 0.0_dp, gVg_h(1:na, 0:lmax12), na)
687 73959811 : DO ir = 1, nr
688 3709498151 : vx(1:na, ir) = vxc_s(1:na, ir, ispin)
689 : END DO
690 : CALL dgemm('N', 'N', na, ld, nr, 1.0_dp, vx(1:na, 1:nr), na, &
691 1225601 : gg(1:nr, 0:lmax12), nr, 0.0_dp, gVg_s(1:na, 0:lmax12), na)
692 :
693 1225601 : matso_h = 0.0_dp
694 1225601 : matso_s = 0.0_dp
695 9503890 : DO iso = 1, max_iso_not0_local
696 25997818 : DO icg = 1, cg_n_list(iso)
697 16493928 : iso1 = cg_list(1, icg, iso)
698 16493928 : iso2 = cg_list(2, icg, iso)
699 16493928 : l = indso(1, iso1) + indso(1, iso2)
700 :
701 16493928 : CPASSERT(l <= lmax_expansion)
702 849160529 : DO ia = 1, na
703 : matso_h(iso1, iso2) = matso_h(iso1, iso2) + &
704 : gVg_h(ia, l)* &
705 : my_CG(iso1, iso2, iso)* &
706 824388312 : harmonics%slm(ia, iso)
707 : matso_s(iso1, iso2) = matso_s(iso1, iso2) + &
708 : gVg_s(ia, l)* &
709 : my_CG(iso1, iso2, iso)* &
710 840882240 : harmonics%slm(ia, iso)
711 : END DO
712 : END DO
713 : END DO
714 :
715 : ! Write in the global matrix
716 5449457 : DO ic = nsoset(lmin(iset2) - 1) + 1, nsoset(lmax(iset2))
717 3174445 : iso1 = nsoset(lmin(iset1) - 1) + 1
718 3174445 : iso2 = ngau2 + ic
719 : CALL daxpy(size1, 1.0_dp, matso_h(iso1, ic), 1, &
720 3174445 : int_hh(ispin)%r_coef(nngau1 + 1, iso2), 1)
721 : CALL daxpy(size1, 1.0_dp, matso_s(iso1, ic), 1, &
722 4400046 : int_ss(ispin)%r_coef(nngau1 + 1, iso2), 1)
723 : END DO
724 :
725 : END DO ! ispin
726 :
727 : END IF ! lmax_expansion
728 :
729 : END DO ! ipfg2
730 : END DO ! ipfg1
731 287432 : m2 = m2 + maxso
732 : END DO ! iset2
733 56701 : m1 = m1 + maxso
734 : END DO ! iset1
735 :
736 16817 : DEALLOCATE (g1, g2, gg, matso_h, matso_s, gVg_s, gVg_h, vx)
737 :
738 16817 : DEALLOCATE (cg_list, cg_n_list)
739 :
740 16817 : CALL timestop(handle)
741 :
742 16817 : END SUBROUTINE gaVxcgb_noGC
743 :
744 : ! **************************************************************************************************
745 : !> \brief ...
746 : !> \param vxc_h ...
747 : !> \param vxc_s ...
748 : !> \param vxg_h ...
749 : !> \param vxg_s ...
750 : !> \param int_hh ...
751 : !> \param int_ss ...
752 : !> \param grid_atom ...
753 : !> \param basis_1c ...
754 : !> \param harmonics ...
755 : !> \param nspins ...
756 : ! **************************************************************************************************
757 31341 : SUBROUTINE gaVxcgb_GC(vxc_h, vxc_s, vxg_h, vxg_s, int_hh, int_ss, &
758 : grid_atom, basis_1c, harmonics, nspins)
759 :
760 : REAL(dp), DIMENSION(:, :, :), POINTER :: vxc_h, vxc_s
761 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: vxg_h, vxg_s
762 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: int_hh, int_ss
763 : TYPE(grid_atom_type), POINTER :: grid_atom
764 : TYPE(gto_basis_set_type), POINTER :: basis_1c
765 : TYPE(harmonics_atom_type), POINTER :: harmonics
766 : INTEGER, INTENT(IN) :: nspins
767 :
768 : CHARACTER(len=*), PARAMETER :: routineN = 'gaVxcgb_GC'
769 :
770 : INTEGER :: dmax_iso_not0_local, handle, ia, ic, icg, ipgf1, ipgf2, ir, iset1, iset2, iso, &
771 : iso1, iso2, ispin, l, lmax12, lmax_expansion, lmin12, m1, m2, max_iso_not0, &
772 : max_iso_not0_local, max_s_harm, maxl, maxso, n1, n2, na, ngau1, ngau2, nngau1, nr, nset, &
773 : size1
774 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list, dcg_n_list
775 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list, dcg_list
776 31341 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf
777 : REAL(dp) :: urad
778 31341 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: g1, g2
779 31341 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: dgg, gg, gVXCg_h, gVXCg_s, matso_h, &
780 31341 : matso_s
781 31341 : REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: gVXGg_h, gVXGg_s
782 31341 : REAL(dp), DIMENSION(:, :), POINTER :: zet
783 : REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG
784 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: my_CG_dxyz
785 :
786 31341 : CALL timeset(routineN, handle)
787 :
788 31341 : NULLIFY (lmin, lmax, npgf, zet, my_CG, my_CG_dxyz)
789 :
790 : CALL get_gto_basis_set(gto_basis_set=basis_1c, lmax=lmax, lmin=lmin, &
791 : maxso=maxso, maxl=maxl, npgf=npgf, &
792 31341 : nset=nset, zet=zet)
793 :
794 31341 : nr = grid_atom%nr
795 31341 : na = grid_atom%ng_sphere
796 31341 : my_CG => harmonics%my_CG
797 31341 : my_CG_dxyz => harmonics%my_CG_dxyz
798 31341 : max_iso_not0 = harmonics%max_iso_not0
799 31341 : lmax_expansion = indso(1, max_iso_not0)
800 31341 : max_s_harm = harmonics%max_s_harm
801 :
802 282069 : ALLOCATE (g1(nr), g2(nr), gg(nr, 0:2*maxl), dgg(nr, 0:2*maxl))
803 188046 : ALLOCATE (gVXCg_h(na, 0:2*maxl), gVXCg_s(na, 0:2*maxl))
804 188046 : ALLOCATE (gVXGg_h(3, na, 0:2*maxl), gVXGg_s(3, na, 0:2*maxl))
805 : ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm), &
806 282069 : dcg_list(2, nsoset(maxl)**2, max_s_harm), dcg_n_list(max_s_harm))
807 :
808 : ALLOCATE (matso_h(nsoset(maxl), nsoset(maxl)), &
809 188046 : matso_s(nsoset(maxl), nsoset(maxl)))
810 :
811 66113 : DO ispin = 1, nspins
812 :
813 34772 : g1 = 0.0_dp
814 34772 : g2 = 0.0_dp
815 34772 : m1 = 0
816 153754 : DO iset1 = 1, nset
817 87641 : n1 = nsoset(lmax(iset1))
818 87641 : m2 = 0
819 386316 : DO iset2 = 1, nset
820 : CALL get_none0_cg_list(my_CG, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
821 298675 : max_s_harm, lmax_expansion, cg_list, cg_n_list, max_iso_not0_local)
822 298675 : CPASSERT(max_iso_not0_local <= max_iso_not0)
823 : CALL get_none0_cg_list(my_CG_dxyz, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
824 298675 : max_s_harm, lmax_expansion, dcg_list, dcg_n_list, dmax_iso_not0_local)
825 :
826 298675 : n2 = nsoset(lmax(iset2))
827 963645 : DO ipgf1 = 1, npgf(iset1)
828 664970 : ngau1 = n1*(ipgf1 - 1) + m1
829 664970 : size1 = nsoset(lmax(iset1)) - nsoset(lmin(iset1) - 1)
830 664970 : nngau1 = nsoset(lmin(iset1) - 1) + ngau1
831 :
832 34605270 : g1(1:nr) = EXP(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))
833 2638136 : DO ipgf2 = 1, npgf(iset2)
834 1674491 : ngau2 = n2*(ipgf2 - 1) + m2
835 :
836 87328041 : g2(1:nr) = EXP(-zet(ipgf2, iset2)*grid_atom%rad2(1:nr))
837 1674491 : lmin12 = lmin(iset1) + lmin(iset2)
838 1674491 : lmax12 = lmax(iset1) + lmax(iset2)
839 :
840 : !test reduce expansion local densities
841 1674491 : IF (lmin12 <= lmax_expansion) THEN
842 :
843 1673891 : gg = 0.0_dp
844 1673891 : dgg = 0.0_dp
845 :
846 1673891 : IF (lmin12 == 0) THEN
847 52703349 : gg(1:nr, lmin12) = g1(1:nr)*g2(1:nr)
848 : ELSE
849 34594092 : gg(1:nr, lmin12) = grid_atom%rad2l(1:nr, lmin12)*g1(1:nr)*g2(1:nr)
850 : END IF
851 :
852 : !test reduce expansion local densities
853 1673891 : IF (lmax12 > lmax_expansion) lmax12 = lmax_expansion
854 :
855 2602413 : DO l = lmin12 + 1, lmax12
856 48789422 : gg(1:nr, l) = grid_atom%rad(1:nr)*gg(:, l - 1)
857 : dgg(1:nr, l - 1) = dgg(1:nr, l - 1) - 2.0_dp*(zet(ipgf1, iset1) + &
858 50463313 : zet(ipgf2, iset2))*gg(1:nr, l)
859 : END DO
860 : dgg(1:nr, lmax12) = dgg(1:nr, lmax12) - 2.0_dp*(zet(ipgf1, iset1) + &
861 : zet(ipgf2, iset2))*grid_atom%rad(1:nr)* &
862 87297441 : gg(1:nr, lmax12)
863 :
864 1673891 : gVXCg_h = 0.0_dp
865 1673891 : gVXCg_s = 0.0_dp
866 1673891 : gVXGg_h = 0.0_dp
867 1673891 : gVXGg_s = 0.0_dp
868 :
869 : ! Cross Term
870 4276304 : DO l = lmin12, lmax12
871 134202170 : DO ia = 1, na
872 6807372859 : DO ir = 1, nr
873 : gVXCg_h(ia, l) = gVXCg_h(ia, l) + &
874 : gg(ir, l)*vxc_h(ia, ir, ispin) + &
875 : dgg(ir, l)* &
876 : (vxg_h(1, ia, ir, ispin)*harmonics%a(1, ia) + &
877 : vxg_h(2, ia, ir, ispin)*harmonics%a(2, ia) + &
878 6674844580 : vxg_h(3, ia, ir, ispin)*harmonics%a(3, ia))
879 :
880 : gVXCg_s(ia, l) = gVXCg_s(ia, l) + &
881 : gg(ir, l)*vxc_s(ia, ir, ispin) + &
882 : dgg(ir, l)* &
883 : (vxg_s(1, ia, ir, ispin)*harmonics%a(1, ia) + &
884 : vxg_s(2, ia, ir, ispin)*harmonics%a(2, ia) + &
885 6674844580 : vxg_s(3, ia, ir, ispin)*harmonics%a(3, ia))
886 :
887 6674844580 : urad = grid_atom%oorad2l(ir, 1)
888 :
889 : gVXGg_h(1, ia, l) = gVXGg_h(1, ia, l) + &
890 : vxg_h(1, ia, ir, ispin)* &
891 6674844580 : gg(ir, l)*urad
892 :
893 : gVXGg_h(2, ia, l) = gVXGg_h(2, ia, l) + &
894 : vxg_h(2, ia, ir, ispin)* &
895 6674844580 : gg(ir, l)*urad
896 :
897 : gVXGg_h(3, ia, l) = gVXGg_h(3, ia, l) + &
898 : vxg_h(3, ia, ir, ispin)* &
899 6674844580 : gg(ir, l)*urad
900 :
901 : gVXGg_s(1, ia, l) = gVXGg_s(1, ia, l) + &
902 : vxg_s(1, ia, ir, ispin)* &
903 6674844580 : gg(ir, l)*urad
904 :
905 : gVXGg_s(2, ia, l) = gVXGg_s(2, ia, l) + &
906 : vxg_s(2, ia, ir, ispin)* &
907 6674844580 : gg(ir, l)*urad
908 :
909 : gVXGg_s(3, ia, l) = gVXGg_s(3, ia, l) + &
910 : vxg_s(3, ia, ir, ispin)* &
911 6804770446 : gg(ir, l)*urad
912 :
913 : END DO ! ir
914 : END DO ! ia
915 : END DO ! l
916 :
917 1673891 : matso_h = 0.0_dp
918 1673891 : matso_s = 0.0_dp
919 12190332 : DO iso = 1, max_iso_not0_local
920 35011831 : DO icg = 1, cg_n_list(iso)
921 22821499 : iso1 = cg_list(1, icg, iso)
922 22821499 : iso2 = cg_list(2, icg, iso)
923 :
924 22821499 : l = indso(1, iso1) + indso(1, iso2)
925 :
926 : !test reduce expansion local densities
927 22821499 : CPASSERT(l <= lmax_expansion)
928 1173548374 : DO ia = 1, na
929 : matso_h(iso1, iso2) = matso_h(iso1, iso2) + &
930 : gVXCg_h(ia, l)* &
931 : harmonics%slm(ia, iso)* &
932 1140210434 : my_CG(iso1, iso2, iso)
933 : matso_s(iso1, iso2) = matso_s(iso1, iso2) + &
934 : gVXCg_s(ia, l)* &
935 : harmonics%slm(ia, iso)* &
936 1163031933 : my_CG(iso1, iso2, iso)
937 : END DO ! ia
938 :
939 : !test reduce expansion local densities
940 :
941 : END DO
942 :
943 : END DO ! iso
944 :
945 6510409 : DO iso = 1, dmax_iso_not0_local
946 43747572 : DO icg = 1, dcg_n_list(iso)
947 37237163 : iso1 = dcg_list(1, icg, iso)
948 37237163 : iso2 = dcg_list(2, icg, iso)
949 :
950 37237163 : l = indso(1, iso1) + indso(1, iso2)
951 : !test reduce expansion local densities
952 37237163 : CPASSERT(l <= lmax_expansion)
953 1902023363 : DO ia = 1, na
954 : matso_h(iso1, iso2) = matso_h(iso1, iso2) + &
955 : (gVXGg_h(1, ia, l)*my_CG_dxyz(1, iso1, iso2, iso) + &
956 : gVXGg_h(2, ia, l)*my_CG_dxyz(2, iso1, iso2, iso) + &
957 : gVXGg_h(3, ia, l)*my_CG_dxyz(3, iso1, iso2, iso))* &
958 1859949682 : harmonics%slm(ia, iso)
959 :
960 : matso_s(iso1, iso2) = matso_s(iso1, iso2) + &
961 : (gVXGg_s(1, ia, l)*my_CG_dxyz(1, iso1, iso2, iso) + &
962 : gVXGg_s(2, ia, l)*my_CG_dxyz(2, iso1, iso2, iso) + &
963 : gVXGg_s(3, ia, l)*my_CG_dxyz(3, iso1, iso2, iso))* &
964 1897186845 : harmonics%slm(ia, iso)
965 :
966 : END DO ! ia
967 :
968 : !test reduce expansion local densities
969 :
970 : END DO ! icg
971 : END DO ! iso
972 : !test reduce expansion local densities
973 : END IF ! lmax_expansion
974 :
975 : ! Write in the global matrix
976 6650865 : DO ic = nsoset(lmin(iset2) - 1) + 1, nsoset(lmax(iset2))
977 4311404 : iso1 = nsoset(lmin(iset1) - 1) + 1
978 4311404 : iso2 = ngau2 + ic
979 : CALL daxpy(size1, 1.0_dp, matso_h(iso1, ic), 1, &
980 4311404 : int_hh(ispin)%r_coef(nngau1 + 1, iso2), 1)
981 : CALL daxpy(size1, 1.0_dp, matso_s(iso1, ic), 1, &
982 5985895 : int_ss(ispin)%r_coef(nngau1 + 1, iso2), 1)
983 : END DO
984 :
985 : END DO ! ipfg2
986 : END DO ! ipfg1
987 983666 : m2 = m2 + maxso
988 : END DO ! iset2
989 122413 : m1 = m1 + maxso
990 : END DO ! iset1
991 : END DO ! ispin
992 :
993 31341 : DEALLOCATE (g1, g2, gg, dgg, matso_h, matso_s, gVXCg_h, gVXCg_s, gVXGg_h, gVXGg_s)
994 31341 : DEALLOCATE (cg_list, cg_n_list, dcg_list, dcg_n_list)
995 :
996 31341 : CALL timestop(handle)
997 :
998 31341 : END SUBROUTINE gaVxcgb_GC
999 :
1000 : ! **************************************************************************************************
1001 : !> \brief Integrates 0.5 * grad_ga .dot. (V_tau * grad_gb) on the atomic grid for meta-GGA
1002 : !> \param vtau_h the hard tau potential
1003 : !> \param vtau_s the soft tau potential
1004 : !> \param int_hh hard one-center matrix contribution
1005 : !> \param int_ss soft one-center matrix contribution
1006 : !> \param tau_cache precomputed compact one-center gradient basis
1007 : !> \param nspins number of spin channels
1008 : !> \note This is a rewrite to correct meta-GGA GAPW bug. This is more brute force than the original
1009 : !> but makes sure that no corner is cut in terms of accuracy (A. Bussy)
1010 : ! **************************************************************************************************
1011 1064 : SUBROUTINE dgaVtaudgb(vtau_h, vtau_s, int_hh, int_ss, tau_cache, nspins)
1012 :
1013 : REAL(dp), DIMENSION(:, :, :), POINTER :: vtau_h, vtau_s
1014 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: int_hh, int_ss
1015 : TYPE(tau_basis_cache_type), INTENT(IN) :: tau_cache
1016 : INTEGER, INTENT(IN) :: nspins
1017 :
1018 : CHARACTER(len=*), PARAMETER :: routineN = 'dgaVtaudgb'
1019 :
1020 : INTEGER :: dir, handle, ia, ibas, igrid, iold, ir, &
1021 : ispin, jbas, jold, max_old_basis, na, &
1022 : nbas, ngrid, nr
1023 1064 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: int_h, int_s, weighted_grad
1024 :
1025 1064 : CALL timeset(routineN, handle)
1026 :
1027 1064 : CPASSERT(ALLOCATED(tau_cache%grad))
1028 1064 : CPASSERT(ASSOCIATED(tau_cache%n2oindex))
1029 :
1030 1064 : nr = tau_cache%nr
1031 1064 : na = tau_cache%na
1032 1064 : nbas = tau_cache%nsatbas
1033 1064 : ngrid = na*nr
1034 64808 : max_old_basis = MAXVAL(tau_cache%n2oindex)
1035 9576 : ALLOCATE (int_h(nbas, nbas), int_s(nbas, nbas), weighted_grad(ngrid, nbas))
1036 :
1037 2136 : DO ispin = 1, nspins
1038 1072 : CPASSERT(SIZE(int_hh(ispin)%r_coef, 1) >= max_old_basis)
1039 1072 : CPASSERT(SIZE(int_hh(ispin)%r_coef, 2) >= max_old_basis)
1040 1072 : CPASSERT(SIZE(int_ss(ispin)%r_coef, 1) >= max_old_basis)
1041 1072 : CPASSERT(SIZE(int_ss(ispin)%r_coef, 2) >= max_old_basis)
1042 1072 : int_h = 0.0_dp
1043 1072 : int_s = 0.0_dp
1044 4288 : DO dir = 1, 3
1045 : !$OMP PARALLEL DO DEFAULT(NONE) COLLAPSE(3) SCHEDULE(STATIC) &
1046 : !$OMP SHARED(dir, ispin, na, nbas, nr, tau_cache, vtau_h, weighted_grad) &
1047 3216 : !$OMP PRIVATE(ia, ibas, igrid, ir)
1048 : DO ibas = 1, nbas
1049 : DO ir = 1, nr
1050 : DO ia = 1, na
1051 : igrid = ia + (ir - 1)*na
1052 : weighted_grad(igrid, ibas) = vtau_h(ia, ir, ispin)* &
1053 : tau_cache%grad(igrid, ibas, dir)
1054 : END DO
1055 : END DO
1056 : END DO
1057 : !$OMP END PARALLEL DO
1058 : CALL dgemm('T', 'N', nbas, nbas, ngrid, 0.5_dp, tau_cache%grad(:, :, dir), &
1059 3216 : ngrid, weighted_grad, ngrid, 1.0_dp, int_h, nbas)
1060 :
1061 : !$OMP PARALLEL DO DEFAULT(NONE) COLLAPSE(3) SCHEDULE(STATIC) &
1062 : !$OMP SHARED(dir, ispin, na, nbas, nr, tau_cache, vtau_s, weighted_grad) &
1063 3216 : !$OMP PRIVATE(ia, ibas, igrid, ir)
1064 : DO ibas = 1, nbas
1065 : DO ir = 1, nr
1066 : DO ia = 1, na
1067 : igrid = ia + (ir - 1)*na
1068 : weighted_grad(igrid, ibas) = vtau_s(ia, ir, ispin)* &
1069 : tau_cache%grad(igrid, ibas, dir)
1070 : END DO
1071 : END DO
1072 : END DO
1073 : !$OMP END PARALLEL DO
1074 : CALL dgemm('T', 'N', nbas, nbas, ngrid, 0.5_dp, tau_cache%grad(:, :, dir), &
1075 4288 : ngrid, weighted_grad, ngrid, 1.0_dp, int_s, nbas)
1076 : END DO
1077 :
1078 : !$OMP PARALLEL DO DEFAULT(NONE) COLLAPSE(2) SCHEDULE(STATIC) &
1079 : !$OMP SHARED(int_h, int_hh, int_s, int_ss, ispin, nbas, tau_cache) &
1080 2136 : !$OMP PRIVATE(ibas, iold, jbas, jold)
1081 : DO jbas = 1, nbas
1082 : DO ibas = 1, nbas
1083 : jold = tau_cache%n2oindex(jbas)
1084 : iold = tau_cache%n2oindex(ibas)
1085 : int_hh(ispin)%r_coef(iold, jold) = int_hh(ispin)%r_coef(iold, jold) + &
1086 : int_h(ibas, jbas)
1087 : int_ss(ispin)%r_coef(iold, jold) = int_ss(ispin)%r_coef(iold, jold) + &
1088 : int_s(ibas, jbas)
1089 : END DO
1090 : END DO
1091 : !$OMP END PARALLEL DO
1092 : END DO
1093 :
1094 1064 : DEALLOCATE (int_h, int_s, weighted_grad)
1095 :
1096 1064 : CALL timestop(handle)
1097 :
1098 2128 : END SUBROUTINE dgaVtaudgb
1099 :
1100 : ! **************************************************************************************************
1101 : !> \brief Calculates the radial weight function for accurate XC integration
1102 : !> \param fun The weight function
1103 : !> \param r2 Radial values of the grid (squared)
1104 : !> \param order Polynomila order (2n)
1105 : !> \param alpha Gaussian exponent
1106 : ! **************************************************************************************************
1107 536 : SUBROUTINE calc_weight_function(fun, r2, order, alpha)
1108 :
1109 : REAL(dp), DIMENSION(:), INTENT(INOUT) :: fun
1110 : REAL(dp), DIMENSION(:), INTENT(IN) :: r2
1111 : INTEGER, INTENT(IN) :: order
1112 : REAL(dp), INTENT(IN) :: alpha
1113 :
1114 : INTEGER :: i
1115 : REAL(dp) :: pval
1116 :
1117 536 : CPASSERT(SIZE(fun) == SIZE(r2))
1118 :
1119 27156 : fun(:) = 1.0_dp
1120 540 : DO i = 1, order
1121 4 : pval = alpha**i/fac(i)
1122 740 : fun(:) = fun(:) + pval*r2**i
1123 : END DO
1124 27156 : fun(:) = fun(:)*EXP(-alpha*r2(:))
1125 :
1126 536 : END SUBROUTINE calc_weight_function
1127 :
1128 : ! **************************************************************************************************
1129 :
1130 0 : END MODULE qs_vxc_atom_utils
|
