diff --git a/src/atom_kind_orbitals.F b/src/atom_kind_orbitals.F index 99a10fd079..e24fe1b65b 100644 --- a/src/atom_kind_orbitals.F +++ b/src/atom_kind_orbitals.F @@ -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 diff --git a/src/atom_sgp.F b/src/atom_sgp.F index 7507ff14f5..91366ed087 100644 --- a/src/atom_sgp.F +++ b/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 diff --git a/src/core_ppnl.F b/src/core_ppnl.F index ff9d89bb63..d674ea59ba 100644 --- a/src/core_ppnl.F +++ b/src/core_ppnl.F @@ -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) diff --git a/src/qs_kind_types.F b/src/qs_kind_types.F index 4e10064697..2e3a1bfc99 100644 --- a/src/qs_kind_types.F +++ b/src/qs_kind_types.F @@ -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) diff --git a/src/qs_neighbor_lists.F b/src/qs_neighbor_lists.F index feaf8d87f9..b72c423db1 100644 --- a/src/qs_neighbor_lists.F +++ b/src/qs_neighbor_lists.F @@ -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), & diff --git a/src/qs_ot_scf.F b/src/qs_ot_scf.F index 74e2b125fb..90e7b4ff95 100644 --- a/src/qs_ot_scf.F +++ b/src/qs_ot_scf.F @@ -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 diff --git a/src/subsys/external_potential_types.F b/src/subsys/external_potential_types.F index bbffff7f1b..49c1de566f 100644 --- a/src/subsys/external_potential_types.F +++ b/src/subsys/external_potential_types.F @@ -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 diff --git a/tests/QS/regtest-kind/C.pbe-hgh.UPF b/tests/QS/regtest-kind/C.pbe-hgh.UPF index 46f988c4d1..2daca6a6d0 100644 --- a/tests/QS/regtest-kind/C.pbe-hgh.UPF +++ b/tests/QS/regtest-kind/C.pbe-hgh.UPF @@ -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"/> diff --git a/tests/QS/regtest-kind/CH2O_upf.inp b/tests/QS/regtest-kind/CH2O_upf.inp new file mode 100644 index 0000000000..809e8fdf51 --- /dev/null +++ b/tests/QS/regtest-kind/CH2O_upf.inp @@ -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 diff --git a/tests/QS/regtest-kind/H.pbe-hgh.UPF b/tests/QS/regtest-kind/H.pbe-hgh.UPF index b8f90f2acf..51d3302feb 100644 --- a/tests/QS/regtest-kind/H.pbe-hgh.UPF +++ b/tests/QS/regtest-kind/H.pbe-hgh.UPF @@ -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"/> diff --git a/tests/QS/regtest-kind/H2O_upf.inp b/tests/QS/regtest-kind/H2O_upf.inp new file mode 100644 index 0000000000..9ddf31b051 --- /dev/null +++ b/tests/QS/regtest-kind/H2O_upf.inp @@ -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 diff --git a/tests/QS/regtest-kind/O.pbe-hgh.UPF b/tests/QS/regtest-kind/O.pbe-hgh.UPF index 3b694e900f..23f03e41c2 100644 --- a/tests/QS/regtest-kind/O.pbe-hgh.UPF +++ b/tests/QS/regtest-kind/O.pbe-hgh.UPF @@ -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"/> diff --git a/tests/QS/regtest-kind/TEST_FILES b/tests/QS/regtest-kind/TEST_FILES index 82d6eb7d70..038ac8e745 100644 --- a/tests/QS/regtest-kind/TEST_FILES +++ b/tests/QS/regtest-kind/TEST_FILES @@ -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