Skip to content

Commit

Permalink
ED analysis : split PP term according to original proposal for core p… (
Browse files Browse the repository at this point in the history 
cp2k#3452)
  • Loading branch information
@juerghutter
juerghutter authored Jun 4, 2024
1 parent 5f00d0c commit 7d74e28
Show file tree
Hide file tree
Showing 5 changed files with 141 additions and 48 deletions.
39 changes: 29 additions & 10 deletions src/core_ae.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -78,9 +78,11 @@ MODULE core_ae
!> \param sac_ae ...
!> \param nimages ...
!> \param cell_to_index ...
!> \param atcore ...
! **************************************************************************************************
SUBROUTINE build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, nimages, cell_to_index)
qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, &
nimages, cell_to_index, atcore)

TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p
TYPE(qs_force_type), DIMENSION(:), POINTER :: force
Expand All @@ -95,6 +97,8 @@ SUBROUTINE build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, us
POINTER :: sab_orb, sac_ae
INTEGER, INTENT(IN) :: nimages
INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
REAL(KIND=dp), DIMENSION(:), INTENT(INOUT), &
OPTIONAL :: atcore

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

Expand All @@ -106,9 +110,10 @@ SUBROUTINE build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, us
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, dokp, 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
Expand All @@ -119,6 +124,7 @@ SUBROUTINE build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, us
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
REAL(KIND=dp), DIMENSION(:, :), POINTER :: h_block, p_block, rpgfa, rpgfb, sphi_a, &
sphi_b, zeta, zetb
REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
Expand All @@ -142,9 +148,10 @@ SUBROUTINE build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, us
nkind = SIZE(atomic_kind_set)
natom = SIZE(particle_set)

doat = PRESENT(atcore)
dokp = (nimages > 1)

IF (calculate_forces) THEN
IF (calculate_forces .OR. doat) THEN
IF (SIZE(matrix_p, 1) == 2) THEN
DO img = 1, nimages
CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
Expand All @@ -156,6 +163,7 @@ SUBROUTINE build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, us
END IF

force_thread = 0.0_dp
at_thread = 0.0_dp
pv_thread = 0.0_dp

ALLOCATE (basis_set_list(nkind))
Expand Down Expand Up @@ -185,17 +193,17 @@ SUBROUTINE build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, us
!$OMP SHARED (ap_iterator, basis_set_list, calculate_forces, use_virial, &
!$OMP matrix_h, matrix_p, atomic_kind_set, qs_kind_set, particle_set, &
!$OMP sab_orb, sac_ae, nthread, ncoset, nkind, cell_to_index, &
!$OMP slot, ldsab, maxnset, ldai, nder, maxl, maxco, dokp, locks, natom) &
!$OMP slot, ldsab, maxnset, ldai, nder, maxl, maxco, dokp, 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, force_a, force_b, mepos, &
!$OMP habd, f0, katom, cellind, img, nij, ff, &
!$OMP atk0, atk1, habd, f0, katom, cellind, img, nij, ff, &
!$OMP sgp_potential, all_potential, hash, hash1, hash2, iatom8) &
!$OMP REDUCTION (+ : pv_thread, force_thread )
!$OMP REDUCTION (+ : pv_thread, force_thread, at_thread )

!$OMP SINGLE
!$ ALLOCATE (locks(nlock))
Expand All @@ -212,7 +220,7 @@ SUBROUTINE build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, us

ALLOCATE (hab(ldsab, ldsab, maxnset*maxnset), work(ldsab, ldsab))
ALLOCATE (verf(ldai, ldai, 2*maxl + nder + 1), vnuc(ldai, ldai, 2*maxl + nder + 1), ff(0:2*maxl + nder))
IF (calculate_forces) THEN
IF (calculate_forces .OR. doat) THEN
ALLOCATE (pab(maxco, maxco, maxnset*maxnset))
END IF

Expand Down Expand Up @@ -281,7 +289,7 @@ SUBROUTINE build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, us
NULLIFY (h_block)
CALL dbcsr_get_block_p(matrix=matrix_h(1, img)%matrix, &
row=irow, col=icol, BLOCK=h_block, found=found)
IF (calculate_forces) THEN
IF (calculate_forces .OR. doat) THEN
NULLIFY (p_block)
CALL dbcsr_get_block_p(matrix=matrix_p(1, img)%matrix, &
row=irow, col=icol, BLOCK=p_block, found=found)
Expand Down Expand Up @@ -352,6 +360,9 @@ SUBROUTINE build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, us
rbc2 = dbc*dbc
nij = jset + (iset - 1)*maxnset
! *** Calculate the GTH pseudo potential forces ***
IF (doat) THEN
atk0 = f0*SUM(hab(1:ncoa, 1:ncob, nij)*pab(1:ncoa, 1:ncob, nij))
END IF
IF (calculate_forces) THEN
na_plus = npgfa(iset)*ncoset(la_max(iset) + nder)
nb_plus = npgfb(jset)*ncoset(lb_max(jset))
Expand Down Expand Up @@ -396,6 +407,11 @@ SUBROUTINE build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, us
alpha_c, core_radius, zeta_c, core_charge, &
rab, rab2, rac, rac2, rbc2, hab(:, :, nij), verf, vnuc, ff(0:))
END IF
! calculate atomic contributions
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 @@ -472,6 +488,9 @@ SUBROUTINE build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, us
END DO
!$OMP END DO
END IF
IF (doat) THEN
atcore(1:natom) = atcore(1:natom) + at_thread(1:natom)
END IF

IF (calculate_forces .AND. use_virial) THEN
virial%pv_ppl = virial%pv_ppl + pv_thread
Expand Down
43 changes: 32 additions & 11 deletions src/core_ppl.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -88,10 +88,11 @@ MODULE core_ppl
!> \param cell_to_index ...
!> \param basis_type ...
!> \param deltaR Weighting factors of the derivatives wrt. nuclear positions
!> \param atcore ...
! **************************************************************************************************
SUBROUTINE build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, &
nimages, cell_to_index, basis_type, deltaR)
nimages, cell_to_index, basis_type, deltaR, atcore)

TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p
TYPE(qs_force_type), DIMENSION(:), POINTER :: force
Expand All @@ -109,6 +110,8 @@ SUBROUTINE build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, u
CHARACTER(LEN=*), INTENT(IN) :: basis_type
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
OPTIONAL :: deltaR
REAL(KIND=dp), DIMENSION(:), INTENT(INOUT), &
OPTIONAL :: atcore

CHARACTER(LEN=*), PARAMETER :: routineN = 'build_core_ppl'
INTEGER, PARAMETER :: nexp_max = 30
Expand All @@ -126,9 +129,11 @@ SUBROUTINE build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, u
nct_lpot, npgfa, npgfb, nsgfa, nsgfb
INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
INTEGER, DIMENSION(nexp_max) :: nct_ppl
LOGICAL :: do_dR, dokp, ecp_local, ecp_semi_local, &
found, lpotextended, only_gaussians
REAL(KIND=dp) :: alpha, dab, dac, dbc, f0, ppl_radius
LOGICAL :: do_dR, doat, dokp, ecp_local, &
ecp_semi_local, found, lpotextended, &
only_gaussians
REAL(KIND=dp) :: alpha, atk0, atk1, dab, dac, dbc, f0, &
ppl_radius
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: qab, work
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: hab2_w, ppl_fwork, ppl_work
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: hab, pab
Expand All @@ -143,6 +148,7 @@ SUBROUTINE build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, u
TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
TYPE(gth_potential_type), POINTER :: gth_potential
REAL(KIND=dp), DIMENSION(SIZE(particle_set)) :: at_thread
REAL(KIND=dp), DIMENSION(nexp_max) :: alpha_ppl
REAL(KIND=dp), DIMENSION(:, :), POINTER :: cval_lpot, h1_1block, h1_2block, &
h1_3block, h_block, p_block, rpgfa, &
Expand All @@ -160,6 +166,7 @@ SUBROUTINE build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, u
!$ INTEGER, PARAMETER :: nlock = 501

do_dR = PRESENT(deltaR)
doat = PRESENT(atcore)
MARK_USED(int_8)

IF (calculate_forces) THEN
Expand All @@ -177,7 +184,7 @@ SUBROUTINE build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, u
CPASSERT(PRESENT(cell_to_index) .AND. ASSOCIATED(cell_to_index))
END IF

IF (calculate_forces) THEN
IF (calculate_forces .OR. doat) THEN
IF (SIZE(matrix_p, 1) == 2) THEN
DO img = 1, nimages
CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
Expand All @@ -188,6 +195,7 @@ SUBROUTINE build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, u
END IF
END IF
force_thread = 0.0_dp
at_thread = 0.0_dp

maxder = ncoset(nder)

Expand Down Expand Up @@ -224,14 +232,14 @@ SUBROUTINE build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, u
!$OMP SHARED (ap_iterator, basis_set_list, calculate_forces, use_virial, &
!$OMP matrix_h, matrix_p, atomic_kind_set, qs_kind_set, particle_set, &
!$OMP sab_orb, sac_ppl, nthread, ncoset, nkind, cell_to_index, &
!$OMP ldsab, maxnset, maxder, do_dR, deltaR, &
!$OMP ldsab, maxnset, maxder, do_dR, deltaR, doat, &
!$OMP maxlgto, nder, 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, dab, irow, icol, h_block, found, iset, ncoa, &
!$OMP sgfa, jset, ncob, sgfb, nsgfb, p_block, work, pab, hab, hab2, hab2_w, qab, &
!$OMP h1_1block, h1_2block, h1_3block, kkind, nseta, &
!$OMP atk0, atk1, h1_1block, h1_2block, h1_3block, kkind, nseta, &
!$OMP gth_potential, sgp_potential, alpha, cexp_ppl, lpotextended, &
!$OMP ppl_radius, nexp_lpot, nexp_ppl, alpha_ppl, alpha_lpot, nct_ppl, &
!$OMP nct_lpot, cval_ppl, cval_lpot, rac, dac, rbc, dbc, &
Expand All @@ -240,7 +248,7 @@ SUBROUTINE build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, u
!$OMP nloc, nrloc, aloc, bloc, n_local, a_local, c_local, &
!$OMP slmax, npot, nrpot, apot, bpot, only_gaussians, &
!$OMP ldai, hash, hash1, hash2, iatom8) &
!$OMP REDUCTION (+ : pv_thread, force_thread )
!$OMP REDUCTION (+ : pv_thread, force_thread, at_thread )

!$OMP SINGLE
!$ ALLOCATE (locks(nlock))
Expand All @@ -258,7 +266,7 @@ SUBROUTINE build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, u
ALLOCATE (hab(ldsab, ldsab, maxnset, maxnset), work(ldsab, ldsab*maxder))
ldai = ncoset(2*maxlgto + 2*nder)
ALLOCATE (ppl_work(ldai, ldai, MAX(maxder, 2*maxlgto + 2*nder + 1)))
IF (calculate_forces) THEN
IF (calculate_forces .OR. doat) THEN
ALLOCATE (pab(maxco, maxco, maxnset, maxnset))
ldai = ncoset(maxlgto)
ALLOCATE (ppl_fwork(ldai, ldai, maxder))
Expand Down Expand Up @@ -350,7 +358,7 @@ SUBROUTINE build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, u

CALL dbcsr_get_block_p(matrix_h(1, img)%matrix, irow, icol, h_block, found)
CPASSERT(found)
IF (calculate_forces) THEN
IF (calculate_forces .OR. doat) THEN
NULLIFY (p_block)
CALL dbcsr_get_block_p(matrix_p(1, img)%matrix, irow, icol, p_block, found)
IF (ASSOCIATED(p_block)) THEN
Expand Down Expand Up @@ -461,6 +469,10 @@ SUBROUTINE build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, u
sgfb = first_sgfb(1, jset)
IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
! *** Calculate the GTH pseudo potential forces ***
IF (doat) THEN
atk0 = f0*SUM(hab(1:ncoa, 1:ncob, iset, jset)* &
pab(1:ncoa, 1:ncob, iset, jset))
END IF
IF (calculate_forces) THEN

force_a(:) = 0.0_dp
Expand Down Expand Up @@ -579,6 +591,12 @@ SUBROUTINE build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, u
!$OMP END CRITICAL(type2)
END IF
END IF
! calculate atomic contributions
IF (doat) THEN
atk1 = f0*SUM(hab(1:ncoa, 1:ncob, iset, jset)* &
pab(1:ncoa, 1:ncob, iset, jset))
at_thread(katom) = at_thread(katom) + (atk1 - atk0)
END IF
END DO
END DO
END DO
Expand Down Expand Up @@ -690,7 +708,7 @@ SUBROUTINE build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, u

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 All @@ -714,6 +732,9 @@ SUBROUTINE build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, u
!$OMP END DO
DEALLOCATE (atom_of_kind, kind_of)
END IF
IF (doat) THEN
atcore(1:natom) = atcore(1:natom) + at_thread(1:natom)
END IF

IF (calculate_forces .AND. use_virial) THEN
virial%pv_ppl = virial%pv_ppl + pv_thread
Expand Down
Loading

0 comments on commit 7d74e28

Please sign in to comment.