diff --git a/src/atom_electronic_structure.F b/src/atom_electronic_structure.F index 6106e10615..7825339bfa 100644 --- a/src/atom_electronic_structure.F +++ b/src/atom_electronic_structure.F @@ -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) diff --git a/src/atom_energy.F b/src/atom_energy.F index 46accb65ba..2e743a589a 100644 --- a/src/atom_energy.F +++ b/src/atom_energy.F @@ -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)') '' + + ! 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 diff --git a/src/atom_types.F b/src/atom_types.F index d83a43ac0b..d5eaa55b28 100644 --- a/src/atom_types.F +++ b/src/atom_types.F @@ -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 ! ************************************************************************************************** diff --git a/src/atom_utils.F b/src/atom_utils.F index 39b5571023..e28aa42366 100644 --- a/src/atom_utils.F +++ b/src/atom_utils.F @@ -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 diff --git a/src/qs_grid_atom.F b/src/qs_grid_atom.F index b8c45919bc..6cc8619371 100644 --- a/src/qs_grid_atom.F +++ b/src/qs_grid_atom.F @@ -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