(git:123285f)
Loading...
Searching...
No Matches
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-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7! **************************************************************************************************
8MODULE 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_cneo_methods, ONLY: vh_1c_nuc_integrals,&
31 calculate_rhoz_cneo
32 USE qs_cneo_types, ONLY: cneo_potential_type,&
33 rhoz_cneo_type
34 USE qs_environment_types, ONLY: get_qs_env,&
35 qs_environment_type
36 USE qs_grid_atom, ONLY: grid_atom_type
37 USE qs_harmonics_atom, ONLY: get_none0_cg_list,&
38 harmonics_atom_type
39 USE qs_kind_types, ONLY: get_qs_kind,&
40 get_qs_kind_set,&
41 qs_kind_type
42 USE qs_local_rho_types, ONLY: get_local_rho,&
43 local_rho_type,&
44 rhoz_type
45 USE qs_rho0_types, ONLY: get_rho0_mpole,&
46 rho0_atom_type,&
47 rho0_mpole_type
48 USE qs_rho_atom_types, ONLY: get_rho_atom,&
49 rho_atom_coeff,&
50 rho_atom_type
51 USE util, ONLY: get_limit
52 USE virial_types, ONLY: virial_type
53#include "./base/base_uses.f90"
54
55 IMPLICIT NONE
56
57 PRIVATE
58
59 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'hartree_local_methods'
60
61 ! Public Subroutine
62
63 PUBLIC :: init_coulomb_local, calculate_vh_1center, vh_1c_gg_integrals
64 PUBLIC :: gapw_background_potential
65
66CONTAINS
67
68! **************************************************************************************************
69!> \brief ...
70!> \param hartree_local ...
71!> \param natom ...
72! **************************************************************************************************
73 SUBROUTINE init_coulomb_local(hartree_local, natom)
74
75 TYPE(hartree_local_type), POINTER :: hartree_local
76 INTEGER, INTENT(IN) :: natom
77
78 CHARACTER(len=*), PARAMETER :: routinen = 'init_coulomb_local'
79
80 INTEGER :: handle
81 TYPE(ecoul_1center_type), DIMENSION(:), POINTER :: ecoul_1c
82
83 CALL timeset(routinen, handle)
84
85 NULLIFY (ecoul_1c)
86 ! Allocate and Initialize 1-center Potentials and Integrals
87 CALL allocate_ecoul_1center(ecoul_1c, natom)
88 hartree_local%ecoul_1c => ecoul_1c
89
90 CALL timestop(handle)
91
92 END SUBROUTINE init_coulomb_local
93
94! **************************************************************************************************
95!> \brief Calculates Hartree potential for hard and soft densities (including
96!> nuclear charge and compensation charges) using numerical integration
97!> \param vrad_h ...
98!> \param vrad_s ...
99!> \param rrad_h ...
100!> \param rrad_s ...
101!> \param rrad_0 ...
102!> \param rrad_z ...
103!> \param grid_atom ...
104!> \par History
105!> 05.2012 JGH refactoring
106!> \author ??
107! **************************************************************************************************
108 SUBROUTINE calculate_vh_1center(vrad_h, vrad_s, rrad_h, rrad_s, rrad_0, rrad_z, grid_atom)
109
110 REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: vrad_h, vrad_s
111 TYPE(rho_atom_coeff), DIMENSION(:), INTENT(IN) :: rrad_h, rrad_s
112 REAL(dp), DIMENSION(:, :), INTENT(IN) :: rrad_0
113 REAL(dp), DIMENSION(:), INTENT(IN) :: rrad_z
114 TYPE(grid_atom_type), POINTER :: grid_atom
115
116 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_Vh_1center'
117
118 INTEGER :: handle, ir, iso, ispin, l_ang, &
119 max_s_harm, nchannels, nr, nspins
120 REAL(dp) :: i1_down, i1_up, i2_down, i2_up, prefactor
121 REAL(dp), ALLOCATABLE, DIMENSION(:) :: rho_1, rho_2
122 REAL(dp), DIMENSION(:), POINTER :: wr
123 REAL(dp), DIMENSION(:, :), POINTER :: oor2l, r2l
124
125 CALL timeset(routinen, handle)
126
127 nr = grid_atom%nr
128 max_s_harm = SIZE(vrad_h, 2)
129 nspins = SIZE(rrad_h, 1)
130 nchannels = SIZE(rrad_0, 2)
131
132 r2l => grid_atom%rad2l
133 oor2l => grid_atom%oorad2l
134 wr => grid_atom%wr
135
136 ALLOCATE (rho_1(nr), rho_2(nr))
137
138 DO iso = 1, max_s_harm
139 rho_1(:) = 0.0_dp
140 rho_2(:) = 0.0_dp
141 IF (iso == 1) rho_1(:) = rrad_z(:)
142 IF (iso <= nchannels) rho_2(:) = rrad_0(:, iso)
143 DO ispin = 1, nspins
144 rho_1(:) = rho_1(:) + rrad_h(ispin)%r_coef(:, iso)
145 rho_2(:) = rho_2(:) + 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(:) = rho_1(:)*wr(:)
152 rho_2(:) = rho_2(:)*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)
160 i2_up = r2l(nr, l_ang)*rho_2(nr)
161
162 DO ir = nr - 1, 1, -1
163 i1_down = i1_down + oor2l(ir, l_ang + 1)*rho_1(ir)
164 i2_down = i2_down + oor2l(ir, l_ang + 1)*rho_2(ir)
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)
174 i1_down = i1_down - oor2l(ir, l_ang + 1)*rho_1(ir)
175 i2_up = i2_up + r2l(ir, l_ang)*rho_2(ir)
176 i2_down = i2_down - oor2l(ir, l_ang + 1)*rho_2(ir)
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, &
225 llmax_nuc, lmax0, lmax0_2nd, lmax_0, m1, max_iso, max_iso_not0, max_iso_not0_nuc, &
226 max_s_harm, max_s_harm_nuc, maxl, maxl_nuc, maxso, maxso_nuc, mepos, n1, nat, nchan_0, &
227 nkind, nr, nset, nset_nuc, nsotot, nsotot_nuc, nspins, num_pe
228 INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list
229 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list
230 INTEGER, DIMENSION(:), POINTER :: atom_list, lmax, lmax_nuc, lmin, &
231 lmin_nuc, npgf, npgf_nuc
232 LOGICAL :: cneo, core_charge, l_2nd_local_rho, &
233 my_core_2nd, my_periodic, paw_atom
234 REAL(dp) :: back_ch, background_energy, &
235 background_nuc, ecoul_1_z_cneo
236 REAL(dp), ALLOCATABLE, DIMENSION(:) :: background_potential, gexp, sqrtwr
237 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: avh1b_00, avh1b_00_nuc, avh1b_hh, &
238 avh1b_hh_nuc, avh1b_ss, avh1b_ss_nuc, &
239 g0_h_w
240 REAL(dp), DIMENSION(:), POINTER :: rrad_z, vrrad_z
241 REAL(dp), DIMENSION(:, :), POINTER :: g0_h, g0_h_2nd, gsph, gsph_nuc, rrad_0, &
242 vh1_h, vh1_s, vrrad_0, zet, zet_nuc
243 REAL(dp), DIMENSION(:, :, :), POINTER :: my_cg, qlm_gg, qlm_gg_2nd
244 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
245 TYPE(cell_type), POINTER :: cell
246 TYPE(cneo_potential_type), POINTER :: cneo_potential
247 TYPE(dft_control_type), POINTER :: dft_control
248 TYPE(grid_atom_type), POINTER :: grid_atom
249 TYPE(gto_basis_set_type), POINTER :: basis_1c, nuc_basis
250 TYPE(harmonics_atom_type), POINTER :: harmonics
251 TYPE(pw_env_type), POINTER :: pw_env
252 TYPE(pw_poisson_type), POINTER :: poisson_env
253 TYPE(qs_charges_type), POINTER :: qs_charges
254 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
255 TYPE(rho0_atom_type), DIMENSION(:), POINTER :: rho0_atom_set, rho0_atom_set_2nd
256 TYPE(rho0_mpole_type), POINTER :: rho0_mpole, rho0_mpole_2nd
257 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set, rho_atom_set_2nd
258 TYPE(rho_atom_type), POINTER :: rho_atom
259 TYPE(rhoz_cneo_type), DIMENSION(:), POINTER :: rhoz_cneo_set
260 TYPE(rhoz_cneo_type), POINTER :: rhoz_cneo
261 TYPE(rhoz_type), DIMENSION(:), POINTER :: rhoz_set, rhoz_set_2nd
262 TYPE(virial_type), POINTER :: virial
263
264 CALL timeset(routinen, handle)
265
266 NULLIFY (cell, dft_control, poisson_env, pw_env, qs_charges)
267 NULLIFY (atomic_kind_set, qs_kind_set, rho_atom_set, rho0_atom_set)
268 NULLIFY (rho0_mpole, rhoz_set)
269 NULLIFY (atom_list, grid_atom, harmonics)
270 NULLIFY (basis_1c, lmin, lmax, npgf, zet)
271 NULLIFY (gsph)
272 NULLIFY (rhoz_cneo, rhoz_cneo_set)
273
274 CALL get_qs_env(qs_env=qs_env, &
275 cell=cell, dft_control=dft_control, &
276 atomic_kind_set=atomic_kind_set, &
277 qs_kind_set=qs_kind_set, &
278 pw_env=pw_env, qs_charges=qs_charges, virial=virial)
279
280 CALL pw_env_get(pw_env, poisson_env=poisson_env)
281 my_periodic = (poisson_env%method == pw_poisson_periodic)
282
283 ! The compensating uniform charge is fixed by the electron number, not
284 ! by the approximate hard/soft quadrature. A planar countercharge makes
285 ! the full Poisson problem neutral. At fixed charge the correction is
286 ! linear in the density, so it contributes no TDDFT Hartree kernel.
287 back_ch = 0.0_dp
288 IF (my_periodic .AND. .NOT. dft_control%do_pcc .AND. .NOT. tddft) back_ch = dft_control%charge
289
290 ! rhoz_set is not accessed in TDDFT
291 CALL get_local_rho(local_rho_set, rho_atom_set, rho0_atom_set, rho0_mpole, rhoz_set, &
292 rhoz_cneo_set) ! for integral space
293
294 ! for forces we need a second local_rho_set
295 l_2nd_local_rho = .false.
296 IF (PRESENT(local_rho_set_2nd)) THEN ! for potential
297 l_2nd_local_rho = .true.
298 NULLIFY (rho_atom_set_2nd, rho0_atom_set_2nd, rhoz_set_2nd) ! for potential
299 CALL get_local_rho(local_rho_set_2nd, rho_atom_set_2nd, rho0_atom_set_2nd, rho0_mpole_2nd, rhoz_set=rhoz_set_2nd)
300 END IF
301
302 nkind = SIZE(atomic_kind_set, 1)
303 nspins = dft_control%nspins
304
305 core_charge = .NOT. tddft ! for forces mixed version
306 my_core_2nd = .true.
307 IF (PRESENT(core_2nd)) my_core_2nd = .NOT. core_2nd ! if my_core_2nd true, include core charge
308
309 ! The aim of the following code was to return immediately if the subroutine
310 ! was called for triplet excited states in spin-restricted case. This check
311 ! is also performed before invocation of this subroutine. It should be save
312 ! to remove the optional argument 'do_triplet' from the subroutine interface.
313 !IF (tddft) THEN
314 ! CPASSERT(PRESENT(do_triplet))
315 ! IF (nspins == 1 .AND. do_triplet) RETURN
316 !END IF
317
318 CALL get_qs_kind_set(qs_kind_set, maxg_iso_not0=max_iso)
319 CALL get_rho0_mpole(rho0_mpole=rho0_mpole, lmax_0=lmax_0)
320
321 ! Put to 0 the local hartree energy contribution from 1 center integrals
322 energy_hartree_1c = 0.0_dp
323 ! Restore total quantum nuclear density to zero
324 qs_charges%total_rho1_hard_nuc = 0.0_dp
325 rho0_mpole%tot_rhoz_cneo_s = 0.0_dp
326
327 ! Here starts the loop over all the atoms
328 DO ikind = 1, nkind
329
330 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
331 NULLIFY (cneo_potential)
332 CALL get_qs_kind(qs_kind_set(ikind), &
333 grid_atom=grid_atom, &
334 harmonics=harmonics, ngrid_rad=nr, &
335 max_iso_not0=max_iso_not0, paw_atom=paw_atom, &
336 cneo_potential=cneo_potential)
337 CALL get_qs_kind(qs_kind_set(ikind), &
338 basis_set=basis_1c, basis_type="GAPW_1C")
339
340 cneo = ASSOCIATED(cneo_potential)
341 IF (cneo .AND. tddft) THEN
342 cpabort("Electronic TDDFT with CNEO quantum nuclei is not implemented.")
343 END IF
344
345 NULLIFY (nuc_basis)
346 max_iso_not0_nuc = 0
347 IF (cneo) THEN
348 cpassert(paw_atom)
349 CALL get_qs_kind(qs_kind_set(ikind), &
350 basis_set=nuc_basis, basis_type="NUC")
351 max_iso_not0_nuc = cneo_potential%harmonics%max_iso_not0
352 END IF
353
354 IF (paw_atom) THEN
355!=========== PAW ===============
356 CALL get_gto_basis_set(gto_basis_set=basis_1c, lmax=lmax, lmin=lmin, &
357 maxso=maxso, npgf=npgf, maxl=maxl, &
358 nset=nset, zet=zet)
359
360 max_s_harm = harmonics%max_s_harm
361 llmax = harmonics%llmax
362
363 nsotot = maxso*nset
364 ALLOCATE (gsph(nr, nsotot))
365 ALLOCATE (gexp(nr))
366 ALLOCATE (sqrtwr(nr), g0_h_w(nr, 0:lmax_0), background_potential(nr))
367 background_potential(:) = gapw_background_potential(back_ch, cell%deth, grid_atom%rad2(:))
368
369 ALLOCATE (avh1b_hh(nsotot, nsotot))
370 ALLOCATE (avh1b_ss(nsotot, nsotot))
371 ALLOCATE (avh1b_00(nsotot, nsotot))
372
373 NULLIFY (qlm_gg, g0_h)
374 CALL get_rho0_mpole(rho0_mpole=rho0_mpole, ikind=ikind, &
375 l0_ikind=lmax0, &
376 qlm_gg=qlm_gg, g0_h=g0_h) ! Qlm_gg of density
377
378 IF (PRESENT(local_rho_set_2nd)) THEN ! for potential
379 NULLIFY (qlm_gg_2nd, g0_h_2nd)
380 CALL get_rho0_mpole(rho0_mpole=rho0_mpole_2nd, ikind=ikind, &
381 l0_ikind=lmax0_2nd, &
382 qlm_gg=qlm_gg_2nd, g0_h=g0_h_2nd) ! Qlm_gg of density
383 END IF
384 nchan_0 = nsoset(lmax0)
385
386 IF (nchan_0 > max(max_iso_not0, max_iso_not0_nuc)) THEN
387 cpabort("channels for rho0 > # max of spherical harmonics")
388 END IF
389
390 NULLIFY (vh1_h, vh1_s)
391 ALLOCATE (vh1_h(nr, max_iso_not0))
392 ALLOCATE (vh1_s(nr, max(max_iso_not0, max_iso_not0_nuc, nchan_0)))
393
394 NULLIFY (lmax_nuc, lmin_nuc, npgf_nuc, zet_nuc, gsph_nuc)
395 maxso_nuc = 0
396 maxl_nuc = -1
397 nset_nuc = 0
398 max_s_harm_nuc = 0
399 llmax_nuc = -1
400 nsotot_nuc = 0
401 IF (cneo) THEN
402 CALL get_gto_basis_set(gto_basis_set=nuc_basis, lmax=lmax_nuc, &
403 lmin=lmin_nuc, maxso=maxso_nuc, npgf=npgf_nuc, &
404 maxl=maxl_nuc, nset=nset_nuc, zet=zet_nuc)
405
406 max_s_harm_nuc = cneo_potential%harmonics%max_s_harm
407 llmax_nuc = cneo_potential%harmonics%llmax
408 nsotot_nuc = maxso_nuc*nset_nuc
409 ALLOCATE (gsph_nuc(nr, nsotot_nuc))
410 gsph_nuc = 0.0_dp
411 ALLOCATE (avh1b_hh_nuc(nsotot_nuc, nsotot_nuc))
412 ALLOCATE (avh1b_ss_nuc(nsotot_nuc, nsotot_nuc))
413 ALLOCATE (avh1b_00_nuc(nsotot_nuc, nsotot_nuc))
414 END IF
415
416 ALLOCATE (cg_list(2, nsoset(max(maxl, maxl_nuc))**2, &
417 max(max_s_harm, max_s_harm_nuc)), &
418 cg_n_list(max(max_s_harm, max_s_harm_nuc)))
419
420 NULLIFY (rrad_z, my_cg)
421 my_cg => harmonics%my_CG
422
423 ! set to zero temporary arrays
424 sqrtwr = 0.0_dp
425 g0_h_w = 0.0_dp
426 gexp = 0.0_dp
427 gsph = 0.0_dp
428
429 sqrtwr(1:nr) = sqrt(grid_atom%wr(1:nr))
430 DO l_ang = 0, lmax0
431 g0_h_w(1:nr, l_ang) = g0_h(1:nr, l_ang)*grid_atom%wr(1:nr)
432 END DO
433
434 m1 = 0
435 DO iset1 = 1, nset
436 n1 = nsoset(lmax(iset1))
437 DO ipgf1 = 1, npgf(iset1)
438 gexp(1:nr) = exp(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))*sqrtwr(1:nr)
439 DO is1 = nsoset(lmin(iset1) - 1) + 1, nsoset(lmax(iset1))
440 iso = is1 + (ipgf1 - 1)*n1 + m1
441 l_ang = indso(1, is1)
442 gsph(1:nr, iso) = grid_atom%rad2l(1:nr, l_ang)*gexp(1:nr)
443 END DO ! is1
444 END DO ! ipgf1
445 m1 = m1 + maxso
446 END DO ! iset1
447
448 IF (cneo) THEN
449 ! initialize nuclear pmat, cpc and e_core to zero
450 DO iat = 1, nat
451 iatom = atom_list(iat)
452 rhoz_cneo_set(iatom)%pmat = 0.0_dp
453 rhoz_cneo_set(iatom)%cpc_h = 0.0_dp
454 rhoz_cneo_set(iatom)%cpc_s = 0.0_dp
455 rhoz_cneo_set(iatom)%e_core = 0.0_dp
456 rhoz_cneo_set(iatom)%ready = .false.
457 END DO
458
459 ! calculate nuclear gsph
460 m1 = 0
461 DO iset1 = 1, nset_nuc
462 n1 = nsoset(lmax_nuc(iset1))
463 DO ipgf1 = 1, npgf_nuc(iset1)
464 gexp(1:nr) = exp(-zet_nuc(ipgf1, iset1)*grid_atom%rad2(1:nr))*sqrtwr(1:nr)
465 DO is1 = nsoset(lmin_nuc(iset1) - 1) + 1, nsoset(lmax_nuc(iset1))
466 iso = is1 + (ipgf1 - 1)*n1 + m1
467 l_ang = indso(1, is1)
468 gsph_nuc(1:nr, iso) = cneo_potential%rad2l(1:nr, l_ang)*gexp(1:nr)
469 END DO ! is1
470 END DO ! ipgf1
471 m1 = m1 + maxso_nuc
472 END DO ! iset1
473 END IF
474
475 ! Distribute the atoms of this kind
476 num_pe = para_env%num_pe
477 mepos = para_env%mepos
478 bo = get_limit(nat, num_pe, mepos)
479
480 DO iat = bo(1), bo(2) !1,nat
481 iatom = atom_list(iat)
482 rho_atom => rho_atom_set(iatom)
483
484 NULLIFY (rrad_z, vrrad_z, rrad_0, vrrad_0)
485 IF (core_charge .AND. .NOT. cneo) THEN
486 rrad_z => rhoz_set(ikind)%r_coef ! for density
487 END IF
488 IF (my_core_2nd .AND. .NOT. cneo) THEN
489 IF (l_2nd_local_rho) THEN
490 vrrad_z => rhoz_set_2nd(ikind)%vr_coef ! for potential
491 ELSE
492 vrrad_z => rhoz_set(ikind)%vr_coef ! for potential
493 END IF
494 END IF
495 rrad_0 => rho0_atom_set(iatom)%rho0_rad_h%r_coef ! for density
496 vrrad_0 => rho0_atom_set(iatom)%vrho0_rad_h%r_coef
497 IF (l_2nd_local_rho) THEN
498 rho_atom => rho_atom_set_2nd(iatom)
499 vrrad_0 => rho0_atom_set_2nd(iatom)%vrho0_rad_h%r_coef ! for potential
500 END IF
501 CALL vh_1c_atom_potential(rho_atom, vrrad_0, &
502 grid_atom, my_core_2nd .AND. .NOT. cneo, & ! core charge for potential (2nd)
503 vrrad_z, vh1_h, vh1_s, &
504 nchan_0, nspins, max_iso_not0, background_potential)
505
506 IF (l_2nd_local_rho) rho_atom => rho_atom_set(iatom) ! rho_atom for density
507
508 ecoul_1_z_cneo = 0.0_dp
509 background_nuc = 0.0_dp
510 IF (cneo) THEN
511 rhoz_cneo => rhoz_cneo_set(iatom)
512 ! Add the soft tail of nuclear Hartree potential to total Vh1_s first.
513 ! vrho already contains the -zeff factor.
514 DO iso = 1, max_iso_not0_nuc
515 vh1_s(:, iso) = vh1_s(:, iso) + rhoz_cneo%vrho_rad_s(:, iso)
516 END DO
517
518 ! Build nuclear 1c integrals according to Vh1_h from electronic density only
519 ! and Vh1_s from electron density and last step (or initial guess) nuclear density
520 ! Vh1_h_nuc = -Z*Vh1_h, Vh1_s_nuc = -Z*Vh1_s
521 CALL vh_1c_nuc_integrals(rhoz_cneo, cneo_potential%zeff, &
522 avh1b_hh_nuc, avh1b_ss_nuc, avh1b_00_nuc, vh1_h, vh1_s, &
523 max_iso_not0, max_iso_not0_nuc, &
524 max_s_harm_nuc, llmax_nuc, cg_list, cg_n_list, &
525 nset_nuc, npgf_nuc, lmin_nuc, lmax_nuc, nsotot_nuc, maxso_nuc, &
526 nchan_0, gsph_nuc, g0_h_w, cneo_potential%harmonics%my_CG, &
527 cneo_potential%Qlm_gg)
528
529 ! Solve the nuclear 1c problem
530 CALL calculate_rhoz_cneo(rhoz_cneo, cneo_potential, cg_list, cg_n_list, nset_nuc, &
531 npgf_nuc, lmin_nuc, lmax_nuc, maxl_nuc, maxso_nuc)
532 ! slm_int(iso=1) = sqrt(4*Pi), slm_int(iso>1) = 0, without using Lebedev grid
533 ! when printing, nuclear density is positive, thus needing the minus sign
534 qs_charges%total_rho1_hard_nuc = qs_charges%total_rho1_hard_nuc - sqrt(fourpi) &
535 *sum(rhoz_cneo%rho_rad_h(:, 1)*grid_atom%wr(:))
536 rho0_mpole%tot_rhoz_cneo_s = rho0_mpole%tot_rhoz_cneo_s - sqrt(fourpi) &
537 *sum(rhoz_cneo%rho_rad_s(:, 1)*grid_atom%wr(:))
538
539 ! Calculate the contributions to Ecoul coming from Vh1_h*rhoz_cneo_h and Vh1_s*rhoz_cneo_s.
540 ! Self-interaction of the quantum nucleus is already removed in ecoul_1_z_cneo.
541 ! rho already contains the -zeff factor.
542 DO iso = 1, min(max_iso_not0, max_iso_not0_nuc)
543 ecoul_1_z_cneo = ecoul_1_z_cneo + 0.5_dp* &
544 (sum(vh1_h(:, iso)*rhoz_cneo%rho_rad_h(:, iso)*grid_atom%wr(:)) &
545 - sum(vh1_s(:, iso)*rhoz_cneo%rho_rad_s(:, iso)*grid_atom%wr(:)))
546 END DO
547 DO iso = max_iso_not0 + 1, max_iso_not0_nuc
548 ecoul_1_z_cneo = ecoul_1_z_cneo - 0.5_dp* &
549 sum(vh1_s(:, iso)*rhoz_cneo%rho_rad_s(:, iso)*grid_atom%wr(:))
550 END DO
551
552 background_nuc = sum(background_potential*(rhoz_cneo%rho_rad_h(:, 1) - &
553 rhoz_cneo%rho_rad_s(:, 1))*grid_atom%wr(:))
554 ! The background interaction is linear, whereas the generic
555 ! Coulomb contraction above includes only one half of it.
556 ecoul_1_z_cneo = ecoul_1_z_cneo + 0.5_dp*background_nuc
557
558 ! Add nuclear Hartree potential to total Vh1_h after solving the nuclear 1c problem
559 ! to avoid nuclear self-interaction.
560 ! vrho already contains the -zeff factor.
561 ! Here the min of two max_iso_not0's is chosen, because even when the nuclear one
562 ! is larger, it is meaningless to let Vh1_h have higher angular momentum components
563 ! as Vh1_h now is only used for the electronic part.
564 DO iso = 1, min(max_iso_not0, max_iso_not0_nuc)
565 vh1_h(:, iso) = vh1_h(:, iso) + rhoz_cneo%vrho_rad_h(:, iso)
566 END DO
567 END IF
568
569 CALL vh_1c_atom_energy(energy_hartree_1c, ecoul_1c, rho_atom, rrad_0, &
570 grid_atom, iatom, core_charge .AND. .NOT. cneo, & ! core charge for density
571 rrad_z, vh1_h, vh1_s, nchan_0, nspins, max_iso_not0, &
572 background_potential, background_energy)
573
574 IF (l_2nd_local_rho) rho_atom => rho_atom_set_2nd(iatom) ! rho_atom for potential (2nd)
575
576 IF (cneo) THEN
577 CALL set_ecoul_1c(ecoul_1c, iatom, ecoul_1_z=ecoul_1c(iatom)%ecoul_1_z + ecoul_1_z_cneo)
578 energy_hartree_1c = energy_hartree_1c + ecoul_1_z_cneo
579 END IF
580
581 ! Besides the OCE coefficient derivatives, E_background has an
582 ! explicit 1/volume dependence. At fixed local radial functions
583 ! its contribution to -dE/d(strain) is E_background times I.
584 IF (virial%pv_calculate .AND. .NOT. tddft) THEN
585 DO iso = 1, 3
586 virial%pv_virial(iso, iso) = virial%pv_virial(iso, iso) + background_energy + background_nuc
587 virial%pv_ehartree(iso, iso) = virial%pv_ehartree(iso, iso) + background_energy + background_nuc
588 END DO
589 END IF
590
591 CALL vh_1c_atom_integrals(rho_atom, & ! results (int_local_h and int_local_s) written on rho_atom_2nd
592 ! int_local_h and int_local_s are used in update_ks_atom
593 ! on int_local_h mixed core / non-core
594 avh1b_hh, avh1b_ss, avh1b_00, vh1_h, vh1_s, max_iso_not0, &
595 max_s_harm, llmax, cg_list, cg_n_list, &
596 nset, npgf, lmin, lmax, nsotot, maxso, nspins, nchan_0, gsph, &
597 g0_h_w, my_cg, qlm_gg) ! Qlm_gg for density from local_rho_set
598
599 END DO ! iat
600
601 DEALLOCATE (avh1b_hh)
602 DEALLOCATE (avh1b_ss)
603 DEALLOCATE (avh1b_00)
604 IF (ALLOCATED(avh1b_hh_nuc)) DEALLOCATE (avh1b_hh_nuc)
605 IF (ALLOCATED(avh1b_ss_nuc)) DEALLOCATE (avh1b_ss_nuc)
606 IF (ALLOCATED(avh1b_00_nuc)) DEALLOCATE (avh1b_00_nuc)
607 DEALLOCATE (vh1_h, vh1_s)
608 DEALLOCATE (cg_list, cg_n_list)
609 DEALLOCATE (gsph)
610 IF (ASSOCIATED(gsph_nuc)) DEALLOCATE (gsph_nuc)
611 DEALLOCATE (gexp)
612 DEALLOCATE (sqrtwr, g0_h_w, background_potential)
613
614 IF (cneo) THEN
615 ! broadcast nuclear pmat, cpc and e_core
616 DO iat = 1, nat
617 iatom = atom_list(iat)
618 CALL para_env%sum(rhoz_cneo_set(iatom)%pmat)
619 CALL para_env%sum(rhoz_cneo_set(iatom)%cpc_h)
620 CALL para_env%sum(rhoz_cneo_set(iatom)%cpc_s)
621 CALL para_env%sum(rhoz_cneo_set(iatom)%e_core)
622 rhoz_cneo_set(iatom)%ready = .true.
623 END DO
624 END IF
625 ELSE
626!=========== NO PAW ===============
627! This term is taken care of using the core density as in GPW
628 cycle
629 END IF ! paw
630 END DO ! ikind
631
632 CALL para_env%sum(energy_hartree_1c)
633 CALL para_env%sum(qs_charges%total_rho1_hard_nuc)
634 CALL para_env%sum(rho0_mpole%tot_rhoz_cneo_s)
635
636 CALL timestop(handle)
637
638 END SUBROUTINE vh_1c_gg_integrals
639
640! **************************************************************************************************
641
642! **************************************************************************************************
643!> \brief ...
644!> \param rho_atom ...
645!> \param vrrad_0 ...
646!> \param grid_atom ...
647!> \param core_charge ...
648!> \param vrrad_z ...
649!> \param Vh1_h ...
650!> \param Vh1_s ...
651!> \param nchan_0 ...
652!> \param nspins ...
653!> \param max_iso_not0 ...
654!> \param background_potential Unweighted l=0 potential of the uniform background.
655! **************************************************************************************************
656 SUBROUTINE vh_1c_atom_potential(rho_atom, vrrad_0, &
657 grid_atom, core_charge, vrrad_z, Vh1_h, Vh1_s, &
658 nchan_0, nspins, max_iso_not0, background_potential)
659
660 TYPE(rho_atom_type), POINTER :: rho_atom
661 REAL(dp), DIMENSION(:, :), POINTER :: vrrad_0
662 TYPE(grid_atom_type), POINTER :: grid_atom
663 LOGICAL, INTENT(IN) :: core_charge
664 REAL(dp), DIMENSION(:), POINTER :: vrrad_z
665 REAL(dp), DIMENSION(:, :), POINTER :: vh1_h, vh1_s
666 INTEGER, INTENT(IN) :: nchan_0, nspins, max_iso_not0
667 REAL(dp), INTENT(IN) :: background_potential(:)
668
669 INTEGER :: iso, ispin, nr
670 TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: vr_h, vr_s
671
672 nr = grid_atom%nr
673
674 NULLIFY (vr_h, vr_s)
675 CALL get_rho_atom(rho_atom=rho_atom, vrho_rad_h=vr_h, vrho_rad_s=vr_s)
676
677 vh1_h = 0.0_dp
678 vh1_s = 0.0_dp
679
680 IF (core_charge) vh1_h(:, 1) = vrrad_z(:)
681
682 DO iso = 1, nchan_0
683 vh1_s(:, iso) = vrrad_0(:, iso)
684 END DO
685
686 DO ispin = 1, nspins
687 DO iso = 1, max_iso_not0
688 vh1_h(:, iso) = vh1_h(:, iso) + vr_h(ispin)%r_coef(:, iso)
689 vh1_s(:, iso) = vh1_s(:, iso) + vr_s(ispin)%r_coef(:, iso)
690 END DO
691 END DO
692
693 vh1_h(:, 1) = vh1_h(:, 1) + background_potential
694 vh1_s(:, 1) = vh1_s(:, 1) + background_potential
695
696 END SUBROUTINE vh_1c_atom_potential
697
698! **************************************************************************************************
699
700! **************************************************************************************************
701!> \brief ...
702!> \param energy_hartree_1c ...
703!> \param ecoul_1c ...
704!> \param rho_atom ...
705!> \param rrad_0 ...
706!> \param grid_atom ...
707!> \param iatom ...
708!> \param core_charge ...
709!> \param rrad_z ...
710!> \param Vh1_h ...
711!> \param Vh1_s ...
712!> \param nchan_0 ...
713!> \param nspins ...
714!> \param max_iso_not0 ...
715!> \param background_potential Unweighted l=0 potential of the uniform background.
716!> \param background_energy Full linear background interaction for the volume derivative.
717! **************************************************************************************************
718 SUBROUTINE vh_1c_atom_energy(energy_hartree_1c, ecoul_1c, rho_atom, rrad_0, &
719 grid_atom, iatom, core_charge, rrad_z, Vh1_h, Vh1_s, &
720 nchan_0, nspins, max_iso_not0, background_potential, background_energy)
721
722 REAL(dp), INTENT(INOUT) :: energy_hartree_1c
723 TYPE(ecoul_1center_type), DIMENSION(:), POINTER :: ecoul_1c
724 TYPE(rho_atom_type), POINTER :: rho_atom
725 REAL(dp), DIMENSION(:, :), POINTER :: rrad_0
726 TYPE(grid_atom_type), POINTER :: grid_atom
727 INTEGER, INTENT(IN) :: iatom
728 LOGICAL, INTENT(IN) :: core_charge
729 REAL(dp), DIMENSION(:), POINTER :: rrad_z
730 REAL(dp), DIMENSION(:, :), POINTER :: vh1_h, vh1_s
731 INTEGER, INTENT(IN) :: nchan_0, nspins, max_iso_not0
732 REAL(dp), INTENT(IN) :: background_potential(:)
733 REAL(dp), INTENT(OUT) :: background_energy
734
735 INTEGER :: iso, ispin, nr
736 REAL(dp) :: bg_0, bg_h, bg_s, bg_z, ecoul_1_0, &
737 ecoul_1_h, ecoul_1_s, ecoul_1_z
738 TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: r_h, r_s
739
740 nr = grid_atom%nr
741
742 NULLIFY (r_h, r_s)
743 CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, rho_rad_s=r_s)
744
745 ! Calculate the contributions to Ecoul coming from Vh1_h*rhoz
746 ecoul_1_z = 0.0_dp
747 bg_z = 0.0_dp
748 IF (core_charge) THEN
749 ecoul_1_z = 0.5_dp*sum(vh1_h(:, 1)*rrad_z(:)*grid_atom%wr(:))
750 bg_z = sum(background_potential*rrad_z(:)*grid_atom%wr(:))
751 END IF
752
753 ! Calculate the contributions to Ecoul coming from Vh1_s*rho0
754 ecoul_1_0 = 0.0_dp
755 DO iso = 1, nchan_0
756 ecoul_1_0 = ecoul_1_0 + 0.5_dp*sum(vh1_s(:, iso)*rrad_0(:, iso)*grid_atom%wr(:))
757 END DO
758 bg_0 = sum(background_potential*rrad_0(:, 1)*grid_atom%wr(:))
759
760 ! Calculate the contributions to Ecoul coming from Vh1_h*rho1_h and Vh1_s*rho1_s
761 ecoul_1_s = 0.0_dp
762 ecoul_1_h = 0.0_dp
763 bg_h = 0.0_dp
764 bg_s = 0.0_dp
765 DO ispin = 1, nspins
766 bg_h = bg_h + sum(background_potential*r_h(ispin)%r_coef(:, 1)*grid_atom%wr(:))
767 bg_s = bg_s + sum(background_potential*r_s(ispin)%r_coef(:, 1)*grid_atom%wr(:))
768 DO iso = 1, max_iso_not0
769 ecoul_1_s = ecoul_1_s + 0.5_dp*sum(vh1_s(:, iso)*r_s(ispin)%r_coef(:, iso)*grid_atom%wr(:))
770 ecoul_1_h = ecoul_1_h + 0.5_dp*sum(vh1_h(:, iso)*r_h(ispin)%r_coef(:, iso)*grid_atom%wr(:))
771 END DO
772 END DO
773
774 ! Unlike the quadratic Hartree self-energy, the interaction with a fixed
775 ! uniform background has no factor 1/2. Its potential is already in Vh1.
776 ecoul_1_z = ecoul_1_z + 0.5_dp*bg_z
777 ecoul_1_0 = ecoul_1_0 + 0.5_dp*bg_0
778 ecoul_1_h = ecoul_1_h + 0.5_dp*bg_h
779 ecoul_1_s = ecoul_1_s + 0.5_dp*bg_s
780 background_energy = bg_z - bg_0 + bg_h - bg_s
781
782 CALL set_ecoul_1c(ecoul_1c, iatom, ecoul_1_z=ecoul_1_z, ecoul_1_0=ecoul_1_0)
783 CALL set_ecoul_1c(ecoul_1c=ecoul_1c, iatom=iatom, ecoul_1_h=ecoul_1_h, ecoul_1_s=ecoul_1_s)
784
785 energy_hartree_1c = energy_hartree_1c + ecoul_1_z - ecoul_1_0
786 energy_hartree_1c = energy_hartree_1c + ecoul_1_h - ecoul_1_s
787
788 END SUBROUTINE vh_1c_atom_energy
789
790! **************************************************************************************************
791!> \brief l=0 potential for the periodic mean-potential correction.
792!> For a neutral localized difference delta_rho between hard and soft
793!> total charges, integral V(delta_rho) = -2*pi/3 integral r^2 delta_rho.
794!> CP2K stores electron-minus-nuclear densities, opposite to net_charge.
795!> Radial quadrature weights belong in the contraction, not in this potential.
796!> \param net_charge Physical cell charge, positive for missing electrons.
797!> \param volume Cell volume in bohr^3.
798!> \param rad2 Squared atom-centered radial coordinates.
799!> \return Real spherical-harmonic l=0 coefficient of the potential.
800! **************************************************************************************************
801 PURE FUNCTION gapw_background_potential(net_charge, volume, rad2) RESULT(potential)
802 REAL(dp), INTENT(IN) :: net_charge, volume, rad2(:)
803 REAL(dp) :: potential(size(rad2))
804
805 potential = -(2.0_dp*pi/3.0_dp)*sqrt(fourpi)*net_charge/volume*rad2
806 END FUNCTION gapw_background_potential
807
808!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
809
810! **************************************************************************************************
811!> \brief ...
812!> \param rho_atom ...
813!> \param aVh1b_hh ...
814!> \param aVh1b_ss ...
815!> \param aVh1b_00 ...
816!> \param Vh1_h ...
817!> \param Vh1_s ...
818!> \param max_iso_not0 ...
819!> \param max_s_harm ...
820!> \param llmax ...
821!> \param cg_list ...
822!> \param cg_n_list ...
823!> \param nset ...
824!> \param npgf ...
825!> \param lmin ...
826!> \param lmax ...
827!> \param nsotot ...
828!> \param maxso ...
829!> \param nspins ...
830!> \param nchan_0 ...
831!> \param gsph ...
832!> \param g0_h_w ...
833!> \param my_CG ...
834!> \param Qlm_gg ...
835! **************************************************************************************************
836 SUBROUTINE vh_1c_atom_integrals(rho_atom, &
837 aVh1b_hh, aVh1b_ss, aVh1b_00, Vh1_h, Vh1_s, max_iso_not0, &
838 max_s_harm, llmax, cg_list, cg_n_list, &
839 nset, npgf, lmin, lmax, nsotot, maxso, nspins, nchan_0, gsph, &
840 g0_h_w, my_CG, Qlm_gg)
841
842 TYPE(rho_atom_type), POINTER :: rho_atom
843 REAL(dp), DIMENSION(:, :) :: avh1b_hh, avh1b_ss, avh1b_00
844 REAL(dp), DIMENSION(:, :), POINTER :: vh1_h, vh1_s
845 INTEGER, INTENT(IN) :: max_iso_not0, max_s_harm, llmax
846 INTEGER, DIMENSION(:, :, :) :: cg_list
847 INTEGER, DIMENSION(:) :: cg_n_list
848 INTEGER, INTENT(IN) :: nset
849 INTEGER, DIMENSION(:), POINTER :: npgf, lmin, lmax
850 INTEGER, INTENT(IN) :: nsotot, maxso, nspins, nchan_0
851 REAL(dp), DIMENSION(:, :), POINTER :: gsph
852 REAL(dp), DIMENSION(:, 0:) :: g0_h_w
853 REAL(dp), DIMENSION(:, :, :), POINTER :: my_cg, qlm_gg
854
855 INTEGER :: icg, ipgf1, ipgf2, ir, is1, is2, iset1, &
856 iset2, iso, iso1, iso2, ispin, l_ang, &
857 m1, m2, max_iso_not0_local, n1, n2, nr
858 REAL(dp) :: gvg_0, gvg_h, gvg_s
859 TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: int_local_h, int_local_s
860
861 NULLIFY (int_local_h, int_local_s)
862 CALL get_rho_atom(rho_atom=rho_atom, &
863 ga_vlocal_gb_h=int_local_h, &
864 ga_vlocal_gb_s=int_local_s)
865
866 ! Calculate the integrals of the potential with 2 primitives
867 avh1b_hh = 0.0_dp
868 avh1b_ss = 0.0_dp
869 avh1b_00 = 0.0_dp
870
871 nr = SIZE(gsph, 1)
872
873 m1 = 0
874 DO iset1 = 1, nset
875 m2 = 0
876 DO iset2 = 1, nset
877 CALL get_none0_cg_list(my_cg, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
878 max_s_harm, llmax, cg_list, cg_n_list, max_iso_not0_local)
879
880 n1 = nsoset(lmax(iset1))
881 DO ipgf1 = 1, npgf(iset1)
882 n2 = nsoset(lmax(iset2))
883 DO ipgf2 = 1, npgf(iset2)
884 ! with contributions to V1_s*rho0
885 DO iso = 1, min(nchan_0, max_iso_not0)
886 l_ang = indso(1, iso)
887 gvg_0 = sum(vh1_s(:, iso)*g0_h_w(:, l_ang))
888 DO icg = 1, cg_n_list(iso)
889 is1 = cg_list(1, icg, iso)
890 is2 = cg_list(2, icg, iso)
891
892 iso1 = is1 + n1*(ipgf1 - 1) + m1
893 iso2 = is2 + n2*(ipgf2 - 1) + m2
894 gvg_h = 0.0_dp
895 gvg_s = 0.0_dp
896
897 DO ir = 1, nr
898 gvg_h = gvg_h + gsph(ir, iso1)*gsph(ir, iso2)*vh1_h(ir, iso)
899 gvg_s = gvg_s + gsph(ir, iso1)*gsph(ir, iso2)*vh1_s(ir, iso)
900 END DO ! ir
901
902 avh1b_hh(iso1, iso2) = avh1b_hh(iso1, iso2) + gvg_h*my_cg(is1, is2, iso)
903 avh1b_ss(iso1, iso2) = avh1b_ss(iso1, iso2) + gvg_s*my_cg(is1, is2, iso)
904 avh1b_00(iso1, iso2) = avh1b_00(iso1, iso2) + gvg_0*qlm_gg(iso1, iso2, iso)
905
906 END DO !icg
907 END DO ! iso
908 ! without contributions to V1_s*rho0
909 DO iso = nchan_0 + 1, max_iso_not0
910 DO icg = 1, cg_n_list(iso)
911 is1 = cg_list(1, icg, iso)
912 is2 = cg_list(2, icg, iso)
913
914 iso1 = is1 + n1*(ipgf1 - 1) + m1
915 iso2 = is2 + n2*(ipgf2 - 1) + m2
916 gvg_h = 0.0_dp
917 gvg_s = 0.0_dp
918
919 DO ir = 1, nr
920 gvg_h = gvg_h + gsph(ir, iso1)*gsph(ir, iso2)*vh1_h(ir, iso)
921 gvg_s = gvg_s + gsph(ir, iso1)*gsph(ir, iso2)*vh1_s(ir, iso)
922 END DO ! ir
923
924 avh1b_hh(iso1, iso2) = avh1b_hh(iso1, iso2) + gvg_h*my_cg(is1, is2, iso)
925 avh1b_ss(iso1, iso2) = avh1b_ss(iso1, iso2) + gvg_s*my_cg(is1, is2, iso)
926
927 END DO !icg
928 END DO ! iso
929 END DO ! ipgf2
930 END DO ! ipgf1
931 m2 = m2 + maxso
932 END DO ! iset2
933 m1 = m1 + maxso
934 END DO !iset1
935 DO ispin = 1, nspins
936 CALL daxpy(nsotot*nsotot, 1.0_dp, avh1b_hh, 1, int_local_h(ispin)%r_coef, 1)
937 CALL daxpy(nsotot*nsotot, 1.0_dp, avh1b_ss, 1, int_local_s(ispin)%r_coef, 1)
938 CALL daxpy(nsotot*nsotot, -1.0_dp, avh1b_00, 1, int_local_h(ispin)%r_coef, 1)
939 CALL daxpy(nsotot*nsotot, -1.0_dp, avh1b_00, 1, int_local_s(ispin)%r_coef, 1)
940 END DO ! ispin
941
942 END SUBROUTINE vh_1c_atom_integrals
943
944!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
945
946END MODULE hartree_local_methods
qs_harmonics_atom::get_none0_cg_list
Definition qs_harmonics_atom.F:54
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:131
basis_set_types
Definition basis_set_types.F:15
basis_set_types::get_gto_basis_set
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum, ccon)
...
Definition basis_set_types.F:643
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::gapw_background_potential
pure real(dp) function, dimension(size(rad2)), public gapw_background_potential(net_charge, volume, rad2)
l=0 potential for the periodic mean-potential correction. For a neutral localized difference delta_rh...
Definition hartree_local_methods.F:802
hartree_local_methods::init_coulomb_local
subroutine, public init_coulomb_local(hartree_local, natom)
...
Definition hartree_local_methods.F:74
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:109
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:27
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:67
qs_charges_types
container for information about total charges on the grids
Definition qs_charges_types.F:14
qs_cneo_methods
A collection of functions used by CNEO-DFT (see J. Chem. Theory Comput. 2025, 21, 16,...
Definition qs_cneo_methods.F:15
qs_cneo_methods::vh_1c_nuc_integrals
subroutine, public vh_1c_nuc_integrals(rhoz_cneo, zeff, avh1b_hh, avh1b_ss, avh1b_00, vh1_h, vh1_s, max_iso_not0_elec, max_iso_not0_nuc, max_s_harm, llmax, cg_list, cg_n_list, nset, npgf, lmin, lmax, nsotot, maxso, nchan_0, gsph, g0_h_w, my_cg, qlm_gg)
Mostly copied from hartree_local_methods::Vh_1c_atom_integrals.
Definition qs_cneo_methods.F:1311
qs_cneo_methods::calculate_rhoz_cneo
subroutine, public calculate_rhoz_cneo(rho, potential, cg_list, cg_n_list, nset, npgf, lmin, lmax, maxl, maxso)
...
Definition qs_cneo_methods.F:1008
qs_cneo_types
Types used by CNEO-DFT (see J. Chem. Theory Comput. 2025, 21, 16, 7865–7877).
Definition qs_cneo_types.F:15
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, mimic, 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_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, 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, matrix_vhxc, 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, xcint_weights, 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, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_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, harris_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, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:544
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, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, 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_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, hund_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, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:483
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, do_mtlr_present, basis_type, total_zeff_corr, npgf_seg, cneo_potential_present, nkind_q, natom_q)
Get attributes of an atomic kind set.
Definition qs_kind_types.F:976
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, rhoz_cneo_set)
...
Definition qs_local_rho_types.F:179
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, rhoz_cneo_s_rs, rhoz_cneo_s_gs)
...
Definition qs_rho0_types.F:449
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:191
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
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
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_control_types::dft_control_type
Definition cp_control_types.F:786
hartree_local_types::ecoul_1center_type
Definition hartree_local_types.F:25
hartree_local_types::hartree_local_type
Definition hartree_local_types.F:34
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:774
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
pw_poisson_types::pw_poisson_type
environment for the poisson solver
Definition pw_poisson_types.F:123
qs_charges_types::qs_charges_type
Container for information about total charges on the grids.
Definition qs_charges_types.F:41
qs_cneo_types::cneo_potential_type
Definition qs_cneo_types.F:53
qs_cneo_types::rhoz_cneo_type
Definition qs_cneo_types.F:29
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:223
qs_grid_atom::grid_atom_type
Definition qs_grid_atom.F:44
qs_harmonics_atom::harmonics_atom_type
Definition qs_harmonics_atom.F:32
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190
qs_local_rho_types::local_rho_type
Definition qs_local_rho_types.F:45
qs_local_rho_types::rhoz_type
Definition qs_local_rho_types.F:37
qs_rho0_types::rho0_atom_type
Definition qs_rho0_types.F:78
qs_rho0_types::rho0_mpole_type
Definition qs_rho0_types.F:52
qs_rho_atom_types::rho_atom_coeff
Definition qs_rho_atom_types.F:19
qs_rho_atom_types::rho_atom_type
Definition qs_rho_atom_types.F:23
virial_types::virial_type
Definition virial_types.F:26