mirror of
https://github.com/cp2k/cp2k.git
synced 2026-07-28 22:25:32 -04:00
ED analysis : split PP term according to original proposal for core p… (#3452)
This commit is contained in:
Juerg Hutter
GitHub
parent
5f00d0ca4a
commit
7d74e28acb
5 changed files with 141 additions and 48 deletions
|
|
@ -78,9 +78,11 @@ CONTAINS
|
|||
!> \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
|
||||
|
|
@ -95,6 +97,8 @@ CONTAINS
|
|||
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'
|
||||
|
||||
|
|
@ -106,9 +110,10 @@ CONTAINS
|
|||
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
|
||||
|
|
@ -119,6 +124,7 @@ CONTAINS
|
|||
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
|
||||
|
|
@ -142,9 +148,10 @@ CONTAINS
|
|||
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, &
|
||||
|
|
@ -156,6 +163,7 @@ CONTAINS
|
|||
END IF
|
||||
|
||||
force_thread = 0.0_dp
|
||||
at_thread = 0.0_dp
|
||||
pv_thread = 0.0_dp
|
||||
|
||||
ALLOCATE (basis_set_list(nkind))
|
||||
|
|
@ -185,7 +193,7 @@ CONTAINS
|
|||
!$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, &
|
||||
|
|
@ -193,9 +201,9 @@ CONTAINS
|
|||
!$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))
|
||||
|
|
@ -212,7 +220,7 @@ CONTAINS
|
|||
|
||||
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
|
||||
|
||||
|
|
@ -281,7 +289,7 @@ CONTAINS
|
|||
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)
|
||||
|
|
@ -352,6 +360,9 @@ CONTAINS
|
|||
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))
|
||||
|
|
@ -396,6 +407,11 @@ CONTAINS
|
|||
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
|
||||
|
|
@ -472,6 +488,9 @@ CONTAINS
|
|||
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
|
||||
|
|
|
|||
|
|
@ -88,10 +88,11 @@ CONTAINS
|
|||
!> \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
|
||||
|
|
@ -109,6 +110,8 @@ CONTAINS
|
|||
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
|
||||
|
|
@ -126,9 +129,11 @@ CONTAINS
|
|||
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
|
||||
|
|
@ -143,6 +148,7 @@ CONTAINS
|
|||
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, &
|
||||
|
|
@ -160,6 +166,7 @@ CONTAINS
|
|||
!$ INTEGER, PARAMETER :: nlock = 501
|
||||
|
||||
do_dR = PRESENT(deltaR)
|
||||
doat = PRESENT(atcore)
|
||||
MARK_USED(int_8)
|
||||
|
||||
IF (calculate_forces) THEN
|
||||
|
|
@ -177,7 +184,7 @@ CONTAINS
|
|||
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, &
|
||||
|
|
@ -188,6 +195,7 @@ CONTAINS
|
|||
END IF
|
||||
END IF
|
||||
force_thread = 0.0_dp
|
||||
at_thread = 0.0_dp
|
||||
|
||||
maxder = ncoset(nder)
|
||||
|
||||
|
|
@ -224,14 +232,14 @@ CONTAINS
|
|||
!$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, &
|
||||
|
|
@ -240,7 +248,7 @@ CONTAINS
|
|||
!$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))
|
||||
|
|
@ -258,7 +266,7 @@ CONTAINS
|
|||
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))
|
||||
|
|
@ -350,7 +358,7 @@ CONTAINS
|
|||
|
||||
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
|
||||
|
|
@ -461,6 +469,10 @@ CONTAINS
|
|||
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
|
||||
|
|
@ -579,6 +591,12 @@ CONTAINS
|
|||
!$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
|
||||
|
|
@ -690,7 +708,7 @@ CONTAINS
|
|||
|
||||
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
|
||||
|
|
@ -714,6 +732,9 @@ CONTAINS
|
|||
!$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
|
||||
|
|
|
|||
|
|
@ -82,10 +82,11 @@ CONTAINS
|
|||
!> \param basis_type ...
|
||||
!> \param deltaR Weighting factors of the derivatives wrt. nuclear positions
|
||||
!> \param matrix_l ...
|
||||
!> \param atcore ...
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
|
||||
qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, eps_ppnl, &
|
||||
nimages, cell_to_index, basis_type, deltaR, matrix_l)
|
||||
nimages, cell_to_index, basis_type, deltaR, matrix_l, atcore)
|
||||
|
||||
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p
|
||||
TYPE(qs_force_type), DIMENSION(:), POINTER :: force
|
||||
|
|
@ -106,6 +107,8 @@ CONTAINS
|
|||
OPTIONAL :: deltaR
|
||||
TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
|
||||
POINTER :: matrix_l
|
||||
REAL(KIND=dp), DIMENSION(:), INTENT(INOUT), &
|
||||
OPTIONAL :: atcore
|
||||
|
||||
CHARACTER(LEN=*), PARAMETER :: routineN = 'build_core_ppnl'
|
||||
|
||||
|
|
@ -119,9 +122,9 @@ CONTAINS
|
|||
INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nprj_ppnl, &
|
||||
nsgf_seta
|
||||
INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
|
||||
LOGICAL :: do_dR, do_gth, do_kp, do_soc, found, &
|
||||
ppnl_present
|
||||
REAL(KIND=dp) :: dac, f0, ppnl_radius
|
||||
LOGICAL :: do_dR, do_gth, do_kp, do_soc, doat, &
|
||||
found, ppnl_present
|
||||
REAL(KIND=dp) :: atk, dac, f0, ppnl_radius
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: radp
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: sab, work
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ai_work, lab, work_l
|
||||
|
|
@ -134,6 +137,7 @@ CONTAINS
|
|||
TYPE(gth_potential_p_type), DIMENSION(:), POINTER :: gpotential
|
||||
TYPE(clist_type), POINTER :: clist
|
||||
TYPE(alist_type), POINTER :: alist_ac, alist_bc
|
||||
REAL(KIND=dp), DIMENSION(SIZE(particle_set)) :: at_thread
|
||||
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: achint, acint, alkint, bchint, bcint, &
|
||||
blkint
|
||||
REAL(KIND=dp), DIMENSION(:, :), POINTER :: cprj, h_block, l_block_x, l_block_y, &
|
||||
|
|
@ -156,6 +160,8 @@ CONTAINS
|
|||
|
||||
do_dR = .FALSE.
|
||||
IF (PRESENT(deltaR)) do_dR = .TRUE.
|
||||
doat = .FALSE.
|
||||
IF (PRESENT(atcore)) doat = .TRUE.
|
||||
|
||||
IF (calculate_forces) THEN
|
||||
CALL timeset(routineN//"_forces", handle)
|
||||
|
|
@ -178,7 +184,7 @@ CONTAINS
|
|||
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, &
|
||||
|
|
@ -478,23 +484,24 @@ CONTAINS
|
|||
! We now calculate the Hamiltonian matrix elements
|
||||
|
||||
force_thread = 0.0_dp
|
||||
at_thread = 0.0_dp
|
||||
pv_thread = 0.0_dp
|
||||
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP SHARED (do_kp, basis_set, matrix_h, matrix_l, cell_to_index,&
|
||||
!$OMP sab_orb, matrix_p, sap_int, nkind, eps_ppnl, force, &
|
||||
!$OMP do_dR, deltaR, maxder, nder, &
|
||||
!$OMP doat, do_dR, deltaR, maxder, nder, &
|
||||
!$OMP locks, virial, use_virial, calculate_forces, do_soc, natom) &
|
||||
!$OMP PRIVATE (ikind, jkind, iatom, jatom, cell_b, rab, &
|
||||
!$OMP slot, iab, atom_a, f0, irow, icol, h_block, &
|
||||
!$OMP l_block_x, l_block_y, l_block_z, &
|
||||
!$OMP r_2block, r_3block, &
|
||||
!$OMP r_2block, r_3block, atk, &
|
||||
!$OMP found,p_block, iac, ibc, alist_ac, alist_bc, acint, bcint, &
|
||||
!$OMP achint, bchint, alkint, blkint, &
|
||||
!$OMP na, np, nb, katom, j, fa, fb, rbc, rac, &
|
||||
!$OMP kkind, kac, kbc, i, img, hash, iatom8) &
|
||||
!$OMP REDUCTION (+ : pv_thread, force_thread )
|
||||
!$OMP REDUCTION (+ : at_thread, pv_thread, force_thread )
|
||||
|
||||
!$OMP SINGLE
|
||||
!$ ALLOCATE (locks(nlock))
|
||||
|
|
@ -557,7 +564,7 @@ CONTAINS
|
|||
CALL dbcsr_get_block_p(matrix_h(3, img)%matrix, irow, icol, r_3block, found)
|
||||
END IF
|
||||
|
||||
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)
|
||||
END IF
|
||||
|
|
@ -711,6 +718,19 @@ CONTAINS
|
|||
END IF
|
||||
|
||||
END IF
|
||||
IF (doat) THEN
|
||||
IF (ASSOCIATED(p_block)) THEN
|
||||
katom = alist_ac%clist(kac)%catom
|
||||
IF (iatom <= jatom) THEN
|
||||
atk = SUM(p_block(1:na, 1:nb)* &
|
||||
MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 1))))
|
||||
ELSE
|
||||
atk = SUM(p_block(1:nb, 1:na)* &
|
||||
MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 1))))
|
||||
END IF
|
||||
at_thread(katom) = at_thread(katom) + f0*atk
|
||||
END IF
|
||||
END IF
|
||||
EXIT ! We have found a match and there can be only one single match
|
||||
END IF
|
||||
END DO
|
||||
|
|
@ -751,7 +771,11 @@ CONTAINS
|
|||
virial%pv_virial = virial%pv_virial + pv_thread
|
||||
END IF
|
||||
|
||||
IF (calculate_forces) THEN
|
||||
IF (doat) THEN
|
||||
atcore(1:natom) = atcore(1:natom) + at_thread
|
||||
END IF
|
||||
|
||||
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
|
||||
|
|
|
|||
|
|
@ -71,6 +71,7 @@ MODULE ed_analysis
|
|||
USE qs_core_energies, ONLY: calculate_ecore_alpha,&
|
||||
calculate_ecore_overlap,&
|
||||
calculate_ecore_self
|
||||
USE qs_core_hamiltonian, ONLY: core_matrices
|
||||
USE qs_dispersion_pairpot, ONLY: calculate_dispersion_pairpot
|
||||
USE qs_dispersion_types, ONLY: qs_dispersion_type
|
||||
USE qs_energy_types, ONLY: qs_energy_type
|
||||
|
|
@ -137,8 +138,8 @@ CONTAINS
|
|||
REAL(KIND=dp) :: ateps, checksum, e1, e2, e_pot, ealpha, &
|
||||
ecc, egcp, ehfx, ekts, evdw, focc, &
|
||||
sum_energy
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: amval, ate1h, ate1xc, atecc, ateks, &
|
||||
atener, atewald, odiag
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: amval, atcore, 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
|
||||
|
|
@ -151,8 +152,9 @@ CONTAINS
|
|||
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_p, matrix_s, matrix_t
|
||||
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: core_mat, matrix_h, matrix_hfx, &
|
||||
matrix_ks, matrix_p, matrix_s, matrix_t
|
||||
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: math, matp
|
||||
TYPE(dbcsr_type) :: dkmat, dmat
|
||||
TYPE(dbcsr_type), POINTER :: smat
|
||||
TYPE(dft_control_type), POINTER :: dft_control
|
||||
|
|
@ -436,11 +438,12 @@ CONTAINS
|
|||
CALL cp_fm_release(cvec2)
|
||||
|
||||
! energy arrays
|
||||
ALLOCATE (atener(natom), ateks(natom), atecc(natom))
|
||||
ALLOCATE (atener(natom), ateks(natom), atecc(natom), atcore(natom))
|
||||
ALLOCATE (ate1xc(natom), ate1h(natom), atewald(natom))
|
||||
atener = 0.0_dp
|
||||
ateks = 0.0_dp
|
||||
atecc = 0.0_dp
|
||||
atcore = 0.0_dp
|
||||
ate1xc = 0.0_dp
|
||||
ate1h = 0.0_dp
|
||||
atewald = 0.0_dp
|
||||
|
|
@ -455,7 +458,7 @@ CONTAINS
|
|||
CALL dbcsr_set(dkmat, 0.0_dp)
|
||||
! KS matrix + correction
|
||||
CALL get_qs_env(qs_env, matrix_h=matrix_h, matrix_ks=matrix_ks, kinetic=matrix_t)
|
||||
ALLOCATE (ks_mat(nspin), vhxc_mat(nspin), exc_mat(1))
|
||||
ALLOCATE (ks_mat(nspin), core_mat(1), vhxc_mat(nspin), exc_mat(1))
|
||||
DO ispin = 1, nspin
|
||||
ALLOCATE (ks_mat(ispin)%matrix)
|
||||
CALL dbcsr_create(ks_mat(ispin)%matrix, template=matrix_h(1)%matrix)
|
||||
|
|
@ -468,6 +471,12 @@ CONTAINS
|
|||
CALL dbcsr_copy(vhxc_mat(ispin)%matrix, smat)
|
||||
CALL dbcsr_set(vhxc_mat(ispin)%matrix, 0.0_dp)
|
||||
END DO
|
||||
!
|
||||
ALLOCATE (core_mat(1)%matrix)
|
||||
CALL dbcsr_create(core_mat(1)%matrix, template=matrix_h(1)%matrix)
|
||||
CALL dbcsr_copy(core_mat(1)%matrix, matrix_h(1)%matrix)
|
||||
CALL dbcsr_set(core_mat(1)%matrix, 0.0_dp)
|
||||
!
|
||||
ALLOCATE (exc_mat(1)%matrix)
|
||||
CALL dbcsr_create(exc_mat(1)%matrix, template=smat)
|
||||
CALL dbcsr_copy(exc_mat(1)%matrix, smat)
|
||||
|
|
@ -475,6 +484,12 @@ CONTAINS
|
|||
!
|
||||
CALL vhxc_correction(qs_env, vhxc_mat, exc_mat, atecc, ate1xc, ate1h)
|
||||
!
|
||||
CALL get_qs_env(qs_env, rho=rho)
|
||||
CALL qs_rho_get(rho, rho_ao=matrix_p)
|
||||
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)
|
||||
!
|
||||
IF (ewald_correction) THEN
|
||||
ALLOCATE (dve_mat%matrix)
|
||||
CALL dbcsr_create(dve_mat%matrix, template=matrix_h(1)%matrix)
|
||||
|
|
@ -485,6 +500,7 @@ CONTAINS
|
|||
!
|
||||
DO ispin = 1, nspin
|
||||
CALL dbcsr_add(ks_mat(ispin)%matrix, vhxc_mat(ispin)%matrix, 1.0_dp, 1.0_dp)
|
||||
CALL dbcsr_add(ks_mat(ispin)%matrix, core_mat(1)%matrix, 1.0_dp, -0.5_dp)
|
||||
END DO
|
||||
!
|
||||
IF (detailed_ener .AND. do_hfx) THEN
|
||||
|
|
@ -575,9 +591,11 @@ CONTAINS
|
|||
DEALLOCATE (ks_mat(ispin)%matrix, vhxc_mat(ispin)%matrix)
|
||||
CALL deallocate_mo_set(mos_loc(ispin))
|
||||
END DO
|
||||
CALL dbcsr_release(core_mat(1)%matrix)
|
||||
DEALLOCATE (core_mat(1)%matrix)
|
||||
CALL dbcsr_release(exc_mat(1)%matrix)
|
||||
DEALLOCATE (exc_mat(1)%matrix)
|
||||
DEALLOCATE (ks_mat, vhxc_mat, exc_mat)
|
||||
DEALLOCATE (ks_mat, core_mat, vhxc_mat, exc_mat)
|
||||
DEALLOCATE (mos_loc)
|
||||
DEALLOCATE (refbas_blk_sizes)
|
||||
DEALLOCATE (ref_basis_set_list)
|
||||
|
|
@ -604,6 +622,9 @@ CONTAINS
|
|||
group = para_env
|
||||
! KS energy
|
||||
atener(1:natom) = ateks(1:natom)
|
||||
! 1/2 of VPP contribution Tr[VPP(K)*P]
|
||||
CALL group%sum(atcore) ! maybe not needed
|
||||
atener(1:natom) = atener(1:natom) + 0.5_dp*atcore(1:natom)
|
||||
! core energy corrections
|
||||
CALL group%sum(atecc)
|
||||
atener(1:natom) = atener(1:natom) + atecc(1:natom)
|
||||
|
|
@ -651,7 +672,9 @@ CONTAINS
|
|||
CALL get_qs_env(qs_env, energy=energy)
|
||||
ekts = energy%kts/REAL(natom, KIND=dp)
|
||||
atener(1:natom) = atener(1:natom) + ekts
|
||||
|
||||
! 0.5 Vpp(at)*D + 0.5 * Vpp*D(at)
|
||||
IF (detailed_ener) atdet(1:natom, 2) = 0.5_dp*atdet(1:natom, 2) + 0.5_dp*atcore(1:natom)
|
||||
!
|
||||
IF (detailed_ener) THEN
|
||||
IF (unit_nr > 0) THEN
|
||||
WRITE (unit_nr, FMT="(/,T2,A)") "Detailed IAO Atomic Energy Components "
|
||||
|
|
@ -726,7 +749,7 @@ CONTAINS
|
|||
"!--------------------------- END OF ED ANALYSIS ------------------------------!"
|
||||
END IF
|
||||
DEALLOCATE (bcenter)
|
||||
DEALLOCATE (atener, ateks, atecc, ate1xc, ate1h, atewald)
|
||||
DEALLOCATE (atener, ateks, atecc, atcore, ate1xc, ate1h, atewald)
|
||||
|
||||
CALL timestop(handle)
|
||||
|
||||
|
|
|
|||
|
|
@ -140,7 +140,7 @@ MODULE qs_core_hamiltonian
|
|||
|
||||
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_core_hamiltonian'
|
||||
|
||||
PUBLIC :: build_core_hamiltonian_matrix
|
||||
PUBLIC :: build_core_hamiltonian_matrix, core_matrices
|
||||
PUBLIC :: dump_info_core_hamiltonian, qs_matrix_h_allocate_imag_from_real
|
||||
|
||||
CONTAINS
|
||||
|
|
@ -412,15 +412,17 @@ CONTAINS
|
|||
!> \param matrix_p ...
|
||||
!> \param calculate_forces ...
|
||||
!> \param nder ...
|
||||
!> \param atcore ...
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE core_matrices(qs_env, matrix_h, matrix_p, calculate_forces, nder)
|
||||
SUBROUTINE core_matrices(qs_env, matrix_h, matrix_p, calculate_forces, nder, atcore)
|
||||
|
||||
TYPE(qs_environment_type), POINTER :: qs_env
|
||||
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p
|
||||
LOGICAL, INTENT(IN) :: calculate_forces
|
||||
INTEGER, INTENT(IN) :: nder
|
||||
REAL(KIND=dp), DIMENSION(:), OPTIONAL :: atcore
|
||||
|
||||
INTEGER :: nimages
|
||||
INTEGER :: natom, nimages
|
||||
INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
|
||||
LOGICAL :: all_present, my_gt_nl, ppl_present, &
|
||||
ppnl_present, use_virial
|
||||
|
|
@ -438,8 +440,11 @@ CONTAINS
|
|||
TYPE(virial_type), POINTER :: virial
|
||||
|
||||
NULLIFY (dft_control)
|
||||
CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, dft_control=dft_control)
|
||||
CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, dft_control=dft_control, natom=natom)
|
||||
nimages = dft_control%nimages
|
||||
IF (PRESENT(atcore)) THEN
|
||||
CPASSERT(SIZE(atcore) >= natom)
|
||||
END IF
|
||||
|
||||
! check whether a gauge transformed version of the non-local potential part has to be used
|
||||
my_gt_nl = .FALSE.
|
||||
|
|
@ -484,7 +489,8 @@ CONTAINS
|
|||
all_present = ASSOCIATED(sac_ae)
|
||||
IF (all_present) THEN
|
||||
CALL 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=atcore)
|
||||
END IF
|
||||
|
||||
! *** compute the ppl contribution to the core hamiltonian ***
|
||||
|
|
@ -500,12 +506,12 @@ CONTAINS
|
|||
ELSE
|
||||
CALL 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, "ORB")
|
||||
nimages, cell_to_index, "ORB", atcore=atcore)
|
||||
END IF
|
||||
ELSE
|
||||
CALL 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, "ORB")
|
||||
nimages, cell_to_index, "ORB", atcore=atcore)
|
||||
END IF
|
||||
END IF
|
||||
END IF
|
||||
|
|
@ -517,7 +523,7 @@ CONTAINS
|
|||
IF (.NOT. my_gt_nl) THEN
|
||||
CALL build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
|
||||
qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, eps_ppnl, &
|
||||
nimages, cell_to_index, "ORB")
|
||||
nimages, cell_to_index, "ORB", atcore=atcore)
|
||||
END IF
|
||||
END IF
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue