LCPQ
/
qmcchem
mirror of https://gitlab.com/scemama/qmcchem.git
10
1
Fork 1
Code Releases Activity

Fixed DMC

This commit is contained in:
Anthony Scemama 2016-01-14 13:07:44 +01:00
parent dc6d6b1637
commit 75d099a2e8
5 changed files with 192 additions and 276 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
ocaml/Input.ml
View File

@ -630,7 +630,7 @@ end = struct
let of_float x =
if (x >= 100.) then
failwith "DMC Projection time should be < 100.";
if (x <= 0.) then
if (x < 0.) then
failwith "DMC Projection time should be positive.";
x
@ -886,6 +886,16 @@ let validate () =
| _ -> ()
in
(*
(* Check Projection time is greater than time step *)
let () =
match (meth, DMC_projection_time.(read () |> to_float) ) with
| (Method.DMC,p) ->
if (p < ts) then failwith "E_ref should not be zero in DMC"
| _ -> ()
in
*)
(* Set block and total time*)
let () =
if ( (Block_time.read ()) > Stop_time.read ()) then

231
src/SAMPLING/dmc_step.irp.f
View File

@ -10,7 +10,7 @@ t = """
&BEGIN_PROVIDER [ $T, $X_2_dmc_block_walk_kahan $D2 ]
implicit none
BEGIN_DOC
! VMC averages of $X
! DMC averages of $X. Computed in E_loc_dmc_block_walk
END_DOC
$X_dmc_block_walk = 0.d0
$X_dmc_block_walk_kahan = 0.d0
@ -27,14 +27,15 @@ for p in properties:
D1 = ", ("+p[2][1:-1]+")"
D2 = ", ("+p[2][1:-1]+",3)"
print t.replace("$X",p[1]).replace("$T",p[0]).replace("$D1",D1).replace("$D2",D2)
END_SHELL
BEGIN_PROVIDER [ double precision, E_loc_dmc_block_walk ]
&BEGIN_PROVIDER [ double precision, E_loc_2_dmc_block_walk ]
BEGIN_PROVIDER [ double precision, E_loc_dmc_block_walk ]
&BEGIN_PROVIDER [ double precision, E_loc_2_dmc_block_walk ]
&BEGIN_PROVIDER [ double precision, E_loc_dmc_block_walk_kahan, (3) ]
&BEGIN_PROVIDER [ double precision, E_loc_2_dmc_block_walk_kahan, (3)
&BEGIN_PROVIDER [ double precision, E_loc_2_dmc_block_walk_kahan, (3) ]
implicit none
include '../types.F'
BEGIN_DOC
@ -43,19 +44,18 @@ END_SHELL
real, allocatable :: elec_coord_tmp(:,:,:)
integer :: mod_align
double precision, allocatable :: psi_grad_psi_inv_save_tmp(:,:,:)
double precision :: psi_value_save(walk_num)
double precision :: psi_value_save_tmp(walk_num)
integer :: trapped_walk_tmp(walk_num)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: psi_grad_psi_inv_save_tmp
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: psi_value_save
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: psi_value_save_tmp
allocate ( elec_coord_tmp(mod_align(elec_num+1),3,walk_num) )
allocate ( psi_grad_psi_inv_save_tmp(elec_num_8,3,walk_num) )
psi_value_save = 0.d0
! Initialization
if (vmc_algo /= t_Brownian) then
call abrt(irp_here,'DMC should run with Brownian algorithm')
endif
PROVIDE E_loc_vmc_block_walk
integer :: k, i_walk, i_step
@ -63,9 +63,13 @@ BEGIN_SHELL [ /usr/bin/python ]
from properties import *
t = """
if (calc_$X) then
!DIR$ VECTOR ALIGNED
$X_dmc_block_walk = 0.d0
!DIR$ VECTOR ALIGNED
$X_dmc_block_walk_kahan = 0.d0
!DIR$ VECTOR ALIGNED
$X_2_dmc_block_walk = 0.d0
!DIR$ VECTOR ALIGNED
$X_2_dmc_block_walk_kahan = 0.d0
$X_min = huge(1.)
$X_max =-huge(1.)
@ -75,65 +79,49 @@ for p in properties:
print t.replace("$X",p[1])
END_SHELL
double precision :: icount
icount = 0.d0
logical :: loop
integer*8 :: cpu0, cpu1, cpu2, count_rate, count_max
loop = .True.
call system_clock(cpu0, count_rate, count_max)
cpu2 = cpu0
block_weight = 0.d0
do while (loop)
! Move to the next projection step
dmc_projection_step = mod(dmc_projection_step,dmc_projection)+1
! Remove contribution of the old value of the weight at the new
! projection step
pop_weight_mult *= 1.d0/pop_weight(dmc_projection_step)
! Compute the new weight of the population
pop_weight(dmc_projection_step) = 0.d0
do k=1,walk_num
pop_weight(dmc_projection_step) += dmc_weight(k)
enddo
! Normalize the weight of the walkers by the weight of the population
do k=1,walk_num
dmc_weight(k) = dmc_weight(k)/pop_weight(dmc_projection_step)
enddo
! Normalize the weight of the population at the current projection step by
! the number of walkers
pop_weight(dmc_projection_step) = pop_weight(dmc_projection_step)/dble(walk_num)
! Update the running population weight
pop_weight_mult *= pop_weight(dmc_projection_step)
SOFT_TOUCH pop_weight_mult
! Every walker makes a step
do i_walk=1,walk_num
integer :: i,j,l
do l=1,3
do i=1,elec_num+1
elec_coord(i,l) = elec_coord_full(i,l,i_walk)
enddo
enddo
TOUCH elec_coord
BEGIN_SHELL [ /usr/bin/python ]
from properties import *
t = """
if (calc_$X) then
! Kahan's summation algorithm to compute these sums reducing the rounding error:
! $X_dmc_block_walk($D) += $X * pop_weight_mult
! $X_2_dmc_block_walk($D) += $X_2 * pop_weight_mult
! see http://en.wikipedia.org/wiki/Kahan_summation_algorithm
$X_dmc_block_walk_kahan($D2 3) = $X * pop_weight_mult - $X_dmc_block_walk_kahan($D2 1)
$X_dmc_block_walk_kahan($D2 2) = $X_dmc_block_walk $D1 + $X_dmc_block_walk_kahan($D2 3)
$X_dmc_block_walk_kahan($D2 1) = ($X_dmc_block_walk_kahan($D2 2) - $X_dmc_block_walk $D1 ) &
- $X_dmc_block_walk_kahan($D2 3)
$X_dmc_block_walk $D1 = $X_dmc_block_walk_kahan($D2 2)
$X_2_dmc_block_walk_kahan($D2 3) = $X_2 * pop_weight_mult - $X_2_dmc_block_walk_kahan($D2 1)
$X_2_dmc_block_walk_kahan($D2 2) = $X_2_dmc_block_walk $D1 + $X_2_dmc_block_walk_kahan($D2 3)
$X_2_dmc_block_walk_kahan($D2 1) = ($X_2_dmc_block_walk_kahan($D2 2) - $X_2_dmc_block_walk $D1 ) &
- $X_2_dmc_block_walk_kahan($D2 3)
$X_2_dmc_block_walk $D1 = $X_2_dmc_block_walk_kahan($D2 2)
endif
if (calc_$X) then
! Kahan's summation algorithm to compute these sums reducing the rounding error:
! $X_dmc_block_walk += $X * pop_weight_mult
! $X_2_dmc_block_walk += $X_2 * pop_weight_mult
! see http://en.wikipedia.org/wiki/Kahan_summation_algorithm
$X_dmc_block_walk_kahan($D2 3) = $X * pop_weight_mult - $X_dmc_block_walk_kahan($D2 1)
$X_dmc_block_walk_kahan($D2 2) = $X_dmc_block_walk $D1 + $X_dmc_block_walk_kahan($D2 3)
$X_dmc_block_walk_kahan($D2 1) = ($X_dmc_block_walk_kahan($D2 2) - $X_dmc_block_walk $D1 ) &
- $X_dmc_block_walk_kahan($D2 3)
$X_dmc_block_walk $D1 = $X_dmc_block_walk_kahan($D2 2)
$X_2_dmc_block_walk_kahan($D2 3) = $X_2 * pop_weight_mult - $X_2_dmc_block_walk_kahan($D2 1)
$X_2_dmc_block_walk_kahan($D2 2) = $X_2_dmc_block_walk $D1 + $X_2_dmc_block_walk_kahan($D2 3)
$X_2_dmc_block_walk_kahan($D2 1) = ($X_2_dmc_block_walk_kahan($D2 2) - $X_2_dmc_block_walk $D1 ) &
- $X_2_dmc_block_walk_kahan($D2 3)
$X_2_dmc_block_walk $D1 = $X_2_dmc_block_walk_kahan($D2 2)
endif
"""
for p in properties:
if p[2] == "":
@ -143,25 +131,83 @@ for p in properties:
D1 = "("+":"*(p[2].count(',')+1)+")"
D2 = ":"*(p[2].count(',')+1)+","
print t.replace("$X",p[1]).replace("$D1",D1).replace("$D2",D2)
END_SHELL
icount += pop_weight_mult
double precision :: p,q
real :: delta_x
logical :: accepted
call brownian_step(p,q,accepted,delta_x)
if (accepted) then
trapped_walk(i_walk) = 0
else
trapped_walk(i_walk) += 1
endif
if ( (trapped_walk(i_walk) < trapped_walk_max).and. &
(psi_value * psi_value_save(i_walk) >= 0.d0) ) then
dmc_weight(i_walk) = dexp(dtime_step*(E_ref - E_loc))
else
dmc_weight(i_walk) = 0.d0
trapped_walk(i_walk) = 0
endif
elec_coord(elec_num+1,1) += p*time_step
elec_coord(elec_num+1,2) = E_loc
elec_coord(elec_num+1,3) = dmc_weight(i_walk)
do l=1,3
do i=1,elec_num+1
elec_coord_full(i,l,i_walk) = elec_coord(i,l)
enddo
enddo
psi_value_save(i_walk) = psi_value
enddo
! Move to the next projection step
if (dmc_projection > 0) then
dmc_projection_step = mod(dmc_projection_step,dmc_projection)+1
else
dmc_projection_step = 1
endif
! Eventually, recompute the weight of the population
if (dmc_projection_step == 1) then
pop_weight_mult = 1.d0
do k=1,dmc_projection
pop_weight_mult *= pop_weight(k)
enddo
endif
! Remove contribution of the old value of the weight at the new
! projection step
pop_weight_mult *= 1.d0/pop_weight(dmc_projection_step)
! Compute the new weight of the population
double precision :: sum_weight
sum_weight = 0.d0
do k=1,walk_num
sum_weight += dmc_weight(k)
enddo
pop_weight(dmc_projection_step) = sum_weight/dble(walk_num)
! Update the running population weight
pop_weight_mult *= pop_weight(dmc_projection_step)
block_weight += pop_weight_mult * dble(walk_num)
! Reconfiguration
integer :: ipos(walk_num)
call reconfigure(ipos,dmc_weight)
do k=1,walk_num
integer :: i, l
do l=1,3
do i=1,elec_num+1
elec_coord_tmp(i,l,k) = elec_coord_full(i,l,k)
enddo
!DIR$ VECTOR ALIGNED
!DIR$ LOOP COUNT(200)
do i=1,elec_num
psi_grad_psi_inv_save_tmp(i,l,k) = psi_grad_psi_inv_save(i,l,k)
enddo
enddo
psi_value_save_tmp(k) = psi_value_save(k)
trapped_walk_tmp(k) = trapped_walk(k)
@ -174,36 +220,11 @@ END_SHELL
do i=1,elec_num+1
elec_coord_full(i,l,k) = elec_coord_tmp(i,l,ipm)
enddo
!DIR$ VECTOR ALIGNED
!DIR$ LOOP COUNT(200)
do i=1,elec_num
psi_grad_psi_inv_save(i,l,k) = psi_grad_psi_inv_save_tmp(i,l,ipm)
enddo
enddo
psi_value_save(k) = psi_value_save_tmp(ipm)
trapped_walk(k) = trapped_walk_tmp(ipm)
enddo
! Set 1st walker
!DIR$ VECTOR ALIGNED
!DIR$ LOOP COUNT(200)
do i=1,elec_num
psi_grad_psi_inv_x(i) = psi_grad_psi_inv_save(i,1,1)
psi_grad_psi_inv_y(i) = psi_grad_psi_inv_save(i,2,1)
psi_grad_psi_inv_z(i) = psi_grad_psi_inv_save(i,3,1)
enddo
!DIR$ VECTOR UNALIGNED
!DIR$ LOOP COUNT(200)
do i=1,elec_num
elec_coord(i,1) = elec_coord_full(i,1,1)
elec_coord(i,2) = elec_coord_full(i,2,1)
elec_coord(i,3) = elec_coord_full(i,3,1)
enddo
psi_value = psi_value_save(1)
TOUCH elec_coord_full psi_value_save psi_grad_psi_inv_save psi_value psi_grad_psi_inv_x psi_grad_psi_inv_y psi_grad_psi_inv_z elec_coord
call system_clock(cpu1, count_rate, count_max)
if (cpu1 < cpu0) then
cpu1 = cpu1+cpu0
@ -216,14 +237,14 @@ END_SHELL
cpu2 = cpu1
endif
SOFT_TOUCH elec_coord_full psi_value psi_grad_psi_inv_x psi_grad_psi_inv_y psi_grad_psi_inv_z elec_coord pop_weight_mult
enddo
double precision :: factor
factor = 1.d0/icount
block_weight *= icount
factor = 1.d0/block_weight
SOFT_TOUCH block_weight
BEGIN_SHELL [ /usr/bin/python ]
from properties import *
t = """
@ -237,7 +258,6 @@ for p in properties:
END_SHELL
deallocate ( elec_coord_tmp )
deallocate ( psi_grad_psi_inv_save_tmp )
END_PROVIDER
@ -272,12 +292,33 @@ END_PROVIDER
dmc_projection_step = 0
END_PROVIDER
BEGIN_PROVIDER [ double precision, pop_weight, (dmc_projection) ]
BEGIN_PROVIDER [ double precision, pop_weight, (0:dmc_projection+1) ]
implicit none
BEGIN_DOC
! Population weight of DMC
END_DOC
pop_weight = 1.d0
pop_weight(dmc_projection) = 1.d0/dble(dmc_projection)
pop_weight(dmc_projection) = 1.d0/dble(size(pop_weight))
END_PROVIDER
BEGIN_PROVIDER [ integer, trapped_walk, (walk_num_8) ]
&BEGIN_PROVIDER [ integer, trapped_walk_max ]
implicit none
BEGIN_DOC
! Number of steps when the walkers were trapped
END_DOC
trapped_walk = 0
trapped_walk_max = 20
END_PROVIDER
BEGIN_PROVIDER [ double precision, dmc_weight, (walk_num_8) ]
implicit none
BEGIN_DOC
! Weight of the walkers in the DMC algorithm: exp(-time_step*(E_loc-E_ref))
END_DOC
!DIR$ VECTOR ALIGNED
dmc_weight = 1.d0
END_PROVIDER

45
src/SAMPLING/reconfigure.irp.f
View File

@ -3,44 +3,49 @@ subroutine reconfigure(ipos,w)
integer, intent(inout) :: ipos(*)
double precision, intent(in) :: w(*)
integer :: kp, km
double precision :: accup, accum
integer :: k
double precision :: dwalk_num
dwalk_num = dble(walk_num)
integer :: kptab(walk_num), kmtab(walk_num)
double precision :: wp(walk_num), wm(walk_num)
double precision :: tmp
double precision :: dwalk_num
tmp = 0.d0
do k=1,walk_num
ipos(k) = k
tmp = tmp + w(k)
enddo
dwalk_num = dble(walk_num)/tmp
integer :: kp, km
kp=0
km=0
double precision :: accup, accum
accup = 0.d0
accum = 0.d0
integer :: k
do k=1,walk_num
tmp = dwalk_num*w(k)-1.d0
if (tmp >= 0.d0) then
kp += 1
wp(kp) = abs(tmp)
accup += wp(kp)
kp = kp+1
wp(kp) = dabs(tmp)
accup = accup + wp(kp)
kptab(kp) = k
else
km += 1
wm(km) = abs(tmp)
accum += wm(km)
km = km+1
wm(km) = dabs(tmp)
accum = accum + wm(km)
kmtab(km) = k
endif
enddo
if(kp+km /= walk_num) then
print *, kp, km
call abrt(irp_here,'pb in reconfiguration +/-')
endif
if(abs(accup-accum).gt.1.d-11) then
if(dabs(accup-accum) > 1.d-11) then
print *, accup, accum
call abrt(irp_here,'pb in reconfiguration')
endif
@ -59,24 +64,26 @@ subroutine reconfigure(ipos,w)
averageconf = accup
kcp = 1
rand = rando(kcp)
do while (rand < averageconf)
k=1
current=wm(k)
do while (rand > current)
k += 1
current += wm(k)
k = k+1
current = current + wm(k)
enddo
kremove = kmtab(k)
k=1
current=wp(k)
do while (rand > current)
k += 1
current += wp(k)
k = k+1
current = current + wp(k)
enddo
kadd = kptab(k)
ipos(kremove) = kadd
kcp += 1
kcp = kcp + 1
rand = rando(kcp)
enddo

122
src/SAMPLING/vmc_step.irp.f
View File

@ -39,7 +39,6 @@ END_SHELL
PROVIDE time_step
BEGIN_SHELL [ /usr/bin/python ]
from properties import *
t = """
@ -61,7 +60,6 @@ for p in properties:
END_SHELL
double precision :: dnorm
!DIR$ VECTOR ALIGNED
block_weight = 0.d0
do i_walk=1,walk_num
integer :: i,j,l
@ -71,24 +69,7 @@ END_SHELL
elec_coord(i,l) = elec_coord_full(i,l,i_walk)
enddo
enddo
endif
PROVIDE psi_grad_psi_inv_save psi_value_save
if (psi_value_save(walk_num) /= 0.) then
psi_value = psi_value_save(i_walk)
!DIR$ VECTOR ALIGNED
!DIR$ LOOP COUNT(200)
do i=1,elec_num
psi_grad_psi_inv_x(i) = psi_grad_psi_inv_save(i,1,i_walk)
psi_grad_psi_inv_y(i) = psi_grad_psi_inv_save(i,2,i_walk)
psi_grad_psi_inv_z(i) = psi_grad_psi_inv_save(i,3,i_walk)
enddo
TOUCH psi_value psi_grad_psi_inv_x psi_grad_psi_inv_y psi_grad_psi_inv_z elec_coord
else
if (i_walk > 1) then
TOUCH elec_coord
endif
TOUCH elec_coord
endif
logical :: loop
@ -101,7 +82,6 @@ END_SHELL
double precision :: p,q
real :: delta_x
logical :: accepted
double precision :: E_old
if (vmc_algo == t_Brownian) then
call brownian_step(p,q,accepted,delta_x)
else if (vmc_algo == t_Langevin) then
@ -109,15 +89,10 @@ END_SHELL
endif
elec_coord(elec_num+1,1) += p*time_step
elec_coord(elec_num+1,2) = E_loc
elec_coord(elec_num+1,3) += p*time_step
if (accepted) then
trapped_walk(i_walk) = 0
else
trapped_walk(i_walk) += 1
endif
elec_coord(elec_num+1,3) = 1.
block_weight += 1.d0
BEGIN_SHELL [ /usr/bin/python ]
from properties import *
t = """
@ -153,21 +128,17 @@ for p in properties:
END_SHELL
if ( qmc_method == t_VMC ) then
call system_clock(cpu1, count_rate, count_max)
if (cpu1 < cpu0) then
cpu1 = cpu1+cpu0
endif
loop = dble(cpu1-cpu0)*dble(walk_num) < dble(block_time)*dble(count_rate)
if (cpu1-cpu2 > count_rate) then
integer :: do_run
call get_running(do_run)
loop = do_run == t_Running
cpu2 = cpu1
endif
else
loop = .False.
endif
call system_clock(cpu1, count_rate, count_max)
if (cpu1 < cpu0) then
cpu1 = cpu1+cpu0
endif
loop = dble(cpu1-cpu0)*dble(walk_num) < dble(block_time)*dble(count_rate)
if (cpu1-cpu2 > count_rate) then
integer :: do_run
call get_running(do_run)
loop = do_run == t_Running
cpu2 = cpu1
endif
enddo ! while (loop)
@ -177,27 +148,6 @@ END_SHELL
enddo
enddo
if (qmc_method == t_DMC) then
if ( (trapped_walk(i_walk) < trapped_walk_max).and. &
(psi_value * psi_value_save(i_walk) >= 0.d0) ) then
dmc_weight(i_walk) = dexp(dtime_step*(E_ref - E_loc))
else
dmc_weight(i_walk) = 0.d0
trapped_walk(i_walk) = 0
endif
psi_value_save(i_walk) = psi_value
!DIR$ VECTOR ALIGNED
!DIR$ LOOP COUNT (200)
do i=1,elec_num
psi_grad_psi_inv_save(i,1,i_walk) = psi_grad_psi_inv_x(i)
psi_grad_psi_inv_save(i,2,i_walk) = psi_grad_psi_inv_y(i)
psi_grad_psi_inv_save(i,3,i_walk) = psi_grad_psi_inv_z(i)
enddo
endif
enddo
double precision :: factor
@ -220,47 +170,3 @@ END_SHELL
END_PROVIDER
BEGIN_PROVIDER [ double precision, dmc_weight, (walk_num_8) ]
implicit none
BEGIN_DOC
! Weight of the walkers in the DMC algorithm: exp(-time_step*(E_loc-E_ref))
END_DOC
!DIR$ VECTOR ALIGNED
dmc_weight = 1.d0
END_PROVIDER
BEGIN_PROVIDER [ double precision, psi_grad_psi_inv_save, (elec_num_8,3,walk_num) ]
&BEGIN_PROVIDER [ double precision, psi_value_save, (walk_num_8) ]
implicit none
BEGIN_DOC
! psi_grad_psi_inv of the previous step to accelerate DMC
!
! updated in vmc_step
END_DOC
integer, save :: ifirst = 0
if (ifirst == 0) then
psi_grad_psi_inv_save = 0.d0
psi_value_save = 0.d0
ifirst = 1
endif
END_PROVIDER
BEGIN_PROVIDER [ integer, trapped_walk, (walk_num_8) ]
implicit none
BEGIN_DOC
! Number of steps when the walkers were trapped
END_DOC
trapped_walk = 0
END_PROVIDER
BEGIN_PROVIDER [ integer, trapped_walk_max ]
implicit none
BEGIN_DOC
! Max number of trapped MC steps before killing walker
END_DOC
trapped_walk_max = 20
END_PROVIDER

58
src/TOOLS/random.irp.f
View File

@ -89,59 +89,11 @@ end
double precision function gauss()
implicit none
! include 'constants.F'
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
double precision :: u1,u2
u1=qmc_ranf()
u2=qmc_ranf()
gauss=dsqrt(-2.d0*dlog(u1))*dcos(dfour_pi*u2)
end