52#include "./base/base_uses.f90"
58 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'hartree_local_methods'
74 INTEGER,
INTENT(IN) :: natom
76 CHARACTER(len=*),
PARAMETER :: routinen =
'init_coulomb_local'
81 CALL timeset(routinen, handle)
86 hartree_local%ecoul_1c => ecoul_1c
108 REAL(
dp),
DIMENSION(:, :),
INTENT(INOUT) :: vrad_h, vrad_s
110 REAL(
dp),
DIMENSION(:, :),
INTENT(IN) :: rrad_0
111 REAL(
dp),
DIMENSION(:),
INTENT(IN) :: rrad_z
114 CHARACTER(len=*),
PARAMETER :: routinen =
'calculate_Vh_1center'
116 INTEGER :: handle, ir, iso, ispin, l_ang, &
117 max_s_harm, nchannels, nr, nspins
118 REAL(
dp) :: i1_down, i1_up, i2_down, i2_up, prefactor
119 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: rho_1, rho_2
120 REAL(
dp),
DIMENSION(:),
POINTER :: wr
121 REAL(
dp),
DIMENSION(:, :),
POINTER :: oor2l, r2l
123 CALL timeset(routinen, handle)
126 max_s_harm =
SIZE(vrad_h, 2)
127 nspins =
SIZE(rrad_h, 1)
128 nchannels =
SIZE(rrad_0, 2)
130 r2l => grid_atom%rad2l
131 oor2l => grid_atom%oorad2l
134 ALLOCATE (rho_1(nr), rho_2(nr))
136 DO iso = 1, max_s_harm
139 IF (iso == 1) rho_1(:) = rrad_z(:)
140 IF (iso <= nchannels) rho_2(:) = rrad_0(:, iso)
142 rho_1(:) = rho_1(:) + rrad_h(ispin)%r_coef(:, iso)
143 rho_2(:) = rho_2(:) + rrad_s(ispin)%r_coef(:, iso)
146 l_ang =
indso(1, iso)
147 prefactor =
fourpi/(2._dp*l_ang + 1._dp)
149 rho_1(:) = rho_1(:)*wr(:)
150 rho_2(:) = rho_2(:)*wr(:)
157 i1_up = r2l(nr, l_ang)*rho_1(nr)
158 i2_up = r2l(nr, l_ang)*rho_2(nr)
160 DO ir = nr - 1, 1, -1
161 i1_down = i1_down + oor2l(ir, l_ang + 1)*rho_1(ir)
162 i2_down = i2_down + oor2l(ir, l_ang + 1)*rho_2(ir)
165 vrad_h(nr, iso) = vrad_h(nr, iso) + prefactor* &
166 (oor2l(nr, l_ang + 1)*i1_up + r2l(nr, l_ang)*i1_down)
167 vrad_s(nr, iso) = vrad_s(nr, iso) + prefactor* &
168 (oor2l(nr, l_ang + 1)*i2_up + r2l(nr, l_ang)*i2_down)
170 DO ir = nr - 1, 1, -1
171 i1_up = i1_up + r2l(ir, l_ang)*rho_1(ir)
172 i1_down = i1_down - oor2l(ir, l_ang + 1)*rho_1(ir)
173 i2_up = i2_up + r2l(ir, l_ang)*rho_2(ir)
174 i2_down = i2_down - oor2l(ir, l_ang + 1)*rho_2(ir)
176 vrad_h(ir, iso) = vrad_h(ir, iso) + prefactor* &
177 (oor2l(ir, l_ang + 1)*i1_up + r2l(ir, l_ang)*i1_down)
178 vrad_s(ir, iso) = vrad_s(ir, iso) + prefactor* &
179 (oor2l(ir, l_ang + 1)*i2_up + r2l(ir, l_ang)*i2_down)
185 DEALLOCATE (rho_1, rho_2)
187 CALL timestop(handle)
208 SUBROUTINE vh_1c_gg_integrals(qs_env, energy_hartree_1c, ecoul_1c, local_rho_set, para_env, tddft, local_rho_set_2nd, &
212 REAL(kind=
dp),
INTENT(out) :: energy_hartree_1c
216 LOGICAL,
INTENT(IN) :: tddft
218 LOGICAL,
INTENT(IN),
OPTIONAL :: core_2nd
220 CHARACTER(LEN=*),
PARAMETER :: routinen =
'Vh_1c_gg_integrals'
222 INTEGER :: bo(2), handle, iat, iatom, ikind, ipgf1, is1, iset1, iso, l_ang, llmax, &
223 llmax_nuc, lmax0, lmax0_2nd, lmax_0, m1, max_iso, max_iso_not0, max_iso_not0_nuc, &
224 max_s_harm, max_s_harm_nuc, maxl, maxl_nuc, maxso, maxso_nuc, mepos, n1, nat, nchan_0, &
225 nkind, nr, nset, nset_nuc, nsotot, nsotot_nuc, nspins, num_pe
226 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: cg_n_list
227 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :) :: cg_list
228 INTEGER,
DIMENSION(:),
POINTER :: atom_list, lmax, lmax_nuc, lmin, &
229 lmin_nuc, npgf, npgf_nuc
230 LOGICAL :: cneo, core_charge, l_2nd_local_rho, &
231 my_core_2nd, my_periodic, paw_atom
232 REAL(
dp) :: back_ch, ecoul_1_z_cneo, factor
233 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: gexp, sqrtwr
234 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: avh1b_00, avh1b_00_nuc, avh1b_hh, &
235 avh1b_hh_nuc, avh1b_ss, avh1b_ss_nuc, &
237 REAL(
dp),
DIMENSION(:),
POINTER :: rrad_z, vrrad_z
238 REAL(
dp),
DIMENSION(:, :),
POINTER :: g0_h, g0_h_2nd, gsph, gsph_nuc, rrad_0, &
239 vh1_h, vh1_s, vrrad_0, zet, zet_nuc
240 REAL(
dp),
DIMENSION(:, :, :),
POINTER :: my_cg, qlm_gg, qlm_gg_2nd
251 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
252 TYPE(
rho0_atom_type),
DIMENSION(:),
POINTER :: rho0_atom_set, rho0_atom_set_2nd
254 TYPE(
rho_atom_type),
DIMENSION(:),
POINTER :: rho_atom_set, rho_atom_set_2nd
258 TYPE(
rhoz_type),
DIMENSION(:),
POINTER :: rhoz_set, rhoz_set_2nd
260 CALL timeset(routinen, handle)
262 NULLIFY (cell, dft_control, poisson_env, pw_env, qs_charges)
263 NULLIFY (atomic_kind_set, qs_kind_set, rho_atom_set, rho0_atom_set)
264 NULLIFY (rho0_mpole, rhoz_set)
265 NULLIFY (atom_list, grid_atom, harmonics)
266 NULLIFY (basis_1c, lmin, lmax, npgf, zet)
268 NULLIFY (rhoz_cneo, rhoz_cneo_set)
271 cell=cell, dft_control=dft_control, &
272 atomic_kind_set=atomic_kind_set, &
273 qs_kind_set=qs_kind_set, &
274 pw_env=pw_env, qs_charges=qs_charges)
276 CALL pw_env_get(pw_env, poisson_env=poisson_env)
279 back_ch = qs_charges%background*cell%deth
282 CALL get_local_rho(local_rho_set, rho_atom_set, rho0_atom_set, rho0_mpole, rhoz_set, &
286 l_2nd_local_rho = .false.
287 IF (
PRESENT(local_rho_set_2nd))
THEN
288 l_2nd_local_rho = .true.
289 NULLIFY (rho_atom_set_2nd, rho0_atom_set_2nd, rhoz_set_2nd)
290 CALL get_local_rho(local_rho_set_2nd, rho_atom_set_2nd, rho0_atom_set_2nd, rho0_mpole_2nd, rhoz_set=rhoz_set_2nd)
293 nkind =
SIZE(atomic_kind_set, 1)
294 nspins = dft_control%nspins
296 core_charge = .NOT. tddft
298 IF (
PRESENT(core_2nd)) my_core_2nd = .NOT. core_2nd
313 energy_hartree_1c = 0.0_dp
315 qs_charges%total_rho1_hard_nuc = 0.0_dp
316 rho0_mpole%tot_rhoz_cneo_s = 0.0_dp
321 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
322 NULLIFY (cneo_potential)
324 grid_atom=grid_atom, &
325 harmonics=harmonics, ngrid_rad=nr, &
326 max_iso_not0=max_iso_not0, paw_atom=paw_atom, &
327 cneo_potential=cneo_potential)
329 basis_set=basis_1c, basis_type=
"GAPW_1C")
331 cneo =
ASSOCIATED(cneo_potential)
332 IF (cneo .AND. tddft)
THEN
333 cpabort(
"Electronic TDDFT with CNEO quantum nuclei is not implemented.")
341 basis_set=nuc_basis, basis_type=
"NUC")
342 max_iso_not0_nuc = cneo_potential%harmonics%max_iso_not0
348 maxso=maxso, npgf=npgf, maxl=maxl, &
351 max_s_harm = harmonics%max_s_harm
352 llmax = harmonics%llmax
355 ALLOCATE (gsph(nr, nsotot))
357 ALLOCATE (sqrtwr(nr), g0_h_w(nr, 0:lmax_0))
359 ALLOCATE (avh1b_hh(nsotot, nsotot))
360 ALLOCATE (avh1b_ss(nsotot, nsotot))
361 ALLOCATE (avh1b_00(nsotot, nsotot))
363 NULLIFY (qlm_gg, g0_h)
366 qlm_gg=qlm_gg, g0_h=g0_h)
368 IF (
PRESENT(local_rho_set_2nd))
THEN
369 NULLIFY (qlm_gg_2nd, g0_h_2nd)
371 l0_ikind=lmax0_2nd, &
372 qlm_gg=qlm_gg_2nd, g0_h=g0_h_2nd)
376 IF (nchan_0 > max(max_iso_not0, max_iso_not0_nuc))
THEN
377 cpabort(
"channels for rho0 > # max of spherical harmonics")
380 NULLIFY (vh1_h, vh1_s)
381 ALLOCATE (vh1_h(nr, max_iso_not0))
382 ALLOCATE (vh1_s(nr, max(max_iso_not0, max_iso_not0_nuc, nchan_0)))
384 NULLIFY (lmax_nuc, lmin_nuc, npgf_nuc, zet_nuc, gsph_nuc)
393 lmin=lmin_nuc, maxso=maxso_nuc, npgf=npgf_nuc, &
394 maxl=maxl_nuc, nset=nset_nuc, zet=zet_nuc)
396 max_s_harm_nuc = cneo_potential%harmonics%max_s_harm
397 llmax_nuc = cneo_potential%harmonics%llmax
398 nsotot_nuc = maxso_nuc*nset_nuc
399 ALLOCATE (gsph_nuc(nr, nsotot_nuc))
401 ALLOCATE (avh1b_hh_nuc(nsotot_nuc, nsotot_nuc))
402 ALLOCATE (avh1b_ss_nuc(nsotot_nuc, nsotot_nuc))
403 ALLOCATE (avh1b_00_nuc(nsotot_nuc, nsotot_nuc))
406 ALLOCATE (cg_list(2,
nsoset(max(maxl, maxl_nuc))**2, &
407 max(max_s_harm, max_s_harm_nuc)), &
408 cg_n_list(max(max_s_harm, max_s_harm_nuc)))
410 NULLIFY (rrad_z, my_cg)
411 my_cg => harmonics%my_CG
419 sqrtwr(1:nr) = sqrt(grid_atom%wr(1:nr))
421 g0_h_w(1:nr, l_ang) = g0_h(1:nr, l_ang)*grid_atom%wr(1:nr)
427 DO ipgf1 = 1, npgf(iset1)
428 gexp(1:nr) = exp(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))*sqrtwr(1:nr)
429 DO is1 =
nsoset(lmin(iset1) - 1) + 1,
nsoset(lmax(iset1))
430 iso = is1 + (ipgf1 - 1)*n1 + m1
431 l_ang =
indso(1, is1)
432 gsph(1:nr, iso) = grid_atom%rad2l(1:nr, l_ang)*gexp(1:nr)
441 iatom = atom_list(iat)
442 rhoz_cneo_set(iatom)%pmat = 0.0_dp
443 rhoz_cneo_set(iatom)%cpc_h = 0.0_dp
444 rhoz_cneo_set(iatom)%cpc_s = 0.0_dp
445 rhoz_cneo_set(iatom)%e_core = 0.0_dp
446 rhoz_cneo_set(iatom)%ready = .false.
451 DO iset1 = 1, nset_nuc
452 n1 =
nsoset(lmax_nuc(iset1))
453 DO ipgf1 = 1, npgf_nuc(iset1)
454 gexp(1:nr) = exp(-zet_nuc(ipgf1, iset1)*grid_atom%rad2(1:nr))*sqrtwr(1:nr)
455 DO is1 =
nsoset(lmin_nuc(iset1) - 1) + 1,
nsoset(lmax_nuc(iset1))
456 iso = is1 + (ipgf1 - 1)*n1 + m1
457 l_ang =
indso(1, is1)
458 gsph_nuc(1:nr, iso) = cneo_potential%rad2l(1:nr, l_ang)*gexp(1:nr)
466 num_pe = para_env%num_pe
467 mepos = para_env%mepos
470 DO iat = bo(1), bo(2)
471 iatom = atom_list(iat)
472 rho_atom => rho_atom_set(iatom)
474 NULLIFY (rrad_z, vrrad_z, rrad_0, vrrad_0)
475 IF (core_charge .AND. .NOT. cneo)
THEN
476 rrad_z => rhoz_set(ikind)%r_coef
478 IF (my_core_2nd .AND. .NOT. cneo)
THEN
479 IF (l_2nd_local_rho)
THEN
480 vrrad_z => rhoz_set_2nd(ikind)%vr_coef
482 vrrad_z => rhoz_set(ikind)%vr_coef
485 rrad_0 => rho0_atom_set(iatom)%rho0_rad_h%r_coef
486 vrrad_0 => rho0_atom_set(iatom)%vrho0_rad_h%r_coef
487 IF (l_2nd_local_rho)
THEN
488 rho_atom => rho_atom_set_2nd(iatom)
489 vrrad_0 => rho0_atom_set_2nd(iatom)%vrho0_rad_h%r_coef
491 IF (my_periodic .AND. back_ch > 1.e-3_dp)
THEN
492 factor = -2.0_dp*
pi/3.0_dp*sqrt(
fourpi)*qs_charges%background
497 CALL vh_1c_atom_potential(rho_atom, vrrad_0, &
498 grid_atom, my_core_2nd .AND. .NOT. cneo, &
499 vrrad_z, vh1_h, vh1_s, &
500 nchan_0, nspins, max_iso_not0, factor)
502 IF (l_2nd_local_rho) rho_atom => rho_atom_set(iatom)
504 ecoul_1_z_cneo = 0.0_dp
506 rhoz_cneo => rhoz_cneo_set(iatom)
509 DO iso = 1, max_iso_not0_nuc
510 vh1_s(:, iso) = vh1_s(:, iso) + rhoz_cneo%vrho_rad_s(:, iso)
517 avh1b_hh_nuc, avh1b_ss_nuc, avh1b_00_nuc, vh1_h, vh1_s, &
518 max_iso_not0, max_iso_not0_nuc, &
519 max_s_harm_nuc, llmax_nuc, cg_list, cg_n_list, &
520 nset_nuc, npgf_nuc, lmin_nuc, lmax_nuc, nsotot_nuc, maxso_nuc, &
521 nchan_0, gsph_nuc, g0_h_w, cneo_potential%harmonics%my_CG, &
522 cneo_potential%Qlm_gg)
526 npgf_nuc, lmin_nuc, lmax_nuc, maxl_nuc, maxso_nuc)
529 qs_charges%total_rho1_hard_nuc = qs_charges%total_rho1_hard_nuc - sqrt(
fourpi) &
530 *sum(rhoz_cneo%rho_rad_h(:, 1)*grid_atom%wr(:))
531 rho0_mpole%tot_rhoz_cneo_s = rho0_mpole%tot_rhoz_cneo_s - sqrt(
fourpi) &
532 *sum(rhoz_cneo%rho_rad_s(:, 1)*grid_atom%wr(:))
537 DO iso = 1, min(max_iso_not0, max_iso_not0_nuc)
538 ecoul_1_z_cneo = ecoul_1_z_cneo + 0.5_dp* &
539 (sum(vh1_h(:, iso)*rhoz_cneo%rho_rad_h(:, iso)*grid_atom%wr(:)) &
540 - sum(vh1_s(:, iso)*rhoz_cneo%rho_rad_s(:, iso)*grid_atom%wr(:)))
542 DO iso = max_iso_not0 + 1, max_iso_not0_nuc
543 ecoul_1_z_cneo = ecoul_1_z_cneo - 0.5_dp* &
544 sum(vh1_s(:, iso)*rhoz_cneo%rho_rad_s(:, iso)*grid_atom%wr(:))
553 DO iso = 1, min(max_iso_not0, max_iso_not0_nuc)
554 vh1_h(:, iso) = vh1_h(:, iso) + rhoz_cneo%vrho_rad_h(:, iso)
558 CALL vh_1c_atom_energy(energy_hartree_1c, ecoul_1c, rho_atom, rrad_0, &
559 grid_atom, iatom, core_charge .AND. .NOT. cneo, &
560 rrad_z, vh1_h, vh1_s, nchan_0, nspins, max_iso_not0)
562 IF (l_2nd_local_rho) rho_atom => rho_atom_set_2nd(iatom)
565 CALL set_ecoul_1c(ecoul_1c, iatom, ecoul_1_z=ecoul_1c(iatom)%ecoul_1_z + ecoul_1_z_cneo)
566 energy_hartree_1c = energy_hartree_1c + ecoul_1_z_cneo
569 CALL vh_1c_atom_integrals(rho_atom, &
572 avh1b_hh, avh1b_ss, avh1b_00, vh1_h, vh1_s, max_iso_not0, &
573 max_s_harm, llmax, cg_list, cg_n_list, &
574 nset, npgf, lmin, lmax, nsotot, maxso, nspins, nchan_0, gsph, &
575 g0_h_w, my_cg, qlm_gg)
579 DEALLOCATE (avh1b_hh)
580 DEALLOCATE (avh1b_ss)
581 DEALLOCATE (avh1b_00)
582 IF (
ALLOCATED(avh1b_hh_nuc))
DEALLOCATE (avh1b_hh_nuc)
583 IF (
ALLOCATED(avh1b_ss_nuc))
DEALLOCATE (avh1b_ss_nuc)
584 IF (
ALLOCATED(avh1b_00_nuc))
DEALLOCATE (avh1b_00_nuc)
585 DEALLOCATE (vh1_h, vh1_s)
586 DEALLOCATE (cg_list, cg_n_list)
588 IF (
ASSOCIATED(gsph_nuc))
DEALLOCATE (gsph_nuc)
590 DEALLOCATE (sqrtwr, g0_h_w)
595 iatom = atom_list(iat)
596 CALL para_env%sum(rhoz_cneo_set(iatom)%pmat)
597 CALL para_env%sum(rhoz_cneo_set(iatom)%cpc_h)
598 CALL para_env%sum(rhoz_cneo_set(iatom)%cpc_s)
599 CALL para_env%sum(rhoz_cneo_set(iatom)%e_core)
600 rhoz_cneo_set(iatom)%ready = .true.
610 CALL para_env%sum(energy_hartree_1c)
611 CALL para_env%sum(qs_charges%total_rho1_hard_nuc)
612 CALL para_env%sum(rho0_mpole%tot_rhoz_cneo_s)
614 CALL timestop(handle)
634 SUBROUTINE vh_1c_atom_potential(rho_atom, vrrad_0, &
635 grid_atom, core_charge, vrrad_z, Vh1_h, Vh1_s, &
636 nchan_0, nspins, max_iso_not0, bfactor)
639 REAL(
dp),
DIMENSION(:, :),
POINTER :: vrrad_0
641 LOGICAL,
INTENT(IN) :: core_charge
642 REAL(
dp),
DIMENSION(:),
POINTER :: vrrad_z
643 REAL(
dp),
DIMENSION(:, :),
POINTER :: vh1_h, vh1_s
644 INTEGER,
INTENT(IN) :: nchan_0, nspins, max_iso_not0
645 REAL(
dp),
INTENT(IN) :: bfactor
647 INTEGER :: ir, iso, ispin, nr
653 CALL get_rho_atom(rho_atom=rho_atom, vrho_rad_h=vr_h, vrho_rad_s=vr_s)
658 IF (core_charge) vh1_h(:, 1) = vrrad_z(:)
661 vh1_s(:, iso) = vrrad_0(:, iso)
665 DO iso = 1, max_iso_not0
666 vh1_h(:, iso) = vh1_h(:, iso) + vr_h(ispin)%r_coef(:, iso)
667 vh1_s(:, iso) = vh1_s(:, iso) + vr_s(ispin)%r_coef(:, iso)
671 IF (bfactor /= 0._dp)
THEN
673 vh1_h(ir, 1) = vh1_h(ir, 1) + bfactor*grid_atom%rad2(ir)*grid_atom%wr(ir)
674 vh1_s(ir, 1) = vh1_s(ir, 1) + bfactor*grid_atom%rad2(ir)*grid_atom%wr(ir)
678 END SUBROUTINE vh_1c_atom_potential
698 SUBROUTINE vh_1c_atom_energy(energy_hartree_1c, ecoul_1c, rho_atom, rrad_0, &
699 grid_atom, iatom, core_charge, rrad_z, Vh1_h, Vh1_s, &
700 nchan_0, nspins, max_iso_not0)
702 REAL(
dp),
INTENT(INOUT) :: energy_hartree_1c
705 REAL(
dp),
DIMENSION(:, :),
POINTER :: rrad_0
707 INTEGER,
INTENT(IN) :: iatom
708 LOGICAL,
INTENT(IN) :: core_charge
709 REAL(
dp),
DIMENSION(:),
POINTER :: rrad_z
710 REAL(
dp),
DIMENSION(:, :),
POINTER :: vh1_h, vh1_s
711 INTEGER,
INTENT(IN) :: nchan_0, nspins, max_iso_not0
713 INTEGER :: iso, ispin, nr
714 REAL(
dp) :: ecoul_1_0, ecoul_1_h, ecoul_1_s, &
721 CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, rho_rad_s=r_s)
725 IF (core_charge)
THEN
726 ecoul_1_z = 0.5_dp*sum(vh1_h(:, 1)*rrad_z(:)*grid_atom%wr(:))
732 ecoul_1_0 = ecoul_1_0 + 0.5_dp*sum(vh1_s(:, iso)*rrad_0(:, iso)*grid_atom%wr(:))
739 DO iso = 1, max_iso_not0
740 ecoul_1_s = ecoul_1_s + 0.5_dp*sum(vh1_s(:, iso)*r_s(ispin)%r_coef(:, iso)*grid_atom%wr(:))
741 ecoul_1_h = ecoul_1_h + 0.5_dp*sum(vh1_h(:, iso)*r_h(ispin)%r_coef(:, iso)*grid_atom%wr(:))
745 CALL set_ecoul_1c(ecoul_1c, iatom, ecoul_1_z=ecoul_1_z, ecoul_1_0=ecoul_1_0)
746 CALL set_ecoul_1c(ecoul_1c=ecoul_1c, iatom=iatom, ecoul_1_h=ecoul_1_h, ecoul_1_s=ecoul_1_s)
748 energy_hartree_1c = energy_hartree_1c + ecoul_1_z - ecoul_1_0
749 energy_hartree_1c = energy_hartree_1c + ecoul_1_h - ecoul_1_s
751 END SUBROUTINE vh_1c_atom_energy
781 SUBROUTINE vh_1c_atom_integrals(rho_atom, &
782 aVh1b_hh, aVh1b_ss, aVh1b_00, Vh1_h, Vh1_s, max_iso_not0, &
783 max_s_harm, llmax, cg_list, cg_n_list, &
784 nset, npgf, lmin, lmax, nsotot, maxso, nspins, nchan_0, gsph, &
785 g0_h_w, my_CG, Qlm_gg)
788 REAL(
dp),
DIMENSION(:, :) :: avh1b_hh, avh1b_ss, avh1b_00
789 REAL(
dp),
DIMENSION(:, :),
POINTER :: vh1_h, vh1_s
790 INTEGER,
INTENT(IN) :: max_iso_not0, max_s_harm, llmax
791 INTEGER,
DIMENSION(:, :, :) :: cg_list
792 INTEGER,
DIMENSION(:) :: cg_n_list
793 INTEGER,
INTENT(IN) :: nset
794 INTEGER,
DIMENSION(:),
POINTER :: npgf, lmin, lmax
795 INTEGER,
INTENT(IN) :: nsotot, maxso, nspins, nchan_0
796 REAL(
dp),
DIMENSION(:, :),
POINTER :: gsph
797 REAL(
dp),
DIMENSION(:, 0:) :: g0_h_w
798 REAL(
dp),
DIMENSION(:, :, :),
POINTER :: my_cg, qlm_gg
800 INTEGER :: icg, ipgf1, ipgf2, ir, is1, is2, iset1, &
801 iset2, iso, iso1, iso2, ispin, l_ang, &
802 m1, m2, max_iso_not0_local, n1, n2, nr
803 REAL(
dp) :: gvg_0, gvg_h, gvg_s
804 TYPE(
rho_atom_coeff),
DIMENSION(:),
POINTER :: int_local_h, int_local_s
806 NULLIFY (int_local_h, int_local_s)
808 ga_vlocal_gb_h=int_local_h, &
809 ga_vlocal_gb_s=int_local_s)
822 CALL get_none0_cg_list(my_cg, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
823 max_s_harm, llmax, cg_list, cg_n_list, max_iso_not0_local)
826 DO ipgf1 = 1, npgf(iset1)
828 DO ipgf2 = 1, npgf(iset2)
830 DO iso = 1, min(nchan_0, max_iso_not0)
831 l_ang =
indso(1, iso)
832 gvg_0 = sum(vh1_s(:, iso)*g0_h_w(:, l_ang))
833 DO icg = 1, cg_n_list(iso)
834 is1 = cg_list(1, icg, iso)
835 is2 = cg_list(2, icg, iso)
837 iso1 = is1 + n1*(ipgf1 - 1) + m1
838 iso2 = is2 + n2*(ipgf2 - 1) + m2
843 gvg_h = gvg_h + gsph(ir, iso1)*gsph(ir, iso2)*vh1_h(ir, iso)
844 gvg_s = gvg_s + gsph(ir, iso1)*gsph(ir, iso2)*vh1_s(ir, iso)
847 avh1b_hh(iso1, iso2) = avh1b_hh(iso1, iso2) + gvg_h*my_cg(is1, is2, iso)
848 avh1b_ss(iso1, iso2) = avh1b_ss(iso1, iso2) + gvg_s*my_cg(is1, is2, iso)
849 avh1b_00(iso1, iso2) = avh1b_00(iso1, iso2) + gvg_0*qlm_gg(iso1, iso2, iso)
854 DO iso = nchan_0 + 1, max_iso_not0
855 DO icg = 1, cg_n_list(iso)
856 is1 = cg_list(1, icg, iso)
857 is2 = cg_list(2, icg, iso)
859 iso1 = is1 + n1*(ipgf1 - 1) + m1
860 iso2 = is2 + n2*(ipgf2 - 1) + m2
865 gvg_h = gvg_h + gsph(ir, iso1)*gsph(ir, iso2)*vh1_h(ir, iso)
866 gvg_s = gvg_s + gsph(ir, iso1)*gsph(ir, iso2)*vh1_s(ir, iso)
869 avh1b_hh(iso1, iso2) = avh1b_hh(iso1, iso2) + gvg_h*my_cg(is1, is2, iso)
870 avh1b_ss(iso1, iso2) = avh1b_ss(iso1, iso2) + gvg_s*my_cg(is1, is2, iso)
881 CALL daxpy(nsotot*nsotot, 1.0_dp, avh1b_hh, 1, int_local_h(ispin)%r_coef, 1)
882 CALL daxpy(nsotot*nsotot, 1.0_dp, avh1b_ss, 1, int_local_s(ispin)%r_coef, 1)
883 CALL daxpy(nsotot*nsotot, -1.0_dp, avh1b_00, 1, int_local_h(ispin)%r_coef, 1)
884 CALL daxpy(nsotot*nsotot, -1.0_dp, avh1b_00, 1, int_local_s(ispin)%r_coef, 1)
887 END SUBROUTINE vh_1c_atom_integrals
Define the atomic kind types and their sub types.
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.
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)
...
Handles all functions related to the CELL.
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
subroutine, public init_coulomb_local(hartree_local, natom)
...
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...
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...
subroutine, public allocate_ecoul_1center(ecoul_1c, natom)
...
subroutine, public set_ecoul_1c(ecoul_1c, iatom, ecoul_1_h, ecoul_1_s, ecoul_1_z, ecoul_1_0)
...
Defines the basic variable types.
integer, parameter, public dp
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
real(kind=dp), parameter, public fourpi
Interface to the message passing library MPI.
Provides Cartesian and spherical orbital pointers and indices.
integer, dimension(:), allocatable, public nsoset
integer, dimension(:, :), allocatable, public indso
container for various plainwaves related things
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
functions related to the poisson solver on regular grids
integer, parameter, public pw_poisson_periodic
container for information about total charges on the grids
A collection of functions used by CNEO-DFT (see J. Chem. Theory Comput. 2025, 21, 16,...
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.
subroutine, public calculate_rhoz_cneo(rho, potential, cg_list, cg_n_list, nset, npgf, lmin, lmax, maxl, maxso)
...
Types used by CNEO-DFT (see J. Chem. Theory Comput. 2025, 21, 16, 7865–7877)
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.
Define the quickstep kind type and their sub types.
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.
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.
subroutine, public get_local_rho(local_rho_set, rho_atom_set, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set)
...
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)
...
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)
...
All kind of helpful little routines.
pure integer function, dimension(2), public get_limit(m, n, me)
divide m entries into n parts, return size of part me
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
stores all the informations relevant to an mpi environment
contained for different pw related things
environment for the poisson solver
Container for information about total charges on the grids.
Provides all information about a quickstep kind.