Skip to content

Commit

Permalink
Fix detailed energies and some small parallel bugs (cp2k#3486)
Browse files Browse the repository at this point in the history
  • Loading branch information
@juerghutter
juerghutter authored Jun 18, 2024
1 parent b95f8ed commit 14d6109
Show file tree
Hide file tree
Showing 3 changed files with 116 additions and 27 deletions.
101 changes: 87 additions & 14 deletions src/core_ae.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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
Expand Down Expand Up @@ -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 ...
Expand All @@ -599,18 +600,22 @@ 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
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
REAL(KIND=dp), DIMENSION(:), INTENT(INOUT), &
OPTIONAL :: atcore

CHARACTER(LEN=*), PARAMETER :: routineN = 'build_erfc'

Expand All @@ -621,21 +626,23 @@ 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
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(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), &
Expand All @@ -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)
Expand All @@ -676,18 +696,19 @@ 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, &
!$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 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))
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion src/core_ppl.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
40 changes: 28 additions & 12 deletions src/ed_analysis.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand Down Expand Up @@ -489,13 +491,15 @@ 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)
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)
CALL group%sum(atewald)
END IF
!
DO ispin = 1, nspin
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -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))
Expand Down Expand Up @@ -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)
Expand Down

0 comments on commit 14d6109

Please sign in to comment.