LCPQ
/
qp2
mirror of https://github.com/QuantumPackage/qp2.git
9
1
Fork 0
Code Releases Activity
qp2/src/utils/util.irp.f

643 lines
12 KiB
Fortran
Raw Permalink Blame History

double precision function derf_mu_x(mu,x)
implicit none
include 'utils/constants.include.F'
double precision, intent(in) :: mu,x
if(dabs(x).gt.1.d-6)then
derf_mu_x = derf(mu * x)/x
else
derf_mu_x = inv_sq_pi * 2.d0 * mu * (1.d0 - mu*mu*x*x/3.d0)
endif
end
double precision function binom_func(i,j)
implicit none
BEGIN_DOC
!.. math ::
!
! \frac{i!}{j!(i-j)!}
!
END_DOC
integer,intent(in) :: i,j
double precision :: logfact
integer, save :: ifirst
double precision, save :: memo(0:15,0:15)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: memo
integer :: k,l
if (ifirst == 0) then
ifirst = 1
do k=0,15
do l=0,15
memo(k,l) = dexp( logfact(k)-logfact(l)-logfact(k-l) )
enddo
enddo
endif
if ( (i<=15).and.(j<=15) ) then
binom_func = memo(i,j)
else
binom_func = dexp( logfact(i)-logfact(j)-logfact(i-j) )
endif
! To avoid .999999 numbers
binom_func = floor(binom_func + 0.5d0)
end
BEGIN_PROVIDER [ double precision, binom, (0:40,0:40) ]
&BEGIN_PROVIDER [ double precision, binom_transp, (0:40,0:40) ]
implicit none
BEGIN_DOC
! Binomial coefficients
END_DOC
integer :: k,l
double precision :: logfact
do k=0,40
do l=0,40
binom(k,l) = dexp( logfact(k)-logfact(l)-logfact(k-l) )
binom_transp(l,k) = binom(k,l)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer*8, binom_int, (0:40,0:40) ]
&BEGIN_PROVIDER [ integer*8, binom_int_transp, (0:40,0:40) ]
implicit none
BEGIN_DOC
! Binomial coefficients, as integers*8
END_DOC
integer :: k,l
double precision :: logfact
do l=0,40
do k=0,40
binom_int(k,l) = int(binom(k,l)+0.1d0,8)
enddo
enddo
END_PROVIDER
double precision function fact(n)
implicit none
BEGIN_DOC
! n!
END_DOC
integer :: n
double precision, save :: memo(1:100)
integer, save :: memomax = 1
if (n<=memomax) then
if (n<2) then
fact = 1.d0
else
fact = memo(n)
endif
return
endif
integer :: i
memo(1) = 1.d0
do i=memomax+1,min(n,100)
memo(i) = memo(i-1)*dble(i)
enddo
memomax = min(n,100)
double precision :: logfact
fact = dexp(logfact(n))
end function
double precision function logfact(n)
implicit none
BEGIN_DOC
! n!
END_DOC
integer :: n
double precision, save :: memo(1:100)
integer, save :: memomax = 1
if (n<=memomax) then
if (n<2) then
logfact = 0.d0
else
logfact = memo(n)
endif
return
endif
integer :: i
memo(1) = 0.d0
do i=memomax+1,min(n,100)
memo(i) = memo(i-1)+dlog(dble(i))
enddo
memomax = min(n,100)
logfact = memo(memomax)
do i=101,n
logfact += dlog(dble(i))
enddo
end function
! ---
BEGIN_PROVIDER [ double precision, fact_inv, (128) ]
implicit none
BEGIN_DOC
! 1/n!
END_DOC
integer :: i
double precision :: fact
do i=1,size(fact_inv)
fact_inv(i) = 1.d0/fact(i)
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, shiftfact_op5_inv, (128) ]
BEGIN_DOC
!
! 1 / Gamma(n + 0.5)
!
END_DOC
implicit none
integer :: i
double precision :: tmp
do i = 1, size(shiftfact_op5_inv)
!tmp = dgamma(dble(i) + 0.5d0)
tmp = gamma(dble(i) + 0.5d0)
shiftfact_op5_inv(i) = 1.d0 / tmp
enddo
END_PROVIDER
! ---
double precision function dble_fact(n)
implicit none
integer :: n
double precision :: dble_fact_even, dble_fact_odd
dble_fact = 1.d0
if(n.lt.0) return
if(iand(n,1).eq.0)then
dble_fact = dble_fact_even(n)
else
dble_fact= dble_fact_odd(n)
endif
end function
double precision function dble_fact_even(n) result(fact2)
implicit none
BEGIN_DOC
! n!!
END_DOC
integer :: n,k
double precision, save :: memo(0:100)
integer, save :: memomax = 0
double precision :: prod
ASSERT (iand(n,1) /= 1)
! prod=1.d0
! do k=2,n,2
! prod=prod*dfloat(k)
! enddo
! fact2=prod
! return
!
if (n <= memomax) then
if (n < 2) then
fact2 = 1.d0
else
fact2 = memo(n)
endif
return
endif
integer :: i
memo(0)=1.d0
memo(1)=1.d0
do i=memomax+2,min(n,100),2
memo(i) = memo(i-2)* dble(i)
enddo
memomax = min(n,100)
fact2 = memo(memomax)
if (n > 100) then
double precision :: dble_logfact
fact2 = dexp(dble_logfact(n))
endif
end function
double precision function dble_fact_odd(n) result(fact2)
implicit none
BEGIN_DOC
! n!!
END_DOC
integer :: n
double precision, save :: memo(1:100)
integer, save :: memomax = 1
ASSERT (iand(n,1) /= 0)
if (n<=memomax) then
if (n<3) then
fact2 = 1.d0
else
fact2 = memo(n)
endif
return
endif
integer :: i
memo(1) = 1.d0
do i=memomax+2,min(n,99),2
memo(i) = memo(i-2)* dble(i)
enddo
memomax = min(n,99)
fact2 = memo(memomax)
if (n > 99) then
double precision :: dble_logfact
fact2 = dexp(dble_logfact(n))
endif
end function
double precision function dble_logfact(n) result(logfact2)
implicit none
BEGIN_DOC
! n!!
END_DOC
integer :: n
integer :: k
double precision :: prod
prod=0.d0
do k=2,n,2
prod=prod+dlog(dfloat(k))
enddo
logfact2=prod
return
end function
subroutine write_git_log(iunit)
implicit none
BEGIN_DOC
! Write the last git commit in file iunit.
END_DOC
integer, intent(in) :: iunit
write(iunit,*) '----------------'
write(iunit,*) 'Last git commit:'
BEGIN_SHELL [ /bin/bash ]
git log -1 2>/dev/null | sed "s/'//g"| sed "s/^/ write(iunit,*) '/g" | sed "s/$/'/g" || echo "Unknown"
END_SHELL
write(iunit,*) '----------------'
end
BEGIN_PROVIDER [ double precision, inv_int, (128) ]
implicit none
BEGIN_DOC
! 1/i
END_DOC
integer :: i
do i=1,128
inv_int(i) = 1.d0/dble(i)
enddo
END_PROVIDER
subroutine wall_time(t)
implicit none
BEGIN_DOC
! The equivalent of cpu_time, but for the wall time.
END_DOC
double precision, intent(out) :: t
integer*8 :: c
integer*8, save :: rate = 0
if (rate == 0) then
CALL SYSTEM_CLOCK(count_rate=rate)
endif
CALL SYSTEM_CLOCK(count=c)
t = dble(c)/dble(rate)
end
BEGIN_PROVIDER [ integer, nproc ]
use omp_lib
implicit none
BEGIN_DOC
! Number of current OpenMP threads
END_DOC
nproc = 1
!$OMP PARALLEL
!$OMP MASTER
!$ nproc = omp_get_num_threads()
!$OMP END MASTER
!$OMP END PARALLEL
END_PROVIDER
double precision function u_dot_v(u,v,sze)
implicit none
BEGIN_DOC
! Compute <u|v>
END_DOC
integer, intent(in) :: sze
double precision, intent(in) :: u(sze),v(sze)
double precision, external :: ddot
!DIR$ FORCEINLINE
u_dot_v = ddot(sze,u,1,v,1)
end
double precision function u_dot_u(u,sze)
implicit none
BEGIN_DOC
! Compute <u|u>
END_DOC
integer, intent(in) :: sze
double precision, intent(in) :: u(sze)
double precision, external :: ddot
!DIR$ FORCEINLINE
u_dot_u = ddot(sze,u,1,u,1)
end
subroutine normalize(u,sze)
implicit none
BEGIN_DOC
! Normalizes vector u
END_DOC
integer, intent(in) :: sze
double precision, intent(inout):: u(sze)
double precision :: d
double precision, external :: dnrm2
integer :: i
!DIR$ FORCEINLINE
d = dnrm2(sze,u,1)
if (d /= 0.d0) then
d = 1.d0/d
endif
if (d /= 1.d0) then
!DIR$ FORCEINLINE
call dscal(sze,d,u,1)
endif
end
double precision function approx_dble(a,n)
implicit none
integer, intent(in) :: n
double precision, intent(in) :: a
double precision :: f
integer :: i
if (a == 0.d0) then
approx_dble = 0.d0
return
endif
f = 1.d0
do i=1,-int(dlog10(dabs(a)))+n
f = f*.1d0
enddo
do i=1,int(dlog10(dabs(a)))-n
f = f*10.d0
enddo
approx_dble = dnint(a/f)*f
end
subroutine lowercase(txt,n)
implicit none
BEGIN_DOC
! Transform to lower case
END_DOC
character*(*), intent(inout) :: txt
integer, intent(in) :: n
character( * ), PARAMETER :: LOWER_CASE = 'abcdefghijklmnopqrstuvwxyz'
character( * ), PARAMETER :: UPPER_CASE = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
integer :: i, ic
do i=1,n
ic = index( UPPER_CASE, txt(i:i) )
if (ic /= 0) then
txt(i:i) = LOWER_CASE(ic:ic)
endif
enddo
end
subroutine v2_over_x(v,x,res)
!BEGIN_DOC
! Two by two diagonalization to avoid the divergence in v^2/x when x goes to 0
!END_DOC
implicit none
double precision, intent(in) :: v, x
double precision, intent(out) :: res
double precision :: delta_E, tmp, val
res = 0d0
delta_E = x
if (v == 0.d0) return
val = 2d0 * v
tmp = dsqrt(delta_E * delta_E + val * val)
if (delta_E < 0.d0) then
tmp = -tmp
endif
res = 0.5d0 * (tmp - delta_E)
end
! ---
subroutine check_sym(A, n)
implicit none
integer, intent(in) :: n
double precision, intent(in) :: A(n,n)
integer :: i, j
double precision :: dev_sym, norm, tmp
dev_sym = 0.d0
norm = 0.d0
do i = 1, n
do j = i+1, n
tmp = A(j,i) - A(i,j)
dev_sym += tmp * tmp
norm += A(j,i) * A(j,i)
enddo
enddo
print*, ' deviation from sym = ', dev_sym
print*, ' norm = ', norm
end subroutine check_sym
! ---
subroutine sum_A_At(A, N)
!BEGIN_DOC
! add a tensor with its transpose without a temporary tensor
!END_DOC
implicit none
integer, intent(in) :: N
double precision, intent(inout) :: A(N,N)
integer :: i, j
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i, j) &
!$OMP SHARED (A, N)
!$OMP DO
do j = 1, N
do i = j, N
A(i,j) += A(j,i)
enddo
enddo
!$OMP END DO
!$OMP DO
do j = 2, N
do i = 1, j-1
A(i,j) = A(j,i)
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
end
! ---
subroutine sub_A_At(A, N)
!BEGIN_DOC
! substruct a tensor with its transpose without a temporary tensor
!END_DOC
implicit none
integer, intent(in) :: N
double precision, intent(inout) :: A(N,N)
integer :: i, j
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i, j) &
!$OMP SHARED (A, N)
!$OMP DO
do j = 1, N
do i = j, N
A(i,j) -= A(j,i)
enddo
enddo
!$OMP END DO
!$OMP DO
do j = 2, N
do i = 1, j-1
A(i,j) = -A(j,i)
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
end
! ---
logical function is_same_spin(sigma_1, sigma_2)
BEGIN_DOC
!
! true if sgn(sigma_1) = sgn(sigma_2)
!
END_DOC
implicit none
double precision, intent(in) :: sigma_1, sigma_2
if((sigma_1 * sigma_2) .gt. 0.d0) then
is_same_spin = .true.
else
is_same_spin = .false.
endif
end function is_same_spin
! ---
function Kronecker_delta(i, j) result(delta)
BEGIN_DOC
! Kronecker Delta
END_DOC
implicit none
integer, intent(in) :: i, j
double precision :: delta
if(i == j) then
delta = 1.d0
else
delta = 0.d0
endif
end function Kronecker_delta
! ---
subroutine diagonalize_sym_matrix(N, A, e)
BEGIN_DOC
!
! Diagonalize a symmetric matrix
!
END_DOC
implicit none
integer, intent(in) :: N
double precision, intent(inout) :: A(N,N)
double precision, intent(out) :: e(N)
integer :: lwork, info
double precision, allocatable :: work(:)
allocate(work(1))
lwork = -1
call dsyev('V', 'U', N, A, N, e, work, lwork, info)
lwork = int(work(1))
deallocate(work)
allocate(work(lwork))
call dsyev('V', 'U', N, A, N, e, work, lwork, info)
deallocate(work)
if(info /= 0) then
print*,'Problem in diagonalize_sym_matrix (dsyev)!!'
endif
end subroutine diagonalize_sym_matrix
! ---
View Git Blame Copy Permalink