(git:b279b6b)
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 ! **************************************************************************************************
30 MODULE 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 
90 CONTAINS
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 
1178 END MODULE qs_integrate_potential_single
pbc
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
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
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