comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-29 06:35:28 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Atom code: new property calculations for basis sets:

Browse source 
Condition number of overlap matrix in a cubic cell
Completeness of basis wrt Slater orbitals with different exponents
3 new regtests


svn-origin-rev: 16661
This commit is contained in:
Jürg Hutter 2016-03-07 11:15:57 +00:00
parent 76eb4310e0
commit 9a797915cf
10 changed files with 637 additions and 21 deletions
Download patch file Download diff file Expand all files Collapse all files

182
src/aobasis/ai_overlap.F
View file

@ -16,9 +16,13 @@
MODULE ai_overlap
USE ai_os_rr, ONLY: os_rr_ovlp
USE kinds, ONLY: dp
USE mathconstants, ONLY: pi
USE mathconstants, ONLY: pi,&
twopi
USE orbital_pointers, ONLY: coset,&
ncoset
nco,&
ncoset,&
nso
USE orbital_transformation_matrices, ONLY: orbtramat
#include "../base/base_uses.f90"
IMPLICIT NONE
@ -28,7 +32,7 @@ MODULE ai_overlap
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_overlap'
! *** Public subroutines ***
PUBLIC :: overlap, overlap_ab, overlap_aab
PUBLIC :: overlap, overlap_ab, overlap_aab, overlap_ab_s, overlap_ab_sp
CONTAINS
@ -1941,4 +1945,176 @@ CONTAINS
END SUBROUTINE overlap_abbb
! *****************************************************************************
!> \brief Calculation of the two-center overlap integrals [a|b] over
!> Spherical Gaussian-type functions.
!> \param la Max L on center A
!> \param zeta Exponents on center A
!> \param lb Max L on center B
!> \param zetb Exponents on center B
!> \param rab Distance vector A-B
!> \param sab Final overlap integrals
!> \date 01.03.2016
!> \author JGH
! *****************************************************************************
SUBROUTINE overlap_ab_s(la,zeta,lb,zetb,rab,sab)
INTEGER, INTENT(IN) :: la
REAL(KIND=dp), INTENT(IN) :: zeta
INTEGER, INTENT(IN) :: lb
REAL(KIND=dp), INTENT(IN) :: zetb
REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
REAL(KIND=dp), DIMENSION(:, :), &
INTENT(INOUT) :: sab
CHARACTER(len=*), PARAMETER :: routineN = 'overlap_ab_s', &
routineP = moduleN//':'//routineN
INTEGER :: nca, ncb, nsa, nsb
REAL(KIND=dp), ALLOCATABLE, &
DIMENSION(:, :) :: cab
REAL(KIND=dp), DIMENSION(1) :: rpgf, za, zb
REAL(KIND=dp), DIMENSION(:, :), POINTER :: c2sa, c2sb
rpgf(1) = HUGE(1._dp)
za(1) = zeta
zb(1) = zetb
nca = nco(la)
ncb = nco(lb)
ALLOCATE(cab(nca,ncb))
nsa = nso(la)
nsb = nso(lb)
CALL overlap_ab(la,la,1,rpgf,za,lb,lb,1,rpgf,zb,rab,cab)
c2sa => orbtramat(la)%c2s
c2sb => orbtramat(lb)%c2s
sab(1:nsa,1:nsb) = MATMUL(c2sa(1:nsa,1:nca),&
MATMUL(cab(1:nca,1:ncb),TRANSPOSE(c2sb(1:nsb,1:ncb))))
DEALLOCATE(cab)
END SUBROUTINE overlap_ab_s
! *****************************************************************************
!> \brief Calculation of the overlap integrals [a|b] over
!> cubic periodic Spherical Gaussian-type functions.
!> \param la Max L on center A
!> \param zeta Exponents on center A
!> \param lb Max L on center B
!> \param zetb Exponents on center B
!> \param alat Lattice constant
!> \param sab Final overlap integrals
!> \date 01.03.2016
!> \author JGH
! *****************************************************************************
SUBROUTINE overlap_ab_sp(la,zeta,lb,zetb,alat,sab)
INTEGER, INTENT(IN) :: la
REAL(KIND=dp), INTENT(IN) :: zeta
INTEGER, INTENT(IN) :: lb
REAL(KIND=dp), INTENT(IN) :: zetb, alat
REAL(KIND=dp), DIMENSION(:, :), &
INTENT(INOUT) :: sab
CHARACTER(len=*), PARAMETER :: routineN = 'overlap_ab_sp', &
routineP = moduleN//':'//routineN
COMPLEX(KIND=dp) :: zfg
COMPLEX(KIND=dp), ALLOCATABLE, &
DIMENSION(:, :) :: fun, gun
INTEGER :: ax, ay, az, bx, by, bz, i, &
ia, ib, l, l1, l2, na, nb, &
nca, ncb, nmax, nsa, nsb
REAL(KIND=dp) :: oa, ob, ovol, zm
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: fexp, gexp, gval
REAL(KIND=dp), ALLOCATABLE, &
DIMENSION(:, :) :: cab
REAL(KIND=dp), DIMENSION(0:3, 0:3) :: fgsum
REAL(KIND=dp), DIMENSION(:, :), POINTER :: c2sa, c2sb
nca = nco(la)
ncb = nco(lb)
ALLOCATE(cab(nca,ncb))
cab = 0.0_dp
nsa = nso(la)
nsb = nso(lb)
zm = MIN(zeta,zetb)
nmax = NINT(1.81_dp*alat*SQRT(zm) + 1.0_dp)
ALLOCATE(fun(-nmax:nmax,0:la),gun(-nmax:nmax,0:lb),&
fexp(-nmax:nmax),gexp(-nmax:nmax),gval(-nmax:nmax))
oa = 1._dp/zeta
ob = 1._dp/zetb
DO i=-nmax,nmax
gval(i) = twopi/alat*REAL(i,KIND=dp)
fexp(i) = SQRT(oa*pi)*EXP(-0.25_dp*oa*gval(i)**2)
gexp(i) = SQRT(ob*pi)*EXP(-0.25_dp*ob*gval(i)**2)
END DO
DO l=0,la
IF(l==0) THEN
fun(:,l) = CMPLX(1.0_dp,0.0_dp,KIND=dp)
ELSEIF(l==1) THEN
fun(:,l) = CMPLX(0.0_dp,0.5_dp*oa*gval(:),KIND=dp)
ELSEIF(l==2) THEN
fun(:,l) = CMPLX(-(0.5_dp*oa*gval(:))**2,0.0_dp,KIND=dp)
fun(:,l) = fun(:,l) + CMPLX(0.5_dp*oa,0.0_dp,KIND=dp)
ELSEIF(l==3) THEN
fun(:,l) = CMPLX(0.0_dp,-(0.5_dp*oa*gval(:))**3,KIND=dp)
fun(:,l) = fun(:,l) + CMPLX(0.0_dp,0.75_dp*oa*oa*gval(:),KIND=dp)
ELSE
CPABORT("l value too high")
END IF
END DO
DO l=0,lb
IF(l==0) THEN
gun(:,l) = CMPLX(1.0_dp,0.0_dp,KIND=dp)
ELSEIF(l==1) THEN
gun(:,l) = CMPLX(0.0_dp,0.5_dp*ob*gval(:),KIND=dp)
ELSEIF(l==2) THEN
gun(:,l) = CMPLX(-(0.5_dp*ob*gval(:))**2,0.0_dp,KIND=dp)
gun(:,l) = gun(:,l) + CMPLX(0.5_dp*ob,0.0_dp,KIND=dp)
ELSEIF(l==3) THEN
gun(:,l) = CMPLX(0.0_dp,-(0.5_dp*ob*gval(:))**3,KIND=dp)
gun(:,l) = gun(:,l) + CMPLX(0.0_dp,0.75_dp*ob*ob*gval(:),KIND=dp)
ELSE
CPABORT("l value too high")
END IF
END DO
fgsum = 0.0_dp
DO l1=0,la
DO l2=0,lb
zfg = SUM(CONJG(fun(:,l1))*fexp(:)*gun(:,l2)*gexp(:))
fgsum(l1,l2) = REAL(zfg,KIND=dp)
END DO
END DO
na = ncoset(la-1)
nb = ncoset(lb-1)
DO ax = 0,la
DO ay = 0,la-ax
az = la - ax - ay
ia = coset(ax,ay,az) - na
DO bx = 0,lb
DO by = 0,lb-bx
bz = lb - bx - by
ib = coset(bx,by,bz) - nb
cab(ia,ib) = fgsum(ax,bx)*fgsum(ay,by)*fgsum(az,bz)
END DO
END DO
END DO
END DO
c2sa => orbtramat(la)%c2s
c2sb => orbtramat(lb)%c2s
sab(1:nsa,1:nsb) = MATMUL(c2sa(1:nsa,1:nca),&
MATMUL(cab(1:nca,1:ncb),TRANSPOSE(c2sb(1:nsb,1:ncb))))
ovol = 1._dp/(alat**3)
sab(1:nsa,1:nsb) = ovol*sab(1:nsa,1:nsb)
DEALLOCATE(cab,fun,gun,fexp,gexp,gval)
END SUBROUTINE overlap_ab_sp
END MODULE ai_overlap

149
src/atom_basis.F
View file

@ -5,17 +5,20 @@
! *****************************************************************************
MODULE atom_basis
USE ai_overlap, ONLY: overlap_ab_s,&
overlap_ab_sp
USE ao_util, ONLY: exp_radius
USE atom_fit, ONLY: atom_fit_basis
USE atom_output, ONLY: atom_print_basis,&
atom_print_info,&
atom_print_method,&
atom_print_potential
USE atom_types, ONLY: &
atom_basis_type, atom_integrals, atom_optimization_type, &
atom_orbitals, atom_p_type, atom_potential_type, atom_state, &
create_atom_orbs, create_atom_type, init_atom_basis, &
init_atom_potential, read_atom_opt_section, release_atom_basis, &
release_atom_potential, release_atom_type, set_atom
CGTO_BASIS, GTO_BASIS, atom_basis_type, atom_integrals, &
atom_optimization_type, atom_orbitals, atom_p_type, &
atom_potential_type, atom_state, create_atom_orbs, create_atom_type, &
init_atom_basis, init_atom_potential, read_atom_opt_section, &
release_atom_basis, release_atom_potential, release_atom_type, set_atom
USE atom_utils, ONLY: atom_consistent_method,&
atom_set_occupation,&
get_maxl_occ,&
@ -31,13 +34,14 @@ MODULE atom_basis
section_vals_val_get
USE kinds, ONLY: default_string_length,&
dp
USE lapack, ONLY: lapack_ssyev
USE periodic_table, ONLY: nelem,&
ptable
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
PUBLIC :: atom_basis_opt
PUBLIC :: atom_basis_opt, atom_basis_condnum
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_basis'
@ -279,5 +283,138 @@ CONTAINS
END SUBROUTINE atom_basis_opt
! *****************************************************************************
!> \brief ...
!> \param basis ...
!> \param rad ...
!> \param cnum ...
! *****************************************************************************
SUBROUTINE atom_basis_condnum(basis,rad,cnum)
TYPE(atom_basis_type), POINTER :: basis
REAL(KIND=dp), INTENT(IN) :: rad
REAL(KIND=dp), INTENT(OUT) :: cnum
CHARACTER(len=*), PARAMETER :: routineN = 'atom_basis_condnum', &
routineP = moduleN//':'//routineN
INTEGER :: ia, ib, imax, info, ix, iy, &
iz, ja, jb, ka, kb, l, la, &
lb, lwork, na, nb, nbas, nna, &
nnb
INTEGER, ALLOCATABLE, DIMENSION(:, :) :: ibptr
REAL(KIND=dp) :: r1, r2, reps, rmax
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: weig, work
REAL(KIND=dp), ALLOCATABLE, &
DIMENSION(:, :) :: smat
REAL(KIND=dp), DIMENSION(3) :: rab
REAL(KIND=dp), DIMENSION(7, 7) :: sab
REAL(KIND=dp), DIMENSION(:), POINTER :: zeta, zetb
! total number of basis functions
nbas = 0
DO l=0,3
nbas = nbas + basis%nbas(l) * (2*l+1)
END DO
ALLOCATE(smat(nbas,nbas),ibptr(nbas,0:3))
smat = 0.0_dp
ibptr = 0
na = 0
DO l=0,3
DO ia=1,basis%nbas(l)
ibptr(ia,l) = na + 1
na = na + (2*l+1)
END DO
END DO
reps = 1.e-14_dp
IF(basis%basis_type == GTO_BASIS .OR. &
basis%basis_type == CGTO_BASIS) THEN
DO la=0,3
na = basis%nprim(la)
nna = 2*la+1
IF(na==0) CYCLE
zeta => basis%am(:,la)
DO lb=0,3
nb = basis%nprim(lb)
nnb = 2*lb+1
IF(nb==0) CYCLE
zetb => basis%am(:,lb)
DO ia=1,na
DO ib=1,nb
IF(rad < 0.1_dp) THEN
imax = 0
ELSE
r1 = exp_radius(la,zeta(ia),reps,1.0_dp)
r2 = exp_radius(lb,zetb(ib),reps,1.0_dp)
rmax = MAX(2._dp*r1,2._dp*r2)
imax = INT(rmax/rad) + 1
END IF
IF(imax > 1) THEN
CALL overlap_ab_sp(la,zeta(ia),lb,zetb(ib),rad,sab)
IF(basis%basis_type == GTO_BASIS) THEN
ja = ibptr(ia,la)
jb = ibptr(ib,lb)
smat(ja:ja+nna-1,jb:jb+nnb-1) = smat(ja:ja+nna-1,jb:jb+nnb-1) + sab(1:nna,1:nnb)
ELSEIF(basis%basis_type == CGTO_BASIS) THEN
DO ka=1,basis%nbas(la)
DO kb=1,basis%nbas(lb)
ja = ibptr(ka,la)
jb = ibptr(kb,lb)
smat(ja:ja+nna-1,jb:jb+nnb-1) = smat(ja:ja+nna-1,jb:jb+nnb-1) +&
sab(1:nna,1:nnb)*basis%cm(ia,ka,la)*basis%cm(ib,kb,lb)
END DO
END DO
END IF
ELSE
DO ix = -imax,imax
rab(1) = rad*ix
DO iy = -imax,imax
rab(2) = rad*iy
DO iz = -imax,imax
rab(3) = rad*iz
CALL overlap_ab_s(la,zeta(ia),lb,zetb(ib),rab,sab)
IF(basis%basis_type == GTO_BASIS) THEN
ja = ibptr(ia,la)
jb = ibptr(ib,lb)
smat(ja:ja+nna-1,jb:jb+nnb-1) = smat(ja:ja+nna-1,jb:jb+nnb-1) + sab(1:nna,1:nnb)
ELSEIF(basis%basis_type == CGTO_BASIS) THEN
DO ka=1,basis%nbas(la)
DO kb=1,basis%nbas(lb)
ja = ibptr(ka,la)
jb = ibptr(kb,lb)
smat(ja:ja+nna-1,jb:jb+nnb-1) = smat(ja:ja+nna-1,jb:jb+nnb-1) +&
sab(1:nna,1:nnb)*basis%cm(ia,ka,la)*basis%cm(ib,kb,lb)
END DO
END DO
END IF
END DO
END DO
END DO
END IF
END DO
END DO
END DO
END DO
ELSE
CPABORT("Condition number not available for this basis type")
END IF
info = 0
lwork = nbas*nbas
ALLOCATE(weig(nbas),work(lwork))
CALL lapack_ssyev ( "N", "U", nbas, smat, nbas, weig, work, lwork, info )
CPASSERT(info==0)
IF(weig(1) < 0.0_dp) THEN
cnum = 100._dp
ELSE
cnum = ABS(weig(nbas)/weig(1))
cnum = LOG10(cnum)
END IF
DEALLOCATE(smat,ibptr,weig,work)
END SUBROUTINE atom_basis_condnum
END MODULE atom_basis

221
src/atom_energy.F
View file

@ -5,6 +5,9 @@
! *****************************************************************************
MODULE atom_energy
USE ai_onecenter, ONLY: sg_overlap,&
sto_overlap
USE atom_basis, ONLY: atom_basis_condnum
USE atom_electronic_structure, ONLY: calculate_atom
USE atom_fit, ONLY: atom_fit_density,&
atom_fit_kgpot
@ -13,7 +16,8 @@ MODULE atom_energy
atom_ppint_release,&
atom_ppint_setup,&
atom_relint_release,&
atom_relint_setup
atom_relint_setup,&
contract2
USE atom_output, ONLY: atom_print_basis,&
atom_print_info,&
atom_print_method,&
@ -21,12 +25,12 @@ MODULE atom_energy
atom_print_potential,&
atom_write_pseudo_param
USE atom_types, ONLY: &
atom_basis_type, atom_gthpot_type, atom_integrals, &
atom_optimization_type, atom_orbitals, atom_p_type, &
atom_potential_type, atom_state, atom_type, create_atom_orbs, &
create_atom_type, gth_pseudo, init_atom_basis, init_atom_potential, &
read_atom_opt_section, release_atom_basis, release_atom_potential, &
release_atom_type, set_atom
CGTO_BASIS, GTO_BASIS, NUM_BASIS, STO_BASIS, atom_basis_type, &
atom_gthpot_type, atom_integrals, atom_optimization_type, &
atom_orbitals, atom_p_type, atom_potential_type, atom_state, &
atom_type, create_atom_orbs, create_atom_type, gth_pseudo, &
init_atom_basis, init_atom_potential, read_atom_opt_section, &
release_atom_basis, release_atom_potential, release_atom_type, set_atom
USE atom_utils, ONLY: &
atom_consistent_method, atom_core_density, atom_density, &
atom_local_potential, atom_read_external_density, &
@ -34,6 +38,7 @@ MODULE atom_energy
get_maxl_occ, get_maxn_occ
USE atom_xc, ONLY: calculate_atom_ext_vxc,&
calculate_atom_zmp
USE basis_set_types, ONLY: srules
USE cp_log_handling, ONLY: cp_get_default_logger,&
cp_logger_type
USE cp_output_handling, ONLY: cp_print_key_finished_output,&
@ -46,11 +51,18 @@ MODULE atom_energy
USE kinds, ONLY: default_string_length,&
dp
USE mathconstants, ONLY: dfac,&
fac,&
gamma1,&
pi
USE mathlib, ONLY: invmat_symm
USE orbital_pointers, ONLY: deallocate_orbital_pointers,&
init_orbital_pointers
USE orbital_transformation_matrices, ONLY: deallocate_spherical_harmonics,&
init_spherical_harmonics
USE periodic_table, ONLY: nelem,&
ptable
USE physcon, ONLY: evolt
USE physcon, ONLY: bohr,&
evolt
#include "./base/base_uses.f90"
IMPLICIT NONE
@ -85,9 +97,9 @@ CONTAINS
INTEGER, DIMENSION(0:3) :: maxn
INTEGER, DIMENSION(:), POINTER :: cn
LOGICAL :: dm, do_zmp, doread, eri_c, &
eri_e, had_ae, had_pp, &
pp_calc, read_vxc
REAL(KIND=dp) :: delta, lambda
eri_e, had_ae, had_pp, lcomp, &
lcond, pp_calc, read_vxc
REAL(KIND=dp) :: crad, delta, lambda
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ext_density, ext_vxc
REAL(KIND=dp), DIMENSION(0:3, 10) :: pocc
TYPE(atom_basis_type), POINTER :: ae_basis, pp_basis
@ -437,6 +449,23 @@ CONTAINS
END DO
END DO
! Analyze basis sets
iw = cp_print_key_unit_nr(logger,atom_section,"PRINT%ANALYZE_BASIS",extension=".log")
IF (iw>0) THEN
CALL section_vals_val_get(atom_section,"PRINT%ANALYZE_BASIS%OVERLAP_CONDITION_NUMBER", l_val=lcond)
CALL section_vals_val_get(atom_section,"PRINT%ANALYZE_BASIS%COMPLETENESS", l_val=lcomp)
crad = ptable(zval)%covalent_radius*bohr
IF ( had_ae ) THEN
IF(lcond) CALL atom_condnumber(ae_basis,crad,iw)
IF(lcomp) CALL atom_completeness(ae_basis,zval,iw)
END IF
IF ( had_pp ) THEN
IF(lcond) CALL atom_condnumber(pp_basis,crad,iw)
IF(lcomp) CALL atom_completeness(pp_basis,zval,iw)
END IF
END IF
CALL cp_print_key_finished_output(iw,logger,atom_section,"PRINT%ANALYZE_BASIS")
! clean up
IF ( had_ae ) THEN
CALL atom_int_release(ae_int)
@ -899,5 +928,175 @@ CONTAINS
END IF
END SUBROUTINE compose
! *****************************************************************************
!> \brief ...
!> \param basis ...
!> \param crad ...
!> \param iw ...
! *****************************************************************************
SUBROUTINE atom_condnumber(basis,crad,iw)
TYPE(atom_basis_type), POINTER :: basis
REAL(KIND=dp) :: crad
INTEGER, INTENT(IN) :: iw
INTEGER :: i
REAL(KIND=dp) :: ci
REAL(KIND=dp), DIMENSION(10) :: cnum, rad
WRITE(iw,'(/,A)') " Basis Set Condition Numbers"
CALL init_orbital_pointers(3)
CALL init_spherical_harmonics(3,0)
cnum = 0.0_dp
DO i=1,9
ci = 2.0_dp * (0.8_dp + i*0.1_dp)
rad(i) = crad*ci
CALL atom_basis_condnum(basis,rad(i),cnum(i))
WRITE(iw,'(A,F15.3,T50,A,F14.4)') " Lattice constant:",&
rad(i),"Condition number:",cnum(i)
END DO
rad(10) = 0.01_dp
CALL atom_basis_condnum(basis,rad(10),cnum(10))
WRITE(iw,'(A,A,T50,A,F14.4)') " Lattice constant:",&
" Inf","Condition number:",cnum(i)
CALL deallocate_orbital_pointers
CALL deallocate_spherical_harmonics
END SUBROUTINE atom_condnumber
! *****************************************************************************
!> \brief ...
!> \param basis ...
!> \param zv ...
!> \param iw ...
! *****************************************************************************
SUBROUTINE atom_completeness(basis,zv,iw)
TYPE(atom_basis_type), POINTER :: basis
INTEGER, INTENT(IN) :: zv, iw
INTEGER :: i, j, l, ll, m, n, nbas, nl, &
nr
INTEGER, DIMENSION(0:3) :: nelem, nlmax, nlmin
INTEGER, DIMENSION(4, 7) :: ne
REAL(KIND=dp) :: al, c1, c2, pf
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: sfun
REAL(KIND=dp), ALLOCATABLE, &
DIMENSION(:, :) :: bmat
REAL(KIND=dp), ALLOCATABLE, &
DIMENSION(:, :, :) :: omat
REAL(KIND=dp), ALLOCATABLE, &
DIMENSION(:, :, :, :) :: sint
REAL(KIND=dp), DIMENSION(0:3, 10) :: snl
REAL(KIND=dp), DIMENSION(2) :: sse
REAL(KIND=dp), DIMENSION(4, 7) :: sexp
ne = 0
nelem(0:3) = ptable(zv)%e_conv(0:3)
DO l=0,3
ll = 2*(2*l+1)
DO i=1,7
IF(nelem(l) >= ll) THEN
ne(l+1,i) = ll
nelem(l) = nelem(l) - ll
ELSE IF(nelem(l) > 0) THEN
ne(l+1,i) = nelem(l)
nelem(l) = 0
ELSE
EXIT
END IF
END DO
END DO
nlmin = 1
nlmax = 1
DO l=0,3
nlmin(l) = l+1
DO i=1,7
IF(ne(l+1,i) > 0) THEN
nlmax(l) = i+l
END IF
END DO
nlmax(l) = MAX(nlmax(l),nlmin(l)+1)
END DO
! Slater exponents
sexp = 0.0_dp
DO l=0,3
sse(1) = 0.05_dp
sse(2) = 10.0_dp
DO i=l+1,7
sexp(l+1,i) = srules(zv,ne,i,l)
IF(ne(l+1,i-l) > 0) THEN
sse(1) = MAX(sse(1),sexp(l+1,i))
sse(2) = MIN(sse(2),sexp(l+1,i))
END IF
END DO
DO i=1,10
snl(l,i) = ABS(2._dp*sse(1) - 0.5_dp*sse(2))/9._dp * REAL(i-1,KIND=dp) + 0.5_dp*MIN(sse(1),sse(2))
END DO
END DO
nbas = MAXVAL(basis%nbas)
ALLOCATE(omat(nbas,nbas,0:3))
nr = SIZE(basis%bf,1)
ALLOCATE(sfun(nr),sint(10,2,nbas,0:3))
sint = 0._dp
! calculate overlaps between test functions and basis
DO l=0,3
DO i=1,10
al = snl(l,i)
nl = nlmin(l)
pf = (2._dp*al)**nl * SQRT(2._dp*al/fac(2*nl))
sfun(1:nr) = pf * basis%grid%rad(1:nr)**(nl-1) * EXP(-al*basis%grid%rad(1:nr))
DO j=1,basis%nbas(l)
sint(i,1,j,l) = SUM(sfun(1:nr)*basis%bf(1:nr,j,l)*basis%grid%wr(1:nr))
END DO
nl = nlmax(l)
pf = (2._dp*al)**nl * SQRT(2._dp*al/fac(2*nl))
sfun(1:nr) = pf * basis%grid%rad(1:nr)**(nl-1) * EXP(-al*basis%grid%rad(1:nr))
DO j=1,basis%nbas(l)
sint(i,2,j,l) = SUM(sfun(1:nr)*basis%bf(1:nr,j,l)*basis%grid%wr(1:nr))
END DO
END DO
END DO
DO l=0,3
n = basis%nbas(l)
IF(n <= 0) CYCLE
m = basis%nprim(l)
SELECT CASE (basis%basis_type)
CASE DEFAULT
CPABORT("")
CASE (GTO_BASIS)
CALL sg_overlap ( omat(1:n,1:n,l), l, basis%am(1:n,l), basis%am(1:n,l) )
CASE (CGTO_BASIS)
ALLOCATE (bmat(m,m))
CALL sg_overlap ( bmat(1:m,1:m), l, basis%am(1:m,l), basis%am(1:m,l) )
CALL contract2(omat(1:n,1:n,l),bmat(1:m,1:m),basis%cm(1:m,1:n,l))
DEALLOCATE (bmat)
CASE (STO_BASIS)
CALL sto_overlap ( omat(1:n,1:n,l), basis%ns(1:n,l), basis%as(1:n,l), &
basis%ns(1:n,l), basis%as(1:n,l) )
CASE (NUM_BASIS)
CPABORT("")
END SELECT
CALL invmat_symm(omat(1:n,1:n,l))
END DO
WRITE(iw,'(/,A)') " Basis Set Completeness Estimates"
DO l=0,3
n = basis%nbas(l)
IF(n <= 0) CYCLE
WRITE(iw,'(A,I3)') " L-quantum number: ",l
WRITE(iw,'(A,T31,A,I2,T61,A,I2)') " Slater Exponent","Completeness n-qm=",nlmin(l),&
"Completeness n-qm=",nlmax(l)
DO i=10,1,-1
c1 = DOT_PRODUCT(sint(i,1,1:n,l),MATMUL(omat(1:n,1:n,l),sint(i,1,1:n,l)))
c2 = DOT_PRODUCT(sint(i,2,1:n,l),MATMUL(omat(1:n,1:n,l),sint(i,2,1:n,l)))
WRITE(iw,"(T6,F14.6,T41,F10.6,T71,F10.6)") snl(l,i),c1,c2
END DO
END DO
DEALLOCATE (omat,sfun,sint)
END SUBROUTINE atom_completeness
END MODULE atom_energy

1
src/atom_operators.F
View file

@ -50,6 +50,7 @@ MODULE atom_operators
PUBLIC :: atom_int_setup, atom_ppint_setup, atom_int_release, atom_ppint_release
PUBLIC :: atom_relint_setup, atom_relint_release
PUBLIC :: contract2
! *****************************************************************************

2
src/atom_types.F
View file

@ -2261,7 +2261,7 @@ CONTAINS
ELSE IF(potential%conf_type==barrier_conf) THEN
potential%acon = 200.0_dp
potential%rcon = 4.0_dp
potential%scon = 8.0_dp
potential%scon = 12.0_dp
potential%confinement = .TRUE.
CALL section_vals_val_get(potential_section,"CONFINEMENT",r_vals=convals)
IF ( SIZE (convals) >= 1 ) THEN

16
src/input_cp2k_atom.F
View file

@ -293,6 +293,22 @@ CONTAINS
CALL section_add_subsection(section,print_key)
CALL section_release(print_key)
CALL cp_print_key_section_create(print_key,"ANALYZE_BASIS",&
description="Calculates some basis set analysis data",&
print_level=high_print_level,filename="__STD_OUT__")
CALL keyword_create(keyword, name="OVERLAP_CONDITION_NUMBER",&
description="Condition number of the basis set overlap matrix calculated for a cubic crystal",&
usage="OVERLAP_CONDITION_NUMBER <logical>",type_of_var=logical_t,default_l_val=.FALSE.)
CALL section_add_keyword(print_key,keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, name="COMPLETENESS",&
description="Calculate a completeness estimate for the basis set.",&
usage="COMPLETENESS <logical>",type_of_var=logical_t,default_l_val=.FALSE.)
CALL section_add_keyword(print_key,keyword)
CALL keyword_release(keyword)
CALL section_add_subsection(section,print_key)
CALL section_release(print_key)
CALL cp_print_key_section_create(print_key,"FIT_PSEUDO",&
description="Controls the printing of FIT PSEUDO task",&
print_level=medium_print_level,filename="__STD_OUT__")

24
tests/ATOM/regtest-pseudo/C_basis1.inp Normal file
View file

@ -0,0 +1,24 @@
&GLOBAL
PROJECT C
PROGRAM_NAME ATOM
&END GLOBAL
&ATOM
ELEMENT C
ELECTRON_CONFIGURATION 1s2 2s2 2p2
&METHOD
METHOD_TYPE KOHN-SHAM
&XC
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&END METHOD
&PRINT
&ANALYZE_BASIS
OVERLAP_CONDITION_NUMBER T
COMPLETENESS T
&END ANALYZE_BASIS
&END
&AE_BASIS
BASIS_TYPE GEOMETRICAL_GTO
&END AE_BASIS
&END ATOM

25
tests/ATOM/regtest-pseudo/C_basis2.inp Normal file
View file

@ -0,0 +1,25 @@
&GLOBAL
PROJECT C
PROGRAM_NAME ATOM
&END GLOBAL
&ATOM
ELEMENT C
ELECTRON_CONFIGURATION 1s2 2s2 2p2
&METHOD
METHOD_TYPE KOHN-SHAM
&XC
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&END METHOD
&PRINT
&ANALYZE_BASIS
OVERLAP_CONDITION_NUMBER T
COMPLETENESS T
&END ANALYZE_BASIS
&END
&AE_BASIS
BASIS_TYPE CONTRACTED_GTO
BASIS_SET 6-31G*
&END AE_BASIS
&END ATOM

35
tests/ATOM/regtest-pseudo/C_basis3.inp Normal file
View file

@ -0,0 +1,35 @@
&GLOBAL
PROJECT C
PROGRAM_NAME ATOM
&END GLOBAL
&ATOM
ELEMENT C
ELECTRON_CONFIGURATION CORE 2s2 2p2
CORE [He]
&METHOD
METHOD_TYPE KOHN-SHAM
&XC
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&END METHOD
&PRINT
&ANALYZE_BASIS
OVERLAP_CONDITION_NUMBER T
COMPLETENESS T
&END ANALYZE_BASIS
&END
&PP_BASIS
BASIS_TYPE GEOMETRICAL_GTO
&END PP_BASIS
&POTENTIAL
PSEUDO_TYPE GTH
&GTH_POTENTIAL
2 4 0 0
0.241474531885 2 -16.691805741058 2.494605958440
1
0.220838245905 1 18.355844903927
&END
&END POTENTIAL
&END ATOM

3
tests/ATOM/regtest-pseudo/TEST_FILES
View file

@ -11,4 +11,7 @@ upf1.inp 35 1.e-12
upf2.inp 35 1.e-12 -7.123642611762
ecp1.inp 35 1.e-12 -5.362932801634
ecp2.inp 35 1.e-12 -5.362932801634
C_basis1.inp 35 1.e-12 -37.747721651890
C_basis2.inp 35 1.e-12 -37.728835366125
C_basis3.inp 35 1.e-12 -13.923782735900
#EOF