mirror of
https://gitlab.com/scemama/qmcchem.git
synced 2024-11-08 23:23:39 +01:00
148 lines
3.9 KiB
Fortran
148 lines
3.9 KiB
Fortran
double precision function qmc_ranf()
|
|
! L'Ecuyer, P. (1999) `Tables of maximally equidistributed combined LFSR
|
|
! generators', Math. of Comput., 68, 261-269.
|
|
implicit none
|
|
integer*8 :: b(2)
|
|
b(1) = ISHFT( IEOR( ISHFT(seed(1),1), seed(1)), -53)
|
|
b(2) = ISHFT( IAND(seed(1),-2_8), 10)
|
|
seed(1) = IEOR( b(2), b(1))
|
|
|
|
b(1) = ISHFT( IEOR( ISHFT(seed(2),24), seed(2)), -50)
|
|
b(2) = ISHFT( IAND(seed(2),-512_8), 5)
|
|
seed(2) = IEOR( b(2), b(1))
|
|
|
|
b(1) = ISHFT( IEOR( ISHFT(seed(3),3), seed(3)), -23)
|
|
b(2) = ISHFT( IAND(seed(3),-4096_8), 29)
|
|
seed(3) = IEOR( b(2), b(1))
|
|
|
|
b(1) = ISHFT( IEOR( ISHFT(seed(4),5), seed(4)), -24)
|
|
b(2) = ISHFT( IAND(seed(4),-131072_8), 23)
|
|
seed(4) = IEOR( b(2), b(1))
|
|
|
|
b(1) = ISHFT( IEOR( ISHFT(seed(5),3), seed(5)), -33)
|
|
b(2) = ISHFT( IAND(seed(5),-8388608_8), 8)
|
|
seed(5) = IEOR( b(2), b(1))
|
|
|
|
qmc_ranf = IEOR( IEOR( IEOR( IEOR(seed(1),seed(2)), seed(3)), &
|
|
seed(4)), seed(5)) * 5.4210108624275221D-20 + 0.5D0
|
|
ASSERT ( qmc_ranf >= 0.d0 )
|
|
ASSERT ( qmc_ranf <= 1.d0 )
|
|
|
|
end
|
|
|
|
subroutine ranf_array(isize,res)
|
|
implicit none
|
|
integer :: isize
|
|
double precision :: res(isize)
|
|
integer :: i
|
|
double precision :: qmc_ranf
|
|
|
|
do i=1,isize
|
|
res(i) = qmc_ranf()
|
|
enddo
|
|
end
|
|
|
|
BEGIN_PROVIDER [ integer*8, seed, (5) ]
|
|
implicit none
|
|
BEGIN_DOC
|
|
! Seeds data
|
|
! Initialized by init_random
|
|
END_DOC
|
|
integer :: iargc
|
|
integer*8 :: i,j
|
|
integer*4 :: clock(12)
|
|
double precision :: r
|
|
integer*8 :: pid8
|
|
read(current_PID,*) pid8
|
|
pid8 = iand( ishft(pid8, 32), pid8)
|
|
do i=1,12
|
|
clock(i) = i
|
|
enddo
|
|
call system_clock(count=clock(1))
|
|
call random_seed(put=clock)
|
|
do i=1,5
|
|
call random_number(r)
|
|
seed(i) = (r-0.5d0)*huge(1_8)
|
|
seed(i) = ieor( seed(i), pid8)
|
|
do j=1,16
|
|
seed(i) = ishft(seed(i),1)+1
|
|
enddo
|
|
enddo
|
|
|
|
END_PROVIDER
|
|
|
|
subroutine gauss_array(isize,res)
|
|
implicit none
|
|
include '../constants.F'
|
|
integer isize
|
|
double precision res(isize)
|
|
|
|
double precision u1(isize),u2(isize)
|
|
integer i
|
|
|
|
call ranf_array(isize,u1)
|
|
call ranf_array(isize,u2)
|
|
do i=1,isize
|
|
res(i)=sqrt(-2.d0*log(u1(i)))*cos(dtwo_pi*u2(i))
|
|
enddo
|
|
end
|
|
|
|
double precision function gauss()
|
|
implicit none
|
|
! include 'constants.F'
|
|
double precision :: qmc_ranf
|
|
! double precision :: u1,u2
|
|
! u1=qmc_ranf()
|
|
! u2=qmc_ranf()
|
|
! gauss=sqrt(-2.d0*dlog(u1))*cos(dfour_pi*u2)
|
|
double precision :: inverse_normal_cdf
|
|
gauss = inverse_normal_cdf(qmc_ranf())
|
|
end
|
|
|
|
double precision function inverse_normal_cdf(p)
|
|
implicit none
|
|
double precision, intent(in) :: p
|
|
double precision :: p_low,p_high
|
|
double precision :: a1,a2,a3,a4,a5,a6
|
|
double precision :: b1,b2,b3,b4,b5
|
|
double precision :: c1,c2,c3,c4,c5,c6
|
|
double precision :: d1,d2,d3,d4
|
|
double precision :: z,q,r
|
|
double precision :: qmc_ranf
|
|
a1=-39.6968302866538d0
|
|
a2=220.946098424521d0
|
|
a3=-275.928510446969d0
|
|
a4=138.357751867269d0
|
|
a5=-30.6647980661472d0
|
|
a6=2.50662827745924d0
|
|
b1=-54.4760987982241d0
|
|
b2=161.585836858041d0
|
|
b3=-155.698979859887d0
|
|
b4=66.8013118877197d0
|
|
b5=-13.2806815528857d0
|
|
c1=-0.00778489400243029d0
|
|
c2=-0.322396458041136d0
|
|
c3=-2.40075827716184d0
|
|
c4=-2.54973253934373d0
|
|
c5=4.37466414146497d0
|
|
c6=2.93816398269878d0
|
|
d1=0.00778469570904146d0
|
|
d2=0.32246712907004d0
|
|
d3=2.445134137143d0
|
|
d4=3.75440866190742d0
|
|
p_low=0.02425d0
|
|
p_high=1.d0-0.02425d0
|
|
if(p < p_low) then
|
|
q=dsqrt(-2.d0*dlog(p))
|
|
inverse_normal_cdf=(((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q+c6)/((((d1*q+d2)*q+d3)*q+d4)*q+1.d0)
|
|
else if(p <= p_high) then
|
|
q=p-0.5d0
|
|
r=q*q
|
|
inverse_normal_cdf=(((((a1*r+a2)*r+a3)*r+a4)*r+a5)*r+a6)*q/(((((b1*r+b2)*r+b3)*r+b4)*r+b5)*r+1.d0)
|
|
else
|
|
q=dsqrt(-2.d0*dlog(max(tiny(1.d0),1.d0-p)))
|
|
inverse_normal_cdf=-(((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q+c6)/((((d1*q+d2)*q+d3)*q+d4)*q+1)
|
|
endif
|
|
end
|
|
|