comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-28 22:25:32 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Energy decomposition:wq (#3146)

Browse source
This commit is contained in:
Juerg Hutter 2023-12-08 12:01:04 +01:00 • committed by GitHub
parent b263019182
commit f5d7817782
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
10 changed files with 776 additions and 25 deletions
Download patch file Download diff file Expand all files Collapse all files

276
src/core_ae.F
View file

@ -13,6 +13,7 @@
! **************************************************************************************************
MODULE core_ae
USE ai_verfc, ONLY: verfc
USE ao_util, ONLY: exp_radius
USE atomic_kind_types, ONLY: atomic_kind_type,&
get_atomic_kind_set
USE basis_set_types, ONLY: gto_basis_set_p_type,&
@ -57,7 +58,7 @@ MODULE core_ae
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'core_ae'
PUBLIC :: build_core_ae
PUBLIC :: build_core_ae, build_erfc
CONTAINS
@ -568,6 +569,279 @@ CONTAINS
END SUBROUTINE verfc_force
! **************************************************************************************************
!> \brief Integrals = -Z*erfc(a*r)/r
!> \param matrix_h ...
!> \param qs_kind_set ...
!> \param atomic_kind_set ...
!> \param particle_set ...
!> \param calpha ...
!> \param ccore ...
!> \param eps_core_charge ...
!> \param sab_orb ...
!> \param sac_ae ...
! **************************************************************************************************
SUBROUTINE build_erfc(matrix_h, qs_kind_set, atomic_kind_set, particle_set, &
calpha, ccore, eps_core_charge, sab_orb, sac_ae)
TYPE(dbcsr_p_type) :: matrix_h
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: calpha, ccore
REAL(KIND=dp), INTENT(IN) :: eps_core_charge
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_orb, sac_ae
CHARACTER(LEN=*), PARAMETER :: routineN = 'build_erfc'
INTEGER :: handle, iatom, icol, ikind, img, irow, iset, jatom, jkind, jset, katom, kkind, &
ldai, ldsab, maxco, maxl, maxnset, maxsgf, mepos, na_plus, natom, nb_plus, ncoa, ncob, &
nij, nkind, nseta, nsetb, nthread, sgfa, sgfb, slot
INTEGER, DIMENSION(3) :: cellind
INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
npgfb, nsgfa, nsgfb
INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
LOGICAL :: dokp, found
REAL(KIND=dp) :: alpha_c, core_charge, core_radius, dab, &
dac, dbc, f0, rab2, rac2, rbc2, zeta_c
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ff
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: habd, work
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: hab, pab, verf, vnuc
REAL(KIND=dp), DIMENSION(3) :: rab, rac, rbc
REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
REAL(KIND=dp), DIMENSION(:, :), POINTER :: h_block, p_block, rpgfa, rpgfb, sphi_a, &
sphi_b, zeta, zetb
TYPE(all_potential_type), POINTER :: all_potential
TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
TYPE(neighbor_list_iterator_p_type), &
DIMENSION(:), POINTER :: ap_iterator
TYPE(sgp_potential_type), POINTER :: sgp_potential
!$ INTEGER(kind=omp_lock_kind), &
!$ ALLOCATABLE, DIMENSION(:) :: locks
!$ INTEGER :: lock_num, hash, hash1, hash2
!$ INTEGER(KIND=int_8) :: iatom8
!$ INTEGER, PARAMETER :: nlock = 501
MARK_USED(int_8)
CALL timeset(routineN, handle)
nkind = SIZE(atomic_kind_set)
natom = SIZE(particle_set)
ALLOCATE (basis_set_list(nkind))
DO ikind = 1, nkind
CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set_a)
IF (ASSOCIATED(basis_set_a)) THEN
basis_set_list(ikind)%gto_basis_set => basis_set_a
ELSE
NULLIFY (basis_set_list(ikind)%gto_basis_set)
END IF
END DO
CALL get_qs_kind_set(qs_kind_set, &
maxco=maxco, maxlgto=maxl, maxsgf=maxsgf, maxnset=maxnset)
CALL init_orbital_pointers(maxl + 1)
ldsab = MAX(maxco, maxsgf)
ldai = ncoset(maxl + 1)
nthread = 1
!$ nthread = omp_get_max_threads()
! iterator for basis/potential list
CALL neighbor_list_iterator_create(ap_iterator, sac_ae, search=.TRUE., nthread=nthread)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP SHARED (ap_iterator, basis_set_list, &
!$OMP matrix_h, atomic_kind_set, qs_kind_set, particle_set, &
!$OMP sab_orb, sac_ae, nthread, ncoset, nkind, calpha, ccore, eps_core_charge, &
!$OMP slot, ldsab, maxnset, ldai, maxl, maxco, dokp, locks, natom) &
!$OMP PRIVATE (ikind, jkind, iatom, jatom, rab, basis_set_a, basis_set_b, &
!$OMP first_sgfa, la_max, la_min, npgfa, nsgfa, sphi_a, &
!$OMP zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, sphi_b, &
!$OMP zetb, zeta_c, alpha_c, core_charge, dab, irow, icol, h_block, found, iset, ncoa, &
!$OMP sgfa, jset, ncob, sgfb, nsgfb, p_block, work, pab, hab, kkind, nseta, &
!$OMP rac, dac, rbc, rab2, rac2, rbc2, dbc, na_plus, nb_plus, verf, vnuc, &
!$OMP set_radius_a, core_radius, rpgfa, mepos, &
!$OMP habd, f0, katom, cellind, img, nij, ff, &
!$OMP sgp_potential, all_potential, hash, hash1, hash2, iatom8)
!$OMP SINGLE
!$ ALLOCATE (locks(nlock))
!$OMP END SINGLE
!$OMP DO
!$ DO lock_num = 1, nlock
!$ call omp_init_lock(locks(lock_num))
!$ END DO
!$OMP END DO
mepos = 0
!$ mepos = omp_get_thread_num()
ALLOCATE (hab(ldsab, ldsab, maxnset*maxnset), work(ldsab, ldsab))
ALLOCATE (verf(ldai, ldai, 2*maxl + 1), vnuc(ldai, ldai, 2*maxl + 1), ff(0:2*maxl))
!$OMP DO SCHEDULE(GUIDED)
DO slot = 1, sab_orb(1)%nl_size
ikind = sab_orb(1)%nlist_task(slot)%ikind
jkind = sab_orb(1)%nlist_task(slot)%jkind
iatom = sab_orb(1)%nlist_task(slot)%iatom
jatom = sab_orb(1)%nlist_task(slot)%jatom
cellind(:) = sab_orb(1)%nlist_task(slot)%cell(:)
rab(1:3) = sab_orb(1)%nlist_task(slot)%r(1:3)
basis_set_a => basis_set_list(ikind)%gto_basis_set
IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
basis_set_b => basis_set_list(jkind)%gto_basis_set
IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
!$ iatom8 = INT(iatom - 1, int_8)*INT(natom, int_8) + INT(jatom, int_8)
!$ hash1 = INT(MOD(iatom8, INT(nlock, int_8)) + 1)
! basis ikind
first_sgfa => basis_set_a%first_sgf
la_max => basis_set_a%lmax
la_min => basis_set_a%lmin
npgfa => basis_set_a%npgf
nseta = basis_set_a%nset
nsgfa => basis_set_a%nsgf_set
rpgfa => basis_set_a%pgf_radius
set_radius_a => basis_set_a%set_radius
sphi_a => basis_set_a%sphi
zeta => basis_set_a%zet
! basis jkind
first_sgfb => basis_set_b%first_sgf
lb_max => basis_set_b%lmax
lb_min => basis_set_b%lmin
npgfb => basis_set_b%npgf
nsetb = basis_set_b%nset
nsgfb => basis_set_b%nsgf_set
rpgfb => basis_set_b%pgf_radius
set_radius_b => basis_set_b%set_radius
sphi_b => basis_set_b%sphi
zetb => basis_set_b%zet
dab = SQRT(SUM(rab*rab))
img = 1
! *** Use the symmetry of the first derivatives ***
IF (iatom == jatom) THEN
f0 = 1.0_dp
ELSE
f0 = 2.0_dp
END IF
! *** Create matrix blocks for a new matrix block column ***
IF (iatom <= jatom) THEN
irow = iatom
icol = jatom
ELSE
irow = jatom
icol = iatom
END IF
NULLIFY (h_block)
CALL dbcsr_get_block_p(matrix=matrix_h%matrix, &
row=irow, col=icol, BLOCK=h_block, found=found)
! loop over all kinds of atoms
hab = 0._dp
DO kkind = 1, nkind
CALL get_qs_kind(qs_kind_set(kkind), zeff=zeta_c)
alpha_c = calpha(kkind)
core_charge = ccore(kkind)
core_radius = exp_radius(0, alpha_c, eps_core_charge, core_charge)
core_radius = MAX(core_radius, 10.0_dp)
CALL nl_set_sub_iterator(ap_iterator, ikind, kkind, iatom, mepos=mepos)
DO WHILE (nl_sub_iterate(ap_iterator, mepos=mepos) == 0)
CALL get_iterator_info(ap_iterator, jatom=katom, r=rac, mepos=mepos)
dac = SQRT(SUM(rac*rac))
rbc(:) = rac(:) - rab(:)
dbc = SQRT(SUM(rbc*rbc))
IF ((MAXVAL(set_radius_a(:)) + core_radius < dac) .OR. &
(MAXVAL(set_radius_b(:)) + core_radius < dbc)) THEN
CYCLE
END IF
DO iset = 1, nseta
IF (set_radius_a(iset) + core_radius < dac) CYCLE
ncoa = npgfa(iset)*ncoset(la_max(iset))
sgfa = first_sgfa(1, iset)
DO jset = 1, nsetb
IF (set_radius_b(jset) + core_radius < dbc) CYCLE
ncob = npgfb(jset)*ncoset(lb_max(jset))
sgfb = first_sgfb(1, jset)
IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
rab2 = dab*dab
rac2 = dac*dac
rbc2 = dbc*dbc
nij = jset + (iset - 1)*maxnset
!
CALL verfc( &
la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
alpha_c, core_radius, zeta_c, core_charge, &
rab, rab2, rac, rac2, rbc2, hab(:, :, nij), verf, vnuc, ff(0:))
END DO
END DO
END DO
END DO
! *** Contract nuclear attraction integrals
DO iset = 1, nseta
ncoa = npgfa(iset)*ncoset(la_max(iset))
sgfa = first_sgfa(1, iset)
DO jset = 1, nsetb
ncob = npgfb(jset)*ncoset(lb_max(jset))
sgfb = first_sgfb(1, jset)
nij = jset + (iset - 1)*maxnset
!$ hash2 = MOD((iset - 1)*nsetb + jset, nlock) + 1
!$ hash = MOD(hash1 + hash2, nlock) + 1
work(1:ncoa, 1:nsgfb(jset)) = MATMUL(hab(1:ncoa, 1:ncob, nij), &
sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1))
!$ CALL omp_set_lock(locks(hash))
IF (iatom <= jatom) THEN
h_block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) = &
h_block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) + &
MATMUL(TRANSPOSE(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1)), work(1:ncoa, 1:nsgfb(jset)))
ELSE
h_block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) = &
h_block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) + &
MATMUL(TRANSPOSE(work(1:ncoa, 1:nsgfb(jset))), sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1))
END IF
!$ CALL omp_unset_lock(locks(hash))
END DO
END DO
END DO
DEALLOCATE (hab, work, verf, vnuc, ff)
!$OMP DO
!$ DO lock_num = 1, nlock
!$ call omp_destroy_lock(locks(lock_num))
!$ END DO
!$OMP END DO
!$OMP SINGLE
!$ DEALLOCATE (locks)
!$OMP END SINGLE NOWAIT
!$OMP END PARALLEL
CALL neighbor_list_iterator_release(ap_iterator)
DEALLOCATE (basis_set_list)
CALL timestop(handle)
END SUBROUTINE build_erfc
! **************************************************************************************************
END MODULE core_ae

255
src/ed_analysis.F
View file

@ -19,6 +19,7 @@ MODULE ed_analysis
USE bibliography, ONLY: Eriksen2020,&
cite_reference
USE cell_types, ONLY: cell_type
USE core_ae, ONLY: build_erfc
USE cp_control_types, ONLY: dft_control_type
USE cp_dbcsr_operations, ONLY: cp_dbcsr_plus_fm_fm_t,&
cp_dbcsr_sm_fm_multiply
@ -48,8 +49,10 @@ MODULE ed_analysis
section_vals_type,&
section_vals_val_get
USE kinds, ONLY: dp
USE mathconstants, ONLY: pi
USE message_passing, ONLY: mp_comm_type,&
mp_para_env_type
USE mulliken, ONLY: mulliken_charges
USE parallel_gemm_api, ONLY: parallel_gemm
USE particle_methods, ONLY: get_particle_set
USE particle_types, ONLY: particle_type
@ -57,12 +60,19 @@ MODULE ed_analysis
USE pw_env_types, ONLY: pw_env_get,&
pw_env_type
USE pw_methods, ONLY: pw_axpy,&
pw_scale
pw_scale,&
pw_transfer
USE pw_poisson_methods, ONLY: pw_poisson_solve
USE pw_poisson_types, ONLY: pw_poisson_type
USE pw_pool_types, ONLY: pw_pool_type
USE pw_types, ONLY: REALDATA3D,&
USE pw_types, ONLY: COMPLEXDATA1D,&
REALDATA3D,&
REALSPACE,&
RECIPROCALSPACE,&
pw_type
USE qs_core_energies, ONLY: calculate_ecore_overlap,&
USE qs_collocate_density, ONLY: calculate_rho_core
USE qs_core_energies, ONLY: calculate_ecore_alpha,&
calculate_ecore_overlap,&
calculate_ecore_self
USE qs_dispersion_pairpot, ONLY: calculate_dispersion_pairpot
USE qs_dispersion_types, ONLY: qs_dispersion_type
@ -72,6 +82,7 @@ MODULE ed_analysis
USE qs_gcp_method, ONLY: calculate_gcp_pairpot
USE qs_gcp_types, ONLY: qs_gcp_type
USE qs_integrate_potential, ONLY: integrate_v_core_rspace,&
integrate_v_gaussian_rspace,&
integrate_v_rspace
USE qs_kind_types, ONLY: get_qs_kind,&
qs_kind_type
@ -83,6 +94,7 @@ MODULE ed_analysis
duplicate_mo_set,&
get_mo_set,&
mo_set_type
USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
USE qs_rho0_ggrid, ONLY: integrate_vhg0_rspace
USE qs_rho_atom_types, ONLY: rho_atom_type,&
zero_rho_atom_integrals
@ -123,15 +135,14 @@ CONTAINS
nref, nspin
INTEGER, DIMENSION(2) :: nocc
INTEGER, DIMENSION(:), POINTER :: refbas_blk_sizes
LOGICAL :: detailed_ener, do_hfx, explicit, gapw, &
gapw_xc, ref_orb_canonical, &
skip_localize, uniform_occupation, uocc
REAL(KIND=dp) :: ateps, checksum, e1, e2, e_pot, ecc, &
egcp, ehfx, ekts, evdw, focc, &
LOGICAL :: detailed_ener, do_hfx, ewald_correction, explicit, gapw, gapw_xc, &
ref_orb_canonical, skip_localize, uniform_occupation, uocc
REAL(KIND=dp) :: ateps, checksum, e1, e2, e_pot, ealpha, &
ecc, egcp, ehfx, ekts, evdw, focc, &
sum_energy
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ate1h, ate1xc, atecc, ateks, atener, &
odiag
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: atdet, mcharge, mweight
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: amval, ate1h, ate1xc, atecc, ateks, &
atener, atewald, odiag
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: atdet, atmul, mcharge, mweight
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: bcenter
REAL(KIND=dp), DIMENSION(:), POINTER :: occupation_numbers
TYPE(cell_type), POINTER :: cell
@ -141,9 +152,10 @@ CONTAINS
rotmat
TYPE(cp_fm_type), POINTER :: cloc, moref
TYPE(dbcsr_distribution_type) :: dbcsr_dist
TYPE(dbcsr_p_type) :: dve_mat
TYPE(dbcsr_p_type), ALLOCATABLE, DIMENSION(:) :: exc_mat, ks_mat, vhxc_mat
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_h, matrix_hfx, matrix_ks, &
matrix_s, matrix_t
matrix_p, matrix_s, matrix_t
TYPE(dbcsr_type) :: dkmat, dmat
TYPE(dbcsr_type), POINTER :: smat
TYPE(dft_control_type), POINTER :: dft_control
@ -159,6 +171,7 @@ CONTAINS
TYPE(qs_gcp_type), POINTER :: gcp_env
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_kind_type), POINTER :: qs_kind
TYPE(qs_rho_type), POINTER :: rho
TYPE(section_vals_type), POINTER :: hfx_sections, input
CALL section_vals_get(input_section, explicit=explicit)
@ -179,6 +192,8 @@ CONTAINS
CALL section_vals_val_get(input_section, "REFERENCE_ORB_CANONICAL", l_val=ref_orb_canonical)
CALL section_vals_val_get(input_section, "DETAILED_ENERGY", l_val=detailed_ener)
CALL section_vals_val_get(input_section, "SKIP_LOCALIZATION", l_val=skip_localize)
CALL section_vals_val_get(input_section, "EWALD_ALPHA_PARAMETER", r_val=ealpha)
ewald_correction = (ealpha /= 0.0_dp)
! k-points?
CALL get_qs_env(qs_env, dft_control=dft_control)
nimages = dft_control%nimages
@ -190,6 +205,12 @@ CONTAINS
END IF
gapw = dft_control%qs_control%gapw
gapw_xc = dft_control%qs_control%gapw_xc
IF (ewald_correction) THEN
IF (gapw .OR. gapw_xc) THEN
CALL cp_warn(__LOCATION__, "Ewald Correction for GAPW and GAPW_XC not available.")
CPABORT("Ewald Correction")
END IF
END IF
CALL get_qs_env(qs_env, input=input)
hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF")
CALL section_vals_get(hfx_sections, explicit=do_hfx)
@ -419,15 +440,17 @@ CONTAINS
! energy arrays
ALLOCATE (atener(natom), ateks(natom), atecc(natom))
ALLOCATE (ate1xc(natom), ate1h(natom))
ALLOCATE (ate1xc(natom), ate1h(natom), atewald(natom))
atener = 0.0_dp
ateks = 0.0_dp
atecc = 0.0_dp
ate1xc = 0.0_dp
ate1h = 0.0_dp
atewald = 0.0_dp
IF (detailed_ener) THEN
ALLOCATE (atdet(natom, 10))
ALLOCATE (atdet(natom, 10), atmul(natom, 10), amval(natom))
atdet = 0.0_dp
atmul = 0.0_dp
END IF
! atom dependent density matrix
CALL dbcsr_create(dkmat, template=smat)
@ -455,6 +478,14 @@ CONTAINS
!
CALL vhxc_correction(qs_env, vhxc_mat, exc_mat, atecc, ate1xc, ate1h)
!
IF (ewald_correction) THEN
ALLOCATE (dve_mat%matrix)
CALL dbcsr_create(dve_mat%matrix, template=matrix_h(1)%matrix)
CALL dbcsr_copy(dve_mat%matrix, matrix_h(1)%matrix)
CALL dbcsr_set(dve_mat%matrix, 0.0_dp)
CALL vh_ewald_correction(qs_env, ealpha, dve_mat, atewald)
END IF
!
DO ispin = 1, nspin
CALL dbcsr_add(ks_mat(ispin)%matrix, vhxc_mat(ispin)%matrix, 1.0_dp, 1.0_dp)
END DO
@ -493,6 +524,7 @@ CONTAINS
CALL dbcsr_scale(matrix_hfx(ispin)%matrix, 0.5_dp)
CALL dbcsr_dot(dkmat, matrix_hfx(ispin)%matrix, ehfx)
atdet(iatom, 5) = atdet(iatom, 5) + ehfx
CALL dbcsr_scale(matrix_hfx(ispin)%matrix, 2.0_dp)
ELSE
ehfx = 0.0_dp
END IF
@ -500,11 +532,45 @@ CONTAINS
atdet(iatom, 3) = atdet(iatom, 3) + ecc - e2 - e1 - ehfx
atdet(iatom, 4) = atdet(iatom, 4) + e1
END IF
IF (ewald_correction) THEN
CALL dbcsr_dot(dkmat, dve_mat%matrix, ecc)
atewald(iatom) = atewald(iatom) + ecc
END IF
!
END DO
DEALLOCATE (mweight)
END DO
!
IF (detailed_ener) THEN
! Mulliken
CALL get_qs_env(qs_env, rho=rho, para_env=para_env)
CALL qs_rho_get(rho, rho_ao=matrix_p)
! kinetic energy
DO ispin = 1, nspin
CALL mulliken_charges(matrix_p(ispin)%matrix, matrix_t(1)%matrix, &
para_env, amval)
atmul(1:natom, 1) = atmul(1:natom, 1) + amval(1:natom)
CALL mulliken_charges(matrix_p(ispin)%matrix, matrix_h(1)%matrix, &
para_env, amval)
atmul(1:natom, 2) = atmul(1:natom, 2) + amval(1:natom)
atmul(1:natom, 3) = atmul(1:natom, 3) - amval(1:natom)
CALL mulliken_charges(matrix_p(ispin)%matrix, ks_mat(ispin)%matrix, &
para_env, amval)
atmul(1:natom, 3) = atmul(1:natom, 3) + amval(1:natom)
IF (do_hfx) THEN
CALL mulliken_charges(matrix_p(ispin)%matrix, matrix_hfx(ispin)%matrix, &
para_env, amval)
atmul(1:natom, 5) = atmul(1:natom, 5) + 0.5_dp*amval(1:natom)
atmul(1:natom, 3) = atmul(1:natom, 3) - 0.5_dp*amval(1:natom)
END IF
CALL mulliken_charges(matrix_p(ispin)%matrix, exc_mat(1)%matrix, &
para_env, amval)
atmul(1:natom, 4) = atmul(1:natom, 4) + amval(1:natom)
atmul(1:natom, 3) = atmul(1:natom, 3) - amval(1:natom)
END DO
atmul(1:natom, 2) = atmul(1:natom, 2) - atmul(1:natom, 1)
END IF
!
CALL dbcsr_release(dkmat)
DO ispin = 1, nspin
CALL dbcsr_release(ks_mat(ispin)%matrix)
@ -525,6 +591,10 @@ CONTAINS
END DO
DEALLOCATE (matrix_hfx)
END IF
IF (ewald_correction) THEN
CALL dbcsr_release(dve_mat%matrix)
DEALLOCATE (dve_mat%matrix)
END IF
CALL cp_fm_release(orb_weight)
CALL cp_fm_release(ciao)
CALL cp_fm_release(rotmat)
@ -555,11 +625,19 @@ CONTAINS
CALL group%sum(atecc)
atener(1:natom) = atener(1:natom) + atecc(1:natom)
IF (detailed_ener) atdet(1:natom, 7) = atdet(1:natom, 7) + atecc(1:natom)
IF (ewald_correction) THEN
atewald(1:natom) = atewald(1:natom) - atecc(1:natom)
END IF
! e self
atecc(1:natom) = 0.0_dp
CALL calculate_ecore_self(qs_env, atecc=atecc)
CALL group%sum(atecc)
atener(1:natom) = atener(1:natom) + atecc(1:natom)
IF (detailed_ener) atdet(1:natom, 7) = atdet(1:natom, 7) + atecc(1:natom)
IF (ewald_correction) THEN
atewald(1:natom) = atewald(1:natom) - atecc(1:natom)
CALL calculate_ecore_alpha(qs_env, ealpha, atewald)
END IF
! vdW pair-potentials
atecc(1:natom) = 0.0_dp
CALL get_qs_env(qs_env=qs_env, dispersion_env=dispersion_env)
@ -579,7 +657,7 @@ CONTAINS
IF (detailed_ener) THEN
IF (unit_nr > 0) THEN
WRITE (unit_nr, FMT="(/,T2,A)") "Detailed Atomic Energy Components "
WRITE (unit_nr, FMT="(/,T2,A)") "Detailed IAO Atomic Energy Components "
CALL write_atdet(unit_nr, atdet(:, 1), "Kinetic Energy")
CALL write_atdet(unit_nr, atdet(:, 2), "Short-Range Core and/or PP Energy")
CALL write_atdet(unit_nr, atdet(:, 3), "Hartree Energy (long-ranged part)")
@ -600,6 +678,24 @@ CONTAINS
atecc(iatom) = SUM(atdet(iatom, 1:10)) - atener(iatom)
END DO
CALL write_atdet(unit_nr, atecc(:), "Missing Energy")
!
WRITE (unit_nr, FMT="(/,T2,A)") "Detailed Mulliken Atomic Energy Components "
CALL write_atdet(unit_nr, atmul(:, 1), "Kinetic Energy")
CALL write_atdet(unit_nr, atmul(:, 2), "Short-Range Core and/or PP Energy")
CALL write_atdet(unit_nr, atmul(:, 3), "Hartree Energy (long-ranged part)")
CALL write_atdet(unit_nr, atmul(:, 5), "Exact Exchange Energy")
CALL write_atdet(unit_nr, atmul(:, 4), "Exchange-Correlation Energy")
CALL write_atdet(unit_nr, atdet(:, 6), "Atomic Core Hartree: Int(nc V(n+nc) dx")
CALL write_atdet(unit_nr, atdet(:, 7), "Core Self Energy Corrections")
IF (gapw) THEN
CALL write_atdet(unit_nr, atdet(:, 9), "GAPW: 1-center Hartree Energy")
END IF
IF (gapw .OR. gapw_xc) THEN
CALL write_atdet(unit_nr, atdet(:, 8), "GAPW: 1-center XC Energy")
END IF
IF (ABS(evdw) > 1.E-14_dp) THEN
CALL write_atdet(unit_nr, atdet(:, 10), "vdW Pairpotential Energy")
END IF
END IF
END IF
@ -609,15 +705,23 @@ CONTAINS
CALL write_atener(unit_nr, particle_set, atener, "Atomic Energy Decomposition")
sum_energy = SUM(atener(:))
checksum = ABS(e_pot - sum_energy)
WRITE (UNIT=unit_nr, FMT="(/,(T2,A,T56,F25.13))") &
WRITE (UNIT=unit_nr, FMT="((T2,A,T56,F25.13))") &
"Potential energy (Atomic):", sum_energy, &
"Potential energy (Total) :", e_pot, &
"Difference :", checksum
CPASSERT((checksum < ateps*ABS(e_pot)))
!
IF (ewald_correction) THEN
WRITE (UNIT=unit_nr, FMT="(/,(T2,A,F10.3,A))") "Universal Ewald Parameter: ", ealpha, " [Angstrom]"
CALL write_atener(unit_nr, particle_set, atewald, "Change in Atomic Energy Decomposition")
sum_energy = SUM(atewald(:))
WRITE (UNIT=unit_nr, FMT="((T2,A,T56,F25.13))") &
"Total Change in Potential energy:", sum_energy
END IF
END IF
IF (detailed_ener) THEN
DEALLOCATE (atdet)
DEALLOCATE (atdet, atmul, amval)
END IF
IF (unit_nr > 0) THEN
@ -625,7 +729,7 @@ CONTAINS
"!--------------------------- END OF ED ANALYSIS ------------------------------!"
END IF
DEALLOCATE (bcenter)
DEALLOCATE (atener, ateks, atecc, ate1xc, ate1h)
DEALLOCATE (atener, ateks, atecc, ate1xc, ate1h, atewald)
CALL timestop(handle)
@ -789,6 +893,119 @@ CONTAINS
END SUBROUTINE vhxc_correction
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param ealpha ...
!> \param vh_mat ...
!> \param atewd ...
! **************************************************************************************************
SUBROUTINE vh_ewald_correction(qs_env, ealpha, vh_mat, atewd)
TYPE(qs_environment_type), POINTER :: qs_env
REAL(KIND=dp), INTENT(IN) :: ealpha
TYPE(dbcsr_p_type) :: vh_mat
REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: atewd
CHARACTER(len=*), PARAMETER :: routineN = 'vh_ewald_correction'
INTEGER :: handle, ikind, ispin, natom, nkind
REAL(KIND=dp) :: eps_core_charge, rhotot, zeff
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: alpha, atecc, ccore
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(dft_control_type), POINTER :: dft_control
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_orb, sac_ae
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(pw_env_type), POINTER :: pw_env
TYPE(pw_poisson_type), POINTER :: poisson_env
TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
TYPE(pw_type) :: rho_tot_ewd_g, rho_tot_ewd_r, &
v_hewd_gspace, v_hewd_rspace
TYPE(pw_type), DIMENSION(:), POINTER :: rho_r
TYPE(pw_type), POINTER :: v_hartree_rspace
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_rho_type), POINTER :: rho
CALL timeset(routineN, handle)
natom = SIZE(atewd)
ALLOCATE (atecc(natom))
CALL get_qs_env(qs_env=qs_env, nkind=nkind, qs_kind_set=qs_kind_set, &
dft_control=dft_control)
ALLOCATE (alpha(nkind), ccore(nkind))
DO ikind = 1, nkind
CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
alpha(ikind) = ealpha
ccore(ikind) = zeff*(ealpha/pi)**1.5_dp
END DO
!
CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
poisson_env=poisson_env)
!
CALL auxbas_pw_pool%create_pw(v_hewd_gspace, &
use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE)
CALL auxbas_pw_pool%create_pw(v_hewd_rspace, &
use_data=REALDATA3D, in_space=REALSPACE)
CALL auxbas_pw_pool%create_pw(rho_tot_ewd_g, &
use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE)
CALL auxbas_pw_pool%create_pw(rho_tot_ewd_r, &
use_data=REALDATA3D, in_space=REALSPACE)
rhotot = 0.0_dp
CALL calculate_rho_core(rho_tot_ewd_r, rhotot, qs_env, calpha=alpha, ccore=ccore)
! Get the total density in g-space [ions + electrons]
CALL get_qs_env(qs_env=qs_env, rho=rho)
CALL qs_rho_get(rho, rho_r=rho_r)
DO ispin = 1, dft_control%nspins
CALL pw_axpy(rho_r(ispin), rho_tot_ewd_r)
END DO
CALL pw_transfer(rho_tot_ewd_r, rho_tot_ewd_g)
!
CALL pw_poisson_solve(poisson_env, density=rho_tot_ewd_g, vhartree=v_hewd_gspace)
CALL pw_transfer(v_hewd_gspace, v_hewd_rspace)
CALL pw_scale(v_hewd_rspace, v_hewd_rspace%pw_grid%dvol)
atecc = 0.0_dp
CALL integrate_v_gaussian_rspace(v_hewd_rspace, qs_env, alpha, ccore, atecc)
atewd(1:natom) = atewd(1:natom) + atecc(1:natom)
!
atecc = 0.0_dp
CALL get_qs_env(qs_env, v_hartree_rspace=v_hartree_rspace)
CALL integrate_v_core_rspace(v_hartree_rspace, qs_env, atecc=atecc)
atewd(1:natom) = atewd(1:natom) - atecc(1:natom)
!
CALL pw_axpy(v_hartree_rspace, v_hewd_rspace, -1.0_dp)
CALL dbcsr_set(vh_mat%matrix, 0.0_dp)
CALL integrate_v_rspace(qs_env=qs_env, v_rspace=v_hewd_rspace, hmat=vh_mat, &
calculate_forces=.FALSE.)
!
CALL dbcsr_scale(vh_mat%matrix, 0.5_dp)
!
! delta erfc
eps_core_charge = dft_control%qs_control%eps_core_charge
CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, &
particle_set=particle_set, sab_orb=sab_orb, sac_ppl=sac_ae)
CALL build_erfc(vh_mat, qs_kind_set, atomic_kind_set, particle_set, &
alpha, ccore, eps_core_charge, sab_orb, sac_ae)
CALL dbcsr_scale(vh_mat%matrix, -1.0_dp)
DO ikind = 1, nkind
CALL get_qs_kind(qs_kind_set(ikind), alpha_core_charge=alpha(ikind), &
ccore_charge=ccore(ikind))
END DO
CALL build_erfc(vh_mat, qs_kind_set, atomic_kind_set, particle_set, &
alpha, ccore, eps_core_charge, sab_orb, sac_ae)
CALL dbcsr_scale(vh_mat%matrix, -1.0_dp)
!
CALL auxbas_pw_pool%give_back_pw(v_hewd_gspace)
CALL auxbas_pw_pool%give_back_pw(v_hewd_rspace)
CALL auxbas_pw_pool%give_back_pw(rho_tot_ewd_g)
CALL auxbas_pw_pool%give_back_pw(rho_tot_ewd_r)
DEALLOCATE (atecc)
DEALLOCATE (alpha, ccore)
CALL timestop(handle)
END SUBROUTINE vh_ewald_correction
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
@ -866,7 +1083,7 @@ CONTAINS
IF (iounit > 0) THEN
natom = SIZE(atener)
WRITE (UNIT=iounit, FMT="(/,T2,A)") TRIM(label)
WRITE (UNIT=iounit, FMT="(T2,A)") TRIM(label)
WRITE (UNIT=iounit, FMT="(5G16.8)") atener(1:natom)
END IF

2
src/fm/cp_fm_basic_linalg.F
View file

@ -2304,11 +2304,11 @@ CONTAINS
CALL pdorgqr(nrow, nrow, nrow, a, first_row, first_col, descc, tau, work, lwork, info)
#else
CPASSERT(first_row == 1 .AND. first_col == 1)
lda = SIZE(a, 1)
CALL dorgqr(nrow, nrow, nrow, a, lda, tau, work, lwork, info)
lwork = INT(work(1))
DEALLOCATE (work)
ALLOCATE (work(lwork))
lda = SIZE(a, 1)
CALL dorgqr(nrow, nrow, nrow, a, lda, tau, work, lwork, info)
#endif
CPASSERT(INFO == 0)

7
src/input_cp2k_dft.F
View file

@ -1828,6 +1828,13 @@ CONTAINS
default_l_val=.FALSE., lone_keyword_l_val=.TRUE., type_of_var=logical_t)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="EWALD_ALPHA_PARAMETER", &
description="Calculate Energy Decomposition for a specific alpha value. "// &
"alpha = 1/(2*rc**2), see GTH pseudopotentials.", &
usage="EWALD_ALPHA_PARAMETER alpha", repeats=.FALSE., n_var=1, &
default_r_val=0.0_dp, type_of_var=real_t)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL section_add_subsection(section, print_key)
CALL section_release(print_key)

17
src/qs_collocate_density.F
View file

@ -782,13 +782,16 @@ CONTAINS
!> \param rho_core ...
!> \param total_rho ...
!> \param qs_env ...
!> \param calpha ...
!> \param ccore ...
!> \param only_nopaw ...
! **************************************************************************************************
SUBROUTINE calculate_rho_core(rho_core, total_rho, qs_env, only_nopaw)
SUBROUTINE calculate_rho_core(rho_core, total_rho, qs_env, calpha, ccore, only_nopaw)
TYPE(pw_type), INTENT(INOUT) :: rho_core
REAL(KIND=dp), INTENT(OUT) :: total_rho
TYPE(qs_environment_type), POINTER :: qs_env
REAL(KIND=dp), DIMENSION(:), OPTIONAL :: calpha, ccore
LOGICAL, INTENT(IN), OPTIONAL :: only_nopaw
CHARACTER(len=*), PARAMETER :: routineN = 'calculate_rho_core'
@ -818,6 +821,9 @@ CONTAINS
my_only_nopaw = .FALSE.
IF (PRESENT(only_nopaw)) my_only_nopaw = only_nopaw
IF (PRESENT(calpha)) THEN
CPASSERT(PRESENT(ccore))
END IF
CALL get_qs_env(qs_env=qs_env, &
atomic_kind_set=atomic_kind_set, &
@ -835,8 +841,13 @@ CONTAINS
DO ikind = 1, SIZE(atomic_kind_set)
CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom, &
alpha_core_charge=alpha, ccore_charge=pab(1, 1))
IF (PRESENT(calpha)) THEN
alpha = calpha(ikind)
pab(1, 1) = ccore(ikind)
ELSE
CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom, &
alpha_core_charge=alpha, ccore_charge=pab(1, 1))
END IF
IF (my_only_nopaw .AND. paw_atom) CYCLE
IF (alpha == 0.0_dp .OR. pab(1, 1) == 0.0_dp) CYCLE

71
src/qs_core_energies.F
View file

@ -17,6 +17,8 @@ MODULE qs_core_energies
get_atomic_kind_set
USE atprop_types, ONLY: atprop_array_init,&
atprop_type
USE cell_types, ONLY: cell_type,&
pbc
USE dbcsr_api, ONLY: dbcsr_dot,&
dbcsr_p_type,&
dbcsr_type
@ -53,7 +55,7 @@ MODULE qs_core_energies
PUBLIC :: calculate_ptrace, &
calculate_ecore_overlap, &
calculate_ecore_self
calculate_ecore_self, calculate_ecore_alpha
INTERFACE calculate_ptrace
MODULE PROCEDURE calculate_ptrace_1, calculate_ptrace_gamma, calculate_ptrace_kp
@ -428,4 +430,71 @@ CONTAINS
END SUBROUTINE calculate_ecore_self
! **************************************************************************************************
!> \brief Calculate the overlap and self energy of the core charge distribution for a given alpha
!> Use a minimum image convention and double loop over all atoms
!> \param qs_env ...
!> \param alpha ...
!> \param atecc ...
!> \author JGH
!> \version 1.0
! **************************************************************************************************
SUBROUTINE calculate_ecore_alpha(qs_env, alpha, atecc)
TYPE(qs_environment_type), POINTER :: qs_env
REAL(KIND=dp), INTENT(IN) :: alpha
REAL(KIND=dp), DIMENSION(:) :: atecc
CHARACTER(len=*), PARAMETER :: routineN = 'calculate_ecore_alpha'
INTEGER :: handle, iatom, ikind, jatom, jkind, &
natom, nkind
INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
REAL(KIND=dp) :: dab, eab, fab, rootaab, zab
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: zeff
REAL(KIND=dp), DIMENSION(3) :: rab
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cell_type), POINTER :: cell
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
CALL timeset(routineN, handle)
CALL get_qs_env(qs_env=qs_env, &
cell=cell, &
atomic_kind_set=atomic_kind_set, &
qs_kind_set=qs_kind_set, &
particle_set=particle_set)
CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
!
nkind = SIZE(atomic_kind_set)
natom = SIZE(particle_set)
ALLOCATE (zeff(nkind))
zeff(:) = 0.0_dp
DO ikind = 1, nkind
CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff(ikind))
END DO
rootaab = SQRT(0.5_dp*alpha)
DO iatom = 1, natom
ikind = kind_of(iatom)
atecc(iatom) = atecc(iatom) - zeff(ikind)**2*SQRT(alpha/twopi)
DO jatom = iatom + 1, natom
jkind = kind_of(jatom)
zab = zeff(ikind)*zeff(jkind)
fab = 2.0_dp*oorootpi*zab*rootaab
rab = particle_set(iatom)%r - particle_set(jatom)%r
rab = pbc(rab, cell)
dab = SQRT(SUM(rab(:)**2))
eab = zab*erfc(rootaab*dab)/dab
atecc(iatom) = atecc(iatom) + 0.5_dp*eab
atecc(jatom) = atecc(jatom) + 0.5_dp*eab
END DO
END DO
DEALLOCATE (zeff)
CALL timestop(handle)
END SUBROUTINE calculate_ecore_alpha
END MODULE qs_core_energies

2
src/qs_integrate_potential.F
View file

@ -26,6 +26,7 @@ MODULE qs_integrate_potential
USE qs_integrate_potential_single, ONLY: integrate_ppl_rspace,&
integrate_rho_nlcc,&
integrate_v_core_rspace,&
integrate_v_gaussian_rspace,&
integrate_v_rspace_diagonal,&
integrate_v_rspace_one_center
#include "./base/base_uses.f90"
@ -46,6 +47,7 @@ MODULE qs_integrate_potential
PUBLIC :: integrate_v_rspace_one_center, &
integrate_v_rspace_diagonal, &
integrate_v_core_rspace, &
integrate_v_gaussian_rspace, &
integrate_ppl_rspace, &
integrate_rho_nlcc

115
src/qs_integrate_potential_single.F
View file

@ -83,6 +83,7 @@ MODULE qs_integrate_potential_single
PUBLIC :: integrate_v_rspace_one_center, &
integrate_v_rspace_diagonal, &
integrate_v_core_rspace, &
integrate_v_gaussian_rspace, &
integrate_ppl_rspace, &
integrate_rho_nlcc
@ -633,6 +634,120 @@ CONTAINS
END SUBROUTINE integrate_v_core_rspace
! **************************************************************************************************
!> \brief computes the overlap of a set of Gaussians with a potential on grid
!> \param v_rspace ...
!> \param qs_env ...
!> \param alpha ...
!> \param ccore ...
!> \param atecc ...
! **************************************************************************************************
SUBROUTINE integrate_v_gaussian_rspace(v_rspace, qs_env, alpha, ccore, atecc)
TYPE(pw_type), INTENT(IN) :: v_rspace
TYPE(qs_environment_type), POINTER :: qs_env
REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: alpha, ccore
REAL(KIND=dp), DIMENSION(:) :: atecc
CHARACTER(len=*), PARAMETER :: routineN = 'integrate_v_gaussian_rspace'
INTEGER :: atom_a, handle, iatom, ikind, j, natom, &
natom_of_kind, npme
INTEGER, DIMENSION(:), POINTER :: atom_list, cores
REAL(KIND=dp) :: alpha_core_charge, eps_rho_rspace, radius
REAL(KIND=dp), DIMENSION(3) :: ra
REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, pab
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cell_type), POINTER :: cell
TYPE(dft_control_type), POINTER :: dft_control
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(pw_env_type), POINTER :: pw_env
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(realspace_grid_type), POINTER :: rs_v
CALL timeset(routineN, handle)
CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
!If gapw, check for gpw kinds
CPASSERT(.NOT. dft_control%qs_control%gapw)
NULLIFY (pw_env)
ALLOCATE (cores(1))
ALLOCATE (hab(1, 1))
ALLOCATE (pab(1, 1))
CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
CALL pw_env_get(pw_env=pw_env, auxbas_rs_grid=rs_v)
CALL transfer_pw2rs(rs_v, v_rspace)
CALL get_qs_env(qs_env=qs_env, &
atomic_kind_set=atomic_kind_set, &
qs_kind_set=qs_kind_set, &
cell=cell, &
dft_control=dft_control, &
particle_set=particle_set, &
pw_env=pw_env)
! atomic energy contributions
natom = SIZE(particle_set)
eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
DO ikind = 1, SIZE(atomic_kind_set)
CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_of_kind, atom_list=atom_list)
pab(1, 1) = -ccore(ikind)
alpha_core_charge = alpha(ikind)
IF (alpha_core_charge == 0.0_dp .OR. pab(1, 1) == 0.0_dp) CYCLE
CALL reallocate(cores, 1, natom_of_kind)
npme = 0
cores = 0
DO iatom = 1, natom_of_kind
atom_a = atom_list(iatom)
ra(:) = pbc(particle_set(atom_a)%r, cell)
IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN
! replicated realspace grid, split the atoms up between procs
IF (MODULO(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos) THEN
npme = npme + 1
cores(npme) = iatom
END IF
ELSE
npme = npme + 1
cores(npme) = iatom
END IF
END DO
DO j = 1, npme
iatom = cores(j)
atom_a = atom_list(iatom)
ra(:) = pbc(particle_set(atom_a)%r, cell)
hab(1, 1) = 0.0_dp
radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
ra=ra, rb=ra, rp=ra, &
zetp=alpha_core_charge, eps=eps_rho_rspace, &
pab=pab, o1=0, o2=0, & ! without map_consistent
prefactor=1.0_dp, cutoff=1.0_dp)
CALL integrate_pgf_product(0, alpha_core_charge, 0, &
0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), &
rs_v, hab, pab=pab, o1=0, o2=0, &
radius=radius, calculate_forces=.FALSE., &
use_subpatch=.TRUE., subpatch_pattern=0)
atecc(atom_a) = atecc(atom_a) + 0.5_dp*hab(1, 1)*pab(1, 1)
END DO
END DO
DEALLOCATE (hab, pab, cores)
CALL timestop(handle)
END SUBROUTINE integrate_v_gaussian_rspace
! **************************************************************************************************
!> \brief computes integrals of product of v_rspace times a one-center function
!> required for LRIGPW
!> \param v_rspace ...

55
tests/QS/regtest-ed/H2O_ewald.inp Normal file
View file

@ -0,0 +1,55 @@
&FORCE_EVAL
METHOD Quickstep
&PROPERTIES
&ATOMIC
ENERGY T
&END
&END PROPERTIES
&DFT
BASIS_SET_FILE_NAME BASIS_SET
POTENTIAL_FILE_NAME GTH_POTENTIALS
&MGRID
CUTOFF 200
&END MGRID
&QS
EPS_DEFAULT 1.E-10
&END QS
&SCF
EPS_SCF 1.0E-8
SCF_GUESS ATOMIC
&END
&XC
&XC_FUNCTIONAL PBE
&END
&END XC
&PRINT
&ENERGY_DECOMPOSITION_ANALYSIS
DETAILED_ENERGY T
EWALD_ALPHA_PARAMETER 5.0
&END ENERGY_DECOMPOSITION_ANALYSIS
&END PRINT
&END DFT
&SUBSYS
&CELL
ABC 5.0 5.0 5.0
&END CELL
&COORD
O 0.000000 0.000000 -0.065587
H 0.000000 -0.757136 0.520545
H 0.000000 0.757136 0.520545
&END COORD
&KIND H
BASIS_SET ORB DZVP-GTH-BLYP
POTENTIAL GTH-PADE-q1
&END KIND
&KIND O
BASIS_SET ORB DZVP-GTH-BLYP
POTENTIAL GTH-PADE-q6
&END KIND
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PROJECT wat
RUN_TYPE ENERGY
PRINT_LEVEL LOW
&END GLOBAL

1
tests/QS/regtest-ed/TEST_FILES
View file

@ -22,4 +22,5 @@ si8_pmix_at.inp 0
H2O_detail.inp 0
H2O_canon.inp 0
H2O_hyb_detail.inp 0
H2O_ewald.inp 0
#EOF