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