From 75d099a2e80de84628d33375813143a08d984193 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 14 Jan 2016 13:07:44 +0100 Subject: [PATCH] Fixed DMC --- ocaml/Input.ml | 12 +- src/SAMPLING/dmc_step.irp.f | 231 +++++++++++++++++++-------------- src/SAMPLING/reconfigure.irp.f | 45 ++++--- src/SAMPLING/vmc_step.irp.f | 122 ++--------------- src/TOOLS/random.irp.f | 58 +-------- 5 files changed, 192 insertions(+), 276 deletions(-) diff --git a/ocaml/Input.ml b/ocaml/Input.ml index 652659d..34a8058 100644 --- a/ocaml/Input.ml +++ b/ocaml/Input.ml @@ -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 diff --git a/src/SAMPLING/dmc_step.irp.f b/src/SAMPLING/dmc_step.irp.f index 706e3e0..3e3c01b 100644 --- a/src/SAMPLING/dmc_step.irp.f +++ b/src/SAMPLING/dmc_step.irp.f @@ -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 + + diff --git a/src/SAMPLING/reconfigure.irp.f b/src/SAMPLING/reconfigure.irp.f index c4e9002..c9434fa 100644 --- a/src/SAMPLING/reconfigure.irp.f +++ b/src/SAMPLING/reconfigure.irp.f @@ -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 diff --git a/src/SAMPLING/vmc_step.irp.f b/src/SAMPLING/vmc_step.irp.f index c9e3c62..7aedf42 100644 --- a/src/SAMPLING/vmc_step.irp.f +++ b/src/SAMPLING/vmc_step.irp.f @@ -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 - diff --git a/src/TOOLS/random.irp.f b/src/TOOLS/random.irp.f index 6939e5c..65bed3e 100644 --- a/src/TOOLS/random.irp.f +++ b/src/TOOLS/random.irp.f @@ -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