From 14d6109598386049dc043f01eb4bd5fdc751d8e4 Mon Sep 17 00:00:00 2001 From: Juerg Hutter Date: Tue, 18 Jun 2024 16:13:51 +0200 Subject: [PATCH] Fix detailed energies and some small parallel bugs (#3486) --- src/core_ae.F | 101 +++++++++++++++++++++++++++++++++++++++------- src/core_ppl.F | 2 +- src/ed_analysis.F | 40 ++++++++++++------ 3 files changed, 116 insertions(+), 27 deletions(-) diff --git a/src/core_ae.F b/src/core_ae.F index 61034f3b9d..aa69a2d8ea 100644 --- a/src/core_ae.F +++ b/src/core_ae.F @@ -445,7 +445,7 @@ SUBROUTINE build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, us END DO DEALLOCATE (hab, work, verf, vnuc, ff) - IF (calculate_forces) THEN + IF (calculate_forces .OR. doat) THEN DEALLOCATE (pab) END IF @@ -465,7 +465,7 @@ SUBROUTINE build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, us DEALLOCATE (basis_set_list) - IF (calculate_forces) THEN + IF (calculate_forces .OR. doat) THEN ! *** If LSD, then recover alpha density and beta density *** ! *** from the total density (1) and the spin density (2) *** IF (SIZE(matrix_p, 1) == 2) THEN @@ -591,6 +591,7 @@ END SUBROUTINE verfc_force ! ************************************************************************************************** !> \brief Integrals = -Z*erfc(a*r)/r !> \param matrix_h ... +!> \param matrix_p ... !> \param qs_kind_set ... !> \param atomic_kind_set ... !> \param particle_set ... @@ -599,11 +600,13 @@ END SUBROUTINE verfc_force !> \param eps_core_charge ... !> \param sab_orb ... !> \param sac_ae ... +!> \param atcore ... ! ************************************************************************************************** - SUBROUTINE build_erfc(matrix_h, qs_kind_set, atomic_kind_set, particle_set, & - calpha, ccore, eps_core_charge, sab_orb, sac_ae) + SUBROUTINE build_erfc(matrix_h, matrix_p, qs_kind_set, atomic_kind_set, particle_set, & + calpha, ccore, eps_core_charge, sab_orb, sac_ae, atcore) TYPE(dbcsr_p_type) :: matrix_h + TYPE(dbcsr_p_type), DIMENSION(:) :: matrix_p 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 @@ -611,6 +614,8 @@ SUBROUTINE build_erfc(matrix_h, qs_kind_set, atomic_kind_set, particle_set, & REAL(KIND=dp), INTENT(IN) :: eps_core_charge TYPE(neighbor_list_set_p_type), DIMENSION(:), & POINTER :: sab_orb, sac_ae + REAL(KIND=dp), DIMENSION(:), INTENT(INOUT), & + OPTIONAL :: atcore CHARACTER(LEN=*), PARAMETER :: routineN = 'build_erfc' @@ -621,9 +626,10 @@ SUBROUTINE build_erfc(matrix_h, qs_kind_set, atomic_kind_set, particle_set, & 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 + LOGICAL :: doat, found + REAL(KIND=dp) :: alpha_c, atk0, atk1, 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 @@ -631,11 +637,12 @@ SUBROUTINE build_erfc(matrix_h, qs_kind_set, atomic_kind_set, particle_set, & 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(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b + TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list + TYPE(all_potential_type), POINTER :: all_potential + REAL(KIND=dp), DIMENSION(SIZE(particle_set)) :: at_thread TYPE(sgp_potential_type), POINTER :: sgp_potential !$ INTEGER(kind=omp_lock_kind), & @@ -651,6 +658,19 @@ SUBROUTINE build_erfc(matrix_h, qs_kind_set, atomic_kind_set, particle_set, & nkind = SIZE(atomic_kind_set) natom = SIZE(particle_set) + doat = PRESENT(atcore) + + IF (doat) THEN + IF (SIZE(matrix_p, 1) == 2) THEN + CALL dbcsr_add(matrix_p(1)%matrix, matrix_p(2)%matrix, & + alpha_scalar=1.0_dp, beta_scalar=1.0_dp) + CALL dbcsr_add(matrix_p(2)%matrix, matrix_p(1)%matrix, & + alpha_scalar=-2.0_dp, beta_scalar=1.0_dp) + END IF + END IF + + at_thread = 0.0_dp + ALLOCATE (basis_set_list(nkind)) DO ikind = 1, nkind CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set_a) @@ -676,9 +696,9 @@ SUBROUTINE build_erfc(matrix_h, qs_kind_set, atomic_kind_set, particle_set, & !$OMP PARALLEL & !$OMP DEFAULT (NONE) & !$OMP SHARED (ap_iterator, basis_set_list, & -!$OMP matrix_h, atomic_kind_set, qs_kind_set, particle_set, & +!$OMP matrix_h, matrix_p, 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 slot, ldsab, maxnset, ldai, maxl, maxco, doat, 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, & @@ -686,8 +706,9 @@ SUBROUTINE build_erfc(matrix_h, qs_kind_set, atomic_kind_set, particle_set, & !$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 atk0, atk1, habd, f0, katom, cellind, img, nij, ff, & +!$OMP sgp_potential, all_potential, hash, hash1, hash2, iatom8) & +!$OMP REDUCTION (+ : at_thread ) !$OMP SINGLE !$ ALLOCATE (locks(nlock)) @@ -704,6 +725,9 @@ SUBROUTINE build_erfc(matrix_h, qs_kind_set, atomic_kind_set, particle_set, & 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)) + IF (doat) THEN + ALLOCATE (pab(maxco, maxco, maxnset*maxnset)) + END IF !$OMP DO SCHEDULE(GUIDED) DO slot = 1, sab_orb(1)%nl_size @@ -765,6 +789,32 @@ SUBROUTINE build_erfc(matrix_h, qs_kind_set, atomic_kind_set, particle_set, & NULLIFY (h_block) CALL dbcsr_get_block_p(matrix=matrix_h%matrix, & row=irow, col=icol, BLOCK=h_block, found=found) + IF (doat) THEN + NULLIFY (p_block) + CALL dbcsr_get_block_p(matrix=matrix_p(1)%matrix, & + row=irow, col=icol, BLOCK=p_block, found=found) + CPASSERT(ASSOCIATED(p_block)) + ! *** Decontract density matrix block *** + 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 + ! *** Decontract density matrix block *** + IF (iatom <= jatom) THEN + work(1:ncoa, 1:nsgfb(jset)) = MATMUL(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1), & + p_block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1)) + ELSE + work(1:ncoa, 1:nsgfb(jset)) = MATMUL(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1), & + TRANSPOSE(p_block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1))) + END IF + pab(1:ncoa, 1:ncob, nij) = MATMUL(work(1:ncoa, 1:nsgfb(jset)), & + TRANSPOSE(sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1))) + END DO + END DO + END IF ! loop over all kinds of atoms hab = 0._dp @@ -801,12 +851,20 @@ SUBROUTINE build_erfc(matrix_h, qs_kind_set, atomic_kind_set, particle_set, & rac2 = dac*dac rbc2 = dbc*dbc nij = jset + (iset - 1)*maxnset + IF (doat) THEN + atk0 = f0*SUM(hab(1:ncoa, 1:ncob, nij)*pab(1:ncoa, 1:ncob, nij)) + END IF ! 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:)) + ! + IF (doat) THEN + atk1 = f0*SUM(hab(1:ncoa, 1:ncob, nij)*pab(1:ncoa, 1:ncob, nij)) + at_thread(katom) = at_thread(katom) + (atk1 - atk0) + END IF END DO END DO END DO @@ -857,6 +915,21 @@ SUBROUTINE build_erfc(matrix_h, qs_kind_set, atomic_kind_set, particle_set, & DEALLOCATE (basis_set_list) + IF (doat) THEN + ! *** If LSD, then recover alpha density and beta density *** + ! *** from the total density (1) and the spin density (2) *** + IF (SIZE(matrix_p, 1) == 2) THEN + CALL dbcsr_add(matrix_p(1)%matrix, matrix_p(2)%matrix, & + alpha_scalar=0.5_dp, beta_scalar=0.5_dp) + CALL dbcsr_add(matrix_p(2)%matrix, matrix_p(1)%matrix, & + alpha_scalar=-1.0_dp, beta_scalar=1.0_dp) + END IF + END IF + + IF (doat) THEN + atcore(1:natom) = atcore(1:natom) + at_thread(1:natom) + END IF + CALL timestop(handle) END SUBROUTINE build_erfc diff --git a/src/core_ppl.F b/src/core_ppl.F index 061df7b0e3..00671df530 100644 --- a/src/core_ppl.F +++ b/src/core_ppl.F @@ -688,7 +688,7 @@ SUBROUTINE build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, u END DO ! slot DEALLOCATE (hab, work, ppl_work) - IF (calculate_forces) THEN + IF (calculate_forces .OR. doat) THEN DEALLOCATE (pab, ppl_fwork) END IF diff --git a/src/ed_analysis.F b/src/ed_analysis.F index a5eb85b736..162c7932f1 100644 --- a/src/ed_analysis.F +++ b/src/ed_analysis.F @@ -437,6 +437,8 @@ SUBROUTINE edmf_analysis(qs_env, input_section, unit_nr) CALL cp_fm_release(cvec) CALL cp_fm_release(cvec2) + CALL get_qs_env(qs_env, para_env=para_env) + group = para_env ! energy arrays ALLOCATE (atener(natom), ateks(natom), atecc(natom), atcore(natom)) ALLOCATE (ate1xc(natom), ate1h(natom), atewald(natom)) @@ -489,6 +491,7 @@ SUBROUTINE edmf_analysis(qs_env, input_section, unit_nr) math(1:1, 1:1) => core_mat(1:1) matp(1:nspin, 1:1) => matrix_p(1:nspin) CALL core_matrices(qs_env, math, matp, .FALSE., 0, atcore=atcore) + CALL group%sum(atcore) ! IF (ewald_correction) THEN ALLOCATE (dve_mat%matrix) @@ -496,6 +499,7 @@ SUBROUTINE edmf_analysis(qs_env, input_section, unit_nr) 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) + CALL group%sum(atewald) END IF ! DO ispin = 1, nspin @@ -620,12 +624,9 @@ SUBROUTINE edmf_analysis(qs_env, input_section, unit_nr) CALL cp_fm_release(fij_mat) END IF - CALL get_qs_env(qs_env, para_env=para_env) - group = para_env ! KS energy atener(1:natom) = ateks(1:natom) ! 1/2 of VPP contribution Tr[VPP(K)*P] - CALL group%sum(atcore) atener(1:natom) = atener(1:natom) + 0.5_dp*atcore(1:natom) ! core energy corrections CALL group%sum(atecc) @@ -932,8 +933,10 @@ SUBROUTINE vh_ewald_correction(qs_env, ealpha, vh_mat, atewd) INTEGER :: handle, ikind, ispin, natom, nkind REAL(KIND=dp) :: eps_core_charge, rhotot, zeff - REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: alpha, atecc, ccore + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: alpha, atcore, atecc, ccore TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set + TYPE(dbcsr_p_type) :: e_mat + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p TYPE(dft_control_type), POINTER :: dft_control TYPE(neighbor_list_set_p_type), DIMENSION(:), & POINTER :: sab_orb, sac_ae @@ -951,7 +954,7 @@ SUBROUTINE vh_ewald_correction(qs_env, ealpha, vh_mat, atewd) CALL timeset(routineN, handle) natom = SIZE(atewd) - ALLOCATE (atecc(natom)) + ALLOCATE (atecc(natom), atcore(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)) @@ -999,25 +1002,38 @@ SUBROUTINE vh_ewald_correction(qs_env, ealpha, vh_mat, atewd) CALL dbcsr_scale(vh_mat%matrix, 0.5_dp) ! ! delta erfc + CALL qs_rho_get(rho, rho_ao=matrix_p) eps_core_charge = dft_control%qs_control%eps_core_charge + ALLOCATE (e_mat%matrix) + CALL dbcsr_create(e_mat%matrix, template=vh_mat%matrix) + CALL dbcsr_copy(e_mat%matrix, vh_mat%matrix) + ! 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) + CALL dbcsr_set(e_mat%matrix, 0.0_dp) + atcore = 0.0_dp + CALL build_erfc(e_mat, matrix_p, qs_kind_set, atomic_kind_set, particle_set, & + alpha, ccore, eps_core_charge, sab_orb, sac_ae, atcore) + atewd(1:natom) = atewd(1:natom) + 0.5_dp*atcore(1:natom) + CALL dbcsr_add(vh_mat%matrix, e_mat%matrix, 1.0_dp, 0.5_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 dbcsr_set(e_mat%matrix, 0.0_dp) + atcore = 0.0_dp + CALL build_erfc(e_mat, matrix_p, qs_kind_set, atomic_kind_set, particle_set, & + alpha, ccore, eps_core_charge, sab_orb, sac_ae, atcore) + atewd(1:natom) = atewd(1:natom) - 0.5_dp*atcore(1:natom) + CALL dbcsr_add(vh_mat%matrix, e_mat%matrix, 1.0_dp, -0.5_dp) ! + CALL dbcsr_release(e_mat%matrix) + DEALLOCATE (e_mat%matrix) 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 (atecc, atcore) DEALLOCATE (alpha, ccore) CALL timestop(handle)