From f5d7817782ccd5e60ecd1b01248fcc43fdeed73e Mon Sep 17 00:00:00 2001 From: Juerg Hutter Date: Fri, 8 Dec 2023 12:01:04 +0100 Subject: [PATCH] Energy decomposition:wq (#3146) --- src/core_ae.F | 276 +++++++++++++++++++++++++++- src/ed_analysis.F | 255 +++++++++++++++++++++++-- src/fm/cp_fm_basic_linalg.F | 2 +- src/input_cp2k_dft.F | 7 + src/qs_collocate_density.F | 17 +- src/qs_core_energies.F | 71 ++++++- src/qs_integrate_potential.F | 2 + src/qs_integrate_potential_single.F | 115 ++++++++++++ tests/QS/regtest-ed/H2O_ewald.inp | 55 ++++++ tests/QS/regtest-ed/TEST_FILES | 1 + 10 files changed, 776 insertions(+), 25 deletions(-) create mode 100644 tests/QS/regtest-ed/H2O_ewald.inp diff --git a/src/core_ae.F b/src/core_ae.F index 89618ed9e8..c3191805f1 100644 --- a/src/core_ae.F +++ b/src/core_ae.F @@ -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 diff --git a/src/ed_analysis.F b/src/ed_analysis.F index abd8b09608..c7e2d24214 100644 --- a/src/ed_analysis.F +++ b/src/ed_analysis.F @@ -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 diff --git a/src/fm/cp_fm_basic_linalg.F b/src/fm/cp_fm_basic_linalg.F index 6e6be57d56..bf053a86b6 100644 --- a/src/fm/cp_fm_basic_linalg.F +++ b/src/fm/cp_fm_basic_linalg.F @@ -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) diff --git a/src/input_cp2k_dft.F b/src/input_cp2k_dft.F index 3c53df61b9..1425f658ea 100644 --- a/src/input_cp2k_dft.F +++ b/src/input_cp2k_dft.F @@ -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) diff --git a/src/qs_collocate_density.F b/src/qs_collocate_density.F index 70f060ba69..2cacf3afad 100644 --- a/src/qs_collocate_density.F +++ b/src/qs_collocate_density.F @@ -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 diff --git a/src/qs_core_energies.F b/src/qs_core_energies.F index ad23868242..1921368734 100644 --- a/src/qs_core_energies.F +++ b/src/qs_core_energies.F @@ -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 diff --git a/src/qs_integrate_potential.F b/src/qs_integrate_potential.F index dfcf085a3a..6c1cbe6e41 100644 --- a/src/qs_integrate_potential.F +++ b/src/qs_integrate_potential.F @@ -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 diff --git a/src/qs_integrate_potential_single.F b/src/qs_integrate_potential_single.F index 818837d64b..115aaea3d8 100644 --- a/src/qs_integrate_potential_single.F +++ b/src/qs_integrate_potential_single.F @@ -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 ... diff --git a/tests/QS/regtest-ed/H2O_ewald.inp b/tests/QS/regtest-ed/H2O_ewald.inp new file mode 100644 index 0000000000..049d746e2c --- /dev/null +++ b/tests/QS/regtest-ed/H2O_ewald.inp @@ -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 diff --git a/tests/QS/regtest-ed/TEST_FILES b/tests/QS/regtest-ed/TEST_FILES index 62a7b6077c..3d9fb16748 100644 --- a/tests/QS/regtest-ed/TEST_FILES +++ b/tests/QS/regtest-ed/TEST_FILES @@ -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