mirror of
https://github.com/cp2k/cp2k.git
synced 2026-07-28 22:25:32 -04:00
Current (#335)
* UPF pseudo HGH type energy correct, forces stress missing * Forces and Stress for UPF HGH potentials * Fit local PP directly * Add UPF regtests * Regtests are not very stable (compiler) * Adjust accuracy
This commit is contained in:
Juerg Hutter
GitHub
parent
0c1cc072f3
commit
efb4dd7a66
13 changed files with 379 additions and 101 deletions
|
|
@ -1308,6 +1308,7 @@ CONTAINS
|
|||
is_nonlocal=sgp_atompot%is_nonlocal, &
|
||||
n_nonlocal=n, a_nonlocal=ap, h_nonlocal=hhp, c_nonlocal=ccp)
|
||||
! nonlocal
|
||||
sgp_atompot%has_nonlocal = ANY(sgp_atompot%is_nonlocal)
|
||||
sgp_atompot%lmax = lm
|
||||
CPASSERT(n <= SIZE(sgp_atompot%a_nonlocal))
|
||||
sgp_atompot%n_nonlocal = n
|
||||
|
|
|
|||
205
src/atom_sgp.F
205
src/atom_sgp.F
|
|
@ -22,7 +22,9 @@ MODULE atom_sgp
|
|||
USE kahan_sum, ONLY: accurate_dot_product
|
||||
USE kinds, ONLY: dp
|
||||
USE mathconstants, ONLY: dfac,&
|
||||
rootpi
|
||||
fourpi,&
|
||||
rootpi,&
|
||||
sqrt2
|
||||
USE mathlib, ONLY: diamat_all,&
|
||||
get_pseudo_inverse_diag
|
||||
USE powell, ONLY: opt_state_type,&
|
||||
|
|
@ -80,8 +82,8 @@ CONTAINS
|
|||
|
||||
INTEGER :: i, n
|
||||
LOGICAL :: is_ecp, is_upf
|
||||
REAL(KIND=dp) :: errcc
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cgauss
|
||||
REAL(KIND=dp) :: errcc, rcut
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cgauss, cutpots, cutpotu
|
||||
TYPE(atom_basis_type) :: basis
|
||||
TYPE(opmat_type), POINTER :: core, hnl, score, shnl
|
||||
|
||||
|
|
@ -94,6 +96,21 @@ CONTAINS
|
|||
IF (PRESENT(upf_pot)) is_upf = .TRUE.
|
||||
CPASSERT(.NOT. (is_ecp .AND. is_upf))
|
||||
|
||||
! upf has often very small grids, use a smooth cutoff function
|
||||
IF (is_upf) THEN
|
||||
n = SIZE(upf_pot%r)
|
||||
ALLOCATE (cutpotu(n))
|
||||
rcut = MAXVAL(upf_pot%r)
|
||||
CALL cutpot(cutpotu, upf_pot%r, rcut, 2.5_dp)
|
||||
n = basis%grid%nr
|
||||
ALLOCATE (cutpots(n))
|
||||
CALL cutpot(cutpots, basis%grid%rad, rcut, 2.5_dp)
|
||||
ELSE
|
||||
n = basis%grid%nr
|
||||
ALLOCATE (cutpots(n))
|
||||
cutpots = 1.0_dp
|
||||
END IF
|
||||
|
||||
! generate the transformed potentials
|
||||
IF (is_ecp) THEN
|
||||
CALL ecp_sgp_constr(ecp_pot, sgp_pot, basis)
|
||||
|
|
@ -113,19 +130,21 @@ CONTAINS
|
|||
IF (is_ecp) THEN
|
||||
CALL ecpints(hnl%op, basis, ecp_pot)
|
||||
ELSEIF (is_upf) THEN
|
||||
CALL upfints(core%op, hnl%op, basis, upf_pot, sgp_pot%ac_local)
|
||||
CALL upfints(core%op, hnl%op, basis, upf_pot, cutpotu, sgp_pot%ac_local)
|
||||
ELSE
|
||||
CPABORT("")
|
||||
END IF
|
||||
!
|
||||
CALL sgpints(score%op, shnl%op, basis, sgp_pot)
|
||||
CALL sgpints(score%op, shnl%op, basis, sgp_pot, cutpots)
|
||||
!
|
||||
error = 0.0_dp
|
||||
IF (sgp_pot%has_local) THEN
|
||||
error(1) = MAXVAL(ABS(core%op-score%op))
|
||||
n = MIN(2, UBOUND(core%op, 3))
|
||||
error(1) = MAXVAL(ABS(core%op(:, :, 0:n)-score%op(:, :, 0:n)))
|
||||
END IF
|
||||
IF (sgp_pot%has_nonlocal) THEN
|
||||
error(2) = MAXVAL(ABS(hnl%op-shnl%op))
|
||||
n = MIN(3, UBOUND(hnl%op, 3))
|
||||
error(2) = MAXVAL(ABS(hnl%op(:, :, 0:n)-shnl%op(:, :, 0:n)))
|
||||
END IF
|
||||
IF (sgp_pot%has_nlcc) THEN
|
||||
IF (is_upf) THEN
|
||||
|
|
@ -144,6 +163,13 @@ CONTAINS
|
|||
error(3) = errcc
|
||||
END IF
|
||||
!
|
||||
IF (is_upf) THEN
|
||||
DEALLOCATE (cutpotu)
|
||||
DEALLOCATE (cutpots)
|
||||
ELSE
|
||||
DEALLOCATE (cutpots)
|
||||
END IF
|
||||
!
|
||||
CALL release_opmat(score)
|
||||
CALL release_opmat(shnl)
|
||||
CALL release_opmat(core)
|
||||
|
|
@ -170,8 +196,8 @@ CONTAINS
|
|||
|
||||
INTEGER :: i, n, ppot_type
|
||||
LOGICAL :: do_transform, explicit, is_ecp, is_upf
|
||||
REAL(KIND=dp) :: errcc
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cgauss
|
||||
REAL(KIND=dp) :: errcc, rcut
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cgauss, cutpots, cutpotu
|
||||
TYPE(atom_ecppot_type), POINTER :: ecp_pot
|
||||
TYPE(atom_sgp_potential_type) :: sgp_pot
|
||||
TYPE(atom_type), POINTER :: atom_ref
|
||||
|
|
@ -227,15 +253,30 @@ CONTAINS
|
|||
CALL create_opmat(score, atom_ref%basis%nbas)
|
||||
CALL create_opmat(shnl, atom_ref%basis%nbas, 5)
|
||||
!
|
||||
! upf has often very small grids, use a smooth cutoff function
|
||||
IF (is_upf) THEN
|
||||
n = SIZE(upf_pot%r)
|
||||
ALLOCATE (cutpotu(n))
|
||||
rcut = MAXVAL(upf_pot%r)
|
||||
CALL cutpot(cutpotu, upf_pot%r, rcut, 2.5_dp)
|
||||
n = atom_ref%basis%grid%nr
|
||||
ALLOCATE (cutpots(n))
|
||||
CALL cutpot(cutpots, atom_ref%basis%grid%rad, rcut, 2.5_dp)
|
||||
ELSE
|
||||
n = atom_ref%basis%grid%nr
|
||||
ALLOCATE (cutpots(n))
|
||||
cutpots = 1.0_dp
|
||||
END IF
|
||||
!
|
||||
IF (is_ecp) THEN
|
||||
CALL ecpints(hnl%op, atom_ref%basis, ecp_pot)
|
||||
ELSEIF (is_upf) THEN
|
||||
CALL upfints(core%op, hnl%op, atom_ref%basis, upf_pot, sgp_pot%ac_local)
|
||||
CALL upfints(core%op, hnl%op, atom_ref%basis, upf_pot, cutpotu, sgp_pot%ac_local)
|
||||
ELSE
|
||||
CPABORT("")
|
||||
END IF
|
||||
!
|
||||
CALL sgpints(score%op, shnl%op, atom_ref%basis, sgp_pot)
|
||||
CALL sgpints(score%op, shnl%op, atom_ref%basis, sgp_pot, cutpots)
|
||||
!
|
||||
IF (sgp_pot%has_local) THEN
|
||||
n = MIN(3, UBOUND(core%op, 3))
|
||||
|
|
@ -276,6 +317,13 @@ CONTAINS
|
|||
END IF
|
||||
END IF
|
||||
!
|
||||
IF (is_upf) THEN
|
||||
DEALLOCATE (cutpotu)
|
||||
DEALLOCATE (cutpots)
|
||||
ELSE
|
||||
DEALLOCATE (cutpots)
|
||||
END IF
|
||||
!
|
||||
CALL release_opmat(score)
|
||||
CALL release_opmat(shnl)
|
||||
CALL release_opmat(core)
|
||||
|
|
@ -454,7 +502,7 @@ CONTAINS
|
|||
ostate%nvar = 2
|
||||
x(1) = 1.00_dp !starting point of geometric series
|
||||
x(2) = 1.20_dp !factor of series
|
||||
ostate%rhoend = 1.e-8_dp
|
||||
ostate%rhoend = 1.e-14_dp
|
||||
ostate%rhobeg = 5.e-2_dp
|
||||
ostate%maxfun = 1000
|
||||
ostate%iprint = 1
|
||||
|
|
@ -465,7 +513,7 @@ CONTAINS
|
|||
DO ir = 1, nl
|
||||
al(ir) = x(1)*x(2)**(ir-1)
|
||||
END DO
|
||||
CALL pplocal_error(nl, al, cl, vloc, vgauss, gbasis, upf_pot%r, ww, ostate%f)
|
||||
CALL pplocal_error(nl, al, cl, vloc, vgauss, gbasis, upf_pot%r, ww, 1, ostate%f)
|
||||
END IF
|
||||
IF (ostate%state == -1) EXIT
|
||||
CALL powell_optimize(ostate%nvar, x, ostate)
|
||||
|
|
@ -475,7 +523,7 @@ CONTAINS
|
|||
DO ir = 1, nl
|
||||
al(ir) = x(1)*x(2)**(ir-1)
|
||||
END DO
|
||||
CALL pplocal_error(nl, al, cl, vloc, vgauss, gbasis, upf_pot%r, ww, errloc)
|
||||
CALL pplocal_error(nl, al, cl, vloc, vgauss, gbasis, upf_pot%r, ww, 1, errloc)
|
||||
!
|
||||
ALLOCATE (sgp_pot%a_local(nl), sgp_pot%c_local(nl))
|
||||
sgp_pot%n_local = nl
|
||||
|
|
@ -566,6 +614,8 @@ CONTAINS
|
|||
sgp_pot%is_nonlocal(la) = .TRUE.
|
||||
DEALLOCATE (score, qmat)
|
||||
END DO
|
||||
! SQRT(4PI)
|
||||
sgp_pot%c_nonlocal = sgp_pot%c_nonlocal/SQRT(fourpi)
|
||||
CALL release_atom_basis(gbasis)
|
||||
DEALLOCATE (pgauss, vloc)
|
||||
DEALLOCATE (al, cl, smat, sinv, tmat, cmat)
|
||||
|
|
@ -645,6 +695,8 @@ CONTAINS
|
|||
sgp_pot%is_nonlocal(la) = .TRUE.
|
||||
DEALLOCATE (score, qmat)
|
||||
END DO
|
||||
! SQRT(4PI)
|
||||
sgp_pot%c_nonlocal = sgp_pot%c_nonlocal/SQRT(fourpi)
|
||||
CALL release_atom_basis(gbasis)
|
||||
DEALLOCATE (pgauss, pproa, pprob)
|
||||
DEALLOCATE (al, cl, smat, sinv, tmat, cmat)
|
||||
|
|
@ -727,12 +779,14 @@ CONTAINS
|
|||
!> \param hnl ...
|
||||
!> \param basis ...
|
||||
!> \param upf_pot ...
|
||||
!> \param cutpotu ...
|
||||
!> \param ac_local ...
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE upfints(core, hnl, basis, upf_pot, ac_local)
|
||||
SUBROUTINE upfints(core, hnl, basis, upf_pot, cutpotu, ac_local)
|
||||
REAL(KIND=dp), DIMENSION(:, :, 0:) :: core, hnl
|
||||
TYPE(atom_basis_type), INTENT(INOUT) :: basis
|
||||
TYPE(atom_upfpot_type) :: upf_pot
|
||||
REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: cutpotu
|
||||
REAL(KIND=dp), INTENT(IN) :: ac_local
|
||||
|
||||
CHARACTER(len=*), PARAMETER :: routineN = 'upfints', routineP = moduleN//':'//routineN
|
||||
|
|
@ -755,12 +809,13 @@ CONTAINS
|
|||
zval = upf_pot%zion
|
||||
DO i = 1, n
|
||||
IF (upf_pot%r(i) < 1.e-12_dp) THEN
|
||||
spot(i) = spot(i)+2.0_dp*zval/rootpi/ac_local
|
||||
spot(i) = spot(i)+sqrt2*zval/rootpi/ac_local
|
||||
ELSE
|
||||
rc = upf_pot%r(i)/ac_local
|
||||
rc = upf_pot%r(i)/ac_local/sqrt2
|
||||
spot(i) = spot(i)+zval/upf_pot%r(i)*erf(rc)
|
||||
END IF
|
||||
END DO
|
||||
spot(:) = spot(:)*cutpotu(:)
|
||||
|
||||
CALL numpot_matrix(core, spot, gbasis, 0)
|
||||
DEALLOCATE (spot)
|
||||
|
|
@ -800,6 +855,7 @@ CONTAINS
|
|||
m = SIZE(upf_pot%vsemi(:, la+1))
|
||||
ALLOCATE (spot(m))
|
||||
spot(:) = upf_pot%vsemi(:, la+1)-upf_pot%vlocal(:)
|
||||
spot(:) = spot(:)*cutpotu(:)
|
||||
n = basis%nbas(la)
|
||||
DO i = 1, n
|
||||
DO j = i, n
|
||||
|
|
@ -868,11 +924,13 @@ CONTAINS
|
|||
!> \param hnl ...
|
||||
!> \param basis ...
|
||||
!> \param sgp_pot ...
|
||||
!> \param cutpots ...
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE sgpints(core, hnl, basis, sgp_pot)
|
||||
SUBROUTINE sgpints(core, hnl, basis, sgp_pot, cutpots)
|
||||
REAL(KIND=dp), DIMENSION(:, :, 0:) :: core, hnl
|
||||
TYPE(atom_basis_type), INTENT(INOUT) :: basis
|
||||
TYPE(atom_sgp_potential_type) :: sgp_pot
|
||||
REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: cutpots
|
||||
|
||||
CHARACTER(len=*), PARAMETER :: routineN = 'sgpints', routineP = moduleN//':'//routineN
|
||||
|
||||
|
|
@ -894,6 +952,7 @@ CONTAINS
|
|||
DO i = 1, sgp_pot%n_local
|
||||
cpot(:) = cpot(:)+sgp_pot%c_local(i)*EXP(-sgp_pot%a_local(i)*rad(:)**2)
|
||||
END DO
|
||||
cpot(:) = cpot(:)*cutpots(:)
|
||||
CALL numpot_matrix(core, cpot, basis, 0)
|
||||
END IF
|
||||
DEALLOCATE (cpot)
|
||||
|
|
@ -917,10 +976,12 @@ CONTAINS
|
|||
c = sgp_pot%c_nonlocal(j, i, l)
|
||||
pgauss(:) = pgauss(:)+c*EXP(-a*rad(:)**2)*rad(:)**l
|
||||
END DO
|
||||
pgauss(:) = pgauss(:)*cutpots(:)
|
||||
DO ia = 1, na
|
||||
qmat(ia, i) = SUM(basis%bf(:, ia, l)*pgauss(:)*basis%grid%wr(:))
|
||||
END DO
|
||||
END DO
|
||||
qmat = SQRT(fourpi)*qmat
|
||||
DO i = 1, na
|
||||
DO j = i, na
|
||||
DO ia = 1, n
|
||||
|
|
@ -1007,13 +1068,15 @@ CONTAINS
|
|||
!> \param gbasis ...
|
||||
!> \param rad ...
|
||||
!> \param ww ...
|
||||
!> \param method ...
|
||||
!> \param errloc ...
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE pplocal_error(nl, al, cl, vloc, vgauss, gbasis, rad, ww, errloc)
|
||||
SUBROUTINE pplocal_error(nl, al, cl, vloc, vgauss, gbasis, rad, ww, method, errloc)
|
||||
INTEGER :: nl
|
||||
REAL(KIND=dp), DIMENSION(:) :: al, cl, vloc, vgauss
|
||||
TYPE(atom_basis_type) :: gbasis
|
||||
REAL(KIND=dp), DIMENSION(:) :: rad, ww
|
||||
INTEGER, INTENT(IN) :: method
|
||||
REAL(KIND=dp) :: errloc
|
||||
|
||||
CHARACTER(len=*), PARAMETER :: routineN = 'pplocal_error', routineP = moduleN//':'//routineN
|
||||
|
|
@ -1024,43 +1087,57 @@ CONTAINS
|
|||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: gmat
|
||||
|
||||
cl = 0.0_dp
|
||||
!
|
||||
ALLOCATE (tv(nl), smat(nl, nl), sinv(nl, nl))
|
||||
!
|
||||
smat = 0.0_dp
|
||||
tv = 0.0_dp
|
||||
DO la = 0, MIN(UBOUND(gbasis%nbas, 1), 3)
|
||||
na = gbasis%nbas(la)
|
||||
ALLOCATE (rmat(na, na), gmat(na, na, nl))
|
||||
rmat = 0.0_dp
|
||||
gmat = 0.0_dp
|
||||
DO ia = 1, na
|
||||
DO ib = ia, na
|
||||
rmat(ia, ib) = SUM(gbasis%bf(:, ia, la)*gbasis%bf(:, ib, la)*vloc(:)*ww(:))
|
||||
rmat(ib, ia) = rmat(ia, ib)
|
||||
END DO
|
||||
END DO
|
||||
IF (method == 1) THEN
|
||||
ALLOCATE (tv(nl), smat(nl, nl), sinv(nl, nl))
|
||||
DO ir = 1, nl
|
||||
vgauss(:) = EXP(-al(ir)*rad(:)**2)
|
||||
DO ix = 1, nl
|
||||
smat(ir, ix) = SUM(vgauss(:)*EXP(-al(ix)*rad(:)**2)*ww(:))
|
||||
END DO
|
||||
tv(ir) = SUM(vloc(:)*vgauss(:)*ww(:))
|
||||
END DO
|
||||
CALL get_pseudo_inverse_diag(smat(1:nl, 1:nl), sinv(1:nl, 1:nl), 1.e-12_dp)
|
||||
cl(1:nl) = MATMUL(sinv(1:nl, 1:nl), tv(1:nl))
|
||||
ELSE
|
||||
!
|
||||
ALLOCATE (tv(nl), smat(nl, nl), sinv(nl, nl))
|
||||
!
|
||||
smat = 0.0_dp
|
||||
tv = 0.0_dp
|
||||
DO la = 0, MIN(UBOUND(gbasis%nbas, 1), 3)
|
||||
na = gbasis%nbas(la)
|
||||
ALLOCATE (rmat(na, na), gmat(na, na, nl))
|
||||
rmat = 0.0_dp
|
||||
gmat = 0.0_dp
|
||||
DO ia = 1, na
|
||||
DO ib = ia, na
|
||||
gmat(ia, ib, ir) = SUM(gbasis%bf(:, ia, la)*gbasis%bf(:, ib, la)*vgauss(:)*ww(:))
|
||||
gmat(ib, ia, ir) = gmat(ia, ib, ir)
|
||||
rmat(ia, ib) = SUM(gbasis%bf(:, ia, la)*gbasis%bf(:, ib, la)*vloc(:)*ww(:))
|
||||
rmat(ib, ia) = rmat(ia, ib)
|
||||
END DO
|
||||
END DO
|
||||
END DO
|
||||
DO ir = 1, nl
|
||||
tv(ir) = tv(ir)+accurate_dot_product(rmat, gmat(:, :, ir))
|
||||
DO ix = ir, nl
|
||||
smat(ir, ix) = smat(ir, ix)+accurate_dot_product(gmat(:, :, ix), gmat(:, :, ir))
|
||||
smat(ix, ir) = smat(ir, ix)
|
||||
DO ir = 1, nl
|
||||
vgauss(:) = EXP(-al(ir)*rad(:)**2)
|
||||
DO ia = 1, na
|
||||
DO ib = ia, na
|
||||
gmat(ia, ib, ir) = SUM(gbasis%bf(:, ia, la)*gbasis%bf(:, ib, la)*vgauss(:)*ww(:))
|
||||
gmat(ib, ia, ir) = gmat(ia, ib, ir)
|
||||
END DO
|
||||
END DO
|
||||
END DO
|
||||
DO ir = 1, nl
|
||||
tv(ir) = tv(ir)+accurate_dot_product(rmat, gmat(:, :, ir))
|
||||
DO ix = ir, nl
|
||||
smat(ir, ix) = smat(ir, ix)+accurate_dot_product(gmat(:, :, ix), gmat(:, :, ir))
|
||||
smat(ix, ir) = smat(ir, ix)
|
||||
END DO
|
||||
END DO
|
||||
DEALLOCATE (rmat, gmat)
|
||||
END DO
|
||||
DEALLOCATE (rmat, gmat)
|
||||
END DO
|
||||
sinv = 0.0_dp
|
||||
CALL get_pseudo_inverse_diag(smat(1:nl, 1:nl), sinv(1:nl, 1:nl), 1.e-12_dp)
|
||||
cl(1:nl) = MATMUL(sinv(1:nl, 1:nl), tv(1:nl))
|
||||
sinv = 0.0_dp
|
||||
CALL get_pseudo_inverse_diag(smat(1:nl, 1:nl), sinv(1:nl, 1:nl), 1.e-12_dp)
|
||||
cl(1:nl) = MATMUL(sinv(1:nl, 1:nl), tv(1:nl))
|
||||
ENDIF
|
||||
!
|
||||
vgauss = 0.0_dp
|
||||
DO ir = 1, nl
|
||||
vgauss(:) = vgauss(:)+cl(ir)*EXP(-al(ir)*rad(:)**2)
|
||||
|
|
@ -1071,4 +1148,36 @@ CONTAINS
|
|||
!
|
||||
END SUBROUTINE pplocal_error
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief ...
|
||||
!> \param pot ...
|
||||
!> \param r ...
|
||||
!> \param rcut ...
|
||||
!> \param rsmooth ...
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE cutpot(pot, r, rcut, rsmooth)
|
||||
REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: pot
|
||||
REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: r
|
||||
REAL(KIND=dp), INTENT(IN) :: rcut, rsmooth
|
||||
|
||||
INTEGER :: i, n
|
||||
REAL(KIND=dp) :: rab, rx, x
|
||||
|
||||
n = SIZE(pot)
|
||||
CPASSERT(n <= SIZE(r))
|
||||
|
||||
pot(:) = 1.0_dp
|
||||
DO i = 1, n
|
||||
rab = r(i)
|
||||
IF (rab > rcut) THEN
|
||||
pot(i) = 0.0_dp
|
||||
ELSE IF (rab > rcut-rsmooth) THEN
|
||||
rx = rab-(rcut-rsmooth)
|
||||
x = rx/rsmooth
|
||||
pot(i) = -6._dp*x**5+15._dp*x**4-10._dp*x**3+1._dp
|
||||
END IF
|
||||
END DO
|
||||
|
||||
END SUBROUTINE cutpot
|
||||
|
||||
END MODULE atom_sgp
|
||||
|
|
|
|||
|
|
@ -263,8 +263,8 @@ CONTAINS
|
|||
ELSE IF (ASSOCIATED(spotential(kkind)%sgp_potential)) THEN
|
||||
! SGP potential
|
||||
dogth = .FALSE.
|
||||
nppnl = spotential(kkind)%sgp_potential%nppnl
|
||||
IF (nppnl == 0) CYCLE
|
||||
nprjc = spotential(kkind)%sgp_potential%nppnl
|
||||
IF (nprjc == 0) CYCLE
|
||||
nnl = spotential(kkind)%sgp_potential%n_nonlocal
|
||||
lppnl = spotential(kkind)%sgp_potential%lmax
|
||||
is_nonlocal = .FALSE.
|
||||
|
|
@ -277,7 +277,7 @@ CONTAINS
|
|||
radp(:) = ppnl_radius
|
||||
cprj => spotential(kkind)%sgp_potential%cprj_ppnl
|
||||
hprj => spotential(kkind)%sgp_potential%vprj_ppnl
|
||||
nprjc = SIZE(cprj, 2)
|
||||
nppnl = SIZE(cprj, 2)
|
||||
ELSE
|
||||
CYCLE
|
||||
END IF
|
||||
|
|
@ -336,38 +336,37 @@ CONTAINS
|
|||
lc_max, lc_min, 1, rprjc, zetc, rac, dac, sab, nder, .TRUE., ai_work, ldai)
|
||||
! *** Transformation step projector functions (cartesian->spherical) ***
|
||||
DO i = 1, maxder
|
||||
first_col = (i-1)*SIZE(work, 1)
|
||||
first_col = (i-1)*ldsab
|
||||
CALL dgemm("N", "N", ncoa, nprjc, ncoc, 1.0_dp, sab(1, first_col+1), SIZE(sab, 1), &
|
||||
cprj(1, prjc), SIZE(cprj, 1), 0.0_dp, work(1, first_col+prjc), SIZE(work, 1))
|
||||
cprj(1, prjc), SIZE(cprj, 1), 0.0_dp, work(1, first_col+prjc), ldsab)
|
||||
END DO
|
||||
prjc = prjc+nprjc
|
||||
END DO
|
||||
DO i = 1, maxder
|
||||
first_col = (i-1)*SIZE(work, 1)+1
|
||||
first_col = (i-1)*ldsab+1
|
||||
! *** Contraction step (basis functions) ***
|
||||
CALL dgemm("T", "N", nsgf_seta(iset), nppnl, ncoa, 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
|
||||
work(1, first_col), SIZE(work, 1), 0.0_dp, clist%acint(sgfa, 1, i), nsgfa)
|
||||
work(1, first_col), ldsab, 0.0_dp, clist%acint(sgfa, 1, i), nsgfa)
|
||||
! *** Multiply with interaction matrix(h) ***
|
||||
CALL dgemm("N", "N", nsgf_seta(iset), nppnl, nppnl, 1.0_dp, clist%acint(sgfa, 1, i), nsgfa, &
|
||||
vprj_ppnl(1, 1), SIZE(vprj_ppnl, 1), 0.0_dp, clist%achint(sgfa, 1, i), nsgfa)
|
||||
END DO
|
||||
ELSE
|
||||
! SGP potential
|
||||
work = 0._dp
|
||||
! *** Calculate the primitive overlap integrals ***
|
||||
CALL overlap(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
|
||||
lppnl, 0, nnl, radp, a_nl, rac, dac, sab, nder, .TRUE., ai_work, ldai)
|
||||
DO i = 1, maxder
|
||||
first_col = (i-1)*nppnl
|
||||
first_col = (i-1)*ldsab+1
|
||||
! *** Transformation step projector functions (cartesian->spherical) ***
|
||||
CALL dgemm("N", "N", ncoa, nppnl, nprjc, 1.0_dp, sab(1, first_col+1), SIZE(sab, 1), &
|
||||
cprj(1, 1), SIZE(cprj, 1), 0.0_dp, work(1, first_col+1), SIZE(work, 1))
|
||||
CALL dgemm("N", "N", ncoa, nppnl, nprjc, 1.0_dp, sab(1, first_col), ldsab, &
|
||||
cprj(1, 1), SIZE(cprj, 1), 0.0_dp, work(1, 1), ldsab)
|
||||
! *** Contraction step (basis functions) ***
|
||||
CALL dgemm("T", "N", nsgf_seta(iset), nprjc, ncoa, 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
|
||||
work(1, first_col), SIZE(work, 1), 0.0_dp, clist%acint(sgfa, 1, i), nsgfa)
|
||||
CALL dgemm("T", "N", nsgf_seta(iset), nppnl, ncoa, 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
|
||||
work(1, 1), ldsab, 0.0_dp, clist%acint(sgfa, 1, i), nsgfa)
|
||||
! *** Multiply with interaction matrix(h) ***
|
||||
ncoc = sgfa+nsgf_seta(iset)-1
|
||||
DO j = 1, nprjc
|
||||
DO j = 1, nppnl
|
||||
clist%achint(sgfa:ncoc, j, i) = clist%acint(sgfa:ncoc, j, i)*hprj(j)
|
||||
END DO
|
||||
END DO
|
||||
|
|
@ -375,9 +374,7 @@ CONTAINS
|
|||
END DO
|
||||
clist%maxac = MAXVAL(ABS(clist%acint(:, :, 1)))
|
||||
clist%maxach = MAXVAL(ABS(clist%achint(:, :, 1)))
|
||||
IF (.NOT. dogth) THEN
|
||||
DEALLOCATE (radp)
|
||||
END IF
|
||||
IF (.NOT. dogth) DEALLOCATE (radp)
|
||||
END DO
|
||||
|
||||
DEALLOCATE (sab, ai_work, work)
|
||||
|
|
|
|||
|
|
@ -43,6 +43,7 @@ MODULE qs_kind_types
|
|||
USE cp_control_types, ONLY: dft_control_type,&
|
||||
qs_control_type
|
||||
USE cp_log_handling, ONLY: cp_get_default_logger,&
|
||||
cp_logger_get_default_io_unit,&
|
||||
cp_logger_type
|
||||
USE cp_output_handling, ONLY: cp_p_file,&
|
||||
cp_print_key_finished_output,&
|
||||
|
|
@ -1404,8 +1405,8 @@ CONTAINS
|
|||
CHARACTER(LEN=default_string_length), &
|
||||
DIMENSION(maxbas) :: basis_set_form, basis_set_name, &
|
||||
basis_set_type
|
||||
INTEGER :: handle, i, i_rep, ipaodesc, ipaopot, ipos, j, jj, k_rep, l, m, n_rep, nb_rep, &
|
||||
nexp, ngauss, nlcc, nloc, nnl, norbitals, npaodesc, npaopot, nppnl, nspin, z
|
||||
INTEGER :: handle, i, i_rep, iounit, ipaodesc, ipaopot, ipos, j, jj, k_rep, l, m, n_rep, &
|
||||
nb_rep, nexp, ngauss, nlcc, nloc, nnl, norbitals, npaodesc, npaopot, nppnl, nspin, z
|
||||
INTEGER, DIMENSION(:), POINTER :: add_el, elec_conf, orbitals
|
||||
LOGICAL :: check, explicit, explicit_basis, explicit_kgpot, explicit_potential, nobasis, &
|
||||
section_enabled, subsection_enabled, update_input
|
||||
|
|
@ -1417,6 +1418,7 @@ CONTAINS
|
|||
TYPE(atom_ecppot_type) :: ecppot
|
||||
TYPE(atom_sgp_potential_type) :: sgppot
|
||||
TYPE(atom_upfpot_type) :: upfpot
|
||||
TYPE(cp_logger_type), POINTER :: logger
|
||||
TYPE(gto_basis_set_type), POINTER :: sup_basis_set, tmp_basis_set
|
||||
TYPE(section_vals_type), POINTER :: basis_section, bs_section, dft_plus_u_section, &
|
||||
dft_section, enforce_occupation_section, kgpot_section, pao_desc_section, &
|
||||
|
|
@ -1425,6 +1427,10 @@ CONTAINS
|
|||
|
||||
CALL timeset(routineN, handle)
|
||||
|
||||
NULLIFY (logger)
|
||||
logger => cp_get_default_logger()
|
||||
iounit = cp_logger_get_default_io_unit(logger)
|
||||
|
||||
NULLIFY (elec_conf)
|
||||
|
||||
update_input = .TRUE.
|
||||
|
|
@ -2219,6 +2225,18 @@ CONTAINS
|
|||
has_nlcc=.FALSE.)
|
||||
! convert PP
|
||||
CALL sgp_construction(sgp_pot=sgppot, ecp_pot=ecppot, error=error)
|
||||
IF (iounit > 0) THEN
|
||||
WRITE (iounit, "(/,T2,'PP Transformation for ',A)") TRIM(ecppot%pname)
|
||||
IF (sgppot%has_local) THEN
|
||||
WRITE (iounit, "(T8,'Accuracy for local part:',T61,F20.12)") error(1)
|
||||
END IF
|
||||
IF (sgppot%has_nonlocal) THEN
|
||||
WRITE (iounit, "(T8,'Accuracy for nonlocal part:',T61,F20.12)") error(2)
|
||||
END IF
|
||||
IF (sgppot%has_nlcc) THEN
|
||||
WRITE (iounit, "(T8,'Accuracy for NLCC density:',T61,F20.12)") error(3)
|
||||
END IF
|
||||
END IF
|
||||
IF (sgppot%has_nonlocal) THEN
|
||||
CALL set_potential(qs_kind%sgp_potential, n_nonlocal=sgppot%n_nonlocal, lmax=sgppot%lmax, &
|
||||
is_nonlocal=sgppot%is_nonlocal)
|
||||
|
|
@ -2264,6 +2282,18 @@ CONTAINS
|
|||
zeff=upfpot%zion, z=z, has_nlcc=upfpot%core_correction)
|
||||
! convert pp
|
||||
CALL sgp_construction(sgp_pot=sgppot, upf_pot=upfpot, error=error)
|
||||
IF (iounit > 0) THEN
|
||||
WRITE (iounit, "(/,T2,'PP Transformation for ',A)") TRIM(upfpot%pname)
|
||||
IF (sgppot%has_local) THEN
|
||||
WRITE (iounit, "(T8,'Accuracy for local part:',T61,F20.12)") error(1)
|
||||
END IF
|
||||
IF (sgppot%has_nonlocal) THEN
|
||||
WRITE (iounit, "(T8,'Accuracy for nonlocal part:',T61,F20.12)") error(2)
|
||||
END IF
|
||||
IF (sgppot%has_nlcc) THEN
|
||||
WRITE (iounit, "(T8,'Accuracy for NLCC density:',T61,F20.12)") error(3)
|
||||
END IF
|
||||
END IF
|
||||
IF (sgppot%has_nonlocal) THEN
|
||||
CALL set_potential(qs_kind%sgp_potential, n_nonlocal=sgppot%n_nonlocal, lmax=sgppot%lmax, &
|
||||
is_nonlocal=sgppot%is_nonlocal)
|
||||
|
|
|
|||
|
|
@ -527,7 +527,7 @@ CONTAINS
|
|||
END IF
|
||||
|
||||
! Pseudopotentials
|
||||
IF (gth_potential_present) THEN
|
||||
IF (gth_potential_present .OR. sgp_potential_present) THEN
|
||||
IF (ASSOCIATED(gth_potential)) THEN
|
||||
CALL get_potential(potential=gth_potential, &
|
||||
ppl_present=ppl_present(ikind), &
|
||||
|
|
|
|||
|
|
@ -395,6 +395,7 @@ CONTAINS
|
|||
|
||||
DO ispin = 1, nspin
|
||||
|
||||
NULLIFY (mo_coeff)
|
||||
CALL get_mo_set(mo_set=mo_array(ispin)%mo_set, mo_coeff_b=mo_coeff, mo_coeff=mo_coeff_fm)
|
||||
CALL copy_fm_to_dbcsr(mo_coeff_fm, mo_coeff) !fm -> dbcsr
|
||||
|
||||
|
|
|
|||
|
|
@ -1316,9 +1316,10 @@ CONTAINS
|
|||
CHARACTER(len=*), PARAMETER :: routineN = 'init_sgp_potential', &
|
||||
routineP = moduleN//':'//routineN
|
||||
|
||||
INTEGER :: i1, i2, i3, j1, j2, j3, l, la, lb, lc, &
|
||||
n1, n2, n3, nnl, nprj
|
||||
REAL(KIND=dp), DIMENSION(:, :), POINTER :: cprj
|
||||
INTEGER :: i1, i2, j1, j2, l, la, lb, n1, n2, nnl, &
|
||||
nprj
|
||||
INTEGER, ALLOCATABLE, DIMENSION(:, :) :: ind1, ind2
|
||||
REAL(KIND=dp), DIMENSION(:, :), POINTER :: cprj, hnl
|
||||
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: cn
|
||||
|
||||
IF (ASSOCIATED(potential)) THEN
|
||||
|
|
@ -1336,46 +1337,72 @@ CONTAINS
|
|||
cprj => potential%cprj_ppnl
|
||||
cprj = 0.0_dp
|
||||
cn => potential%c_nonlocal
|
||||
!
|
||||
nnl = potential%n_nonlocal
|
||||
!
|
||||
n1 = 0
|
||||
DO la = 0, potential%lmax
|
||||
n1 = n1+nnl*nco(la)
|
||||
END DO
|
||||
ALLOCATE (ind1(n1, 3))
|
||||
n1 = 0
|
||||
DO i1 = 1, nnl
|
||||
DO la = 0, potential%lmax
|
||||
DO j1 = 1, nco(la)
|
||||
n1 = n1+1
|
||||
!
|
||||
n2 = 0
|
||||
DO i2 = 1, nnl
|
||||
DO lb = 0, potential%lmax
|
||||
IF (.NOT. potential%is_nonlocal(lb)) CYCLE
|
||||
DO j2 = 1, nso(lb)
|
||||
n2 = n2+1
|
||||
IF (la /= lb) CYCLE
|
||||
!
|
||||
n3 = 0
|
||||
DO i3 = 1, nnl
|
||||
DO lc = 0, potential%lmax
|
||||
DO j3 = 1, nso(lc)
|
||||
n3 = n3+1
|
||||
!
|
||||
IF (i1 /= i3) CYCLE
|
||||
IF (la /= lc) CYCLE
|
||||
cprj(n1, n2) = cprj(n1, n2)+orbtramat(la)%c2s(j3, j1)*cn(n3, n2, la)
|
||||
END DO
|
||||
END DO
|
||||
END DO
|
||||
END DO
|
||||
END DO
|
||||
END DO
|
||||
ind1(n1, 1) = la
|
||||
ind1(n1, 2) = j1
|
||||
ind1(n1, 3) = i1
|
||||
END DO
|
||||
END DO
|
||||
END DO
|
||||
!
|
||||
n2 = 0
|
||||
DO i2 = 1, nnl
|
||||
DO lb = 0, potential%lmax
|
||||
IF (.NOT. potential%is_nonlocal(lb)) CYCLE
|
||||
n2 = n2+nso(lb)
|
||||
END DO
|
||||
END DO
|
||||
ALLOCATE (ind2(n2, 3))
|
||||
n2 = 0
|
||||
DO i2 = 1, nnl
|
||||
DO lb = 0, potential%lmax
|
||||
IF (.NOT. potential%is_nonlocal(lb)) CYCLE
|
||||
DO j2 = 1, nso(lb)
|
||||
n2 = n2+1
|
||||
ind2(n2, 1) = lb
|
||||
ind2(n2, 2) = j2
|
||||
ind2(n2, 3) = i2
|
||||
END DO
|
||||
END DO
|
||||
END DO
|
||||
!
|
||||
DO n1 = 1, SIZE(ind1, 1)
|
||||
la = ind1(n1, 1)
|
||||
j1 = ind1(n1, 2)
|
||||
i1 = ind1(n1, 3)
|
||||
DO n2 = 1, SIZE(ind2, 1)
|
||||
lb = ind2(n2, 1)
|
||||
IF (la /= lb) CYCLE
|
||||
j2 = ind2(n2, 2)
|
||||
i2 = ind2(n2, 3)
|
||||
cprj(n1, n2) = orbtramat(la)%c2s(j2, j1)*cn(i1, i2, la)
|
||||
END DO
|
||||
END DO
|
||||
!
|
||||
hnl => potential%h_nonlocal
|
||||
IF (ASSOCIATED(potential%vprj_ppnl)) THEN
|
||||
DEALLOCATE (potential%vprj_ppnl)
|
||||
END IF
|
||||
ALLOCATE (potential%vprj_ppnl(nprj))
|
||||
potential%vprj_ppnl = 0.0_dp
|
||||
DO n2 = 1, SIZE(ind2, 1)
|
||||
lb = ind2(n2, 1)
|
||||
i2 = ind2(n2, 3)
|
||||
potential%vprj_ppnl(n2) = hnl(i2, lb)
|
||||
END DO
|
||||
!
|
||||
DEALLOCATE (ind1, ind2)
|
||||
END IF
|
||||
END IF
|
||||
|
||||
|
|
|
|||
|
|
@ -46,7 +46,7 @@
|
|||
l_max_rho="0"
|
||||
l_local="-3"
|
||||
mesh_size="1073"
|
||||
number_of_wfc=" 2"
|
||||
number_of_wfc="2"
|
||||
number_of_proj="1"/>
|
||||
<PP_MESH dx="1.250000000000000E-002" mesh="1073" xmin="-7.000000000000000E+000" rmax="1.003075063120137E+002"
|
||||
zmesh="6.000000000000000E+000">
|
||||
|
|
|
|||
58
tests/QS/regtest-kind/CH2O_upf.inp
Normal file
58
tests/QS/regtest-kind/CH2O_upf.inp
Normal file
|
|
@ -0,0 +1,58 @@
|
|||
&GLOBAL
|
||||
PROJECT H2
|
||||
RUN_TYPE energy
|
||||
PRINT_LEVEL LOW
|
||||
&END GLOBAL
|
||||
|
||||
&FORCE_EVAL
|
||||
METHOD Quickstep
|
||||
|
||||
&DFT
|
||||
BASIS_SET_FILE_NAME GTH_BASIS_SETS
|
||||
|
||||
&MGRID
|
||||
CUTOFF 400
|
||||
&END MGRID
|
||||
|
||||
&QS
|
||||
METHOD GAPW
|
||||
&END QS
|
||||
|
||||
&SCF
|
||||
MAX_SCF 20
|
||||
EPS_SCF 1e-7
|
||||
SCF_GUESS atomic
|
||||
&END SCF
|
||||
|
||||
&XC
|
||||
&XC_FUNCTIONAL PADE
|
||||
&END XC_FUNCTIONAL
|
||||
&END XC
|
||||
&END DFT
|
||||
|
||||
&SUBSYS
|
||||
&CELL
|
||||
ABC 8.0 8.0 8.0
|
||||
&END CELL
|
||||
|
||||
&COORD
|
||||
O 0.9588431900 1.1234806613 1.8643358699
|
||||
C 1.0045827842 1.0372747429 0.6713062328
|
||||
H 1.0304990091 1.9328340670 0.0202209074
|
||||
H 1.0234151411 0.0574091066 0.1554705073
|
||||
&END COORD
|
||||
|
||||
&KIND O
|
||||
BASIS_SET aug-TZV2P-GTH
|
||||
POTENTIAL UPF O.pbe-hgh.UPF
|
||||
&END KIND
|
||||
&KIND C
|
||||
BASIS_SET aug-TZV2P-GTH
|
||||
POTENTIAL UPF C.pbe-hgh.UPF
|
||||
&END KIND
|
||||
&KIND H
|
||||
BASIS_SET aug-TZV2P-GTH
|
||||
POTENTIAL UPF H.pbe-hgh.UPF
|
||||
&END KIND
|
||||
&END SUBSYS
|
||||
&END FORCE_EVAL
|
||||
|
|
@ -46,7 +46,7 @@
|
|||
l_max_rho="0"
|
||||
l_local="-3"
|
||||
mesh_size="929"
|
||||
number_of_wfc=" 1"
|
||||
number_of_wfc="1"
|
||||
number_of_proj="0"/>
|
||||
<PP_MESH dx="1.250000000000000E-002" mesh="929" xmin="-7.000000000000000E+000" rmax="9.948431564193386E+001"
|
||||
zmesh="1.000000000000000E+000">
|
||||
|
|
|
|||
53
tests/QS/regtest-kind/H2O_upf.inp
Normal file
53
tests/QS/regtest-kind/H2O_upf.inp
Normal file
|
|
@ -0,0 +1,53 @@
|
|||
&GLOBAL
|
||||
PROJECT H2
|
||||
RUN_TYPE energy
|
||||
PRINT_LEVEL LOW
|
||||
&END GLOBAL
|
||||
|
||||
&FORCE_EVAL
|
||||
METHOD Quickstep
|
||||
|
||||
&DFT
|
||||
BASIS_SET_FILE_NAME GTH_BASIS_SETS
|
||||
|
||||
&MGRID
|
||||
CUTOFF 400
|
||||
&END MGRID
|
||||
|
||||
&QS
|
||||
METHOD GAPW
|
||||
&END QS
|
||||
|
||||
&SCF
|
||||
MAX_SCF 20
|
||||
EPS_SCF 1e-7
|
||||
SCF_GUESS atomic
|
||||
&END SCF
|
||||
|
||||
&XC
|
||||
&XC_FUNCTIONAL PADE
|
||||
&END XC_FUNCTIONAL
|
||||
&END XC
|
||||
&END DFT
|
||||
|
||||
&SUBSYS
|
||||
&CELL
|
||||
ABC 8.0 8.0 8.0
|
||||
&END CELL
|
||||
|
||||
&COORD
|
||||
O 4.000000 4.000000 4.000000
|
||||
H 4.000000 3.250000 4.520000
|
||||
H 4.000000 4.750000 4.520000
|
||||
&END COORD
|
||||
|
||||
&KIND O
|
||||
BASIS_SET aug-TZV2P-GTH
|
||||
POTENTIAL UPF O.pbe-hgh.UPF
|
||||
&END KIND
|
||||
&KIND H
|
||||
BASIS_SET aug-TZV2P-GTH
|
||||
POTENTIAL UPF H.pbe-hgh.UPF
|
||||
&END KIND
|
||||
&END SUBSYS
|
||||
&END FORCE_EVAL
|
||||
|
|
@ -46,7 +46,7 @@
|
|||
l_max_rho="0"
|
||||
l_local="-3"
|
||||
mesh_size="1095"
|
||||
number_of_wfc=" 2"
|
||||
number_of_wfc="2"
|
||||
number_of_proj="1"/>
|
||||
<PP_MESH dx="1.250000000000000E-002" mesh="1095" xmin="-7.000000000000000E+000" rmax="9.904343173443276E+001"
|
||||
zmesh="8.000000000000000E+000">
|
||||
|
|
|
|||
|
|
@ -12,6 +12,8 @@ ch2o_6.inp 1 3e-13
|
|||
ch2o_7.inp 1 3e-13 -22.56644338703015
|
||||
ch2o_8.inp 1 3e-13 -22.56644338703015
|
||||
ch2o_9.inp 1 3e-13 -22.56644338703015
|
||||
# disable unstable test
|
||||
#H2_upf.inp 1 3e-13 -1.135018770048071
|
||||
# tests maybe unstable
|
||||
H2_upf.inp 1 1e-06 -1.13501864985782
|
||||
H2O_upf.inp 1 5e-05 -17.11484046900296
|
||||
CH2O_upf.inp 1 5e-05 -22.74205305294306
|
||||
#EOF
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue