mirror of
https://github.com/cp2k/cp2k.git
synced 2026-07-21 06:25:15 -04:00
Remove unused beta routines and move digamma to mathlib (#5568)
This commit is contained in:
SY Wang
GitHub
parent
d4bfb39614
commit
7de47fd987
5 changed files with 132 additions and 2086 deletions
|
|
@ -70,7 +70,6 @@ list(
|
|||
soc_pseudopotential_methods.F
|
||||
soc_pseudopotential_utils.F
|
||||
basis_set_output.F
|
||||
beta_gamma_psi.F
|
||||
block_p_types.F
|
||||
bse_main.F
|
||||
bse_full_diag.F
|
||||
|
|
|
|||
2081
src/beta_gamma_psi.F
2081
src/beta_gamma_psi.F
File diff suppressed because it is too large
Load diff
|
|
@ -55,6 +55,7 @@ MODULE mathlib
|
|||
multinomial, &
|
||||
det_3x3, &
|
||||
dihedral_angle, &
|
||||
digamma, &
|
||||
gcd, &
|
||||
inv_3x3, &
|
||||
lcm, &
|
||||
|
|
@ -1385,6 +1386,134 @@ CONTAINS
|
|||
|
||||
END FUNCTION ei
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief Computes the digamma function, the logarithmic derivative of the gamma function.
|
||||
!> \param xx Argument of the digamma function.
|
||||
!> \return The digamma function at xx. Returns zero at poles and for arguments outside the
|
||||
!> supported numerical range.
|
||||
!> \par History
|
||||
!> Adapted from the FUNPACK routine PSI, written at Argonne National Laboratory and modified
|
||||
!> by A. H. Morris (NSWC).
|
||||
!> \note The implementation uses the rational Chebyshev approximations of Cody, Strecok, and
|
||||
!> Thacher, Math. Comp. 27, 123-127 (1973), together with the reflection formula.
|
||||
! **************************************************************************************************
|
||||
PURE FUNCTION digamma(xx) RESULT(fn_val)
|
||||
REAL(dp), INTENT(IN) :: xx
|
||||
REAL(dp) :: fn_val
|
||||
|
||||
REAL(dp), PARAMETER :: dx0 = 1.461632144968362341262659542325721325e0_dp, p1(7) = [ &
|
||||
.895385022981970e-02_dp, .477762828042627e+01_dp, .142441585084029e+03_dp, &
|
||||
.118645200713425e+04_dp, .363351846806499e+04_dp, .413810161269013e+04_dp, &
|
||||
.130560269827897e+04_dp], p2(4) = [-.212940445131011e+01_dp, -.701677227766759e+01_dp, &
|
||||
-.448616543918019e+01_dp, -.648157123766197e+00_dp], piov4 = .785398163397448e0_dp, q1(6) &
|
||||
= [.448452573429826e+02_dp, .520752771467162e+03_dp, .221000799247830e+04_dp, &
|
||||
.364127349079381e+04_dp, .190831076596300e+04_dp, .691091682714533e-05_dp]
|
||||
REAL(dp), PARAMETER :: q2(4) = [.322703493791143e+02_dp, .892920700481861e+02_dp, &
|
||||
.546117738103215e+02_dp, .777788548522962e+01_dp]
|
||||
|
||||
INTEGER :: i, m, n, nq
|
||||
REAL(dp) :: aug, den, sgn, upper, w, x, xmax1, xmx0, &
|
||||
xsmall, z
|
||||
|
||||
! piov4 is pi/4 and dx0 is the positive zero of the digamma function. The coefficient sets
|
||||
! p1/q1 approximate digamma(x)/(x - dx0) for 0.5 <= x <= 3, while p2/q2 approximate
|
||||
! digamma(x) - log(x) + 1/(2*x) for x > 3.
|
||||
|
||||
! xmax1 limits integer argument reduction for large negative values and marks the point beyond
|
||||
! which the asymptotic value log(x) is sufficient. For very small x, pi*cot(pi*x) is
|
||||
! approximated by 1/x.
|
||||
xmax1 = MIN(REAL(HUGE(0), KIND=dp), 1.0e0_dp/EPSILON(1.0e0_dp))
|
||||
xsmall = 1.e-9_dp
|
||||
|
||||
x = xx
|
||||
aug = 0.0e0_dp
|
||||
IF (x < 0.5e0_dp) THEN
|
||||
! Use digamma(x) = digamma(1 - x) - pi*cot(pi*x).
|
||||
IF (ABS(x) <= xsmall) THEN
|
||||
IF (x == 0.0e0_dp) THEN
|
||||
! The digamma function has a pole at zero.
|
||||
fn_val = 0.0e0_dp
|
||||
RETURN
|
||||
END IF
|
||||
! For small nonzero x, use -1/x for -pi*cot(pi*x).
|
||||
aug = -1.0e0_dp/x
|
||||
x = 1.0e0_dp - x
|
||||
ELSE
|
||||
! Reduce the cotangent argument to the first quadrant and keep track of its sign.
|
||||
w = -x
|
||||
sgn = piov4
|
||||
IF (w <= 0.0e0_dp) THEN
|
||||
w = -w
|
||||
sgn = -sgn
|
||||
END IF
|
||||
IF (w >= xmax1) THEN
|
||||
! Argument reduction is no longer reliable.
|
||||
fn_val = 0.0e0_dp
|
||||
RETURN
|
||||
END IF
|
||||
nq = INT(w)
|
||||
w = w - nq
|
||||
nq = INT(w*4.0e0_dp)
|
||||
w = 4.0e0_dp*(w - nq*.25e0_dp)
|
||||
|
||||
n = nq/2
|
||||
IF ((n + n) /= nq) w = 1.0e0_dp - w
|
||||
z = piov4*w
|
||||
m = n/2
|
||||
IF ((m + m) /= n) sgn = -sgn
|
||||
|
||||
! Evaluate -pi*cot(pi*x), using either tan(z) or cot(z) after quadrant reduction.
|
||||
n = (nq + 1)/2
|
||||
m = n/2
|
||||
m = m + m
|
||||
IF (m /= n) THEN
|
||||
aug = sgn*((SIN(z)/COS(z))*4.0e0_dp)
|
||||
ELSE
|
||||
IF (z == 0.0e0_dp) THEN
|
||||
! Negative integers are poles of the digamma function.
|
||||
fn_val = 0.0e0_dp
|
||||
RETURN
|
||||
END IF
|
||||
aug = sgn*((COS(z)/SIN(z))*4.0e0_dp)
|
||||
END IF
|
||||
x = 1.0e0_dp - x
|
||||
END IF
|
||||
END IF
|
||||
|
||||
IF (x <= 3.0e0_dp) THEN
|
||||
! Rational approximation on 0.5 <= x <= 3.
|
||||
den = x
|
||||
upper = p1(1)*x
|
||||
|
||||
DO i = 1, 5
|
||||
den = (den + q1(i))*x
|
||||
upper = (upper + p1(i + 1))*x
|
||||
END DO
|
||||
|
||||
den = (upper + p1(7))/(den + q1(6))
|
||||
xmx0 = x - dx0
|
||||
fn_val = den*xmx0 + aug
|
||||
RETURN
|
||||
END IF
|
||||
|
||||
IF (x < xmax1) THEN
|
||||
! Asymptotic rational correction for 3 < x < xmax1.
|
||||
w = 1.0e0_dp/(x*x)
|
||||
den = w
|
||||
upper = p2(1)*w
|
||||
|
||||
DO i = 1, 3
|
||||
den = (den + q2(i))*w
|
||||
upper = (upper + p2(i + 1))*w
|
||||
END DO
|
||||
|
||||
aug = upper/(den + q2(4)) - 0.5e0_dp/x + aug
|
||||
END IF
|
||||
! For x >= xmax1, the correction is negligible and this reduces to log(x).
|
||||
fn_val = aug + LOG(x)
|
||||
|
||||
END FUNCTION digamma
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief computes the exponential integral
|
||||
!> En(x) = Int(exp(-x*t)/t^n,t=1..infinity) x>0, n=0,1,..
|
||||
|
|
|
|||
|
|
@ -13,7 +13,6 @@
|
|||
! **************************************************************************************************
|
||||
MODULE fp_methods
|
||||
|
||||
USE beta_gamma_psi, ONLY: psi
|
||||
USE cell_types, ONLY: cell_type,&
|
||||
pbc
|
||||
USE cp_log_handling, ONLY: cp_get_default_logger,&
|
||||
|
|
@ -29,6 +28,7 @@ MODULE fp_methods
|
|||
USE mathconstants, ONLY: fac,&
|
||||
maxfac,&
|
||||
oorootpi
|
||||
USE mathlib, ONLY: digamma
|
||||
USE particle_list_types, ONLY: particle_list_type
|
||||
USE particle_types, ONLY: particle_type
|
||||
#include "./base/base_uses.f90"
|
||||
|
|
@ -162,7 +162,7 @@ CONTAINS
|
|||
rab = pbc(rab, cell)
|
||||
r = SQRT(SUM(rab**2))
|
||||
CALL smooth_count(r, fp_env%inner_radius, fp_env%smooth_width, c, dcdr)
|
||||
sf = kT*(psi(fp_env%ro1 + fp_env%ri2 + 1) - psi(fp_env%ri2 + 1))*(-dcdr)/r
|
||||
sf = kT*(digamma(fp_env%ro1 + fp_env%ri2 + 1) - digamma(fp_env%ri2 + 1))*(-dcdr)/r
|
||||
particles(iparticle)%f = particles(iparticle)%f - sf*rab
|
||||
particles(icenter)%f = particles(icenter)%f + sf*rab
|
||||
END DO
|
||||
|
|
@ -173,7 +173,7 @@ CONTAINS
|
|||
rab = pbc(rab, cell)
|
||||
r = SQRT(SUM(rab**2))
|
||||
CALL smooth_count(r, fp_env%outer_radius, fp_env%smooth_width, c, dcdr)
|
||||
sf = kT*(psi(fp_env%ro1 + fp_env%ri2 + 1) - psi(fp_env%ro1 + 1))*(dcdr)/r
|
||||
sf = kT*(digamma(fp_env%ro1 + fp_env%ri2 + 1) - digamma(fp_env%ro1 + 1))*(dcdr)/r
|
||||
particles(iparticle)%f = particles(iparticle)%f - sf*rab
|
||||
particles(icenter)%f = particles(icenter)%f + sf*rab
|
||||
END DO
|
||||
|
|
|
|||
|
|
@ -6,7 +6,6 @@ base_hooks.F: Found STOP statement in procedure "cp__b" https://cp2k.org/conv#c2
|
|||
base_hooks.F: Found WRITE statement with hardcoded unit in "cp_abort" https://cp2k.org/conv#c012
|
||||
base_hooks.F: Found WRITE statement with hardcoded unit in "cp_hint" https://cp2k.org/conv#c012
|
||||
base_hooks.F: Found WRITE statement with hardcoded unit in "cp_warn" https://cp2k.org/conv#c012
|
||||
beta_gamma_psi.F: Found GOTO statement in procedure "bratio" https://cp2k.org/conv#c201
|
||||
bfgs_optimizer.F: Found GOTO statement in procedure "rat_fun_opt" https://cp2k.org/conv#c201
|
||||
cg_test.F: Found WRITE statement with hardcoded unit in "clebsch_gordon_test" https://cp2k.org/conv#c012
|
||||
cp_dbcsr_operations.F: Found CALL cp_fm_gemm in procedure "cp_dbcsr_plus_fm_fm_t" https://cp2k.org/conv#c101
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue