mirror of
https://github.com/cp2k/cp2k.git
synced 2026-07-26 13:15:21 -04:00
Detail abort messages and add warning for bad electron count
This commit is contained in:
Matthias Krack
parent
4a2d255c13
commit
e4d6f4841e
5 changed files with 41 additions and 34 deletions
|
|
@ -98,9 +98,9 @@ CONTAINS
|
|||
CASE (do_uks_atom, do_uhf_atom)
|
||||
CALL calculate_atom_unrestricted(atom, iw, noguess, converged)
|
||||
CASE (do_rohf_atom)
|
||||
CPABORT("")
|
||||
CPABORT("The ROHF method is not yet implemented in the ATOM code.")
|
||||
CASE DEFAULT
|
||||
CPABORT("")
|
||||
CPABORT("Invalid method specified for ATOM code. Check the code!")
|
||||
END SELECT
|
||||
|
||||
CALL timestop(handle)
|
||||
|
|
|
|||
|
|
@ -51,8 +51,8 @@ MODULE atom_energy
|
|||
USE kinds, ONLY: default_string_length,&
|
||||
dp
|
||||
USE mathconstants, ONLY: dfac,&
|
||||
fourpi,&
|
||||
gamma1,&
|
||||
pi,&
|
||||
rootpi
|
||||
USE periodic_table, ONLY: nelem,&
|
||||
ptable
|
||||
|
|
@ -657,7 +657,8 @@ CONTAINS
|
|||
INTEGER :: i, ibeta, ii, j, k, kk, l, lmax, n0, &
|
||||
nbeta, nr, nwfc, nwfn
|
||||
LOGICAL :: soc, up
|
||||
REAL(KIND=dp) :: occa, occb, pf, rl, rlp, rmax, vnlm, vnlp
|
||||
REAL(KIND=dp) :: occa, occb, pf, rhoatom, rl, rlp, rmax, &
|
||||
vnlm, vnlp
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: beta, corden, dens, ef, locpot, rp
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dij
|
||||
TYPE(atom_gthpot_type), POINTER :: pot
|
||||
|
|
@ -963,11 +964,18 @@ CONTAINS
|
|||
CALL atom_density(dens, atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, &
|
||||
"RHO", atom%basis%grid%rad)
|
||||
IF (up) THEN
|
||||
WRITE (iw, '(T8,4ES25.12E3)') (4._dp*pi*dens(j)*atom%basis%grid%rad2(j), j=1, nr)
|
||||
WRITE (iw, '(T8,4ES25.12E3)') (fourpi*dens(j)*atom%basis%grid%rad2(j), j=1, nr)
|
||||
ELSE
|
||||
WRITE (iw, '(T8,4ES25.12E3)') (4._dp*pi*dens(j)*atom%basis%grid%rad2(j), j=nr, 1, -1)
|
||||
WRITE (iw, '(T8,4ES25.12E3)') (fourpi*dens(j)*atom%basis%grid%rad2(j), j=nr, 1, -1)
|
||||
END IF
|
||||
WRITE (iw, '(T4,A)') '</PP_RHOATOM>'
|
||||
|
||||
! Check atomic electron count
|
||||
rhoatom = SUM(fourpi*atom%basis%grid%wr(:)*dens(:))
|
||||
IF (ABS(pot%zion - rhoatom) > 0.01_dp) THEN
|
||||
CALL cp_warn(__LOCATION__, "Valence charge and electron count differ by more than 0.01")
|
||||
END IF
|
||||
|
||||
DEALLOCATE (dens)
|
||||
|
||||
! PP SOC information
|
||||
|
|
|
|||
|
|
@ -380,6 +380,7 @@ CONTAINS
|
|||
INTEGER, INTENT(IN) :: zval
|
||||
CHARACTER(LEN=2) :: btyp
|
||||
|
||||
CHARACTER(LEN=*), PARAMETER :: routineN = 'init_atom_basis'
|
||||
INTEGER, PARAMETER :: nua = 40, nup = 16
|
||||
REAL(KIND=dp), DIMENSION(nua), PARAMETER :: ugbs = (/0.007299_dp, 0.013705_dp, 0.025733_dp, &
|
||||
0.048316_dp, 0.090718_dp, 0.170333_dp, 0.319819_dp, 0.600496_dp, 1.127497_dp, 2.117000_dp,&
|
||||
|
|
@ -392,14 +393,16 @@ CONTAINS
|
|||
341890134.751331_dp/)
|
||||
|
||||
CHARACTER(LEN=default_string_length) :: basis_fn, basis_name
|
||||
INTEGER :: basistype, i, j, k, l, ll, m, ngp, nl, &
|
||||
nr, nu, quadtype
|
||||
INTEGER :: basistype, handle, i, j, k, l, ll, m, &
|
||||
ngp, nl, nr, nu, quadtype
|
||||
INTEGER, DIMENSION(0:lmat) :: starti
|
||||
INTEGER, DIMENSION(:), POINTER :: nqm, num_gto, num_slater, sindex
|
||||
REAL(KIND=dp) :: al, amax, aval, cval, ear, pf, rk
|
||||
REAL(KIND=dp), DIMENSION(:), POINTER :: expo
|
||||
TYPE(section_vals_type), POINTER :: gto_basis_section
|
||||
|
||||
CALL timeset(routineN, handle)
|
||||
|
||||
! btyp = AE : standard all-electron basis
|
||||
! btyp = PP : standard pseudopotential basis
|
||||
! btyp = AA : high accuracy all-electron basis
|
||||
|
|
@ -700,6 +703,8 @@ CONTAINS
|
|||
CPABORT("")
|
||||
END SELECT
|
||||
|
||||
CALL timestop(handle)
|
||||
|
||||
END SUBROUTINE init_atom_basis
|
||||
|
||||
! **************************************************************************************************
|
||||
|
|
@ -709,6 +714,7 @@ CONTAINS
|
|||
SUBROUTINE init_atom_basis_default_pp(basis)
|
||||
TYPE(atom_basis_type), INTENT(INOUT) :: basis
|
||||
|
||||
CHARACTER(LEN=*), PARAMETER :: routineN = 'init_atom_basis_default_pp'
|
||||
INTEGER, PARAMETER :: nua = 40, nup = 20
|
||||
REAL(KIND=dp), DIMENSION(nua), PARAMETER :: ugbs = (/0.007299_dp, 0.013705_dp, 0.025733_dp, &
|
||||
0.048316_dp, 0.090718_dp, 0.170333_dp, 0.319819_dp, 0.600496_dp, 1.127497_dp, 2.117000_dp,&
|
||||
|
|
@ -720,11 +726,14 @@ CONTAINS
|
|||
27508345.793637_dp, 51649961.080194_dp, 96978513.342764_dp, 182087882.613702_dp, &
|
||||
341890134.751331_dp/)
|
||||
|
||||
INTEGER :: i, k, l, m, ngp, nr, nu, quadtype
|
||||
INTEGER :: handle, i, k, l, m, ngp, nr, nu, quadtype
|
||||
REAL(KIND=dp) :: al, ear, rk
|
||||
|
||||
CALL timeset(routineN, handle)
|
||||
|
||||
NULLIFY (basis%am, basis%cm, basis%as, basis%ns, basis%bf, basis%dbf, basis%ddbf)
|
||||
! allocate and initialize the atomic grid
|
||||
|
||||
! Allocate and initialize the atomic grid
|
||||
NULLIFY (basis%grid)
|
||||
CALL allocate_grid_atom(basis%grid)
|
||||
quadtype = do_gapw_log
|
||||
|
|
@ -768,6 +777,8 @@ CONTAINS
|
|||
END DO
|
||||
END DO
|
||||
|
||||
CALL timestop(handle)
|
||||
|
||||
END SUBROUTINE init_atom_basis_default_pp
|
||||
|
||||
! **************************************************************************************************
|
||||
|
|
|
|||
|
|
@ -402,7 +402,7 @@ CONTAINS
|
|||
+ 2._dp*ff*basis%dbf(:, i, l)*basis%bf(:, j, l)/rr(:) &
|
||||
+ 2._dp*ff*basis%bf(:, i, l)*basis%dbf(:, j, l)/rr(:)
|
||||
ELSE
|
||||
CPABORT("")
|
||||
CPABORT("Unknown matrix type specified. Check the code!")
|
||||
END IF
|
||||
END DO
|
||||
END DO
|
||||
|
|
@ -1148,8 +1148,6 @@ CONTAINS
|
|||
END DO
|
||||
expa(:) = EXP(-z(:)**2)*ff/(a + b)**(l + 1)
|
||||
SELECT CASE (l)
|
||||
CASE DEFAULT
|
||||
CPABORT("")
|
||||
CASE (0)
|
||||
cpot(:) = cpot(:) + 0.25_dp*erfa(:)
|
||||
CASE (1)
|
||||
|
|
@ -1158,6 +1156,8 @@ CONTAINS
|
|||
cpot(:) = cpot(:) + 0.9375_dp*erfa(:) - expa(:)*(0.875_dp + 0.25_dp*z(:)**2)
|
||||
CASE (3)
|
||||
cpot(:) = cpot(:) + 3.28125_dp*erfa(:) - expa(:)*(3.5625_dp + 1.375_dp*z(:)**2 + 0.25*z(:)**4)
|
||||
CASE DEFAULT
|
||||
CPABORT("Invalid l number for GTO specified. Check the code!")
|
||||
END SELECT
|
||||
END DO
|
||||
END DO
|
||||
|
|
@ -1183,8 +1183,6 @@ CONTAINS
|
|||
END DO
|
||||
expa(:) = EXP(-z(:)**2)*ff/(a + b)**(l + 1)
|
||||
SELECT CASE (l)
|
||||
CASE DEFAULT
|
||||
CPABORT("")
|
||||
CASE (0)
|
||||
cpot(:) = cpot(:) + 0.25_dp*erfa(:)
|
||||
CASE (1)
|
||||
|
|
@ -1193,6 +1191,8 @@ CONTAINS
|
|||
cpot(:) = cpot(:) + 0.9375_dp*erfa(:) - expa(:)*(0.875_dp + 0.25_dp*z(:)**2)
|
||||
CASE (3)
|
||||
cpot(:) = cpot(:) + 3.28125_dp*erfa(:) - expa(:)*(3.5625_dp + 1.375_dp*z(:)**2 + 0.25*z(:)**4)
|
||||
CASE DEFAULT
|
||||
CPABORT("Invalid l number for CGTO specified. Check the code!")
|
||||
END SELECT
|
||||
END DO
|
||||
END DO
|
||||
|
|
|
|||
|
|
@ -602,15 +602,12 @@ CONTAINS
|
|||
INTEGER :: i
|
||||
REAL(dp) :: cost, f, sint, sint2, t, w, x
|
||||
|
||||
f = pi/REAL(n + 1, dp)
|
||||
f = pi/REAL(n + 1, KIND=dp)
|
||||
|
||||
SELECT CASE (radial_quadrature)
|
||||
|
||||
CASE (do_gapw_gcs)
|
||||
|
||||
! *** Gauss-Chebyshev quadrature formula of the second kind ***
|
||||
! *** u [-1,+1] -> r [0,infinity] => r = (1 + u)/(1 - u) ***
|
||||
|
||||
! Gauss-Chebyshev quadrature formula of the second kind
|
||||
! u [-1,+1] -> r [0,infinity] => r = (1 + u)/(1 - u)
|
||||
DO i = 1, n
|
||||
t = REAL(i, dp)*f
|
||||
x = COS(t)
|
||||
|
|
@ -620,12 +617,9 @@ CONTAINS
|
|||
wr(i) = w/SQRT(1.0_dp - x**2)
|
||||
wr(i) = 2.0_dp*wr(i)*r2(i)/(1.0_dp - x)**2
|
||||
END DO
|
||||
|
||||
CASE (do_gapw_gct)
|
||||
|
||||
! *** Transformed Gauss-Chebyshev quadrature formula of the second kind ***
|
||||
! *** u [-1,+1] -> r [0,infinity] => r = (1 + u)/(1 - u) ***
|
||||
|
||||
! Transformed Gauss-Chebyshev quadrature formula of the second kind
|
||||
! u [-1,+1] -> r [0,infinity] => r = (1 + u)/(1 - u)
|
||||
DO i = 1, n
|
||||
t = REAL(i, dp)*f
|
||||
cost = COS(t)
|
||||
|
|
@ -638,12 +632,9 @@ CONTAINS
|
|||
r2(n + 1 - i) = r(n + 1 - i)**2
|
||||
wr(n + 1 - i) = 2.0_dp*w*r2(n + 1 - i)/(1.0_dp - x)**2
|
||||
END DO
|
||||
|
||||
CASE (do_gapw_log)
|
||||
|
||||
! *** Transformed Gauss-Chebyshev quadrature formula of the second kind ***
|
||||
! *** u [-1,+1] -> r [0,infinity] => r = ln(2/(1 - u))/ln(2) ***
|
||||
|
||||
! Logarithmic transformed Gauss-Chebyshev quadrature formula of the second kind
|
||||
! u [-1,+1] -> r [0,infinity] => r = ln(2/(1 - u))/ln(2)
|
||||
DO i = 1, n
|
||||
t = REAL(i, dp)*f
|
||||
cost = COS(t)
|
||||
|
|
@ -656,11 +647,8 @@ CONTAINS
|
|||
r2(n + 1 - i) = r(n + 1 - i)**2
|
||||
wr(n + 1 - i) = w*r2(n + 1 - i)/(LOG(2.0_dp)*(1.0_dp - x))
|
||||
END DO
|
||||
|
||||
CASE DEFAULT
|
||||
|
||||
CPABORT("Invalid radial quadrature type specified")
|
||||
|
||||
END SELECT
|
||||
|
||||
END SUBROUTINE radial_grid
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue