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

Better leverage mathconstants module (#3359)

Browse source
This commit is contained in:
Hans Pabst 2024-04-22 14:18:18 +02:00 • committed by GitHub
parent b195825777
commit 1f285aab9e
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
10 changed files with 38 additions and 51 deletions
Download patch file Download diff file Expand all files Collapse all files

7
src/aobasis/ai_onecenter.F
View file

@ -32,7 +32,8 @@ MODULE ai_onecenter
fac,&
gamma0,&
gamma1,&
pi
pi,&
rootpi
#include "../base/base_uses.f90"
IMPLICIT NONE
@ -74,7 +75,7 @@ CONTAINS
CPASSERT(.NOT. (n > SIZE(smat, 1) .OR. m > SIZE(smat, 2)))
spi = SQRT(pi)/2.0_dp**(l + 2)*dfac(2*l + 1)
spi = rootpi/2.0_dp**(l + 2)*dfac(2*l + 1)
el = REAL(l, dp) + 1.5_dp
DO iq = 1, m
@ -111,7 +112,7 @@ CONTAINS
CPASSERT(.NOT. (n > SIZE(kmat, 1) .OR. m > SIZE(kmat, 2)))
spi = dfac(2*l + 3)*SQRT(pi)/2.0_dp**(l + 2)
spi = dfac(2*l + 3)*rootpi/2.0_dp**(l + 2)
DO iq = 1, m
DO ip = 1, n
kmat(ip, iq) = spi*pa(ip)*pb(iq)/(pa(ip) + pb(iq))**(l + 2.5_dp)

5
src/atom_energy.F
View file

@ -49,7 +49,8 @@ MODULE atom_energy
dp
USE mathconstants, ONLY: dfac,&
gamma1,&
pi
pi,&
rootpi
USE periodic_table, ONLY: nelem,&
ptable
USE physcon, ONLY: bohr
@ -583,7 +584,7 @@ CONTAINS
WRITE (iw, '(/,A,I0,A,I0,A)') " Exponent Coef.(Quickstep Normalization), first ", &
n - nder, " valence ", nder, " response"
expzet = 0.25_dp*REAL(2*l + 3, dp)
prefac = SQRT(SQRT(pi)/2._dp**(l + 2)*dfac(2*l + 1))
prefac = SQRT(rootpi/2._dp**(l + 2)*dfac(2*l + 1))
DO i = 1, m
zeta = (2._dp*atom%basis%am(i, l))**expzet
WRITE (iw, '(4X,F20.10,4X,15ES20.6)') atom%basis%am(i, l), ((prefac*rbasis(i, k, l)/zeta), k=1, n)

8
src/atom_grb.F
View file

@ -39,7 +39,7 @@ MODULE atom_grb
dp
USE lapack, ONLY: lapack_ssygv
USE mathconstants, ONLY: dfac,&
pi
rootpi
USE orbital_pointers, ONLY: deallocate_orbital_pointers,&
init_orbital_pointers
USE orbital_transformation_matrices, ONLY: deallocate_spherical_harmonics,&
@ -322,7 +322,7 @@ CONTAINS
! Quickstep normalization
expzet = 0.25_dp*REAL(2*l + 3, dp)
prefac = SQRT(SQRT(pi)/2._dp**(l + 2)*dfac(2*l + 1))
prefac = SQRT(rootpi/2._dp**(l + 2)*dfac(2*l + 1))
DO i = 1, m
zeta = (2._dp*atom%basis%am(i, l))**expzet
qbasis(i, 1:n, l) = rbasis(i, 1:n, l)*prefac/zeta
@ -630,7 +630,7 @@ CONTAINS
! Quickstep normalization polarization basis
DO l = 0, 7
expzet = 0.25_dp*REAL(2*l + 3, dp)
prefac = SQRT(SQRT(pi)/2._dp**(l + 2)*dfac(2*l + 1))
prefac = SQRT(rootpi/2._dp**(l + 2)*dfac(2*l + 1))
DO i = 1, num_pol
zeta = (2._dp*alp(i))**expzet
pbasis(i, 1:num_pol, l) = pbasis(i, 1:num_pol, l)*prefac/zeta
@ -639,7 +639,7 @@ CONTAINS
! Quickstep normalization extended basis
DO l = 0, lmat
expzet = 0.25_dp*REAL(2*l + 3, dp)
prefac = SQRT(SQRT(pi)/2._dp**(l + 2)*dfac(2*l + 1))
prefac = SQRT(rootpi/2._dp**(l + 2)*dfac(2*l + 1))
DO i = 1, next_prim(l)
zeta = (2._dp*ale(i))**expzet
ebasis(i, 1:next_bas(l), l) = ebasis(i, 1:next_bas(l), l)*prefac/zeta

5
src/atom_output.F
View file

@ -30,7 +30,8 @@ MODULE atom_output
USE kinds, ONLY: default_string_length,&
dp
USE mathconstants, ONLY: dfac,&
pi
pi,&
rootpi
USE periodic_table, ONLY: ptable
USE physcon, ONLY: evolt
USE xc_derivatives, ONLY: xc_functional_get_info
@ -472,7 +473,7 @@ CONTAINS
WRITE (iw, '(T3,A,I3)') "L Quantum Number:", l
! Quickstep normalization
expzet = 0.25_dp*REAL(2*l + 3, dp)
prefac = SQRT(SQRT(pi)/2._dp**(l + 2)*dfac(2*l + 1))
prefac = SQRT(rootpi/2._dp**(l + 2)*dfac(2*l + 1))
DO i = 1, atom_basis%nbas(l)
zeta = (2._dp*atom_basis%am(i, l))**expzet
WRITE (iw, '(T5,F14.8,2x,6F12.8)') atom_basis%am(i, l), wfn(i, 1:im, l)*prefac/zeta

2
src/atom_types.F
View file

@ -2374,7 +2374,7 @@ CONTAINS
! Normalization
DO j = 0, lmat
expzet = 0.25_dp*REAL(2*j + 3, dp)
prefac = SQRT(SQRT(pi)/2._dp**(j + 2)*dfac(2*j + 1))
prefac = SQRT(rootpi/2._dp**(j + 2)*dfac(2*j + 1))
DO ipgf = 1, basis%nprim(j)
DO ii = 1, basis%nbas(j)
gcca = basis%cm(ipgf, ii, j)

35
src/common/mathlib.F
View file

@ -18,7 +18,7 @@ MODULE mathlib
dp
USE mathconstants, ONLY: euler,&
fac,&
pi
oorootpi
#include "../base/base_uses.f90"
IMPLICIT NONE
@ -172,7 +172,6 @@ CONTAINS
!> \version 1.0
! **************************************************************************************************
PURE FUNCTION angle(a, b) RESULT(angle_ab)
REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: a, b
REAL(KIND=dp) :: angle_ab
@ -203,7 +202,6 @@ CONTAINS
!> \version 1.0
! **************************************************************************************************
ELEMENTAL FUNCTION binomial(n, k) RESULT(n_over_k)
INTEGER, INTENT(IN) :: n, k
REAL(KIND=dp) :: n_over_k
@ -228,7 +226,6 @@ CONTAINS
!> \version 1.0
! **************************************************************************************************
ELEMENTAL FUNCTION binomial_gen(z, k) RESULT(z_over_k)
REAL(KIND=dp), INTENT(IN) :: z
INTEGER, INTENT(IN) :: k
REAL(KIND=dp) :: z_over_k
@ -254,7 +251,6 @@ CONTAINS
!> \author Ole Schuett
! **************************************************************************************************
PURE FUNCTION multinomial(n, k) RESULT(res)
INTEGER, INTENT(IN) :: n
INTEGER, DIMENSION(:), INTENT(IN) :: k
REAL(KIND=dp) :: res
@ -517,7 +513,6 @@ CONTAINS
!> \version 1.0
! **************************************************************************************************
PURE FUNCTION inv_3x3(a) RESULT(a_inv)
REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: a
REAL(KIND=dp), DIMENSION(3, 3) :: a_inv
@ -951,7 +946,6 @@ CONTAINS
!> \author Dorothea Golze [02.2015]
! **************************************************************************************************
SUBROUTINE get_pseudo_inverse_svd(a, a_pinverse, rskip, determinant, sval)
REAL(KIND=dp), DIMENSION(:, :) :: a, a_pinverse
REAL(KIND=dp), INTENT(IN) :: rskip
REAL(KIND=dp), INTENT(OUT), OPTIONAL :: determinant
@ -1033,7 +1027,6 @@ CONTAINS
!> \author Dorothea Golze [02.2015]
! **************************************************************************************************
SUBROUTINE get_pseudo_inverse_diag(a, a_pinverse, rskip)
REAL(KIND=dp), DIMENSION(:, :) :: a, a_pinverse
REAL(KIND=dp), INTENT(IN) :: rskip
@ -1101,7 +1094,6 @@ CONTAINS
!> \version 1.0
! **************************************************************************************************
PURE FUNCTION reflect_vector(a, b) RESULT(a_mirror)
REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: a, b
REAL(KIND=dp), DIMENSION(3) :: a_mirror
@ -1469,7 +1461,7 @@ CONTAINS
END DO
CPABORT("series failed in expint")
END IF
RETURN
END FUNCTION expint
! **************************************************************************************************
@ -1512,7 +1504,6 @@ CONTAINS
!> - Creation (20.11.98, Matthias Krack)
! **************************************************************************************************
SUBROUTINE diag(n, a, d, v)
INTEGER, INTENT(IN) :: n
REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: a
REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: d
@ -1736,8 +1727,11 @@ CONTAINS
REAL(dp), INTENT(in) :: r, eps, omega
REAL(dp), INTENT(out) :: fn, df
fn = erfc(omega*r) - r*eps
df = -omega*2*EXP(-omega**2*r**2)/SQRT(pi) - eps
REAL(dp) :: qr
qr = omega*r
fn = ERFC(qr) - r*eps
df = -2.0_dp*oorootpi*omega*EXP(-qr**2) - eps
END SUBROUTINE eval_transc_func
END SUBROUTINE erfc_cutoff
@ -1749,8 +1743,7 @@ CONTAINS
!> \author A. Bussy
! **************************************************************************************************
SUBROUTINE complex_diag(matrix, eigenvectors, eigenvalues)
COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: matrix, eigenvectors
COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: matrix, eigenvectors
REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
COMPLEX(KIND=dp), DIMENSION(:), ALLOCATABLE :: work
@ -1775,20 +1768,14 @@ CONTAINS
liwork = iwork(1)
DEALLOCATE (iwork, rwork, work)
ALLOCATE (iwork(liwork), rwork(lrwork), work(lwork))
ALLOCATE (iwork(liwork))
iwork(:) = 0
ALLOCATE (rwork(lrwork))
rwork(:) = 0.0_dp
ALLOCATE (work(lwork))
work(:) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
!diagonalization proper
! diagonalization proper
CALL ZHEEVD('V', 'U', n, eigenvectors, n, eigenvalues, work, lwork, rwork, lrwork, iwork, liwork, info)
DEALLOCATE (iwork, rwork, work)
IF (info /= 0) &
CPABORT("Diagonalisation of a complex matrix failed")
DEALLOCATE (iwork, rwork, work)
END SUBROUTINE complex_diag

10
src/fm/cp_cfm_diag.F
View file

@ -102,13 +102,7 @@ CONTAINS
#endif
DEALLOCATE (iwork, rwork, work)
ALLOCATE (iwork(liwork))
iwork(:) = 0
ALLOCATE (rwork(lrwork))
rwork(:) = 0.0_dp
ALLOCATE (work(lwork))
work(:) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
ALLOCATE (iwork(liwork), rwork(lrwork), work(lwork))
#if defined(__SCALAPACK)
! Scalapack takes advantage of IEEE754 exceptions for speedup.
@ -130,9 +124,9 @@ CONTAINS
eigenvectors%local_data = matrix%local_data
#endif
DEALLOCATE (iwork, rwork, work)
IF (info /= 0) &
CPABORT("Diagonalisation of a complex matrix failed")
DEALLOCATE (iwork, rwork, work)
CALL timestop(handle)

7
src/lri_environment_init.F
View file

@ -33,7 +33,8 @@ MODULE lri_environment_init
lri_env_create,&
lri_environment_type
USE mathconstants, ONLY: fac,&
pi
pi,&
rootpi
USE mathlib, ONLY: invert_matrix
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
@ -474,7 +475,7 @@ CONTAINS
DO ishell = 1, basis%nshell(iset)
l = basis%l(ishell, iset)
expa = 0.5_dp*REAL(2*l + 3, dp)
ppl = fac(2*l + 2)*SQRT(pi)/2._dp**REAL(2*l + 3, dp)/fac(l + 1)
ppl = fac(2*l + 2)*rootpi/2._dp**REAL(2*l + 3, dp)/fac(l + 1)
DO isgf = basis%first_sgf(ishell, iset), basis%last_sgf(ishell, iset)
DO ipgf = 1, basis%npgf(iset)
cci = basis%gcc(ipgf, ishell, iset)
@ -561,7 +562,7 @@ CONTAINS
IF (li == lj) THEN
l = li
expa = 0.5_dp*REAL(2*l + 3, dp)
ppl = fac(2*l + 2)*SQRT(pi)/2._dp**REAL(2*l + 3, dp)/fac(l + 1)
ppl = fac(2*l + 2)*rootpi/2._dp**REAL(2*l + 3, dp)/fac(l + 1)
DO isgf = basis%first_sgf(ishell, iset), basis%last_sgf(ishell, iset)
m_i = basis%m(isgf)
DO jsgf = basis%first_sgf(jshell, jset), basis%last_sgf(jshell, jset)

5
src/qs_dispersion_nonloc.F
View file

@ -29,7 +29,8 @@ MODULE qs_dispersion_nonloc
xc_vdw_fun_nonloc
USE kinds, ONLY: default_string_length,&
dp
USE mathconstants, ONLY: pi
USE mathconstants, ONLY: pi,&
rootpi
USE message_passing, ONLY: mp_para_env_type
USE pw_grid_types, ONLY: HALFSPACE,&
pw_grid_type
@ -211,7 +212,7 @@ CONTAINS
beta = 0.03125_dp*(3.0_dp/(b_value**2.0_dp))**0.75_dp
nspin = SIZE(rho_r)
const = 1.0_dp/(3.0_dp*SQRT(pi)*b_value**1.5_dp)/(pi**0.75_dp)
const = 1.0_dp/(3.0_dp*rootpi*b_value**1.5_dp)/(pi**0.75_dp)
! temporary arrays for FFT
CALL pw_pool%create_pw(tmp_g)

5
src/xc/xc_xpbe_hole_t_c_lr.F
View file

@ -19,7 +19,8 @@ MODULE xc_xpbe_hole_t_c_lr
section_vals_val_get
USE kinds, ONLY: dp
USE mathconstants, ONLY: euler,&
pi
pi,&
rootpi
USE xc_derivative_desc, ONLY: deriv_norm_drho,&
deriv_norm_drhoa,&
deriv_norm_drhob,&
@ -537,7 +538,7 @@ CONTAINS
t61 = D**2
t64 = D*t21
t65 = t34*B
t69 = SQRT(pi)
t69 = rootpi
t71 = F1*t21
t73 = t71*t34 + F2
t77 = C*(1 + t73*t12*t13)