(git:6d276e9)
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 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
54CONTAINS
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 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 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 REAL(dp), DIMENSION(:), POINTER :: a_nlcc, alpha_nlcc, c_nlcc
83 REAL(dp), DIMENSION(:, :), POINTER :: cval_nlcc
84
85 NULLIFY (a_nlcc, alpha_nlcc, c_nlcc, cval_nlcc, nct_nlcc)
86 rho = 0.0_dp
87 rho_x = 0.0_dp
88 rho_xx = 0.0_dp
89 displacement = point - center
90 r2 = dot_product(displacement, displacement)
91
92 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 nct_nlcc=nct_nlcc, cval_nlcc=cval_nlcc)
96 IF (nlcc_present) THEN
97 DO iexp = 1, nexp_nlcc
98 alpha = alpha_nlcc(iexp)
99 beta = 0.5_dp/(alpha*alpha)
100 scaled_r2 = r2/(alpha*alpha)
101 exponential = exp(-0.5_dp*scaled_r2)
102 DO ic = 1, nct_nlcc(iexp)
103 power = ic - 1
104 poly = cval_nlcc(ic, iexp)*scaled_r2**power
105 dpoly = 0.0_dp
106 IF (power > 0) THEN
107 dpoly = cval_nlcc(ic, iexp)*real(power, dp)* &
108 scaled_r2**(power - 1)/(alpha*alpha)
109 END IF
110 d2poly = 0.0_dp
111 IF (power > 1) THEN
112 d2poly = cval_nlcc(ic, iexp)*real(power*(power - 1), dp)* &
113 scaled_r2**(power - 2)/(alpha**4)
114 END IF
115 rho = rho + exponential*poly
116 rho_x = rho_x + exponential*(dpoly - beta*poly)
117 rho_xx = rho_xx + exponential*(d2poly - 2.0_dp*beta*dpoly + beta*beta*poly)
118 END DO
119 END DO
120 END IF
121 ELSE IF (ASSOCIATED(sgp_potential)) THEN
122 CALL get_potential(sgp_potential, has_nlcc=has_sgp_nlcc, n_nlcc=n_nlcc, &
123 a_nlcc=a_nlcc, c_nlcc=c_nlcc)
124 IF (has_sgp_nlcc) THEN
125 DO iexp = 1, n_nlcc
126 exponential = exp(-a_nlcc(iexp)*r2)
127 rho = rho + c_nlcc(iexp)*exponential
128 rho_x = rho_x - a_nlcc(iexp)*c_nlcc(iexp)*exponential
129 rho_xx = rho_xx + a_nlcc(iexp)**2*c_nlcc(iexp)*exponential
130 END DO
131 END IF
132 END IF
133
134 gradient = 2.0_dp*rho_x*displacement
135 DO idir = 1, 3
136 DO jdir = 1, 3
137 hessian(idir, jdir) = 4.0_dp*rho_xx*displacement(idir)*displacement(jdir)
138 END DO
139 hessian(idir, idir) = hessian(idir, idir) + 2.0_dp*rho_x
140 END DO
141 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 FUNCTION gapw_atom_grid_support_radius( &
155 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 cutoff = 0.0_dp
167 descending = grid_atom%rad(1) > grid_atom%rad(grid_atom%nr)
168 IF (descending) THEN
169 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 maxval(abs(tau_h(:, ir, :) - tau_s(:, ir, :))))
174 IF (amplitude > 1.0e-12_dp) THEN
175 support_ir = max(1, ir - 2)
176 cutoff = grid_atom%rad(support_ir)
177 RETURN
178 END IF
179 END DO
180 ELSE
181 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 maxval(abs(tau_h(:, ir, :) - tau_s(:, ir, :))))
186 IF (amplitude > 1.0e-12_dp) THEN
187 support_ir = min(grid_atom%nr, ir + 2)
188 cutoff = grid_atom%rad(support_ir)
189 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 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 cpassert(ASSOCIATED(r_h))
233 cpassert(ASSOCIATED(r_s))
234 cpassert(ASSOCIATED(rho_h))
235 cpassert(ASSOCIATED(rho_s))
236 IF (grad_func) THEN
237 cpassert(ASSOCIATED(dr_h))
238 cpassert(ASSOCIATED(dr_s))
239 cpassert(ASSOCIATED(r_h_d))
240 cpassert(ASSOCIATED(r_s_d))
241 cpassert(ASSOCIATED(drho_h))
242 cpassert(ASSOCIATED(drho_s))
243 END IF
244
245 na = grid_atom%ng_sphere
246 rad = grid_atom%rad(ir)
247 urad = grid_atom%oorad2l(ir, 1)
248 DO ispin = 1, nspins
249 DO iso = 1, harmonics%max_iso_not0
250 DO ia = 1, na
251 rho_h(ia, ir, ispin) = rho_h(ia, ir, ispin) + &
252 r_h(ispin)%r_coef(ir, iso)*harmonics%slm(ia, iso)
253 rho_s(ia, ir, ispin) = rho_s(ia, ir, ispin) + &
254 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 IF (grad_func) THEN
260 DO ispin = 1, nspins
261 DO iso = 1, harmonics%max_iso_not0
262 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 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 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 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 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 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 harmonics%slm(ia, iso)
301
302 END DO ! ia
303 END DO ! iso
304 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 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 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 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 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 REAL(dp), ALLOCATABLE, DIMENSION(:) :: a1, a2, gexp, r1, r2
337 REAL(dp), DIMENSION(:, :), POINTER :: slm
338 REAL(dp), DIMENSION(:, :, :), POINTER :: dslm_dxyz
339
340 NULLIFY (slm, dslm_dxyz)
341
342 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 zet=tau_cache%zet)
348 CALL get_paw_basis_info(basis_1c, o2nindex=tau_cache%o2nindex, &
349 n2oindex=tau_cache%n2oindex, &
350 nsatbas=tau_cache%nsatbas)
351
352 tau_cache%nr = grid_atom%nr
353 tau_cache%na = grid_atom%ng_sphere
354 slm => harmonics%slm
355 dslm_dxyz => harmonics%dslm_dxyz
356
357 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 r1(tau_cache%nr), r2(tau_cache%nr))
360 tau_cache%grad = 0.0_dp
361
362 DO iset = 1, tau_cache%nset
363 DO ipgf = 1, tau_cache%npgf(iset)
364 starti = (iset - 1)*tau_cache%maxso + &
365 (ipgf - 1)*nsoset(tau_cache%lmax(iset))
366 gexp(1:tau_cache%nr) = exp(-tau_cache%zet(ipgf, iset)* &
367 grid_atom%rad2(1:tau_cache%nr))
368 DO iso = nsoset(tau_cache%lmin(iset) - 1) + 1, nsoset(tau_cache%lmax(iset))
369 ip = tau_cache%o2nindex(starti + iso)
370 IF (ip == 0) cycle
371 l = indso(1, iso)
372
373 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 grid_atom%rad2(1:tau_cache%nr)*r1(1:tau_cache%nr)
376
377 DO dir = 1, 3
378 a1(1:tau_cache%na) = dslm_dxyz(dir, 1:tau_cache%na, iso)
379 a2(1:tau_cache%na) = harmonics%a(dir, 1:tau_cache%na)*slm(1:tau_cache%na, iso)
380 DO ir = 1, tau_cache%nr
381 DO ia = 1, tau_cache%na
382 igrid = ia + (ir - 1)*tau_cache%na
383 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 DEALLOCATE (a1, a2, gexp, r1, r2)
392
393 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 SUBROUTINE release_tau_basis_cache(tau_cache)
400
401 TYPE(tau_basis_cache_type), INTENT(INOUT) :: tau_cache
402
403 IF (ALLOCATED(tau_cache%grad)) DEALLOCATE (tau_cache%grad)
404 IF (ASSOCIATED(tau_cache%n2oindex)) DEALLOCATE (tau_cache%n2oindex)
405 IF (ASSOCIATED(tau_cache%o2nindex)) DEALLOCATE (tau_cache%o2nindex)
406 NULLIFY (tau_cache%lmax, tau_cache%lmin, tau_cache%n2oindex, tau_cache%npgf, &
407 tau_cache%zet, tau_cache%o2nindex)
408 tau_cache%maxso = 0
409 tau_cache%na = 0
410 tau_cache%nr = 0
411 tau_cache%nsatbas = 0
412 tau_cache%nset = 0
413
414 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 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 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: work
440
441 CALL timeset(routinen, handle)
442
443 cpassert(ALLOCATED(tau_cache%grad))
444
445 !zeroing tau, assuming it is already allocated
446 tau_h = 0.0_dp
447 tau_s = 0.0_dp
448
449 nr = tau_cache%nr
450 na = tau_cache%na
451 nbas = tau_cache%nsatbas
452 ngrid = na*nr
453 ALLOCATE (work(ngrid, nbas))
454
455 DO ispin = 1, nspins
456 DO dir = 1, 3
457 CALL dgemm('N', 'T', ngrid, nbas, nbas, 0.5_dp, tau_cache%grad(:, :, dir), &
458 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!$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 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!$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 DEALLOCATE (work)
494
495 CALL timestop(handle)
496
497 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 SUBROUTINE calc_rho_nlcc(grid_atom, nspins, grad_func, &
513 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 cpassert(ASSOCIATED(rho_h))
528 cpassert(ASSOCIATED(rho_s))
529 IF (grad_func) THEN
530 cpassert(ASSOCIATED(drho_h))
531 cpassert(ASSOCIATED(drho_s))
532 END IF
533
534 na = grid_atom%ng_sphere
535 rad = grid_atom%rad(ir)
536 urad = grid_atom%oorad2l(ir, 1)
537
538 xsp = real(nspins, kind=dp)
539 rho = rho_nlcc(ir)/xsp
540 DO ispin = 1, nspins
541 rho_h(1:na, ir, ispin) = rho_h(1:na, ir, ispin) + rho
542 rho_s(1:na, ir, ispin) = rho_s(1:na, ir, ispin) + rho
543 END DO ! ispin
544
545 IF (grad_func) THEN
546 drho = drho_nlcc(ir)/xsp
547 DO ispin = 1, nspins
548 DO ia = 1, na
549 IF (grid_atom%azi(ia) == 0.0_dp) THEN
550 dx = 0.0_dp
551 dy = 0.0_dp
552 ELSE
553 dx = grid_atom%sin_pol(ia)*grid_atom%sin_azi(ia)
554 dy = grid_atom%sin_pol(ia)*grid_atom%cos_azi(ia)
555 END IF
556 dz = grid_atom%cos_pol(ia)
557 ! components of the gradient of rho1 hard
558 drho_h(1, ia, ir, ispin) = drho_h(1, ia, ir, ispin) + drho*dx
559 drho_h(2, ia, ir, ispin) = drho_h(2, ia, ir, ispin) + drho*dy
560 drho_h(3, ia, ir, ispin) = drho_h(3, ia, ir, ispin) + drho*dz
561 ! components of the gradient of rho1 soft
562 drho_s(1, ia, ir, ispin) = drho_s(1, ia, ir, ispin) + drho*dx
563 drho_s(2, ia, ir, ispin) = drho_s(2, ia, ir, ispin) + drho*dy
564 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 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 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 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 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 INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf
609 REAL(dp), ALLOCATABLE, DIMENSION(:) :: g1, g2
610 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: gg, gvg_h, gvg_s, matso_h, matso_s, vx
611 REAL(dp), DIMENSION(:, :), POINTER :: zet
612 REAL(dp), DIMENSION(:, :, :), POINTER :: my_cg
613
614 CALL timeset(routinen, handle)
615
616 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 nset=nset, zet=zet)
621
622 nr = grid_atom%nr
623 na = grid_atom%ng_sphere
624 my_cg => harmonics%my_CG
625 max_iso_not0 = harmonics%max_iso_not0
626 lmax_expansion = indso(1, max_iso_not0)
627 max_s_harm = harmonics%max_s_harm
628
629 ALLOCATE (g1(nr), g2(nr), gg(nr, 0:2*maxl))
630 ALLOCATE (gvg_h(na, 0:2*maxl), gvg_s(na, 0:2*maxl))
631 ALLOCATE (matso_h(nsoset(maxl), nsoset(maxl)), &
632 matso_s(nsoset(maxl), nsoset(maxl)))
633 ALLOCATE (vx(na, nr))
634 ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm))
635
636 g1 = 0.0_dp
637 g2 = 0.0_dp
638 m1 = 0
639 DO iset1 = 1, nset
640 n1 = nsoset(lmax(iset1))
641 m2 = 0
642 DO iset2 = 1, nset
643 CALL get_none0_cg_list(my_cg, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
644 max_s_harm, lmax_expansion, cg_list, cg_n_list, max_iso_not0_local)
645 cpassert(max_iso_not0_local <= max_iso_not0)
646
647 n2 = nsoset(lmax(iset2))
648 DO ipgf1 = 1, npgf(iset1)
649 ngau1 = n1*(ipgf1 - 1) + m1
650 size1 = nsoset(lmax(iset1)) - nsoset(lmin(iset1) - 1)
651 nngau1 = nsoset(lmin(iset1) - 1) + ngau1
652
653 g1(1:nr) = exp(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))
654 DO ipgf2 = 1, npgf(iset2)
655 ngau2 = n2*(ipgf2 - 1) + m2
656
657 g2(1:nr) = exp(-zet(ipgf2, iset2)*grid_atom%rad2(1:nr))
658 lmin12 = lmin(iset1) + lmin(iset2)
659 lmax12 = lmax(iset1) + lmax(iset2)
660
661 ! reduce expansion local densities
662 IF (lmin12 <= lmax_expansion) THEN
663
664 gg = 0.0_dp
665 IF (lmin12 == 0) THEN
666 gg(1:nr, lmin12) = g1(1:nr)*g2(1:nr)
667 ELSE
668 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 IF (lmax12 > lmax_expansion) lmax12 = lmax_expansion
673
674 DO l = lmin12 + 1, lmax12
675 gg(1:nr, l) = grid_atom%rad(1:nr)*gg(:, l - 1)
676 END DO
677
678 DO ispin = 1, nspins
679 ld = lmax12 + 1
680 DO ir = 1, nr
681 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 gg(1:nr, 0:lmax12), nr, 0.0_dp, gvg_h(1:na, 0:lmax12), na)
685 DO ir = 1, nr
686 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 gg(1:nr, 0:lmax12), nr, 0.0_dp, gvg_s(1:na, 0:lmax12), na)
690
691 matso_h = 0.0_dp
692 matso_s = 0.0_dp
693 DO iso = 1, max_iso_not0_local
694 DO icg = 1, cg_n_list(iso)
695 iso1 = cg_list(1, icg, iso)
696 iso2 = cg_list(2, icg, iso)
697 l = indso(1, iso1) + indso(1, iso2)
698
699 cpassert(l <= lmax_expansion)
700 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 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 harmonics%slm(ia, iso)
709 END DO
710 END DO
711 END DO
712
713 ! Write in the global matrix
714 DO ic = nsoset(lmin(iset2) - 1) + 1, nsoset(lmax(iset2))
715 iso1 = nsoset(lmin(iset1) - 1) + 1
716 iso2 = ngau2 + ic
717 CALL daxpy(size1, 1.0_dp, matso_h(iso1, ic), 1, &
718 int_hh(ispin)%r_coef(nngau1 + 1, iso2), 1)
719 CALL daxpy(size1, 1.0_dp, matso_s(iso1, ic), 1, &
720 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 m2 = m2 + maxso
730 END DO ! iset2
731 m1 = m1 + maxso
732 END DO ! iset1
733
734 DEALLOCATE (g1, g2, gg, matso_h, matso_s, gvg_s, gvg_h, vx)
735
736 DEALLOCATE (cg_list, cg_n_list)
737
738 CALL timestop(handle)
739
740 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 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 INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf
775 REAL(dp) :: urad
776 REAL(dp), ALLOCATABLE, DIMENSION(:) :: g1, g2
777 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: dgg, gg, gvxcg_h, gvxcg_s, matso_h, &
778 matso_s
779 REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: gvxgg_h, gvxgg_s
780 REAL(dp), DIMENSION(:, :), POINTER :: zet
781 REAL(dp), DIMENSION(:, :, :), POINTER :: my_cg
782 REAL(dp), DIMENSION(:, :, :, :), POINTER :: my_cg_dxyz
783
784 CALL timeset(routinen, handle)
785
786 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 nset=nset, zet=zet)
791
792 nr = grid_atom%nr
793 na = grid_atom%ng_sphere
794 my_cg => harmonics%my_CG
795 my_cg_dxyz => harmonics%my_CG_dxyz
796 max_iso_not0 = harmonics%max_iso_not0
797 lmax_expansion = indso(1, max_iso_not0)
798 max_s_harm = harmonics%max_s_harm
799
800 ALLOCATE (g1(nr), g2(nr), gg(nr, 0:2*maxl), dgg(nr, 0:2*maxl))
801 ALLOCATE (gvxcg_h(na, 0:2*maxl), gvxcg_s(na, 0:2*maxl))
802 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 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 matso_s(nsoset(maxl), nsoset(maxl)))
808
809 DO ispin = 1, nspins
810
811 g1 = 0.0_dp
812 g2 = 0.0_dp
813 m1 = 0
814 DO iset1 = 1, nset
815 n1 = nsoset(lmax(iset1))
816 m2 = 0
817 DO iset2 = 1, nset
818 CALL get_none0_cg_list(my_cg, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
819 max_s_harm, lmax_expansion, cg_list, cg_n_list, max_iso_not0_local)
820 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 max_s_harm, lmax_expansion, dcg_list, dcg_n_list, dmax_iso_not0_local)
823
824 n2 = nsoset(lmax(iset2))
825 DO ipgf1 = 1, npgf(iset1)
826 ngau1 = n1*(ipgf1 - 1) + m1
827 size1 = nsoset(lmax(iset1)) - nsoset(lmin(iset1) - 1)
828 nngau1 = nsoset(lmin(iset1) - 1) + ngau1
829
830 g1(1:nr) = exp(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))
831 DO ipgf2 = 1, npgf(iset2)
832 ngau2 = n2*(ipgf2 - 1) + m2
833
834 g2(1:nr) = exp(-zet(ipgf2, iset2)*grid_atom%rad2(1:nr))
835 lmin12 = lmin(iset1) + lmin(iset2)
836 lmax12 = lmax(iset1) + lmax(iset2)
837
838 !test reduce expansion local densities
839 IF (lmin12 <= lmax_expansion) THEN
840
841 gg = 0.0_dp
842 dgg = 0.0_dp
843
844 IF (lmin12 == 0) THEN
845 gg(1:nr, lmin12) = g1(1:nr)*g2(1:nr)
846 ELSE
847 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 IF (lmax12 > lmax_expansion) lmax12 = lmax_expansion
852
853 DO l = lmin12 + 1, lmax12
854 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 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 gg(1:nr, lmax12)
861
862 gvxcg_h = 0.0_dp
863 gvxcg_s = 0.0_dp
864 gvxgg_h = 0.0_dp
865 gvxgg_s = 0.0_dp
866
867 ! Cross Term
868 DO l = lmin12, lmax12
869 DO ia = 1, na
870 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 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 vxg_s(3, ia, ir, ispin)*harmonics%a(3, ia))
884
885 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 gg(ir, l)*urad
890
891 gvxgg_h(2, ia, l) = gvxgg_h(2, ia, l) + &
892 vxg_h(2, ia, ir, ispin)* &
893 gg(ir, l)*urad
894
895 gvxgg_h(3, ia, l) = gvxgg_h(3, ia, l) + &
896 vxg_h(3, ia, ir, ispin)* &
897 gg(ir, l)*urad
898
899 gvxgg_s(1, ia, l) = gvxgg_s(1, ia, l) + &
900 vxg_s(1, ia, ir, ispin)* &
901 gg(ir, l)*urad
902
903 gvxgg_s(2, ia, l) = gvxgg_s(2, ia, l) + &
904 vxg_s(2, ia, ir, ispin)* &
905 gg(ir, l)*urad
906
907 gvxgg_s(3, ia, l) = gvxgg_s(3, ia, l) + &
908 vxg_s(3, ia, ir, ispin)* &
909 gg(ir, l)*urad
910
911 END DO ! ir
912 END DO ! ia
913 END DO ! l
914
915 matso_h = 0.0_dp
916 matso_s = 0.0_dp
917 DO iso = 1, max_iso_not0_local
918 DO icg = 1, cg_n_list(iso)
919 iso1 = cg_list(1, icg, iso)
920 iso2 = cg_list(2, icg, iso)
921
922 l = indso(1, iso1) + indso(1, iso2)
923
924 !test reduce expansion local densities
925 cpassert(l <= lmax_expansion)
926 DO ia = 1, na
927 matso_h(iso1, iso2) = matso_h(iso1, iso2) + &
928 gvxcg_h(ia, l)* &
929 harmonics%slm(ia, iso)* &
930 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 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 DO iso = 1, dmax_iso_not0_local
944 DO icg = 1, dcg_n_list(iso)
945 iso1 = dcg_list(1, icg, iso)
946 iso2 = dcg_list(2, icg, iso)
947
948 l = indso(1, iso1) + indso(1, iso2)
949 !test reduce expansion local densities
950 cpassert(l <= lmax_expansion)
951 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 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 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 DO ic = nsoset(lmin(iset2) - 1) + 1, nsoset(lmax(iset2))
975 iso1 = nsoset(lmin(iset1) - 1) + 1
976 iso2 = ngau2 + ic
977 CALL daxpy(size1, 1.0_dp, matso_h(iso1, ic), 1, &
978 int_hh(ispin)%r_coef(nngau1 + 1, iso2), 1)
979 CALL daxpy(size1, 1.0_dp, matso_s(iso1, ic), 1, &
980 int_ss(ispin)%r_coef(nngau1 + 1, iso2), 1)
981 END DO
982
983 END DO ! ipfg2
984 END DO ! ipfg1
985 m2 = m2 + maxso
986 END DO ! iset2
987 m1 = m1 + maxso
988 END DO ! iset1
989 END DO ! ispin
990
991 DEALLOCATE (g1, g2, gg, dgg, matso_h, matso_s, gvxcg_h, gvxcg_s, gvxgg_h, gvxgg_s)
992 DEALLOCATE (cg_list, cg_n_list, dcg_list, dcg_n_list)
993
994 CALL timestop(handle)
995
996 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 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 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: int_h, int_s, weighted_grad
1022
1023 CALL timeset(routinen, handle)
1024
1025 cpassert(ALLOCATED(tau_cache%grad))
1026 cpassert(ASSOCIATED(tau_cache%n2oindex))
1027
1028 nr = tau_cache%nr
1029 na = tau_cache%na
1030 nbas = tau_cache%nsatbas
1031 ngrid = na*nr
1032 max_old_basis = maxval(tau_cache%n2oindex)
1033 ALLOCATE (int_h(nbas, nbas), int_s(nbas, nbas), weighted_grad(ngrid, nbas))
1034
1035 DO ispin = 1, nspins
1036 cpassert(SIZE(int_hh(ispin)%r_coef, 1) >= max_old_basis)
1037 cpassert(SIZE(int_hh(ispin)%r_coef, 2) >= max_old_basis)
1038 cpassert(SIZE(int_ss(ispin)%r_coef, 1) >= max_old_basis)
1039 cpassert(SIZE(int_ss(ispin)%r_coef, 2) >= max_old_basis)
1040 int_h = 0.0_dp
1041 int_s = 0.0_dp
1042 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!$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 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!$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 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!$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 DEALLOCATE (int_h, int_s, weighted_grad)
1093
1094 CALL timestop(handle)
1095
1096 END SUBROUTINE dgavtaudgb
1097
1098END 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:278
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
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:4286
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:217
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:400
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:328
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:593
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:757
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:514
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:1010
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:428
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:156
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:68
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
external_potential_types::gth_potential_type
Definition external_potential_types.F:116
external_potential_types::sgp_potential_type
Definition external_potential_types.F:174
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:36