(git:ccc2433)
hartree_local_methods.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 MODULE hartree_local_methods
9  USE atomic_kind_types, ONLY: atomic_kind_type,&
10  get_atomic_kind
11  USE basis_set_types, ONLY: get_gto_basis_set,&
12  gto_basis_set_type
13  USE cell_types, ONLY: cell_type
14  USE cp_control_types, ONLY: dft_control_type
15  USE hartree_local_types, ONLY: allocate_ecoul_1center,&
16  ecoul_1center_type,&
17  hartree_local_type,&
18  set_ecoul_1c
19  USE kinds, ONLY: dp
20  USE mathconstants, ONLY: fourpi,&
21  pi
22  USE message_passing, ONLY: mp_para_env_type
23  USE orbital_pointers, ONLY: indso,&
24  nsoset
25  USE pw_env_types, ONLY: pw_env_get,&
26  pw_env_type
27  USE pw_poisson_types, ONLY: pw_poisson_periodic,&
28  pw_poisson_type
29  USE qs_charges_types, ONLY: qs_charges_type
30  USE qs_environment_types, ONLY: get_qs_env,&
31  qs_environment_type
32  USE qs_grid_atom, ONLY: grid_atom_type
33  USE qs_harmonics_atom, ONLY: get_none0_cg_list,&
34  harmonics_atom_type
35  USE qs_kind_types, ONLY: get_qs_kind,&
36  get_qs_kind_set,&
37  qs_kind_type
38  USE qs_local_rho_types, ONLY: get_local_rho,&
39  local_rho_type,&
40  rhoz_type
41  USE qs_rho0_types, ONLY: get_rho0_mpole,&
42  rho0_atom_type,&
43  rho0_mpole_type
44  USE qs_rho_atom_types, ONLY: get_rho_atom,&
45  rho_atom_coeff,&
46  rho_atom_type
47  USE util, ONLY: get_limit
48 #include "./base/base_uses.f90"
49 
50  IMPLICIT NONE
51 
52  PRIVATE
53 
54  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'hartree_local_methods'
55 
56  ! Public Subroutine
57 
58  PUBLIC :: init_coulomb_local, calculate_vh_1center, vh_1c_gg_integrals
59 
60 CONTAINS
61 
62 ! **************************************************************************************************
63 !> \brief ...
64 !> \param hartree_local ...
65 !> \param natom ...
66 ! **************************************************************************************************
67  SUBROUTINE init_coulomb_local(hartree_local, natom)
68 
69  TYPE(hartree_local_type), POINTER :: hartree_local
70  INTEGER, INTENT(IN) :: natom
71 
72  CHARACTER(len=*), PARAMETER :: routinen = 'init_coulomb_local'
73 
74  INTEGER :: handle
75  TYPE(ecoul_1center_type), DIMENSION(:), POINTER :: ecoul_1c
76 
77  CALL timeset(routinen, handle)
78 
79  NULLIFY (ecoul_1c)
80  ! Allocate and Initialize 1-center Potentials and Integrals
81  CALL allocate_ecoul_1center(ecoul_1c, natom)
82  hartree_local%ecoul_1c => ecoul_1c
83 
84  CALL timestop(handle)
85 
86  END SUBROUTINE init_coulomb_local
87 
88 ! **************************************************************************************************
89 !> \brief Calculates Hartree potential for hard and soft densities (including
90 !> nuclear charge and compensation charges) using numerical integration
91 !> \param vrad_h ...
92 !> \param vrad_s ...
93 !> \param rrad_h ...
94 !> \param rrad_s ...
95 !> \param rrad_0 ...
96 !> \param rrad_z ...
97 !> \param grid_atom ...
98 !> \par History
99 !> 05.2012 JGH refactoring
100 !> \author ??
101 ! **************************************************************************************************
102  SUBROUTINE calculate_vh_1center(vrad_h, vrad_s, rrad_h, rrad_s, rrad_0, rrad_z, grid_atom)
103 
104  REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: vrad_h, vrad_s
105  TYPE(rho_atom_coeff), DIMENSION(:), INTENT(IN) :: rrad_h, rrad_s
106  REAL(dp), DIMENSION(:, :), INTENT(IN) :: rrad_0
107  REAL(dp), DIMENSION(:), INTENT(IN) :: rrad_z
108  TYPE(grid_atom_type), POINTER :: grid_atom
109 
110  CHARACTER(len=*), PARAMETER :: routinen = 'calculate_Vh_1center'
111 
112  INTEGER :: handle, ir, iso, ispin, l_ang, &
113  max_s_harm, nchannels, nr, nspins
114  REAL(dp) :: i1_down, i1_up, i2_down, i2_up, prefactor
115  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rho_1, rho_2
116  REAL(dp), DIMENSION(:), POINTER :: wr
117  REAL(dp), DIMENSION(:, :), POINTER :: oor2l, r2l
118 
119  CALL timeset(routinen, handle)
120 
121  nr = grid_atom%nr
122  max_s_harm = SIZE(vrad_h, 2)
123  nspins = SIZE(rrad_h, 1)
124  nchannels = SIZE(rrad_0, 2)
125 
126  r2l => grid_atom%rad2l
127  oor2l => grid_atom%oorad2l
128  wr => grid_atom%wr
129 
130  ALLOCATE (rho_1(nr, max_s_harm), rho_2(nr, max_s_harm))
131  rho_1 = 0.0_dp
132  rho_2 = 0.0_dp
133 
134  ! Case lm = 0
135  rho_1(:, 1) = rrad_z(:)
136  rho_2(:, 1) = rrad_0(:, 1)
137 
138  DO iso = 2, nchannels
139  rho_2(:, iso) = rrad_0(:, iso)
140  END DO
141 
142  DO iso = 1, max_s_harm
143  DO ispin = 1, nspins
144  rho_1(:, iso) = rho_1(:, iso) + rrad_h(ispin)%r_coef(:, iso)
145  rho_2(:, iso) = rho_2(:, iso) + rrad_s(ispin)%r_coef(:, iso)
146  END DO
147 
148  l_ang = indso(1, iso)
149  prefactor = fourpi/(2._dp*l_ang + 1._dp)
150 
151  rho_1(:, iso) = rho_1(:, iso)*wr(:)
152  rho_2(:, iso) = rho_2(:, iso)*wr(:)
153 
154  i1_up = 0.0_dp
155  i1_down = 0.0_dp
156  i2_up = 0.0_dp
157  i2_down = 0.0_dp
158 
159  i1_up = r2l(nr, l_ang)*rho_1(nr, iso)
160  i2_up = r2l(nr, l_ang)*rho_2(nr, iso)
161 
162  DO ir = nr - 1, 1, -1
163  i1_down = i1_down + oor2l(ir, l_ang + 1)*rho_1(ir, iso)
164  i2_down = i2_down + oor2l(ir, l_ang + 1)*rho_2(ir, iso)
165  END DO
166 
167  vrad_h(nr, iso) = vrad_h(nr, iso) + prefactor* &
168  (oor2l(nr, l_ang + 1)*i1_up + r2l(nr, l_ang)*i1_down)
169  vrad_s(nr, iso) = vrad_s(nr, iso) + prefactor* &
170  (oor2l(nr, l_ang + 1)*i2_up + r2l(nr, l_ang)*i2_down)
171 
172  DO ir = nr - 1, 1, -1
173  i1_up = i1_up + r2l(ir, l_ang)*rho_1(ir, iso)
174  i1_down = i1_down - oor2l(ir, l_ang + 1)*rho_1(ir, iso)
175  i2_up = i2_up + r2l(ir, l_ang)*rho_2(ir, iso)
176  i2_down = i2_down - oor2l(ir, l_ang + 1)*rho_2(ir, iso)
177 
178  vrad_h(ir, iso) = vrad_h(ir, iso) + prefactor* &
179  (oor2l(ir, l_ang + 1)*i1_up + r2l(ir, l_ang)*i1_down)
180  vrad_s(ir, iso) = vrad_s(ir, iso) + prefactor* &
181  (oor2l(ir, l_ang + 1)*i2_up + r2l(ir, l_ang)*i2_down)
182 
183  END DO
184 
185  END DO
186 
187  DEALLOCATE (rho_1, rho_2)
188 
189  CALL timestop(handle)
190 
191  END SUBROUTINE calculate_vh_1center
192 
193 ! **************************************************************************************************
194 !> \brief Calculates one center GAPW Hartree energies and matrix elements
195 !> Hartree potentials are input
196 !> Takes possible background charge into account
197 !> Special case for densities without core charge
198 !> \param qs_env ...
199 !> \param energy_hartree_1c ...
200 !> \param ecoul_1c ...
201 !> \param local_rho_set ...
202 !> \param para_env ...
203 !> \param tddft ...
204 !> \param local_rho_set_2nd ...
205 !> \param core_2nd ...
206 !> \par History
207 !> 05.2012 JGH refactoring
208 !> \author ??
209 ! **************************************************************************************************
210  SUBROUTINE vh_1c_gg_integrals(qs_env, energy_hartree_1c, ecoul_1c, local_rho_set, para_env, tddft, local_rho_set_2nd, &
211  core_2nd)
212 
213  TYPE(qs_environment_type), POINTER :: qs_env
214  REAL(kind=dp), INTENT(out) :: energy_hartree_1c
215  TYPE(ecoul_1center_type), DIMENSION(:), POINTER :: ecoul_1c
216  TYPE(local_rho_type), POINTER :: local_rho_set
217  TYPE(mp_para_env_type), POINTER :: para_env
218  LOGICAL, INTENT(IN) :: tddft
219  TYPE(local_rho_type), OPTIONAL, POINTER :: local_rho_set_2nd
220  LOGICAL, INTENT(IN), OPTIONAL :: core_2nd
221 
222  CHARACTER(LEN=*), PARAMETER :: routinen = 'Vh_1c_gg_integrals'
223 
224  INTEGER :: bo(2), handle, iat, iatom, ikind, ipgf1, is1, iset1, iso, l_ang, llmax, lmax0, &
225  lmax0_2nd, lmax_0, m1, max_iso, max_iso_not0, max_s_harm, maxl, maxso, mepos, n1, nat, &
226  nchan_0, nkind, nr, nset, nsotot, nspins, num_pe
227  INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list
228  INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list
229  INTEGER, DIMENSION(:), POINTER :: atom_list, lmax, lmin, npgf
230  LOGICAL :: core_charge, l_2nd_local_rho, &
231  my_core_2nd, my_periodic, paw_atom
232  REAL(dp) :: back_ch, factor
233  REAL(dp), ALLOCATABLE, DIMENSION(:) :: gexp, sqrtwr
234  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: avh1b_00, avh1b_hh, avh1b_ss, g0_h_w
235  REAL(dp), DIMENSION(:), POINTER :: rrad_z, vrrad_z
236  REAL(dp), DIMENSION(:, :), POINTER :: g0_h, g0_h_2nd, gsph, rrad_0, vh1_h, &
237  vh1_s, vrrad_0, zet
238  REAL(dp), DIMENSION(:, :, :), POINTER :: my_cg, qlm_gg, qlm_gg_2nd
239  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
240  TYPE(cell_type), POINTER :: cell
241  TYPE(dft_control_type), POINTER :: dft_control
242  TYPE(grid_atom_type), POINTER :: grid_atom
243  TYPE(gto_basis_set_type), POINTER :: basis_1c
244  TYPE(harmonics_atom_type), POINTER :: harmonics
245  TYPE(pw_env_type), POINTER :: pw_env
246  TYPE(pw_poisson_type), POINTER :: poisson_env
247  TYPE(qs_charges_type), POINTER :: qs_charges
248  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
249  TYPE(rho0_atom_type), DIMENSION(:), POINTER :: rho0_atom_set, rho0_atom_set_2nd
250  TYPE(rho0_mpole_type), POINTER :: rho0_mpole, rho0_mpole_2nd
251  TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set, rho_atom_set_2nd
252  TYPE(rho_atom_type), POINTER :: rho_atom
253  TYPE(rhoz_type), DIMENSION(:), POINTER :: rhoz_set, rhoz_set_2nd
254 
255  CALL timeset(routinen, handle)
256 
257  NULLIFY (cell, dft_control, poisson_env, pw_env, qs_charges)
258  NULLIFY (atomic_kind_set, qs_kind_set, rho_atom_set, rho0_atom_set)
259  NULLIFY (rho0_mpole, rhoz_set)
260  NULLIFY (atom_list, grid_atom, harmonics)
261  NULLIFY (basis_1c, lmin, lmax, npgf, zet)
262  NULLIFY (gsph)
263 
264  CALL get_qs_env(qs_env=qs_env, &
265  cell=cell, dft_control=dft_control, &
266  atomic_kind_set=atomic_kind_set, &
267  qs_kind_set=qs_kind_set, &
268  pw_env=pw_env, qs_charges=qs_charges)
269 
270  CALL pw_env_get(pw_env, poisson_env=poisson_env)
271  my_periodic = (poisson_env%method == pw_poisson_periodic)
272 
273  back_ch = qs_charges%background*cell%deth
274 
275  ! rhoz_set is not accessed in TDDFT
276  CALL get_local_rho(local_rho_set, rho_atom_set, rho0_atom_set, rho0_mpole, rhoz_set) ! for integral space
277 
278  ! for forces we need a second local_rho_set
279  l_2nd_local_rho = .false.
280  IF (PRESENT(local_rho_set_2nd)) THEN ! for potential
281  l_2nd_local_rho = .true.
282  NULLIFY (rho_atom_set_2nd, rho0_atom_set_2nd, rhoz_set_2nd) ! for potential
283  CALL get_local_rho(local_rho_set_2nd, rho_atom_set_2nd, rho0_atom_set_2nd, rho0_mpole_2nd, rhoz_set=rhoz_set_2nd)
284  END IF
285 
286  nkind = SIZE(atomic_kind_set, 1)
287  nspins = dft_control%nspins
288 
289  core_charge = .NOT. tddft ! for forces mixed version
290  my_core_2nd = .true.
291  IF (PRESENT(core_2nd)) my_core_2nd = .NOT. core_2nd ! if my_core_2nd true, include core charge
292 
293  ! The aim of the following code was to return immediately if the subroutine
294  ! was called for triplet excited states in spin-restricted case. This check
295  ! is also performed before invocation of this subroutine. It should be save
296  ! to remove the optional argument 'do_triplet' from the subroutine interface.
297  !IF (tddft) THEN
298  ! CPASSERT(PRESENT(do_triplet))
299  ! IF (nspins == 1 .AND. do_triplet) RETURN
300  !END IF
301 
302  CALL get_qs_kind_set(qs_kind_set, maxg_iso_not0=max_iso)
303  CALL get_rho0_mpole(rho0_mpole=rho0_mpole, lmax_0=lmax_0)
304 
305  ! Put to 0 the local hartree energy contribution from 1 center integrals
306  energy_hartree_1c = 0.0_dp
307 
308  ! Here starts the loop over all the atoms
309  DO ikind = 1, nkind
310 
311  CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
312  CALL get_qs_kind(qs_kind_set(ikind), &
313  grid_atom=grid_atom, &
314  harmonics=harmonics, ngrid_rad=nr, &
315  max_iso_not0=max_iso_not0, paw_atom=paw_atom)
316  CALL get_qs_kind(qs_kind_set(ikind), &
317  basis_set=basis_1c, basis_type="GAPW_1C")
318 
319  IF (paw_atom) THEN
320 !=========== PAW ===============
321  CALL get_gto_basis_set(gto_basis_set=basis_1c, lmax=lmax, lmin=lmin, &
322  maxso=maxso, npgf=npgf, maxl=maxl, &
323  nset=nset, zet=zet)
324 
325  max_s_harm = harmonics%max_s_harm
326  llmax = harmonics%llmax
327 
328  nsotot = maxso*nset
329  ALLOCATE (gsph(nr, nsotot))
330  ALLOCATE (gexp(nr))
331  ALLOCATE (sqrtwr(nr), g0_h_w(nr, 0:lmax_0))
332 
333  NULLIFY (vh1_h, vh1_s)
334  ALLOCATE (vh1_h(nr, max_iso_not0))
335  ALLOCATE (vh1_s(nr, max_iso_not0))
336 
337  ALLOCATE (avh1b_hh(nsotot, nsotot))
338  ALLOCATE (avh1b_ss(nsotot, nsotot))
339  ALLOCATE (avh1b_00(nsotot, nsotot))
340  ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm))
341 
342  NULLIFY (qlm_gg, g0_h)
343  CALL get_rho0_mpole(rho0_mpole=rho0_mpole, ikind=ikind, &
344  l0_ikind=lmax0, &
345  qlm_gg=qlm_gg, g0_h=g0_h) ! Qlm_gg of density
346 
347  IF (PRESENT(local_rho_set_2nd)) THEN ! for potential
348  NULLIFY (qlm_gg_2nd, g0_h_2nd)
349  CALL get_rho0_mpole(rho0_mpole=rho0_mpole_2nd, ikind=ikind, &
350  l0_ikind=lmax0_2nd, &
351  qlm_gg=qlm_gg_2nd, g0_h=g0_h_2nd) ! Qlm_gg of density
352  END IF
353  nchan_0 = nsoset(lmax0)
354 
355  IF (nchan_0 > max_iso_not0) cpabort("channels for rho0 > # max of spherical harmonics")
356 
357  NULLIFY (rrad_z, my_cg)
358  my_cg => harmonics%my_CG
359 
360  ! set to zero temporary arrays
361  sqrtwr = 0.0_dp
362  g0_h_w = 0.0_dp
363  gexp = 0.0_dp
364  gsph = 0.0_dp
365 
366  sqrtwr(1:nr) = sqrt(grid_atom%wr(1:nr))
367  DO l_ang = 0, lmax0
368  g0_h_w(1:nr, l_ang) = g0_h(1:nr, l_ang)*grid_atom%wr(1:nr)
369  END DO
370 
371  m1 = 0
372  DO iset1 = 1, nset
373  n1 = nsoset(lmax(iset1))
374  DO ipgf1 = 1, npgf(iset1)
375  gexp(1:nr) = exp(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))*sqrtwr(1:nr)
376  DO is1 = nsoset(lmin(iset1) - 1) + 1, nsoset(lmax(iset1))
377  iso = is1 + (ipgf1 - 1)*n1 + m1
378  l_ang = indso(1, is1)
379  gsph(1:nr, iso) = grid_atom%rad2l(1:nr, l_ang)*gexp(1:nr)
380  END DO ! is1
381  END DO ! ipgf1
382  m1 = m1 + maxso
383  END DO ! iset1
384 
385  ! Distribute the atoms of this kind
386  num_pe = para_env%num_pe
387  mepos = para_env%mepos
388  bo = get_limit(nat, num_pe, mepos)
389 
390  DO iat = bo(1), bo(2) !1,nat
391  iatom = atom_list(iat)
392  rho_atom => rho_atom_set(iatom)
393 
394  NULLIFY (rrad_z, vrrad_z, rrad_0, vrrad_0)
395  IF (core_charge) THEN
396  rrad_z => rhoz_set(ikind)%r_coef ! for density
397  END IF
398  IF (my_core_2nd) THEN
399  IF (l_2nd_local_rho) THEN
400  vrrad_z => rhoz_set_2nd(ikind)%vr_coef ! for potential
401  ELSE
402  vrrad_z => rhoz_set(ikind)%vr_coef ! for potential
403  END IF
404  END IF
405  rrad_0 => rho0_atom_set(iatom)%rho0_rad_h%r_coef ! for density
406  vrrad_0 => rho0_atom_set(iatom)%vrho0_rad_h%r_coef
407  IF (l_2nd_local_rho) THEN
408  rho_atom => rho_atom_set_2nd(iatom)
409  vrrad_0 => rho0_atom_set_2nd(iatom)%vrho0_rad_h%r_coef ! for potential
410  END IF
411  IF (my_periodic .AND. back_ch .GT. 1.e-3_dp) THEN
412  factor = -2.0_dp*pi/3.0_dp*sqrt(fourpi)*qs_charges%background
413  ELSE
414  factor = 0._dp
415  END IF
416 
417  CALL vh_1c_atom_potential(rho_atom, vrrad_0, &
418  grid_atom, my_core_2nd, vrrad_z, vh1_h, vh1_s, & ! core charge for potential (2nd)
419  nchan_0, nspins, max_iso_not0, factor)
420 
421  IF (l_2nd_local_rho) rho_atom => rho_atom_set(iatom) ! rho_atom for density
422 
423  CALL vh_1c_atom_energy(energy_hartree_1c, ecoul_1c, rho_atom, rrad_0, &
424  grid_atom, iatom, core_charge, rrad_z, vh1_h, vh1_s, & ! core charge for density
425  nchan_0, nspins, max_iso_not0)
426 
427  IF (l_2nd_local_rho) rho_atom => rho_atom_set_2nd(iatom) ! rho_atom for potential (2nd)
428 
429  CALL vh_1c_atom_integrals(rho_atom, & ! results (int_local_h and int_local_s) written on rho_atom_2nd
430  ! int_local_h and int_local_s are used in update_ks_atom
431  ! on int_local_h mixed core / non-core
432  avh1b_hh, avh1b_ss, avh1b_00, vh1_h, vh1_s, max_iso_not0, &
433  max_s_harm, llmax, cg_list, cg_n_list, &
434  nset, npgf, lmin, lmax, nsotot, maxso, nspins, nchan_0, gsph, &
435  g0_h_w, my_cg, qlm_gg) ! Qlm_gg for density from local_rho_set
436 
437  END DO ! iat
438 
439  DEALLOCATE (avh1b_hh)
440  DEALLOCATE (avh1b_ss)
441  DEALLOCATE (avh1b_00)
442  DEALLOCATE (vh1_h, vh1_s)
443  DEALLOCATE (cg_list, cg_n_list)
444  DEALLOCATE (gsph)
445  DEALLOCATE (gexp)
446  DEALLOCATE (sqrtwr, g0_h_w)
447 
448  ELSE
449 !=========== NO PAW ===============
450 ! This term is taken care of using the core density as in GPW
451  cycle
452  END IF ! paw
453  END DO ! ikind
454 
455  CALL para_env%sum(energy_hartree_1c)
456 
457  CALL timestop(handle)
458 
459  END SUBROUTINE vh_1c_gg_integrals
460 
461 ! **************************************************************************************************
462 
463 ! **************************************************************************************************
464 !> \brief ...
465 !> \param rho_atom ...
466 !> \param vrrad_0 ...
467 !> \param grid_atom ...
468 !> \param core_charge ...
469 !> \param vrrad_z ...
470 !> \param Vh1_h ...
471 !> \param Vh1_s ...
472 !> \param nchan_0 ...
473 !> \param nspins ...
474 !> \param max_iso_not0 ...
475 !> \param bfactor ...
476 ! **************************************************************************************************
477  SUBROUTINE vh_1c_atom_potential(rho_atom, vrrad_0, &
478  grid_atom, core_charge, vrrad_z, Vh1_h, Vh1_s, &
479  nchan_0, nspins, max_iso_not0, bfactor)
480 
481  TYPE(rho_atom_type), POINTER :: rho_atom
482  REAL(dp), DIMENSION(:, :), POINTER :: vrrad_0
483  TYPE(grid_atom_type), POINTER :: grid_atom
484  LOGICAL, INTENT(IN) :: core_charge
485  REAL(dp), DIMENSION(:), POINTER :: vrrad_z
486  REAL(dp), DIMENSION(:, :), POINTER :: vh1_h, vh1_s
487  INTEGER, INTENT(IN) :: nchan_0, nspins, max_iso_not0
488  REAL(dp), INTENT(IN) :: bfactor
489 
490  INTEGER :: ir, iso, ispin, nr
491  TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: vr_h, vr_s
492 
493  nr = grid_atom%nr
494 
495  NULLIFY (vr_h, vr_s)
496  CALL get_rho_atom(rho_atom=rho_atom, vrho_rad_h=vr_h, vrho_rad_s=vr_s)
497 
498  vh1_h = 0.0_dp
499  vh1_s = 0.0_dp
500 
501  IF (core_charge) vh1_h(:, 1) = vrrad_z(:)
502 
503  DO iso = 1, nchan_0
504  vh1_s(:, iso) = vrrad_0(:, iso)
505  END DO
506 
507  DO ispin = 1, nspins
508  DO iso = 1, max_iso_not0
509  vh1_h(:, iso) = vh1_h(:, iso) + vr_h(ispin)%r_coef(:, iso)
510  vh1_s(:, iso) = vh1_s(:, iso) + vr_s(ispin)%r_coef(:, iso)
511  END DO
512  END DO
513 
514  IF (bfactor /= 0._dp) THEN
515  DO ir = 1, nr
516  vh1_h(ir, 1) = vh1_h(ir, 1) + bfactor*grid_atom%rad2(ir)*grid_atom%wr(ir)
517  vh1_s(ir, 1) = vh1_s(ir, 1) + bfactor*grid_atom%rad2(ir)*grid_atom%wr(ir)
518  END DO
519  END IF
520 
521  END SUBROUTINE vh_1c_atom_potential
522 
523 ! **************************************************************************************************
524 
525 ! **************************************************************************************************
526 !> \brief ...
527 !> \param energy_hartree_1c ...
528 !> \param ecoul_1c ...
529 !> \param rho_atom ...
530 !> \param rrad_0 ...
531 !> \param grid_atom ...
532 !> \param iatom ...
533 !> \param core_charge ...
534 !> \param rrad_z ...
535 !> \param Vh1_h ...
536 !> \param Vh1_s ...
537 !> \param nchan_0 ...
538 !> \param nspins ...
539 !> \param max_iso_not0 ...
540 ! **************************************************************************************************
541  SUBROUTINE vh_1c_atom_energy(energy_hartree_1c, ecoul_1c, rho_atom, rrad_0, &
542  grid_atom, iatom, core_charge, rrad_z, Vh1_h, Vh1_s, &
543  nchan_0, nspins, max_iso_not0)
544 
545  REAL(dp), INTENT(INOUT) :: energy_hartree_1c
546  TYPE(ecoul_1center_type), DIMENSION(:), POINTER :: ecoul_1c
547  TYPE(rho_atom_type), POINTER :: rho_atom
548  REAL(dp), DIMENSION(:, :), POINTER :: rrad_0
549  TYPE(grid_atom_type), POINTER :: grid_atom
550  INTEGER, INTENT(IN) :: iatom
551  LOGICAL, INTENT(IN) :: core_charge
552  REAL(dp), DIMENSION(:), POINTER :: rrad_z
553  REAL(dp), DIMENSION(:, :), POINTER :: vh1_h, vh1_s
554  INTEGER, INTENT(IN) :: nchan_0, nspins, max_iso_not0
555 
556  INTEGER :: iso, ispin, nr
557  REAL(dp) :: ecoul_1_0, ecoul_1_h, ecoul_1_s, &
558  ecoul_1_z
559  TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: r_h, r_s
560 
561  nr = grid_atom%nr
562 
563  NULLIFY (r_h, r_s)
564  CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, rho_rad_s=r_s)
565 
566  ! Calculate the contributions to Ecoul coming from Vh1_h*rhoz
567  ecoul_1_z = 0.0_dp
568  IF (core_charge) THEN
569  ecoul_1_z = 0.5_dp*sum(vh1_h(:, 1)*rrad_z(:)*grid_atom%wr(:))
570  END IF
571 
572  ! Calculate the contributions to Ecoul coming from Vh1_s*rho0
573  ecoul_1_0 = 0.0_dp
574  DO iso = 1, nchan_0
575  ecoul_1_0 = ecoul_1_0 + 0.5_dp*sum(vh1_s(:, iso)*rrad_0(:, iso)*grid_atom%wr(:))
576  END DO
577 
578  ! Calculate the contributions to Ecoul coming from Vh1_h*rho1_h and Vh1_s*rho1_s
579  ecoul_1_s = 0.0_dp
580  ecoul_1_h = 0.0_dp
581  DO ispin = 1, nspins
582  DO iso = 1, max_iso_not0
583  ecoul_1_s = ecoul_1_s + 0.5_dp*sum(vh1_s(:, iso)*r_s(ispin)%r_coef(:, iso)*grid_atom%wr(:))
584  ecoul_1_h = ecoul_1_h + 0.5_dp*sum(vh1_h(:, iso)*r_h(ispin)%r_coef(:, iso)*grid_atom%wr(:))
585  END DO
586  END DO
587 
588  CALL set_ecoul_1c(ecoul_1c, iatom, ecoul_1_z=ecoul_1_z, ecoul_1_0=ecoul_1_0)
589  CALL set_ecoul_1c(ecoul_1c=ecoul_1c, iatom=iatom, ecoul_1_h=ecoul_1_h, ecoul_1_s=ecoul_1_s)
590 
591  energy_hartree_1c = energy_hartree_1c + ecoul_1_z - ecoul_1_0
592  energy_hartree_1c = energy_hartree_1c + ecoul_1_h - ecoul_1_s
593 
594  END SUBROUTINE vh_1c_atom_energy
595 
596 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
597 
598 ! **************************************************************************************************
599 !> \brief ...
600 !> \param rho_atom ...
601 !> \param aVh1b_hh ...
602 !> \param aVh1b_ss ...
603 !> \param aVh1b_00 ...
604 !> \param Vh1_h ...
605 !> \param Vh1_s ...
606 !> \param max_iso_not0 ...
607 !> \param max_s_harm ...
608 !> \param llmax ...
609 !> \param cg_list ...
610 !> \param cg_n_list ...
611 !> \param nset ...
612 !> \param npgf ...
613 !> \param lmin ...
614 !> \param lmax ...
615 !> \param nsotot ...
616 !> \param maxso ...
617 !> \param nspins ...
618 !> \param nchan_0 ...
619 !> \param gsph ...
620 !> \param g0_h_w ...
621 !> \param my_CG ...
622 !> \param Qlm_gg ...
623 ! **************************************************************************************************
624  SUBROUTINE vh_1c_atom_integrals(rho_atom, &
625  aVh1b_hh, aVh1b_ss, aVh1b_00, Vh1_h, Vh1_s, max_iso_not0, &
626  max_s_harm, llmax, cg_list, cg_n_list, &
627  nset, npgf, lmin, lmax, nsotot, maxso, nspins, nchan_0, gsph, &
628  g0_h_w, my_CG, Qlm_gg)
629 
630  TYPE(rho_atom_type), POINTER :: rho_atom
631  REAL(dp), DIMENSION(:, :) :: avh1b_hh, avh1b_ss, avh1b_00
632  REAL(dp), DIMENSION(:, :), POINTER :: vh1_h, vh1_s
633  INTEGER, INTENT(IN) :: max_iso_not0, max_s_harm, llmax
634  INTEGER, DIMENSION(:, :, :) :: cg_list
635  INTEGER, DIMENSION(:) :: cg_n_list
636  INTEGER, INTENT(IN) :: nset
637  INTEGER, DIMENSION(:), POINTER :: npgf, lmin, lmax
638  INTEGER, INTENT(IN) :: nsotot, maxso, nspins, nchan_0
639  REAL(dp), DIMENSION(:, :), POINTER :: gsph
640  REAL(dp), DIMENSION(:, 0:) :: g0_h_w
641  REAL(dp), DIMENSION(:, :, :), POINTER :: my_cg, qlm_gg
642 
643  INTEGER :: icg, ipgf1, ipgf2, ir, is1, is2, iset1, &
644  iset2, iso, iso1, iso2, ispin, l_ang, &
645  m1, m2, max_iso_not0_local, n1, n2, nr
646  REAL(dp) :: gvg_0, gvg_h, gvg_s
647  TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: int_local_h, int_local_s
648 
649  NULLIFY (int_local_h, int_local_s)
650  CALL get_rho_atom(rho_atom=rho_atom, &
651  ga_vlocal_gb_h=int_local_h, &
652  ga_vlocal_gb_s=int_local_s)
653 
654  ! Calculate the integrals of the potential with 2 primitives
655  avh1b_hh = 0.0_dp
656  avh1b_ss = 0.0_dp
657  avh1b_00 = 0.0_dp
658 
659  nr = SIZE(gsph, 1)
660 
661  m1 = 0
662  DO iset1 = 1, nset
663  m2 = 0
664  DO iset2 = 1, nset
665  CALL get_none0_cg_list(my_cg, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
666  max_s_harm, llmax, cg_list, cg_n_list, max_iso_not0_local)
667 
668  n1 = nsoset(lmax(iset1))
669  DO ipgf1 = 1, npgf(iset1)
670  n2 = nsoset(lmax(iset2))
671  DO ipgf2 = 1, npgf(iset2)
672  ! with contributions to V1_s*rho0
673  DO iso = 1, nchan_0
674  l_ang = indso(1, iso)
675  gvg_0 = sum(vh1_s(:, iso)*g0_h_w(:, l_ang))
676  DO icg = 1, cg_n_list(iso)
677  is1 = cg_list(1, icg, iso)
678  is2 = cg_list(2, icg, iso)
679 
680  iso1 = is1 + n1*(ipgf1 - 1) + m1
681  iso2 = is2 + n2*(ipgf2 - 1) + m2
682  gvg_h = 0.0_dp
683  gvg_s = 0.0_dp
684 
685  DO ir = 1, nr
686  gvg_h = gvg_h + gsph(ir, iso1)*gsph(ir, iso2)*vh1_h(ir, iso)
687  gvg_s = gvg_s + gsph(ir, iso1)*gsph(ir, iso2)*vh1_s(ir, iso)
688  END DO ! ir
689 
690  avh1b_hh(iso1, iso2) = avh1b_hh(iso1, iso2) + gvg_h*my_cg(is1, is2, iso)
691  avh1b_ss(iso1, iso2) = avh1b_ss(iso1, iso2) + gvg_s*my_cg(is1, is2, iso)
692  avh1b_00(iso1, iso2) = avh1b_00(iso1, iso2) + gvg_0*qlm_gg(iso1, iso2, iso)
693 
694  END DO !icg
695  END DO ! iso
696  ! without contributions to V1_s*rho0
697  DO iso = nchan_0 + 1, max_iso_not0
698  DO icg = 1, cg_n_list(iso)
699  is1 = cg_list(1, icg, iso)
700  is2 = cg_list(2, icg, iso)
701 
702  iso1 = is1 + n1*(ipgf1 - 1) + m1
703  iso2 = is2 + n2*(ipgf2 - 1) + m2
704  gvg_h = 0.0_dp
705  gvg_s = 0.0_dp
706 
707  DO ir = 1, nr
708  gvg_h = gvg_h + gsph(ir, iso1)*gsph(ir, iso2)*vh1_h(ir, iso)
709  gvg_s = gvg_s + gsph(ir, iso1)*gsph(ir, iso2)*vh1_s(ir, iso)
710  END DO ! ir
711 
712  avh1b_hh(iso1, iso2) = avh1b_hh(iso1, iso2) + gvg_h*my_cg(is1, is2, iso)
713  avh1b_ss(iso1, iso2) = avh1b_ss(iso1, iso2) + gvg_s*my_cg(is1, is2, iso)
714 
715  END DO !icg
716  END DO ! iso
717  END DO ! ipgf2
718  END DO ! ipgf1
719  m2 = m2 + maxso
720  END DO ! iset2
721  m1 = m1 + maxso
722  END DO !iset1
723  DO ispin = 1, nspins
724  CALL daxpy(nsotot*nsotot, 1.0_dp, avh1b_hh, 1, int_local_h(ispin)%r_coef, 1)
725  CALL daxpy(nsotot*nsotot, 1.0_dp, avh1b_ss, 1, int_local_s(ispin)%r_coef, 1)
726  CALL daxpy(nsotot*nsotot, -1.0_dp, avh1b_00, 1, int_local_h(ispin)%r_coef, 1)
727  CALL daxpy(nsotot*nsotot, -1.0_dp, avh1b_00, 1, int_local_s(ispin)%r_coef, 1)
728  END DO ! ispin
729 
730  END SUBROUTINE vh_1c_atom_integrals
731 
732 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
733 
734 END MODULE hartree_local_methods
735 
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
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
hartree_local_methods
Definition: hartree_local_methods.F:8
hartree_local_methods::init_coulomb_local
subroutine, public init_coulomb_local(hartree_local, natom)
...
Definition: hartree_local_methods.F:68
hartree_local_methods::vh_1c_gg_integrals
subroutine, public vh_1c_gg_integrals(qs_env, energy_hartree_1c, ecoul_1c, local_rho_set, para_env, tddft, local_rho_set_2nd, core_2nd)
Calculates one center GAPW Hartree energies and matrix elements Hartree potentials are input Takes po...
Definition: hartree_local_methods.F:212
hartree_local_methods::calculate_vh_1center
subroutine, public calculate_vh_1center(vrad_h, vrad_s, rrad_h, rrad_s, rrad_0, rrad_z, grid_atom)
Calculates Hartree potential for hard and soft densities (including nuclear charge and compensation c...
Definition: hartree_local_methods.F:103
hartree_local_types
Definition: hartree_local_types.F:8
hartree_local_types::allocate_ecoul_1center
subroutine, public allocate_ecoul_1center(ecoul_1c, natom)
...
Definition: hartree_local_types.F:58
hartree_local_types::set_ecoul_1c
subroutine, public set_ecoul_1c(ecoul_1c, iatom, ecoul_1_h, ecoul_1_s, ecoul_1_z, ecoul_1_0)
...
Definition: hartree_local_types.F:165
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition: mathconstants.F:110
mathconstants::fourpi
real(kind=dp), parameter, public fourpi
Definition: mathconstants.F:113
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::nsoset
integer, dimension(:), allocatable, public nsoset
Definition: orbital_pointers.F:42
orbital_pointers::indso
integer, dimension(:, :), allocatable, public indso
Definition: orbital_pointers.F:43
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_poisson_types
functions related to the poisson solver on regular grids
Definition: pw_poisson_types.F:22
pw_poisson_types::pw_poisson_periodic
integer, parameter, public pw_poisson_periodic
Definition: pw_poisson_types.F:66
qs_charges_types
container for information about total charges on the grids
Definition: qs_charges_types.F:14
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_grid_atom
Definition: qs_grid_atom.F:8
qs_harmonics_atom
Definition: qs_harmonics_atom.F:8
qs_kind_types
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
qs_kind_types::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
qs_kind_types::get_qs_kind_set
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr)
Get attributes of an atomic kind set.
Definition: qs_kind_types.F:864
qs_local_rho_types
Definition: qs_local_rho_types.F:8
qs_local_rho_types::get_local_rho
subroutine, public get_local_rho(local_rho_set, rho_atom_set, rho0_atom_set, rho0_mpole, rhoz_set)
...
Definition: qs_local_rho_types.F:173
qs_rho0_types
Definition: qs_rho0_types.F:8
qs_rho0_types::get_rho0_mpole
subroutine, public get_rho0_mpole(rho0_mpole, g0_h, vg0_h, iat, ikind, lmax_0, l0_ikind, mp_gau_ikind, mp_rho, norm_g0l_h, Qlm_gg, Qlm_car, Qlm_tot, zet0_h, igrid_zet0_s, rpgf0_h, rpgf0_s, max_rpgf0_s, rho0_s_rs, rho0_s_gs)
...
Definition: qs_rho0_types.F:425
qs_rho_atom_types
Definition: qs_rho_atom_types.F:8
qs_rho_atom_types::get_rho_atom
subroutine, public get_rho_atom(rho_atom, cpc_h, cpc_s, rho_rad_h, rho_rad_s, drho_rad_h, drho_rad_s, vrho_rad_h, vrho_rad_s, rho_rad_h_d, rho_rad_s_d, ga_Vlocal_gb_h, ga_Vlocal_gb_s, int_scr_h, int_scr_s)
...
Definition: qs_rho_atom_types.F:195
util
All kind of helpful little routines.
Definition: util.F:14
util::get_limit
pure integer function, dimension(2), public get_limit(m, n, me)
divide m entries into n parts, return size of part me
Definition: util.F:333