(git:e90df42)
Loading...
Searching...
No Matches
qs_vxc_atom_utils.F
Go to the documentation of this file.
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! **************************************************************************************************
12MODULE 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
56CONTAINS
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 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 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 REAL(dp), DIMENSION(:), POINTER :: a_nlcc, alpha_nlcc, c_nlcc
85 REAL(dp), DIMENSION(:, :), POINTER :: cval_nlcc
86
87 NULLIFY (a_nlcc, alpha_nlcc, c_nlcc, cval_nlcc, nct_nlcc)
88 rho = 0.0_dp
89 rho_x = 0.0_dp
90 rho_xx = 0.0_dp
91 displacement = point - center
92 r2 = dot_product(displacement, displacement)
93
94 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 nct_nlcc=nct_nlcc, cval_nlcc=cval_nlcc)
98 IF (nlcc_present) THEN
99 DO iexp = 1, nexp_nlcc
100 alpha = alpha_nlcc(iexp)
101 beta = 0.5_dp/(alpha*alpha)
102 scaled_r2 = r2/(alpha*alpha)
103 exponential = exp(-0.5_dp*scaled_r2)
104 DO ic = 1, nct_nlcc(iexp)
105 power = ic - 1
106 poly = cval_nlcc(ic, iexp)*scaled_r2**power
107 dpoly = 0.0_dp
108 IF (power > 0) THEN
109 dpoly = cval_nlcc(ic, iexp)*real(power, dp)* &
110 scaled_r2**(power - 1)/(alpha*alpha)
111 END IF
112 d2poly = 0.0_dp
113 IF (power > 1) THEN
114 d2poly = cval_nlcc(ic, iexp)*real(power*(power - 1), dp)* &
115 scaled_r2**(power - 2)/(alpha**4)
116 END IF
117 rho = rho + exponential*poly
118 rho_x = rho_x + exponential*(dpoly - beta*poly)
119 rho_xx = rho_xx + exponential*(d2poly - 2.0_dp*beta*dpoly + beta*beta*poly)
120 END DO
121 END DO
122 END IF
123 ELSE IF (ASSOCIATED(sgp_potential)) THEN
124 CALL get_potential(sgp_potential, has_nlcc=has_sgp_nlcc, n_nlcc=n_nlcc, &
125 a_nlcc=a_nlcc, c_nlcc=c_nlcc)
126 IF (has_sgp_nlcc) THEN
127 DO iexp = 1, n_nlcc
128 exponential = exp(-a_nlcc(iexp)*r2)
129 rho = rho + c_nlcc(iexp)*exponential
130 rho_x = rho_x - a_nlcc(iexp)*c_nlcc(iexp)*exponential
131 rho_xx = rho_xx + a_nlcc(iexp)**2*c_nlcc(iexp)*exponential
132 END DO
133 END IF
134 END IF
135
136 gradient = 2.0_dp*rho_x*displacement
137 DO idir = 1, 3
138 DO jdir = 1, 3
139 hessian(idir, jdir) = 4.0_dp*rho_xx*displacement(idir)*displacement(jdir)
140 END DO
141 hessian(idir, idir) = hessian(idir, idir) + 2.0_dp*rho_x
142 END DO
143 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 FUNCTION gapw_atom_grid_support_radius( &
157 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 cutoff = 0.0_dp
169 descending = grid_atom%rad(1) > grid_atom%rad(grid_atom%nr)
170 IF (descending) THEN
171 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 maxval(abs(tau_h(:, ir, :) - tau_s(:, ir, :))))
176 IF (amplitude > 1.0e-12_dp) THEN
177 support_ir = max(1, ir - 2)
178 cutoff = grid_atom%rad(support_ir)
179 RETURN
180 END IF
181 END DO
182 ELSE
183 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 maxval(abs(tau_h(:, ir, :) - tau_s(:, ir, :))))
188 IF (amplitude > 1.0e-12_dp) THEN
189 support_ir = min(grid_atom%nr, ir + 2)
190 cutoff = grid_atom%rad(support_ir)
191 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 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 cpassert(ASSOCIATED(r_h))
235 cpassert(ASSOCIATED(r_s))
236 cpassert(ASSOCIATED(rho_h))
237 cpassert(ASSOCIATED(rho_s))
238 IF (grad_func) THEN
239 cpassert(ASSOCIATED(dr_h))
240 cpassert(ASSOCIATED(dr_s))
241 cpassert(ASSOCIATED(r_h_d))
242 cpassert(ASSOCIATED(r_s_d))
243 cpassert(ASSOCIATED(drho_h))
244 cpassert(ASSOCIATED(drho_s))
245 END IF
246
247 na = grid_atom%ng_sphere
248 rad = grid_atom%rad(ir)
249 urad = grid_atom%oorad2l(ir, 1)
250 DO ispin = 1, nspins
251 DO iso = 1, harmonics%max_iso_not0
252 DO ia = 1, na
253 rho_h(ia, ir, ispin) = rho_h(ia, ir, ispin) + &
254 r_h(ispin)%r_coef(ir, iso)*harmonics%slm(ia, iso)
255 rho_s(ia, ir, ispin) = rho_s(ia, ir, ispin) + &
256 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 IF (grad_func) THEN
262 DO ispin = 1, nspins
263 DO iso = 1, harmonics%max_iso_not0
264 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 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 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 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 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 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 harmonics%slm(ia, iso)
303
304 END DO ! ia
305 END DO ! iso
306 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 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 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 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 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 REAL(dp), ALLOCATABLE, DIMENSION(:) :: a1, a2, gexp, r1, r2
339 REAL(dp), DIMENSION(:, :), POINTER :: slm
340 REAL(dp), DIMENSION(:, :, :), POINTER :: dslm_dxyz
341
342 NULLIFY (slm, dslm_dxyz)
343
344 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 zet=tau_cache%zet)
350 CALL get_paw_basis_info(basis_1c, o2nindex=tau_cache%o2nindex, &
351 n2oindex=tau_cache%n2oindex, &
352 nsatbas=tau_cache%nsatbas)
353
354 tau_cache%nr = grid_atom%nr
355 tau_cache%na = grid_atom%ng_sphere
356 slm => harmonics%slm
357 dslm_dxyz => harmonics%dslm_dxyz
358
359 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 r1(tau_cache%nr), r2(tau_cache%nr))
362 tau_cache%grad = 0.0_dp
363
364 DO iset = 1, tau_cache%nset
365 DO ipgf = 1, tau_cache%npgf(iset)
366 starti = (iset - 1)*tau_cache%maxso + &
367 (ipgf - 1)*nsoset(tau_cache%lmax(iset))
368 gexp(1:tau_cache%nr) = exp(-tau_cache%zet(ipgf, iset)* &
369 grid_atom%rad2(1:tau_cache%nr))
370 DO iso = nsoset(tau_cache%lmin(iset) - 1) + 1, nsoset(tau_cache%lmax(iset))
371 ip = tau_cache%o2nindex(starti + iso)
372 IF (ip == 0) cycle
373 l = indso(1, iso)
374
375 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 grid_atom%rad2(1:tau_cache%nr)*r1(1:tau_cache%nr)
378
379 DO dir = 1, 3
380 a1(1:tau_cache%na) = dslm_dxyz(dir, 1:tau_cache%na, iso)
381 a2(1:tau_cache%na) = harmonics%a(dir, 1:tau_cache%na)*slm(1:tau_cache%na, iso)
382 DO ir = 1, tau_cache%nr
383 DO ia = 1, tau_cache%na
384 igrid = ia + (ir - 1)*tau_cache%na
385 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 DEALLOCATE (a1, a2, gexp, r1, r2)
394
395 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 SUBROUTINE release_tau_basis_cache(tau_cache)
402
403 TYPE(tau_basis_cache_type), INTENT(INOUT) :: tau_cache
404
405 IF (ALLOCATED(tau_cache%grad)) DEALLOCATE (tau_cache%grad)
406 IF (ASSOCIATED(tau_cache%n2oindex)) DEALLOCATE (tau_cache%n2oindex)
407 IF (ASSOCIATED(tau_cache%o2nindex)) DEALLOCATE (tau_cache%o2nindex)
408 NULLIFY (tau_cache%lmax, tau_cache%lmin, tau_cache%n2oindex, tau_cache%npgf, &
409 tau_cache%zet, tau_cache%o2nindex)
410 tau_cache%maxso = 0
411 tau_cache%na = 0
412 tau_cache%nr = 0
413 tau_cache%nsatbas = 0
414 tau_cache%nset = 0
415
416 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 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 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: work
442
443 CALL timeset(routinen, handle)
444
445 cpassert(ALLOCATED(tau_cache%grad))
446
447 !zeroing tau, assuming it is already allocated
448 tau_h = 0.0_dp
449 tau_s = 0.0_dp
450
451 nr = tau_cache%nr
452 na = tau_cache%na
453 nbas = tau_cache%nsatbas
454 ngrid = na*nr
455 ALLOCATE (work(ngrid, nbas))
456
457 DO ispin = 1, nspins
458 DO dir = 1, 3
459 CALL dgemm('N', 'T', ngrid, nbas, nbas, 0.5_dp, tau_cache%grad(:, :, dir), &
460 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!$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 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!$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 DEALLOCATE (work)
496
497 CALL timestop(handle)
498
499 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 SUBROUTINE calc_rho_nlcc(grid_atom, nspins, grad_func, &
515 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 cpassert(ASSOCIATED(rho_h))
530 cpassert(ASSOCIATED(rho_s))
531 IF (grad_func) THEN
532 cpassert(ASSOCIATED(drho_h))
533 cpassert(ASSOCIATED(drho_s))
534 END IF
535
536 na = grid_atom%ng_sphere
537 rad = grid_atom%rad(ir)
538 urad = grid_atom%oorad2l(ir, 1)
539
540 xsp = real(nspins, kind=dp)
541 rho = rho_nlcc(ir)/xsp
542 DO ispin = 1, nspins
543 rho_h(1:na, ir, ispin) = rho_h(1:na, ir, ispin) + rho
544 rho_s(1:na, ir, ispin) = rho_s(1:na, ir, ispin) + rho
545 END DO ! ispin
546
547 IF (grad_func) THEN
548 drho = drho_nlcc(ir)/xsp
549 DO ispin = 1, nspins
550 DO ia = 1, na
551 IF (grid_atom%azi(ia) == 0.0_dp) THEN
552 dx = 0.0_dp
553 dy = 0.0_dp
554 ELSE
555 dx = grid_atom%sin_pol(ia)*grid_atom%sin_azi(ia)
556 dy = grid_atom%sin_pol(ia)*grid_atom%cos_azi(ia)
557 END IF
558 dz = grid_atom%cos_pol(ia)
559 ! components of the gradient of rho1 hard
560 drho_h(1, ia, ir, ispin) = drho_h(1, ia, ir, ispin) + drho*dx
561 drho_h(2, ia, ir, ispin) = drho_h(2, ia, ir, ispin) + drho*dy
562 drho_h(3, ia, ir, ispin) = drho_h(3, ia, ir, ispin) + drho*dz
563 ! components of the gradient of rho1 soft
564 drho_s(1, ia, ir, ispin) = drho_s(1, ia, ir, ispin) + drho*dx
565 drho_s(2, ia, ir, ispin) = drho_s(2, ia, ir, ispin) + drho*dy
566 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 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 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 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 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 INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf
611 REAL(dp), ALLOCATABLE, DIMENSION(:) :: g1, g2
612 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: gg, gvg_h, gvg_s, matso_h, matso_s, vx
613 REAL(dp), DIMENSION(:, :), POINTER :: zet
614 REAL(dp), DIMENSION(:, :, :), POINTER :: my_cg
615
616 CALL timeset(routinen, handle)
617
618 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 nset=nset, zet=zet)
623
624 nr = grid_atom%nr
625 na = grid_atom%ng_sphere
626 my_cg => harmonics%my_CG
627 max_iso_not0 = harmonics%max_iso_not0
628 lmax_expansion = indso(1, max_iso_not0)
629 max_s_harm = harmonics%max_s_harm
630
631 ALLOCATE (g1(nr), g2(nr), gg(nr, 0:2*maxl))
632 ALLOCATE (gvg_h(na, 0:2*maxl), gvg_s(na, 0:2*maxl))
633 ALLOCATE (matso_h(nsoset(maxl), nsoset(maxl)), &
634 matso_s(nsoset(maxl), nsoset(maxl)))
635 ALLOCATE (vx(na, nr))
636 ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm))
637
638 g1 = 0.0_dp
639 g2 = 0.0_dp
640 m1 = 0
641 DO iset1 = 1, nset
642 n1 = nsoset(lmax(iset1))
643 m2 = 0
644 DO iset2 = 1, nset
645 CALL get_none0_cg_list(my_cg, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
646 max_s_harm, lmax_expansion, cg_list, cg_n_list, max_iso_not0_local)
647 cpassert(max_iso_not0_local <= max_iso_not0)
648
649 n2 = nsoset(lmax(iset2))
650 DO ipgf1 = 1, npgf(iset1)
651 ngau1 = n1*(ipgf1 - 1) + m1
652 size1 = nsoset(lmax(iset1)) - nsoset(lmin(iset1) - 1)
653 nngau1 = nsoset(lmin(iset1) - 1) + ngau1
654
655 g1(1:nr) = exp(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))
656 DO ipgf2 = 1, npgf(iset2)
657 ngau2 = n2*(ipgf2 - 1) + m2
658
659 g2(1:nr) = exp(-zet(ipgf2, iset2)*grid_atom%rad2(1:nr))
660 lmin12 = lmin(iset1) + lmin(iset2)
661 lmax12 = lmax(iset1) + lmax(iset2)
662
663 ! reduce expansion local densities
664 IF (lmin12 <= lmax_expansion) THEN
665
666 gg = 0.0_dp
667 IF (lmin12 == 0) THEN
668 gg(1:nr, lmin12) = g1(1:nr)*g2(1:nr)
669 ELSE
670 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 IF (lmax12 > lmax_expansion) lmax12 = lmax_expansion
675
676 DO l = lmin12 + 1, lmax12
677 gg(1:nr, l) = grid_atom%rad(1:nr)*gg(:, l - 1)
678 END DO
679
680 DO ispin = 1, nspins
681 ld = lmax12 + 1
682 DO ir = 1, nr
683 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 gg(1:nr, 0:lmax12), nr, 0.0_dp, gvg_h(1:na, 0:lmax12), na)
687 DO ir = 1, nr
688 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 gg(1:nr, 0:lmax12), nr, 0.0_dp, gvg_s(1:na, 0:lmax12), na)
692
693 matso_h = 0.0_dp
694 matso_s = 0.0_dp
695 DO iso = 1, max_iso_not0_local
696 DO icg = 1, cg_n_list(iso)
697 iso1 = cg_list(1, icg, iso)
698 iso2 = cg_list(2, icg, iso)
699 l = indso(1, iso1) + indso(1, iso2)
700
701 cpassert(l <= lmax_expansion)
702 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 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 harmonics%slm(ia, iso)
711 END DO
712 END DO
713 END DO
714
715 ! Write in the global matrix
716 DO ic = nsoset(lmin(iset2) - 1) + 1, nsoset(lmax(iset2))
717 iso1 = nsoset(lmin(iset1) - 1) + 1
718 iso2 = ngau2 + ic
719 CALL daxpy(size1, 1.0_dp, matso_h(iso1, ic), 1, &
720 int_hh(ispin)%r_coef(nngau1 + 1, iso2), 1)
721 CALL daxpy(size1, 1.0_dp, matso_s(iso1, ic), 1, &
722 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 m2 = m2 + maxso
732 END DO ! iset2
733 m1 = m1 + maxso
734 END DO ! iset1
735
736 DEALLOCATE (g1, g2, gg, matso_h, matso_s, gvg_s, gvg_h, vx)
737
738 DEALLOCATE (cg_list, cg_n_list)
739
740 CALL timestop(handle)
741
742 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 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 INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf
777 REAL(dp) :: urad
778 REAL(dp), ALLOCATABLE, DIMENSION(:) :: g1, g2
779 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: dgg, gg, gvxcg_h, gvxcg_s, matso_h, &
780 matso_s
781 REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: gvxgg_h, gvxgg_s
782 REAL(dp), DIMENSION(:, :), POINTER :: zet
783 REAL(dp), DIMENSION(:, :, :), POINTER :: my_cg
784 REAL(dp), DIMENSION(:, :, :, :), POINTER :: my_cg_dxyz
785
786 CALL timeset(routinen, handle)
787
788 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 nset=nset, zet=zet)
793
794 nr = grid_atom%nr
795 na = grid_atom%ng_sphere
796 my_cg => harmonics%my_CG
797 my_cg_dxyz => harmonics%my_CG_dxyz
798 max_iso_not0 = harmonics%max_iso_not0
799 lmax_expansion = indso(1, max_iso_not0)
800 max_s_harm = harmonics%max_s_harm
801
802 ALLOCATE (g1(nr), g2(nr), gg(nr, 0:2*maxl), dgg(nr, 0:2*maxl))
803 ALLOCATE (gvxcg_h(na, 0:2*maxl), gvxcg_s(na, 0:2*maxl))
804 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 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 matso_s(nsoset(maxl), nsoset(maxl)))
810
811 DO ispin = 1, nspins
812
813 g1 = 0.0_dp
814 g2 = 0.0_dp
815 m1 = 0
816 DO iset1 = 1, nset
817 n1 = nsoset(lmax(iset1))
818 m2 = 0
819 DO iset2 = 1, nset
820 CALL get_none0_cg_list(my_cg, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
821 max_s_harm, lmax_expansion, cg_list, cg_n_list, max_iso_not0_local)
822 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 max_s_harm, lmax_expansion, dcg_list, dcg_n_list, dmax_iso_not0_local)
825
826 n2 = nsoset(lmax(iset2))
827 DO ipgf1 = 1, npgf(iset1)
828 ngau1 = n1*(ipgf1 - 1) + m1
829 size1 = nsoset(lmax(iset1)) - nsoset(lmin(iset1) - 1)
830 nngau1 = nsoset(lmin(iset1) - 1) + ngau1
831
832 g1(1:nr) = exp(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))
833 DO ipgf2 = 1, npgf(iset2)
834 ngau2 = n2*(ipgf2 - 1) + m2
835
836 g2(1:nr) = exp(-zet(ipgf2, iset2)*grid_atom%rad2(1:nr))
837 lmin12 = lmin(iset1) + lmin(iset2)
838 lmax12 = lmax(iset1) + lmax(iset2)
839
840 !test reduce expansion local densities
841 IF (lmin12 <= lmax_expansion) THEN
842
843 gg = 0.0_dp
844 dgg = 0.0_dp
845
846 IF (lmin12 == 0) THEN
847 gg(1:nr, lmin12) = g1(1:nr)*g2(1:nr)
848 ELSE
849 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 IF (lmax12 > lmax_expansion) lmax12 = lmax_expansion
854
855 DO l = lmin12 + 1, lmax12
856 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 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 gg(1:nr, lmax12)
863
864 gvxcg_h = 0.0_dp
865 gvxcg_s = 0.0_dp
866 gvxgg_h = 0.0_dp
867 gvxgg_s = 0.0_dp
868
869 ! Cross Term
870 DO l = lmin12, lmax12
871 DO ia = 1, na
872 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 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 vxg_s(3, ia, ir, ispin)*harmonics%a(3, ia))
886
887 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 gg(ir, l)*urad
892
893 gvxgg_h(2, ia, l) = gvxgg_h(2, ia, l) + &
894 vxg_h(2, ia, ir, ispin)* &
895 gg(ir, l)*urad
896
897 gvxgg_h(3, ia, l) = gvxgg_h(3, ia, l) + &
898 vxg_h(3, ia, ir, ispin)* &
899 gg(ir, l)*urad
900
901 gvxgg_s(1, ia, l) = gvxgg_s(1, ia, l) + &
902 vxg_s(1, ia, ir, ispin)* &
903 gg(ir, l)*urad
904
905 gvxgg_s(2, ia, l) = gvxgg_s(2, ia, l) + &
906 vxg_s(2, ia, ir, ispin)* &
907 gg(ir, l)*urad
908
909 gvxgg_s(3, ia, l) = gvxgg_s(3, ia, l) + &
910 vxg_s(3, ia, ir, ispin)* &
911 gg(ir, l)*urad
912
913 END DO ! ir
914 END DO ! ia
915 END DO ! l
916
917 matso_h = 0.0_dp
918 matso_s = 0.0_dp
919 DO iso = 1, max_iso_not0_local
920 DO icg = 1, cg_n_list(iso)
921 iso1 = cg_list(1, icg, iso)
922 iso2 = cg_list(2, icg, iso)
923
924 l = indso(1, iso1) + indso(1, iso2)
925
926 !test reduce expansion local densities
927 cpassert(l <= lmax_expansion)
928 DO ia = 1, na
929 matso_h(iso1, iso2) = matso_h(iso1, iso2) + &
930 gvxcg_h(ia, l)* &
931 harmonics%slm(ia, iso)* &
932 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 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 DO iso = 1, dmax_iso_not0_local
946 DO icg = 1, dcg_n_list(iso)
947 iso1 = dcg_list(1, icg, iso)
948 iso2 = dcg_list(2, icg, iso)
949
950 l = indso(1, iso1) + indso(1, iso2)
951 !test reduce expansion local densities
952 cpassert(l <= lmax_expansion)
953 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 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 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 DO ic = nsoset(lmin(iset2) - 1) + 1, nsoset(lmax(iset2))
977 iso1 = nsoset(lmin(iset1) - 1) + 1
978 iso2 = ngau2 + ic
979 CALL daxpy(size1, 1.0_dp, matso_h(iso1, ic), 1, &
980 int_hh(ispin)%r_coef(nngau1 + 1, iso2), 1)
981 CALL daxpy(size1, 1.0_dp, matso_s(iso1, ic), 1, &
982 int_ss(ispin)%r_coef(nngau1 + 1, iso2), 1)
983 END DO
984
985 END DO ! ipfg2
986 END DO ! ipfg1
987 m2 = m2 + maxso
988 END DO ! iset2
989 m1 = m1 + maxso
990 END DO ! iset1
991 END DO ! ispin
992
993 DEALLOCATE (g1, g2, gg, dgg, matso_h, matso_s, gvxcg_h, gvxcg_s, gvxgg_h, gvxgg_s)
994 DEALLOCATE (cg_list, cg_n_list, dcg_list, dcg_n_list)
995
996 CALL timestop(handle)
997
998 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 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 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: int_h, int_s, weighted_grad
1024
1025 CALL timeset(routinen, handle)
1026
1027 cpassert(ALLOCATED(tau_cache%grad))
1028 cpassert(ASSOCIATED(tau_cache%n2oindex))
1029
1030 nr = tau_cache%nr
1031 na = tau_cache%na
1032 nbas = tau_cache%nsatbas
1033 ngrid = na*nr
1034 max_old_basis = maxval(tau_cache%n2oindex)
1035 ALLOCATE (int_h(nbas, nbas), int_s(nbas, nbas), weighted_grad(ngrid, nbas))
1036
1037 DO ispin = 1, nspins
1038 cpassert(SIZE(int_hh(ispin)%r_coef, 1) >= max_old_basis)
1039 cpassert(SIZE(int_hh(ispin)%r_coef, 2) >= max_old_basis)
1040 cpassert(SIZE(int_ss(ispin)%r_coef, 1) >= max_old_basis)
1041 cpassert(SIZE(int_ss(ispin)%r_coef, 2) >= max_old_basis)
1042 int_h = 0.0_dp
1043 int_s = 0.0_dp
1044 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!$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 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!$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 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!$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 DEALLOCATE (int_h, int_s, weighted_grad)
1095
1096 CALL timestop(handle)
1097
1098 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 Polynomial order (2n)
1105!> \param alpha Gaussian exponent
1106! **************************************************************************************************
1107 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 cpassert(SIZE(fun) == SIZE(r2))
1118
1119 fun(:) = 1.0_dp
1120 DO i = 1, order
1121 pval = alpha**i/fac(i)
1122 fun(:) = fun(:) + pval*r2**i
1123 END DO
1124 fun(:) = fun(:)*exp(-alpha*r2(:))
1125
1126 END SUBROUTINE calc_weight_function
1127
1128! **************************************************************************************************
1129
1130END MODULE qs_vxc_atom_utils
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition grid_cpu_task_list.c:225
external_potential_types::get_potential
Definition external_potential_types.F:279
qs_harmonics_atom::get_none0_cg_list
Definition qs_harmonics_atom.F:54
basis_set_types
Definition basis_set_types.F:15
basis_set_types::get_gto_basis_set
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum, ccon)
...
Definition basis_set_types.F:643
external_potential_types
Definition of the atomic potential types.
Definition external_potential_types.F:14
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::fac
real(kind=dp), dimension(0:maxfac), parameter, public fac
Definition mathconstants.F:37
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::nsoset
integer, dimension(:), allocatable, public nsoset
Definition orbital_pointers.F:42
orbital_pointers::indso
integer, dimension(:, :), allocatable, public indso
Definition orbital_pointers.F:43
paw_basis_types
Definition paw_basis_types.F:13
paw_basis_types::get_paw_basis_info
subroutine, public get_paw_basis_info(basis_1c, o2nindex, n2oindex, nsatbas)
Return some info on the PAW basis derived from a GTO basis set.
Definition paw_basis_types.F:40
qs_grid_atom
Definition qs_grid_atom.F:8
qs_harmonics_atom
Definition qs_harmonics_atom.F:8
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::has_nlcc
logical function, public has_nlcc(qs_kind_set)
finds if a given qs run needs to use nlcc
Definition qs_kind_types.F:4283
qs_rho_atom_types
Definition qs_rho_atom_types.F:8
qs_vxc_atom_utils
Support routines for integrals of the Vxc/Fxc/Gxc potentials calculated for the atomic density in the...
Definition qs_vxc_atom_utils.F:12
qs_vxc_atom_utils::calc_rho_angular
subroutine, public calc_rho_angular(grid_atom, harmonics, nspins, grad_func, ir, r_h, r_s, rho_h, rho_s, dr_h, dr_s, r_h_d, r_s_d, drho_h, drho_s)
...
Definition qs_vxc_atom_utils.F:219
qs_vxc_atom_utils::release_tau_basis_cache
subroutine, public release_tau_basis_cache(tau_cache)
Release precomputed GAPW meta-GGA tau factors.
Definition qs_vxc_atom_utils.F:402
qs_vxc_atom_utils::create_tau_basis_cache
subroutine, public create_tau_basis_cache(tau_cache, grid_atom, basis_1c, harmonics)
Precompute radial and angular factors for GAPW meta-GGA tau contractions.
Definition qs_vxc_atom_utils.F:330
qs_vxc_atom_utils::gavxcgb_nogc
subroutine, public gavxcgb_nogc(vxc_h, vxc_s, int_hh, int_ss, grid_atom, basis_1c, harmonics, nspins)
...
Definition qs_vxc_atom_utils.F:595
qs_vxc_atom_utils::gavxcgb_gc
subroutine, public gavxcgb_gc(vxc_h, vxc_s, vxg_h, vxg_s, int_hh, int_ss, grid_atom, basis_1c, harmonics, nspins)
...
Definition qs_vxc_atom_utils.F:759
qs_vxc_atom_utils::calc_weight_function
subroutine, public calc_weight_function(fun, r2, order, alpha)
Calculates the radial weight function for accurate XC integration.
Definition qs_vxc_atom_utils.F:1108
qs_vxc_atom_utils::calc_rho_nlcc
subroutine, public calc_rho_nlcc(grid_atom, nspins, grad_func, ir, rho_nlcc, rho_h, rho_s, drho_nlcc, drho_h, drho_s)
...
Definition qs_vxc_atom_utils.F:516
qs_vxc_atom_utils::dgavtaudgb
subroutine, public dgavtaudgb(vtau_h, vtau_s, int_hh, int_ss, tau_cache, nspins)
Integrates 0.5 * grad_ga .dot. (V_tau * grad_gb) on the atomic grid for meta-GGA.
Definition qs_vxc_atom_utils.F:1012
qs_vxc_atom_utils::calc_tau_atom
subroutine, public calc_tau_atom(tau_h, tau_s, rho_atom, tau_cache, nspins)
Computes tau hard and soft on the atomic grids for meta-GGA calculations.
Definition qs_vxc_atom_utils.F:430
qs_vxc_atom_utils::gapw_atom_grid_support_radius
real(dp) function, public gapw_atom_grid_support_radius(grid_atom, rho_h, rho_s, drho_h, drho_s, tau_h, tau_s)
Compact support radius of one hard-minus-soft atom-grid field.
Definition qs_vxc_atom_utils.F:158
qs_vxc_atom_utils::evaluate_nlcc_primitive_fields
subroutine, public evaluate_nlcc_primitive_fields(point, center, gth_potential, sgp_potential, rho, gradient, hessian)
Evaluate an NLCC density and its first two Cartesian derivatives.
Definition qs_vxc_atom_utils.F:70
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
external_potential_types::gth_potential_type
Definition external_potential_types.F:121
external_potential_types::sgp_potential_type
Definition external_potential_types.F:179
qs_grid_atom::grid_atom_type
Definition qs_grid_atom.F:44
qs_harmonics_atom::harmonics_atom_type
Definition qs_harmonics_atom.F:32
qs_rho_atom_types::rho_atom_coeff
Definition qs_rho_atom_types.F:19
qs_rho_atom_types::rho_atom_type
Definition qs_rho_atom_types.F:23
qs_vxc_atom_utils::tau_basis_cache_type
Definition qs_vxc_atom_utils.F:37