From 1f285aab9e80d5f88cf770a48ae4166d0bfccb77 Mon Sep 17 00:00:00 2001 From: Hans Pabst Date: Mon, 22 Apr 2024 14:18:18 +0200 Subject: [PATCH] Better leverage mathconstants module (#3359) --- src/aobasis/ai_onecenter.F | 7 ++++--- src/atom_energy.F | 5 +++-- src/atom_grb.F | 8 ++++---- src/atom_output.F | 5 +++-- src/atom_types.F | 2 +- src/common/mathlib.F | 35 +++++++++++------------------------ src/fm/cp_cfm_diag.F | 10 ++-------- src/lri_environment_init.F | 7 ++++--- src/qs_dispersion_nonloc.F | 5 +++-- src/xc/xc_xpbe_hole_t_c_lr.F | 5 +++-- 10 files changed, 38 insertions(+), 51 deletions(-) diff --git a/src/aobasis/ai_onecenter.F b/src/aobasis/ai_onecenter.F index ae27c96797..42e2188194 100644 --- a/src/aobasis/ai_onecenter.F +++ b/src/aobasis/ai_onecenter.F @@ -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) diff --git a/src/atom_energy.F b/src/atom_energy.F index b977bedef1..bc290c7bf7 100644 --- a/src/atom_energy.F +++ b/src/atom_energy.F @@ -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) diff --git a/src/atom_grb.F b/src/atom_grb.F index 551f19004d..ebc89ddb6e 100644 --- a/src/atom_grb.F +++ b/src/atom_grb.F @@ -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 diff --git a/src/atom_output.F b/src/atom_output.F index 595d2f990c..14a304726c 100644 --- a/src/atom_output.F +++ b/src/atom_output.F @@ -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 diff --git a/src/atom_types.F b/src/atom_types.F index 5e782c9848..d1a6e1effb 100644 --- a/src/atom_types.F +++ b/src/atom_types.F @@ -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) diff --git a/src/common/mathlib.F b/src/common/mathlib.F index 7c2505e794..8f263081e4 100644 --- a/src/common/mathlib.F +++ b/src/common/mathlib.F @@ -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 diff --git a/src/fm/cp_cfm_diag.F b/src/fm/cp_cfm_diag.F index 531cb828f0..ac430e2645 100644 --- a/src/fm/cp_cfm_diag.F +++ b/src/fm/cp_cfm_diag.F @@ -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) diff --git a/src/lri_environment_init.F b/src/lri_environment_init.F index 993fe3f777..c1e7a8a0c8 100644 --- a/src/lri_environment_init.F +++ b/src/lri_environment_init.F @@ -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) diff --git a/src/qs_dispersion_nonloc.F b/src/qs_dispersion_nonloc.F index 2abc54a4b2..6f4cd68963 100644 --- a/src/qs_dispersion_nonloc.F +++ b/src/qs_dispersion_nonloc.F @@ -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) diff --git a/src/xc/xc_xpbe_hole_t_c_lr.F b/src/xc/xc_xpbe_hole_t_c_lr.F index 5f89e2c508..03e7c1788d 100644 --- a/src/xc/xc_xpbe_hole_t_c_lr.F +++ b/src/xc/xc_xpbe_hole_t_c_lr.F @@ -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)