(git:374b731)
Loading...
Searching...
No Matches
qs_integrate_potential_single.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Build up the plane wave density by collocating the primitive Gaussian
10!> functions (pgf).
11!> \par History
12!> Joost VandeVondele (02.2002)
13!> 1) rewrote collocate_pgf for increased accuracy and speed
14!> 2) collocate_core hack for PGI compiler
15!> 3) added multiple grid feature
16!> 4) new way to go over the grid
17!> Joost VandeVondele (05.2002)
18!> 1) prelim. introduction of the real space grid type
19!> JGH [30.08.02] multigrid arrays independent from potential
20!> JGH [17.07.03] distributed real space code
21!> JGH [23.11.03] refactoring and new loop ordering
22!> JGH [04.12.03] OpneMP parallelization of main loops
23!> Joost VandeVondele (12.2003)
24!> 1) modified to compute tau
25!> Joost removed incremental build feature
26!> Joost introduced map consistent
27!> Rewrote grid integration/collocation routines, [Joost VandeVondele,03.2007]
28!> \author Matthias Krack (03.04.2001)
29! **************************************************************************************************
30MODULE qs_integrate_potential_single
31 USE ao_util, ONLY: exp_radius_very_extended
32 USE atomic_kind_types, ONLY: atomic_kind_type,&
33 get_atomic_kind
34 USE atprop_types, ONLY: atprop_array_init,&
35 atprop_type
36 USE basis_set_types, ONLY: get_gto_basis_set,&
37 gto_basis_set_type
38 USE cell_types, ONLY: cell_type,&
39 pbc
40 USE cp_control_types, ONLY: dft_control_type
41 USE dbcsr_api, ONLY: dbcsr_get_block_p,&
42 dbcsr_type
43 USE external_potential_types, ONLY: get_potential,&
44 gth_potential_type
45 USE gaussian_gridlevels, ONLY: gaussian_gridlevel,&
46 gridlevel_info_type
47 USE grid_api, ONLY: integrate_pgf_product
48 USE kinds, ONLY: dp
49 USE lri_environment_types, ONLY: lri_kind_type
50 USE memory_utilities, ONLY: reallocate
51 USE message_passing, ONLY: mp_para_env_type
52 USE orbital_pointers, ONLY: coset,&
53 ncoset
54 USE particle_types, ONLY: particle_type
55 USE pw_env_types, ONLY: pw_env_get,&
56 pw_env_type
57 USE pw_types, ONLY: pw_r3d_rs_type
58 USE qs_environment_types, ONLY: get_qs_env,&
59 qs_environment_type
60 USE qs_force_types, ONLY: qs_force_type
61 USE qs_kind_types, ONLY: get_qs_kind,&
62 qs_kind_type
63 USE realspace_grid_types, ONLY: map_gaussian_here,&
64 realspace_grid_type,&
65 rs_grid_zero,&
66 transfer_pw2rs
67 USE rs_pw_interface, ONLY: potential_pw2rs
68 USE virial_types, ONLY: virial_type
69#include "./base/base_uses.f90"
70
71 IMPLICIT NONE
72
73 PRIVATE
74
75 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .false.
76
77 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_integrate_potential_single'
78
79! *** Public subroutines ***
80! *** Don't include this routines directly, use the interface to
81! *** qs_integrate_potential
82
83 PUBLIC :: integrate_v_rspace_one_center, &
84 integrate_v_rspace_diagonal, &
85 integrate_v_core_rspace, &
86 integrate_v_gaussian_rspace, &
87 integrate_ppl_rspace, &
88 integrate_rho_nlcc
89
90CONTAINS
91
92! **************************************************************************************************
93!> \brief computes the forces/virial due to the local pseudopotential
94!> \param rho_rspace ...
95!> \param qs_env ...
96! **************************************************************************************************
97 SUBROUTINE integrate_ppl_rspace(rho_rspace, qs_env)
98 TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_rspace
99 TYPE(qs_environment_type), POINTER :: qs_env
100
101 CHARACTER(len=*), PARAMETER :: routinen = 'integrate_ppl_rspace'
102
103 INTEGER :: atom_a, handle, iatom, ikind, j, lppl, &
104 n, natom_of_kind, ni, npme
105 INTEGER, DIMENSION(:), POINTER :: atom_list, cores
106 LOGICAL :: use_virial
107 REAL(kind=dp) :: alpha, eps_rho_rspace, radius
108 REAL(kind=dp), DIMENSION(3) :: force_a, force_b, ra
109 REAL(kind=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
110 REAL(kind=dp), DIMENSION(:), POINTER :: cexp_ppl
111 REAL(kind=dp), DIMENSION(:, :), POINTER :: hab, pab
112 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
113 TYPE(cell_type), POINTER :: cell
114 TYPE(dft_control_type), POINTER :: dft_control
115 TYPE(gth_potential_type), POINTER :: gth_potential
116 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
117 TYPE(pw_env_type), POINTER :: pw_env
118 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
119 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
120 TYPE(realspace_grid_type), POINTER :: rs_v
121 TYPE(virial_type), POINTER :: virial
122
123 CALL timeset(routinen, handle)
124
125 NULLIFY (pw_env, cores)
126
127 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
128 CALL pw_env_get(pw_env=pw_env, auxbas_rs_grid=rs_v)
129
130 CALL transfer_pw2rs(rs_v, rho_rspace)
131
132 CALL get_qs_env(qs_env=qs_env, &
133 atomic_kind_set=atomic_kind_set, &
134 qs_kind_set=qs_kind_set, &
135 cell=cell, &
136 dft_control=dft_control, &
137 particle_set=particle_set, &
138 pw_env=pw_env, &
139 force=force, virial=virial)
140
141 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
142
143 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
144
145 DO ikind = 1, SIZE(atomic_kind_set)
146
147 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_of_kind, atom_list=atom_list)
148 CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential)
149
150 IF (.NOT. ASSOCIATED(gth_potential)) cycle
151 CALL get_potential(potential=gth_potential, alpha_ppl=alpha, nexp_ppl=lppl, cexp_ppl=cexp_ppl)
152
153 IF (lppl <= 0) cycle
154
155 ni = ncoset(2*lppl - 2)
156 ALLOCATE (hab(ni, 1), pab(ni, 1))
157 pab = 0._dp
158
159 CALL reallocate(cores, 1, natom_of_kind)
160 npme = 0
161 cores = 0
162
163 ! prepare core function
164 DO j = 1, lppl
165 SELECT CASE (j)
166 CASE (1)
167 pab(1, 1) = cexp_ppl(1)
168 CASE (2)
169 n = coset(2, 0, 0)
170 pab(n, 1) = cexp_ppl(2)
171 n = coset(0, 2, 0)
172 pab(n, 1) = cexp_ppl(2)
173 n = coset(0, 0, 2)
174 pab(n, 1) = cexp_ppl(2)
175 CASE (3)
176 n = coset(4, 0, 0)
177 pab(n, 1) = cexp_ppl(3)
178 n = coset(0, 4, 0)
179 pab(n, 1) = cexp_ppl(3)
180 n = coset(0, 0, 4)
181 pab(n, 1) = cexp_ppl(3)
182 n = coset(2, 2, 0)
183 pab(n, 1) = 2._dp*cexp_ppl(3)
184 n = coset(2, 0, 2)
185 pab(n, 1) = 2._dp*cexp_ppl(3)
186 n = coset(0, 2, 2)
187 pab(n, 1) = 2._dp*cexp_ppl(3)
188 CASE (4)
189 n = coset(6, 0, 0)
190 pab(n, 1) = cexp_ppl(4)
191 n = coset(0, 6, 0)
192 pab(n, 1) = cexp_ppl(4)
193 n = coset(0, 0, 6)
194 pab(n, 1) = cexp_ppl(4)
195 n = coset(4, 2, 0)
196 pab(n, 1) = 3._dp*cexp_ppl(4)
197 n = coset(4, 0, 2)
198 pab(n, 1) = 3._dp*cexp_ppl(4)
199 n = coset(2, 4, 0)
200 pab(n, 1) = 3._dp*cexp_ppl(4)
201 n = coset(2, 0, 4)
202 pab(n, 1) = 3._dp*cexp_ppl(4)
203 n = coset(0, 4, 2)
204 pab(n, 1) = 3._dp*cexp_ppl(4)
205 n = coset(0, 2, 4)
206 pab(n, 1) = 3._dp*cexp_ppl(4)
207 n = coset(2, 2, 2)
208 pab(n, 1) = 6._dp*cexp_ppl(4)
209 CASE DEFAULT
210 cpabort("")
211 END SELECT
212 END DO
213
214 DO iatom = 1, natom_of_kind
215 atom_a = atom_list(iatom)
216 ra(:) = pbc(particle_set(atom_a)%r, cell)
217 IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN
218 ! replicated realspace grid, split the atoms up between procs
219 IF (modulo(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos) THEN
220 npme = npme + 1
221 cores(npme) = iatom
222 END IF
223 ELSE
224 npme = npme + 1
225 cores(npme) = iatom
226 END IF
227 END DO
228
229 DO j = 1, npme
230
231 iatom = cores(j)
232 atom_a = atom_list(iatom)
233 ra(:) = pbc(particle_set(atom_a)%r, cell)
234 hab(:, 1) = 0.0_dp
235 force_a(:) = 0.0_dp
236 force_b(:) = 0.0_dp
237 IF (use_virial) THEN
238 my_virial_a = 0.0_dp
239 my_virial_b = 0.0_dp
240 END IF
241 ni = 2*lppl - 2
242
243 radius = exp_radius_very_extended(la_min=0, la_max=ni, lb_min=0, lb_max=0, &
244 ra=ra, rb=ra, rp=ra, &
245 zetp=alpha, eps=eps_rho_rspace, &
246 pab=pab, o1=0, o2=0, & ! without map_consistent
247 prefactor=1.0_dp, cutoff=1.0_dp)
248
249 CALL integrate_pgf_product(ni, alpha, 0, &
250 0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), &
251 rs_v, hab, pab=pab, o1=0, o2=0, &
252 radius=radius, &
253 calculate_forces=.true., force_a=force_a, &
254 force_b=force_b, use_virial=use_virial, my_virial_a=my_virial_a, &
255 my_virial_b=my_virial_b, use_subpatch=.true., subpatch_pattern=0)
256
257 force(ikind)%gth_ppl(:, iatom) = &
258 force(ikind)%gth_ppl(:, iatom) + force_a(:)*rho_rspace%pw_grid%dvol
259
260 IF (use_virial) THEN
261 virial%pv_ppl = virial%pv_ppl + my_virial_a*rho_rspace%pw_grid%dvol
262 virial%pv_virial = virial%pv_virial + my_virial_a*rho_rspace%pw_grid%dvol
263 cpabort("Virial not debuged for CORE_PPL")
264 END IF
265 END DO
266
267 DEALLOCATE (hab, pab)
268
269 END DO
270
271 DEALLOCATE (cores)
272
273 CALL timestop(handle)
274
275 END SUBROUTINE integrate_ppl_rspace
276
277! **************************************************************************************************
278!> \brief computes the forces/virial due to the nlcc pseudopotential
279!> \param rho_rspace ...
280!> \param qs_env ...
281! **************************************************************************************************
282 SUBROUTINE integrate_rho_nlcc(rho_rspace, qs_env)
283 TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_rspace
284 TYPE(qs_environment_type), POINTER :: qs_env
285
286 CHARACTER(len=*), PARAMETER :: routinen = 'integrate_rho_nlcc'
287
288 INTEGER :: atom_a, handle, iatom, iexp_nlcc, ikind, &
289 ithread, j, n, natom, nc, nexp_nlcc, &
290 ni, npme, nthread
291 INTEGER, DIMENSION(:), POINTER :: atom_list, cores, nct_nlcc
292 LOGICAL :: nlcc, use_virial
293 REAL(kind=dp) :: alpha, eps_rho_rspace, radius
294 REAL(kind=dp), DIMENSION(3) :: force_a, force_b, ra
295 REAL(kind=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
296 REAL(kind=dp), DIMENSION(:), POINTER :: alpha_nlcc
297 REAL(kind=dp), DIMENSION(:, :), POINTER :: cval_nlcc, hab, pab
298 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
299 TYPE(cell_type), POINTER :: cell
300 TYPE(dft_control_type), POINTER :: dft_control
301 TYPE(gth_potential_type), POINTER :: gth_potential
302 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
303 TYPE(pw_env_type), POINTER :: pw_env
304 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
305 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
306 TYPE(realspace_grid_type), POINTER :: rs_v
307 TYPE(virial_type), POINTER :: virial
308
309 CALL timeset(routinen, handle)
310
311 NULLIFY (pw_env, cores)
312
313 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
314 CALL pw_env_get(pw_env=pw_env, auxbas_rs_grid=rs_v)
315
316 CALL transfer_pw2rs(rs_v, rho_rspace)
317
318 CALL get_qs_env(qs_env=qs_env, &
319 atomic_kind_set=atomic_kind_set, &
320 qs_kind_set=qs_kind_set, &
321 cell=cell, &
322 dft_control=dft_control, &
323 particle_set=particle_set, &
324 pw_env=pw_env, &
325 force=force, virial=virial)
326
327 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
328
329 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
330
331 DO ikind = 1, SIZE(atomic_kind_set)
332
333 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
334 CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential)
335
336 IF (.NOT. ASSOCIATED(gth_potential)) cycle
337 CALL get_potential(potential=gth_potential, nlcc_present=nlcc, nexp_nlcc=nexp_nlcc, &
338 alpha_nlcc=alpha_nlcc, nct_nlcc=nct_nlcc, cval_nlcc=cval_nlcc)
339
340 IF (.NOT. nlcc) cycle
341
342 DO iexp_nlcc = 1, nexp_nlcc
343
344 alpha = alpha_nlcc(iexp_nlcc)
345 nc = nct_nlcc(iexp_nlcc)
346
347 ni = ncoset(2*nc - 2)
348
349 nthread = 1
350 ithread = 0
351
352 ALLOCATE (hab(ni, 1), pab(ni, 1))
353 pab = 0._dp
354
355 CALL reallocate(cores, 1, natom)
356 npme = 0
357 cores = 0
358
359 ! prepare core function
360 DO j = 1, nc
361 SELECT CASE (j)
362 CASE (1)
363 pab(1, 1) = cval_nlcc(1, iexp_nlcc)
364 CASE (2)
365 n = coset(2, 0, 0)
366 pab(n, 1) = cval_nlcc(2, iexp_nlcc)/alpha**2
367 n = coset(0, 2, 0)
368 pab(n, 1) = cval_nlcc(2, iexp_nlcc)/alpha**2
369 n = coset(0, 0, 2)
370 pab(n, 1) = cval_nlcc(2, iexp_nlcc)/alpha**2
371 CASE (3)
372 n = coset(4, 0, 0)
373 pab(n, 1) = cval_nlcc(3, iexp_nlcc)/alpha**4
374 n = coset(0, 4, 0)
375 pab(n, 1) = cval_nlcc(3, iexp_nlcc)/alpha**4
376 n = coset(0, 0, 4)
377 pab(n, 1) = cval_nlcc(3, iexp_nlcc)/alpha**4
378 n = coset(2, 2, 0)
379 pab(n, 1) = 2._dp*cval_nlcc(3, iexp_nlcc)/alpha**4
380 n = coset(2, 0, 2)
381 pab(n, 1) = 2._dp*cval_nlcc(3, iexp_nlcc)/alpha**4
382 n = coset(0, 2, 2)
383 pab(n, 1) = 2._dp*cval_nlcc(3, iexp_nlcc)/alpha**4
384 CASE (4)
385 n = coset(6, 0, 0)
386 pab(n, 1) = cval_nlcc(4, iexp_nlcc)/alpha**6
387 n = coset(0, 6, 0)
388 pab(n, 1) = cval_nlcc(4, iexp_nlcc)/alpha**6
389 n = coset(0, 0, 6)
390 pab(n, 1) = cval_nlcc(4, iexp_nlcc)/alpha**6
391 n = coset(4, 2, 0)
392 pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
393 n = coset(4, 0, 2)
394 pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
395 n = coset(2, 4, 0)
396 pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
397 n = coset(2, 0, 4)
398 pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
399 n = coset(0, 4, 2)
400 pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
401 n = coset(0, 2, 4)
402 pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
403 n = coset(2, 2, 2)
404 pab(n, 1) = 6._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
405 CASE DEFAULT
406 cpabort("")
407 END SELECT
408 END DO
409 IF (dft_control%nspins == 2) pab = pab*0.5_dp
410
411 DO iatom = 1, natom
412 atom_a = atom_list(iatom)
413 ra(:) = pbc(particle_set(atom_a)%r, cell)
414 IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN
415 ! replicated realspace grid, split the atoms up between procs
416 IF (modulo(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos) THEN
417 npme = npme + 1
418 cores(npme) = iatom
419 END IF
420 ELSE
421 npme = npme + 1
422 cores(npme) = iatom
423 END IF
424 END DO
425
426 DO j = 1, npme
427
428 iatom = cores(j)
429 atom_a = atom_list(iatom)
430 ra(:) = pbc(particle_set(atom_a)%r, cell)
431 hab(:, 1) = 0.0_dp
432 force_a(:) = 0.0_dp
433 force_b(:) = 0.0_dp
434 IF (use_virial) THEN
435 my_virial_a = 0.0_dp
436 my_virial_b = 0.0_dp
437 END IF
438 ni = 2*nc - 2
439
440 radius = exp_radius_very_extended(la_min=0, la_max=ni, lb_min=0, lb_max=0, &
441 ra=ra, rb=ra, rp=ra, &
442 zetp=1/(2*alpha**2), eps=eps_rho_rspace, &
443 pab=pab, o1=0, o2=0, & ! without map_consistent
444 prefactor=1.0_dp, cutoff=1.0_dp)
445
446 CALL integrate_pgf_product(ni, 1/(2*alpha**2), 0, &
447 0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), &
448 rs_v, hab, pab=pab, o1=0, o2=0, &
449 radius=radius, &
450 calculate_forces=.true., force_a=force_a, &
451 force_b=force_b, use_virial=use_virial, my_virial_a=my_virial_a, &
452 my_virial_b=my_virial_b, use_subpatch=.true., subpatch_pattern=0)
453
454 force(ikind)%gth_nlcc(:, iatom) = &
455 force(ikind)%gth_nlcc(:, iatom) + force_a(:)*rho_rspace%pw_grid%dvol
456
457 IF (use_virial) THEN
458 virial%pv_nlcc = virial%pv_nlcc + my_virial_a*rho_rspace%pw_grid%dvol
459 virial%pv_virial = virial%pv_virial + my_virial_a*rho_rspace%pw_grid%dvol
460 END IF
461 END DO
462
463 DEALLOCATE (hab, pab)
464
465 END DO
466
467 END DO
468
469 DEALLOCATE (cores)
470
471 CALL timestop(handle)
472
473 END SUBROUTINE integrate_rho_nlcc
474
475! **************************************************************************************************
476!> \brief computes the forces/virial due to the ionic cores with a potential on
477!> grid
478!> \param v_rspace ...
479!> \param qs_env ...
480!> \param atecc ...
481! **************************************************************************************************
482 SUBROUTINE integrate_v_core_rspace(v_rspace, qs_env, atecc)
483 TYPE(pw_r3d_rs_type), INTENT(IN) :: v_rspace
484 TYPE(qs_environment_type), POINTER :: qs_env
485 REAL(kind=dp), DIMENSION(:), OPTIONAL :: atecc
486
487 CHARACTER(len=*), PARAMETER :: routinen = 'integrate_v_core_rspace'
488
489 INTEGER :: atom_a, handle, iatom, ikind, j, natom, &
490 natom_of_kind, npme
491 INTEGER, DIMENSION(:), POINTER :: atom_list, cores
492 LOGICAL :: paw_atom, skip_fcore, use_virial
493 REAL(kind=dp) :: alpha_core_charge, ccore_charge, &
494 eps_rho_rspace, radius
495 REAL(kind=dp), DIMENSION(3) :: force_a, force_b, ra
496 REAL(kind=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
497 REAL(kind=dp), DIMENSION(:, :), POINTER :: hab, pab
498 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
499 TYPE(atprop_type), POINTER :: atprop
500 TYPE(cell_type), POINTER :: cell
501 TYPE(dft_control_type), POINTER :: dft_control
502 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
503 TYPE(pw_env_type), POINTER :: pw_env
504 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
505 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
506 TYPE(realspace_grid_type), POINTER :: rs_v
507 TYPE(virial_type), POINTER :: virial
508
509 CALL timeset(routinen, handle)
510 NULLIFY (virial, force, atprop, dft_control)
511
512 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
513
514 !If gapw, check for gpw kinds
515 skip_fcore = .false.
516 IF (dft_control%qs_control%gapw) THEN
517 IF (.NOT. dft_control%qs_control%gapw_control%nopaw_as_gpw) skip_fcore = .true.
518 END IF
519
520 IF (.NOT. skip_fcore) THEN
521 NULLIFY (pw_env)
522 ALLOCATE (cores(1))
523 ALLOCATE (hab(1, 1))
524 ALLOCATE (pab(1, 1))
525
526 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
527 CALL pw_env_get(pw_env=pw_env, auxbas_rs_grid=rs_v)
528
529 CALL transfer_pw2rs(rs_v, v_rspace)
530
531 CALL get_qs_env(qs_env=qs_env, &
532 atomic_kind_set=atomic_kind_set, &
533 qs_kind_set=qs_kind_set, &
534 cell=cell, &
535 dft_control=dft_control, &
536 particle_set=particle_set, &
537 pw_env=pw_env, &
538 force=force, &
539 virial=virial, &
540 atprop=atprop)
541
542 ! atomic energy contributions
543 natom = SIZE(particle_set)
544 IF (ASSOCIATED(atprop)) THEN
545 CALL atprop_array_init(atprop%ateb, natom)
546 END IF
547
548 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
549
550 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
551
552 DO ikind = 1, SIZE(atomic_kind_set)
553
554 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_of_kind, atom_list=atom_list)
555 CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom, &
556 alpha_core_charge=alpha_core_charge, &
557 ccore_charge=ccore_charge)
558
559 IF (dft_control%qs_control%gapw .AND. paw_atom) cycle
560
561 pab(1, 1) = -ccore_charge
562 IF (alpha_core_charge == 0.0_dp .OR. pab(1, 1) == 0.0_dp) cycle
563
564 CALL reallocate(cores, 1, natom_of_kind)
565 npme = 0
566 cores = 0
567
568 DO iatom = 1, natom_of_kind
569 atom_a = atom_list(iatom)
570 ra(:) = pbc(particle_set(atom_a)%r, cell)
571 IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN
572 ! replicated realspace grid, split the atoms up between procs
573 IF (modulo(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos) THEN
574 npme = npme + 1
575 cores(npme) = iatom
576 END IF
577 ELSE
578 npme = npme + 1
579 cores(npme) = iatom
580 END IF
581 END DO
582
583 DO j = 1, npme
584
585 iatom = cores(j)
586 atom_a = atom_list(iatom)
587 ra(:) = pbc(particle_set(atom_a)%r, cell)
588 hab(1, 1) = 0.0_dp
589 force_a(:) = 0.0_dp
590 force_b(:) = 0.0_dp
591 IF (use_virial) THEN
592 my_virial_a = 0.0_dp
593 my_virial_b = 0.0_dp
594 END IF
595
596 radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
597 ra=ra, rb=ra, rp=ra, &
598 zetp=alpha_core_charge, eps=eps_rho_rspace, &
599 pab=pab, o1=0, o2=0, & ! without map_consistent
600 prefactor=1.0_dp, cutoff=1.0_dp)
601
602 CALL integrate_pgf_product(0, alpha_core_charge, 0, &
603 0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), &
604 rs_v, hab, pab=pab, o1=0, o2=0, &
605 radius=radius, &
606 calculate_forces=.true., force_a=force_a, &
607 force_b=force_b, use_virial=use_virial, my_virial_a=my_virial_a, &
608 my_virial_b=my_virial_b, use_subpatch=.true., subpatch_pattern=0)
609
610 IF (ASSOCIATED(force)) THEN
611 force(ikind)%rho_core(:, iatom) = force(ikind)%rho_core(:, iatom) + force_a(:)
612 END IF
613 IF (use_virial) THEN
614 virial%pv_ehartree = virial%pv_ehartree + my_virial_a
615 virial%pv_virial = virial%pv_virial + my_virial_a
616 END IF
617 IF (ASSOCIATED(atprop)) THEN
618 atprop%ateb(atom_a) = atprop%ateb(atom_a) + 0.5_dp*hab(1, 1)*pab(1, 1)
619 END IF
620 IF (PRESENT(atecc)) THEN
621 atecc(atom_a) = atecc(atom_a) + 0.5_dp*hab(1, 1)*pab(1, 1)
622 END IF
623
624 END DO
625
626 END DO
627
628 DEALLOCATE (hab, pab, cores)
629
630 END IF
631
632 CALL timestop(handle)
633
634 END SUBROUTINE integrate_v_core_rspace
635
636! **************************************************************************************************
637!> \brief computes the overlap of a set of Gaussians with a potential on grid
638!> \param v_rspace ...
639!> \param qs_env ...
640!> \param alpha ...
641!> \param ccore ...
642!> \param atecc ...
643! **************************************************************************************************
644 SUBROUTINE integrate_v_gaussian_rspace(v_rspace, qs_env, alpha, ccore, atecc)
645 TYPE(pw_r3d_rs_type), INTENT(IN) :: v_rspace
646 TYPE(qs_environment_type), POINTER :: qs_env
647 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: alpha, ccore
648 REAL(kind=dp), DIMENSION(:) :: atecc
649
650 CHARACTER(len=*), PARAMETER :: routinen = 'integrate_v_gaussian_rspace'
651
652 INTEGER :: atom_a, handle, iatom, ikind, j, natom, &
653 natom_of_kind, npme
654 INTEGER, DIMENSION(:), POINTER :: atom_list, cores
655 REAL(kind=dp) :: alpha_core_charge, eps_rho_rspace, radius
656 REAL(kind=dp), DIMENSION(3) :: ra
657 REAL(kind=dp), DIMENSION(:, :), POINTER :: hab, pab
658 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
659 TYPE(cell_type), POINTER :: cell
660 TYPE(dft_control_type), POINTER :: dft_control
661 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
662 TYPE(pw_env_type), POINTER :: pw_env
663 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
664 TYPE(realspace_grid_type), POINTER :: rs_v
665
666 CALL timeset(routinen, handle)
667
668 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
669
670 !If gapw, check for gpw kinds
671 cpassert(.NOT. dft_control%qs_control%gapw)
672
673 NULLIFY (pw_env)
674 ALLOCATE (cores(1))
675 ALLOCATE (hab(1, 1))
676 ALLOCATE (pab(1, 1))
677
678 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
679 CALL pw_env_get(pw_env=pw_env, auxbas_rs_grid=rs_v)
680
681 CALL transfer_pw2rs(rs_v, v_rspace)
682
683 CALL get_qs_env(qs_env=qs_env, &
684 atomic_kind_set=atomic_kind_set, &
685 qs_kind_set=qs_kind_set, &
686 cell=cell, &
687 dft_control=dft_control, &
688 particle_set=particle_set, &
689 pw_env=pw_env)
690
691 ! atomic energy contributions
692 natom = SIZE(particle_set)
693 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
694
695 DO ikind = 1, SIZE(atomic_kind_set)
696
697 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_of_kind, atom_list=atom_list)
698 pab(1, 1) = -ccore(ikind)
699 alpha_core_charge = alpha(ikind)
700 IF (alpha_core_charge == 0.0_dp .OR. pab(1, 1) == 0.0_dp) cycle
701
702 CALL reallocate(cores, 1, natom_of_kind)
703 npme = 0
704 cores = 0
705
706 DO iatom = 1, natom_of_kind
707 atom_a = atom_list(iatom)
708 ra(:) = pbc(particle_set(atom_a)%r, cell)
709 IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN
710 ! replicated realspace grid, split the atoms up between procs
711 IF (modulo(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos) THEN
712 npme = npme + 1
713 cores(npme) = iatom
714 END IF
715 ELSE
716 npme = npme + 1
717 cores(npme) = iatom
718 END IF
719 END DO
720
721 DO j = 1, npme
722
723 iatom = cores(j)
724 atom_a = atom_list(iatom)
725 ra(:) = pbc(particle_set(atom_a)%r, cell)
726 hab(1, 1) = 0.0_dp
727
728 radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
729 ra=ra, rb=ra, rp=ra, &
730 zetp=alpha_core_charge, eps=eps_rho_rspace, &
731 pab=pab, o1=0, o2=0, & ! without map_consistent
732 prefactor=1.0_dp, cutoff=1.0_dp)
733
734 CALL integrate_pgf_product(0, alpha_core_charge, 0, &
735 0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), &
736 rs_v, hab, pab=pab, o1=0, o2=0, &
737 radius=radius, calculate_forces=.false., &
738 use_subpatch=.true., subpatch_pattern=0)
739 atecc(atom_a) = atecc(atom_a) + 0.5_dp*hab(1, 1)*pab(1, 1)
740
741 END DO
742
743 END DO
744
745 DEALLOCATE (hab, pab, cores)
746
747 CALL timestop(handle)
748
749 END SUBROUTINE integrate_v_gaussian_rspace
750! **************************************************************************************************
751!> \brief computes integrals of product of v_rspace times a one-center function
752!> required for LRIGPW
753!> \param v_rspace ...
754!> \param qs_env ...
755!> \param int_res ...
756!> \param calculate_forces ...
757!> \param basis_type ...
758!> \param atomlist ...
759!> \author Dorothea Golze
760! **************************************************************************************************
761 SUBROUTINE integrate_v_rspace_one_center(v_rspace, qs_env, int_res, &
762 calculate_forces, basis_type, atomlist)
763 TYPE(pw_r3d_rs_type), INTENT(IN) :: v_rspace
764 TYPE(qs_environment_type), POINTER :: qs_env
765 TYPE(lri_kind_type), DIMENSION(:), POINTER :: int_res
766 LOGICAL, INTENT(IN) :: calculate_forces
767 CHARACTER(len=*), INTENT(IN) :: basis_type
768 INTEGER, DIMENSION(:), OPTIONAL :: atomlist
769
770 CHARACTER(len=*), PARAMETER :: routinen = 'integrate_v_rspace_one_center'
771
772 INTEGER :: atom_a, group_size, handle, i, iatom, igrid_level, ikind, ipgf, iset, m1, maxco, &
773 maxsgf_set, my_pos, na1, natom_of_kind, ncoa, nkind, nseta, offset, sgfa
774 INTEGER, DIMENSION(:), POINTER :: atom_list, la_max, la_min, npgfa, &
775 nsgf_seta
776 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
777 LOGICAL :: use_virial
778 LOGICAL, ALLOCATABLE, DIMENSION(:) :: map_it
779 REAL(kind=dp) :: eps_rho_rspace, radius
780 REAL(kind=dp), DIMENSION(3) :: force_a, force_b, ra
781 REAL(kind=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
782 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a
783 REAL(kind=dp), DIMENSION(:, :), POINTER :: hab, pab, rpgfa, sphi_a, work_f, work_i, &
784 zeta
785 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
786 TYPE(cell_type), POINTER :: cell
787 TYPE(dft_control_type), POINTER :: dft_control
788 TYPE(gridlevel_info_type), POINTER :: gridlevel_info
789 TYPE(gto_basis_set_type), POINTER :: lri_basis_set
790 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
791 TYPE(pw_env_type), POINTER :: pw_env
792 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
793 TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_v
794 TYPE(realspace_grid_type), POINTER :: rs_grid
795 TYPE(virial_type), POINTER :: virial
796
797 CALL timeset(routinen, handle)
798
799 NULLIFY (atomic_kind_set, qs_kind_set, atom_list, cell, dft_control, &
800 first_sgfa, gridlevel_info, hab, la_max, la_min, lri_basis_set, &
801 npgfa, nsgf_seta, pab, particle_set, pw_env, rpgfa, &
802 rs_grid, rs_v, virial, set_radius_a, sphi_a, work_f, &
803 work_i, zeta)
804
805 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
806
807 CALL pw_env_get(pw_env, rs_grids=rs_v)
808 DO i = 1, SIZE(rs_v)
809 CALL rs_grid_zero(rs_v(i))
810 END DO
811
812 gridlevel_info => pw_env%gridlevel_info
813
814 CALL potential_pw2rs(rs_v, v_rspace, pw_env)
815
816 CALL get_qs_env(qs_env=qs_env, &
817 atomic_kind_set=atomic_kind_set, &
818 qs_kind_set=qs_kind_set, &
819 cell=cell, &
820 dft_control=dft_control, &
821 nkind=nkind, &
822 particle_set=particle_set, &
823 pw_env=pw_env, &
824 virial=virial)
825
826 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
827
828 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
829
830 offset = 0
831 my_pos = v_rspace%pw_grid%para%my_pos
832 group_size = v_rspace%pw_grid%para%group_size
833
834 DO ikind = 1, nkind
835
836 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_of_kind, atom_list=atom_list)
837 CALL get_qs_kind(qs_kind_set(ikind), basis_set=lri_basis_set, basis_type=basis_type)
838 CALL get_gto_basis_set(gto_basis_set=lri_basis_set, &
839 first_sgf=first_sgfa, &
840 lmax=la_max, &
841 lmin=la_min, &
842 maxco=maxco, &
843 maxsgf_set=maxsgf_set, &
844 npgf=npgfa, &
845 nset=nseta, &
846 nsgf_set=nsgf_seta, &
847 pgf_radius=rpgfa, &
848 set_radius=set_radius_a, &
849 sphi=sphi_a, &
850 zet=zeta)
851
852 ALLOCATE (hab(maxco, 1), pab(maxco, 1))
853 hab = 0._dp
854 pab(:, 1) = 0._dp
855
856 DO iatom = 1, natom_of_kind
857
858 atom_a = atom_list(iatom)
859 IF (PRESENT(atomlist)) THEN
860 IF (atomlist(atom_a) == 0) cycle
861 END IF
862 ra(:) = pbc(particle_set(atom_a)%r, cell)
863 force_a(:) = 0._dp
864 force_b(:) = 0._dp
865 my_virial_a(:, :) = 0._dp
866 my_virial_b(:, :) = 0._dp
867
868 m1 = maxval(npgfa(1:nseta))
869 ALLOCATE (map_it(m1))
870
871 DO iset = 1, nseta
872 !
873 map_it = .false.
874 DO ipgf = 1, npgfa(iset)
875 igrid_level = gaussian_gridlevel(gridlevel_info, zeta(ipgf, iset))
876 rs_grid => rs_v(igrid_level)
877 map_it(ipgf) = map_gaussian_here(rs_grid, cell%h_inv, ra, offset, group_size, my_pos)
878 END DO
879 offset = offset + 1
880 !
881 IF (any(map_it(1:npgfa(iset)))) THEN
882 sgfa = first_sgfa(1, iset)
883 ncoa = npgfa(iset)*ncoset(la_max(iset))
884 hab(:, 1) = 0._dp
885 ALLOCATE (work_i(nsgf_seta(iset), 1))
886 work_i = 0.0_dp
887
888 ! get fit coefficients for forces
889 IF (calculate_forces) THEN
890 m1 = sgfa + nsgf_seta(iset) - 1
891 ALLOCATE (work_f(nsgf_seta(iset), 1))
892 work_f(1:nsgf_seta(iset), 1) = int_res(ikind)%acoef(iatom, sgfa:m1)
893 CALL dgemm("N", "N", ncoa, 1, nsgf_seta(iset), 1.0_dp, sphi_a(1, sgfa), &
894 SIZE(sphi_a, 1), work_f(1, 1), SIZE(work_f, 1), 0.0_dp, pab(1, 1), &
895 SIZE(pab, 1))
896 DEALLOCATE (work_f)
897 END IF
898
899 DO ipgf = 1, npgfa(iset)
900 na1 = (ipgf - 1)*ncoset(la_max(iset))
901 igrid_level = gaussian_gridlevel(gridlevel_info, zeta(ipgf, iset))
902 rs_grid => rs_v(igrid_level)
903
904 radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
905 lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, &
906 zetp=zeta(ipgf, iset), eps=eps_rho_rspace, &
907 prefactor=1.0_dp, cutoff=1.0_dp)
908
909 IF (map_it(ipgf)) THEN
910 IF (.NOT. calculate_forces) THEN
911 CALL integrate_pgf_product(la_max=la_max(iset), &
912 zeta=zeta(ipgf, iset), la_min=la_min(iset), &
913 lb_max=0, zetb=0.0_dp, lb_min=0, &
914 ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
915 rsgrid=rs_grid, &
916 hab=hab, o1=na1, o2=0, radius=radius, &
917 calculate_forces=calculate_forces)
918 ELSE
919 CALL integrate_pgf_product(la_max=la_max(iset), &
920 zeta=zeta(ipgf, iset), la_min=la_min(iset), &
921 lb_max=0, zetb=0.0_dp, lb_min=0, &
922 ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
923 rsgrid=rs_grid, &
924 hab=hab, pab=pab, o1=na1, o2=0, radius=radius, &
925 calculate_forces=calculate_forces, &
926 force_a=force_a, force_b=force_b, &
927 use_virial=use_virial, &
928 my_virial_a=my_virial_a, my_virial_b=my_virial_b)
929 END IF
930 END IF
931 END DO
932 ! contract hab
933 CALL dgemm("T", "N", nsgf_seta(iset), 1, ncoa, 1.0_dp, sphi_a(1, sgfa), &
934 SIZE(sphi_a, 1), hab(1, 1), SIZE(hab, 1), 0.0_dp, work_i(1, 1), SIZE(work_i, 1))
935
936 int_res(ikind)%v_int(iatom, sgfa:sgfa - 1 + nsgf_seta(iset)) = &
937 int_res(ikind)%v_int(iatom, sgfa:sgfa - 1 + nsgf_seta(iset)) + work_i(1:nsgf_seta(iset), 1)
938 DEALLOCATE (work_i)
939 END IF
940 END DO
941 !
942 IF (calculate_forces) THEN
943 int_res(ikind)%v_dfdr(iatom, :) = int_res(ikind)%v_dfdr(iatom, :) + force_a(:)
944 IF (use_virial) THEN
945 virial%pv_lrigpw = virial%pv_lrigpw + my_virial_a
946 virial%pv_virial = virial%pv_virial + my_virial_a
947 END IF
948 END IF
949
950 DEALLOCATE (map_it)
951
952 END DO
953
954 DEALLOCATE (hab, pab)
955 END DO
956
957 CALL timestop(handle)
958
959 END SUBROUTINE integrate_v_rspace_one_center
960
961! **************************************************************************************************
962!> \brief computes integrals of product of v_rspace times the diagonal block basis functions
963!> required for LRIGPW with exact 1c terms
964!> \param v_rspace ...
965!> \param ksmat ...
966!> \param pmat ...
967!> \param qs_env ...
968!> \param calculate_forces ...
969!> \param basis_type ...
970!> \author JGH
971! **************************************************************************************************
972 SUBROUTINE integrate_v_rspace_diagonal(v_rspace, ksmat, pmat, qs_env, calculate_forces, basis_type)
973 TYPE(pw_r3d_rs_type), INTENT(IN) :: v_rspace
974 TYPE(dbcsr_type), INTENT(INOUT) :: ksmat, pmat
975 TYPE(qs_environment_type), POINTER :: qs_env
976 LOGICAL, INTENT(IN) :: calculate_forces
977 CHARACTER(len=*), INTENT(IN) :: basis_type
978
979 CHARACTER(len=*), PARAMETER :: routinen = 'integrate_v_rspace_diagonal'
980
981 INTEGER :: atom_a, group_size, handle, iatom, igrid_level, ikind, ipgf, iset, jpgf, jset, &
982 m1, maxco, maxsgf_set, my_pos, na1, na2, natom_of_kind, nb1, nb2, ncoa, ncob, nkind, &
983 nseta, nsgfa, offset, sgfa, sgfb
984 INTEGER, DIMENSION(:), POINTER :: atom_list, la_max, la_min, npgfa, &
985 nsgf_seta
986 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
987 LOGICAL :: found, use_virial
988 LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: map_it2
989 REAL(kind=dp) :: eps_rho_rspace, radius, zetp
990 REAL(kind=dp), DIMENSION(3) :: force_a, force_b, ra
991 REAL(kind=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
992 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a
993 REAL(kind=dp), DIMENSION(:, :), POINTER :: h_block, hab, hmat, p_block, pab, pblk, &
994 rpgfa, sphi_a, work, zeta
995 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
996 TYPE(cell_type), POINTER :: cell
997 TYPE(dft_control_type), POINTER :: dft_control
998 TYPE(gridlevel_info_type), POINTER :: gridlevel_info
999 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1000 TYPE(mp_para_env_type), POINTER :: para_env
1001 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1002 TYPE(pw_env_type), POINTER :: pw_env
1003 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1004 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1005 TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_v
1006 TYPE(realspace_grid_type), POINTER :: rs_grid
1007 TYPE(virial_type), POINTER :: virial
1008
1009 CALL timeset(routinen, handle)
1010
1011 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
1012 CALL pw_env_get(pw_env, rs_grids=rs_v)
1013 CALL potential_pw2rs(rs_v, v_rspace, pw_env)
1014
1015 gridlevel_info => pw_env%gridlevel_info
1016
1017 CALL get_qs_env(qs_env=qs_env, &
1018 atomic_kind_set=atomic_kind_set, &
1019 qs_kind_set=qs_kind_set, &
1020 cell=cell, &
1021 dft_control=dft_control, &
1022 nkind=nkind, &
1023 particle_set=particle_set, &
1024 force=force, &
1025 virial=virial, &
1026 para_env=para_env)
1027
1028 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
1029 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1030
1031 offset = 0
1032 my_pos = v_rspace%pw_grid%para%my_pos
1033 group_size = v_rspace%pw_grid%para%group_size
1034
1035 DO ikind = 1, nkind
1036
1037 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_of_kind, atom_list=atom_list)
1038 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type=basis_type)
1039 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
1040 lmax=la_max, lmin=la_min, maxco=maxco, maxsgf_set=maxsgf_set, &
1041 npgf=npgfa, nset=nseta, nsgf_set=nsgf_seta, nsgf=nsgfa, &
1042 first_sgf=first_sgfa, pgf_radius=rpgfa, set_radius=set_radius_a, &
1043 sphi=sphi_a, zet=zeta)
1044
1045 ALLOCATE (hab(maxco, maxco), work(maxco, maxsgf_set), hmat(nsgfa, nsgfa))
1046 IF (calculate_forces) ALLOCATE (pab(maxco, maxco), pblk(nsgfa, nsgfa))
1047
1048 DO iatom = 1, natom_of_kind
1049 atom_a = atom_list(iatom)
1050 ra(:) = pbc(particle_set(atom_a)%r, cell)
1051 hmat = 0.0_dp
1052 IF (calculate_forces) THEN
1053 CALL dbcsr_get_block_p(matrix=pmat, row=atom_a, col=atom_a, block=p_block, found=found)
1054 IF (found) THEN
1055 pblk(1:nsgfa, 1:nsgfa) = p_block(1:nsgfa, 1:nsgfa)
1056 ELSE
1057 pblk = 0.0_dp
1058 END IF
1059 CALL para_env%sum(pblk)
1060 force_a(:) = 0._dp
1061 force_b(:) = 0._dp
1062 IF (use_virial) THEN
1063 my_virial_a = 0.0_dp
1064 my_virial_b = 0.0_dp
1065 END IF
1066 END IF
1067 m1 = maxval(npgfa(1:nseta))
1068 ALLOCATE (map_it2(m1, m1))
1069 DO iset = 1, nseta
1070 sgfa = first_sgfa(1, iset)
1071 ncoa = npgfa(iset)*ncoset(la_max(iset))
1072 DO jset = 1, nseta
1073 sgfb = first_sgfa(1, jset)
1074 ncob = npgfa(jset)*ncoset(la_max(jset))
1075 !
1076 map_it2 = .false.
1077 DO ipgf = 1, npgfa(iset)
1078 DO jpgf = 1, npgfa(jset)
1079 zetp = zeta(ipgf, iset) + zeta(jpgf, jset)
1080 igrid_level = gaussian_gridlevel(gridlevel_info, zetp)
1081 rs_grid => rs_v(igrid_level)
1082 map_it2(ipgf, jpgf) = map_gaussian_here(rs_grid, cell%h_inv, ra, offset, group_size, my_pos)
1083 END DO
1084 END DO
1085 offset = offset + 1
1086 !
1087 IF (any(map_it2(1:npgfa(iset), 1:npgfa(jset)))) THEN
1088 hab(:, :) = 0._dp
1089 IF (calculate_forces) THEN
1090 CALL dgemm("N", "N", ncoa, nsgf_seta(jset), nsgf_seta(iset), &
1091 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1092 pblk(sgfa, sgfb), SIZE(pblk, 1), &
1093 0.0_dp, work(1, 1), SIZE(work, 1))
1094 CALL dgemm("N", "T", ncoa, ncob, nsgf_seta(jset), &
1095 1.0_dp, work(1, 1), SIZE(work, 1), &
1096 sphi_a(1, sgfb), SIZE(sphi_a, 1), &
1097 0.0_dp, pab(1, 1), SIZE(pab, 1))
1098 END IF
1099
1100 DO ipgf = 1, npgfa(iset)
1101 na1 = (ipgf - 1)*ncoset(la_max(iset))
1102 na2 = ipgf*ncoset(la_max(iset))
1103 DO jpgf = 1, npgfa(jset)
1104 nb1 = (jpgf - 1)*ncoset(la_max(jset))
1105 nb2 = jpgf*ncoset(la_max(jset))
1106 zetp = zeta(ipgf, iset) + zeta(jpgf, jset)
1107 igrid_level = gaussian_gridlevel(gridlevel_info, zetp)
1108 rs_grid => rs_v(igrid_level)
1109
1110 radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
1111 lb_min=la_min(jset), lb_max=la_max(jset), &
1112 ra=ra, rb=ra, rp=ra, &
1113 zetp=zetp, eps=eps_rho_rspace, &
1114 prefactor=1.0_dp, cutoff=1.0_dp)
1115
1116 IF (map_it2(ipgf, jpgf)) THEN
1117 IF (calculate_forces) THEN
1118 CALL integrate_pgf_product( &
1119 la_max(iset), zeta(ipgf, iset), la_min(iset), &
1120 la_max(jset), zeta(jpgf, jset), la_min(jset), &
1121 ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
1122 rsgrid=rs_v(igrid_level), &
1123 hab=hab, pab=pab, o1=na1, o2=nb1, &
1124 radius=radius, &
1125 calculate_forces=.true., &
1126 force_a=force_a, force_b=force_b, &
1127 use_virial=use_virial, my_virial_a=my_virial_a, my_virial_b=my_virial_b)
1128 ELSE
1129 CALL integrate_pgf_product( &
1130 la_max(iset), zeta(ipgf, iset), la_min(iset), &
1131 la_max(jset), zeta(jpgf, jset), la_min(jset), &
1132 ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
1133 rsgrid=rs_v(igrid_level), &
1134 hab=hab, o1=na1, o2=nb1, &
1135 radius=radius, &
1136 calculate_forces=.false.)
1137 END IF
1138 END IF
1139 END DO
1140 END DO
1141 ! contract hab
1142 CALL dgemm("N", "N", ncoa, nsgf_seta(jset), ncob, &
1143 1.0_dp, hab(1, 1), SIZE(hab, 1), &
1144 sphi_a(1, sgfb), SIZE(sphi_a, 1), &
1145 0.0_dp, work(1, 1), SIZE(work, 1))
1146 CALL dgemm("T", "N", nsgf_seta(iset), nsgf_seta(jset), ncoa, &
1147 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1148 work(1, 1), SIZE(work, 1), &
1149 1.0_dp, hmat(sgfa, sgfb), SIZE(hmat, 1))
1150 END IF
1151 END DO
1152 END DO
1153 DEALLOCATE (map_it2)
1154 ! update KS matrix
1155 CALL para_env%sum(hmat)
1156 CALL dbcsr_get_block_p(matrix=ksmat, row=atom_a, col=atom_a, block=h_block, found=found)
1157 IF (found) THEN
1158 h_block(1:nsgfa, 1:nsgfa) = h_block(1:nsgfa, 1:nsgfa) + hmat(1:nsgfa, 1:nsgfa)
1159 END IF
1160 IF (calculate_forces) THEN
1161 force(ikind)%rho_elec(:, iatom) = force(ikind)%rho_elec(:, iatom) + 2.0_dp*force_a(:)
1162 IF (use_virial) THEN
1163 IF (use_virial .AND. calculate_forces) THEN
1164 virial%pv_lrigpw = virial%pv_lrigpw + 2.0_dp*my_virial_a
1165 virial%pv_virial = virial%pv_virial + 2.0_dp*my_virial_a
1166 END IF
1167 END IF
1168 END IF
1169 END DO
1170 DEALLOCATE (hab, work, hmat)
1171 IF (calculate_forces) DEALLOCATE (pab, pblk)
1172 END DO
1173
1174 CALL timestop(handle)
1175
1176 END SUBROUTINE integrate_v_rspace_diagonal
1177
1178END MODULE qs_integrate_potential_single
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:117
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:195
cell_types::pbc
Definition cell_types.F:92
external_potential_types::get_potential
Definition external_potential_types.F:276
memory_utilities::reallocate
Definition memory_utilities.F:27
ao_util
All kind of helpful little routines.
Definition ao_util.F:14
ao_util::exp_radius_very_extended
real(kind=dp) function, public exp_radius_very_extended(la_min, la_max, lb_min, lb_max, pab, o1, o2, ra, rb, rp, zetp, eps, prefactor, cutoff, epsabs)
computes the radius of the Gaussian outside of which it is smaller than eps
Definition ao_util.F:209
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
atprop_types
Holds information on atomic properties.
Definition atprop_types.F:14
atprop_types::atprop_array_init
subroutine, public atprop_array_init(atarray, natom)
...
Definition atprop_types.F:98
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)
...
Definition basis_set_types.F:615
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
external_potential_types
Definition of the atomic potential types.
Definition external_potential_types.F:14
gaussian_gridlevels
Definition gaussian_gridlevels.F:13
gaussian_gridlevels::gaussian_gridlevel
integer function, public gaussian_gridlevel(gridlevel_info, exponent)
...
Definition gaussian_gridlevels.F:148
grid_api
Fortran API for the grid package, which is written in C.
Definition grid_api.F:12
grid_api::integrate_pgf_product
subroutine, public integrate_pgf_product(la_max, zeta, la_min, lb_max, zetb, lb_min, ra, rab, rsgrid, hab, pab, o1, o2, radius, calculate_forces, force_a, force_b, compute_tau, use_virial, my_virial_a, my_virial_b, hdab, hadb, a_hdab, use_subpatch, subpatch_pattern)
low level function to compute matrix elements of primitive gaussian functions
Definition grid_api.F:279
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
lri_environment_types
contains the types and subroutines for dealing with the lri_env lri : local resolution of the identit...
Definition lri_environment_types.F:18
memory_utilities
Utility routines for the memory handling.
Definition memory_utilities.F:14
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
orbital_pointers::coset
integer, dimension(:, :, :), allocatable, public coset
Definition orbital_pointers.F:45
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
pw_env_types
container for various plainwaves related things
Definition pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition pw_env_types.F:113
pw_types
Definition pw_types.F:24
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:502
qs_force_types
Definition qs_force_types.F:13
qs_integrate_potential_single
Build up the plane wave density by collocating the primitive Gaussian functions (pgf).
Definition qs_integrate_potential_single.F:30
qs_integrate_potential_single::integrate_v_gaussian_rspace
subroutine, public integrate_v_gaussian_rspace(v_rspace, qs_env, alpha, ccore, atecc)
computes the overlap of a set of Gaussians with a potential on grid
Definition qs_integrate_potential_single.F:645
qs_integrate_potential_single::integrate_v_rspace_diagonal
subroutine, public integrate_v_rspace_diagonal(v_rspace, ksmat, pmat, qs_env, calculate_forces, basis_type)
computes integrals of product of v_rspace times the diagonal block basis functions required for LRIGP...
Definition qs_integrate_potential_single.F:973
qs_integrate_potential_single::integrate_v_core_rspace
subroutine, public integrate_v_core_rspace(v_rspace, qs_env, atecc)
computes the forces/virial due to the ionic cores with a potential on grid
Definition qs_integrate_potential_single.F:483
qs_integrate_potential_single::integrate_v_rspace_one_center
subroutine, public integrate_v_rspace_one_center(v_rspace, qs_env, int_res, calculate_forces, basis_type, atomlist)
computes integrals of product of v_rspace times a one-center function required for LRIGPW
Definition qs_integrate_potential_single.F:763
qs_integrate_potential_single::integrate_ppl_rspace
subroutine, public integrate_ppl_rspace(rho_rspace, qs_env)
computes the forces/virial due to the local pseudopotential
Definition qs_integrate_potential_single.F:98
qs_integrate_potential_single::integrate_rho_nlcc
subroutine, public integrate_rho_nlcc(rho_rspace, qs_env)
computes the forces/virial due to the nlcc pseudopotential
Definition qs_integrate_potential_single.F:283
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:451
realspace_grid_types
Definition realspace_grid_types.F:18
realspace_grid_types::transfer_pw2rs
subroutine, public transfer_pw2rs(rs, pw)
...
Definition realspace_grid_types.F:698
realspace_grid_types::map_gaussian_here
pure logical function, public map_gaussian_here(rs_grid, h_inv, ra, offset, group_size, my_pos)
...
Definition realspace_grid_types.F:2148
realspace_grid_types::rs_grid_zero
subroutine, public rs_grid_zero(rs)
Initialize grid to zero.
Definition realspace_grid_types.F:1909
rs_pw_interface
Transfers densities from PW to RS grids and potentials from PW to RS.
Definition rs_pw_interface.F:14
rs_pw_interface::potential_pw2rs
subroutine, public potential_pw2rs(rs_v, v_rspace, pw_env)
transfers a potential from a pw_grid to a vector of realspace multigrids
Definition rs_pw_interface.F:152
virial_types
Definition virial_types.F:13
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
atprop_types::atprop_type
type for the atomic properties
Definition atprop_types.F:31
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:71
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_control_types::dft_control_type
Definition cp_control_types.F:559
external_potential_types::gth_potential_type
Definition external_potential_types.F:116
gaussian_gridlevels::gridlevel_info_type
Definition gaussian_gridlevels.F:34
lri_environment_types::lri_kind_type
Definition lri_environment_types.F:330
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
particle_types::particle_type
Definition particle_types.F:35
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:215
qs_force_types::qs_force_type
Definition qs_force_types.F:28
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171
realspace_grid_types::realspace_grid_type
Definition realspace_grid_types.F:137
virial_types::virial_type
Definition virial_types.F:26