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

ZV energy

This commit is contained in:
Anthony Scemama 2016-05-02 21:19:36 +02:00
parent 8cc2c6a24b
commit 571df84d9d
11 changed files with 416 additions and 52 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

32
ocaml/Input.ml
View File

@ -387,7 +387,7 @@ end
module Method : sig
type t = VMC | DMC | SRMC | FKMC
type t = VMC | DMC | SRMC | FKMC | PDMC
val doc : string
val read : unit -> t
val write : t -> unit
@ -396,22 +396,24 @@ module Method : sig
end = struct
type t = VMC | DMC | SRMC | FKMC
type t = VMC | DMC | SRMC | FKMC | PDMC
let doc = "QMC Method : [ VMC | DMC | SRMC | FKMC ]"
let doc = "QMC Method : [ VMC | DMC | SRMC | FKMC | PDMC ]"
let of_string = function
| "VMC" | "vmc" -> VMC
| "DMC" | "dmc" -> DMC
| "SRMC" | "srmc" -> SRMC
| "PDMC" | "pdmc" -> PDMC
| "FKMC" | "fkmc" -> FKMC
| x -> failwith ("Method should be [ VMC | DMC | SRMC | FKMC ], not "^x^".")
| x -> failwith ("Method should be [ VMC | DMC | SRMC | FKMC | PDMC ], not "^x^".")
let to_string = function
| VMC -> "VMC"
| DMC -> "DMC"
| SRMC -> "SRMC"
| PDMC -> "PDMC"
| FKMC -> "FKMC"
@ -852,25 +854,19 @@ let validate () =
(* Check sampling and time steps *)
let () =
match (sampling, meth, Pseudo.to_bool do_pseudo) with
| (Sampling.Brownian, Method.DMC, true)
| (Sampling.Brownian, Method.FKMC, true)
| (Sampling.Brownian, Method.SRMC, true) ->
if ( (Time_step.to_float ts) >= 0.5 ) then
warn ( "Time step seems large for "^(Method.to_string meth) )
| (Sampling.Brownian, Method.SRMC, false)
| (Sampling.Brownian, Method.FKMC, false)
| (Sampling.Brownian, Method.DMC, false) ->
if ( (Time_step.to_float ts) >= 0.01 ) then
warn ( "Time step seems large for "^(Method.to_string meth) )
| (Sampling.Brownian, Method.VMC, _) ->
if ( (Time_step.to_float ts) >= 10. ) then
warn "Time step seems large for VMC."
| (Sampling.Langevin, Method.VMC, _) ->
if ( (Time_step.to_float ts) <= 0.01 ) then
warn "Time step seems small for Langevin sampling."
| (Sampling.Langevin, Method.SRMC, _)
| (Sampling.Langevin, Method.FKMC, _)
| (Sampling.Langevin, Method.DMC, _) ->
| (Sampling.Brownian, _, true) ->
if ( (Time_step.to_float ts) >= 0.5 ) then
warn ( "Time step seems large for "^(Method.to_string meth) )
| (Sampling.Brownian, _, false) ->
if ( (Time_step.to_float ts) >= 0.01 ) then
warn ( "Time step seems large for "^(Method.to_string meth) )
| (Sampling.Langevin, _, _) ->
failwith "Lanvegin sampling is incompatible with DMC"
in
@ -879,6 +875,7 @@ let validate () =
let () =
match (meth, Ref_energy.(read () |> to_float) ) with
| (Method.SRMC,0.)
| (Method.PDMC,0.)
| (Method.FKMC,0.)
| (Method.DMC,0.) -> failwith ("E_ref should not be zero in "^(Method.to_string meth) )
| _ -> ()
@ -894,6 +891,7 @@ let validate () =
let () =
match (meth, Property.(calc E_loc)) with
| (Method.SRMC, false)
| (Method.PDMC, false)
| (Method.FKMC, false)
| (Method.DMC, false) -> failwith ( "E_loc should be sampled in "^(Method.to_string meth) )
| (Method.VMC, false) -> warn "Sampling of E_loc is not activated in input"

9
ocaml/Qmcchem_result.ml
View File

@ -145,6 +145,15 @@ let display_summary ~range =
in
List.iter properties ~f:print_property ;
(** TODO *)
let open Random_variable in
let p = (of_raw_data ~range Property.E_loc_zv)
+! (of_raw_data ~range Property.E_loc)
in
Printf.printf "%20s : %s\n"
("E_loc_zv(+)")
(Random_variable.to_string p);
(** TODO *)
let cpu =
Random_variable.of_raw_data ~range Property.Cpu

33
ocaml/Random_variable.ml
View File

@ -162,6 +162,8 @@ let average { property ; data } =
(** Compute sum (for CPU/Wall time) *)
let sum { property ; data } =
List.fold data ~init:0. ~f:(fun accu x ->
@ -342,7 +344,9 @@ let max_block =
(** Create a hash table for merging *)
let create_hash ~hashable ~create_key ?(update_block_id=(fun x->x)) t =
let create_hash ~hashable ~create_key ?(update_block_id=(fun x->x))
?(update_value=(fun wc vc wb vb sw -> (wc *. vc +. wb *. vb) /. sw) )
?(update_weight=(fun wc wb -> wc +. wb) ) t =
let table = Hashtbl.create ~hashable:hashable ()
in
List.iter t.data ~f:(fun block ->
@ -356,7 +360,7 @@ let create_hash ~hashable ~create_key ?(update_block_id=(fun x->x)) t =
Weight.to_float block.weight
in
let sw =
wc +. wb
update_weight wc wb
in
if (Property.is_scalar current.property) then
let vc, vb =
@ -365,7 +369,7 @@ let create_hash ~hashable ~create_key ?(update_block_id=(fun x->x)) t =
in Some
{ property = current.property ;
weight = Weight.of_float sw ;
value = Sample.of_float ((wc *. vc +. wb *. vb) /. sw);
value = Sample.of_float (update_value wc vc wb vb sw);
block_id = update_block_id block.block_id;
pid = block.pid ;
compute_node = block.compute_node;
@ -380,7 +384,7 @@ let create_hash ~hashable ~create_key ?(update_block_id=(fun x->x)) t =
{ property = current.property ;
weight = Weight.of_float sw ;
value =
Array.init dim ~f:(fun i -> ((wc *. vc.(i) +. wb *. vb.(i)) /. sw))
Array.init dim ~f:(fun i -> update_value wc vc.(i) wb vb.(i) sw)
|> Sample.of_float_array ~dim ;
block_id = update_block_id block.block_id;
pid = block.pid ;
@ -401,9 +405,8 @@ let create_hash ~hashable ~create_key ?(update_block_id=(fun x->x)) t =
(** Genergic merge function *)
let merge ~hashable ~create_key ?update_block_id t =
let table = create_hash ~hashable:hashable ~create_key:create_key
?update_block_id:update_block_id t
let merge ~hashable ~create_key ?update_block_id ?update_value ?update_weight t =
let table = create_hash ~hashable ~create_key ?update_block_id ?update_value ?update_weight t
in
{ property = t.property ;
data = Hashtbl.to_alist table
@ -454,6 +457,20 @@ let merge_per_compute_node_and_block_id =
(Block_id.to_int block.Block.block_id) )
(** Create a new random variable which is a sum of 2 *)
let (+!) p1 p2 =
merge
~hashable:String.hashable
~create_key:(fun block ->
Printf.sprintf "%s %10.10d %10.10d"
(Compute_node.to_string block.Block.compute_node)
(Block_id.to_int block.Block.block_id)
(Pid.to_int block.Block.pid) )
~update_value: (fun wc vc wb vb sw -> (vc +. vb) )
~update_weight:(fun wc wb -> wc )
{ property = p1.property ;
data = List.concat [ p1.data ; p2.data ] }
(** Merge two consecutive blocks *)
@ -467,7 +484,7 @@ let compress =
((Block_id.to_int block_id)+1)/2
|> Block_id.of_int )
(** Last value on each compute node (for wall_time) *)

2
src/PROPERTIES/properties_energy.irp.f
View File

@ -253,7 +253,7 @@ BEGIN_PROVIDER [ double precision, E_loc_zv ]
BEGIN_DOC
! Zero-variance parameter on E_loc
END_DOC
E_loc_zv = E_ref - E_loc
E_loc_zv = 0.d0
END_PROVIDER

6
src/SAMPLING/block.irp.f
View File

@ -30,6 +30,12 @@ t = """
$X_block_walk = $X_srmc_block_walk
$X_2_block_walk = $X_2_srmc_block_walk
endif
else if (qmc_method == t_PDMC) then
PROVIDE E_loc_pdmc_block_walk
if (calc_$X) then
$X_block_walk = $X_pdmc_block_walk
$X_2_block_walk = $X_2_pdmc_block_walk
endif
else if (qmc_method == t_FKMC) then
PROVIDE E_loc_fkmc_block_walk
if (calc_$X) then

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

@ -150,11 +150,6 @@ END_SHELL
psi_value_save(i_walk) = psi_value
E_loc_save(i_walk) = E_loc
if (calc_E_loc_zv) then
E_loc_zv = (E_ref - E_loc)/dmc_weight(i_walk)
TOUCH E_loc_zv
endif
BEGIN_SHELL [ /usr/bin/python ]
from properties import *
t = """

325
src/SAMPLING/pdmc_step.irp.f Normal file
View File

@ -0,0 +1,325 @@
! Providers of *_pdmc_block_walk
!==============================
BEGIN_SHELL [ /usr/bin/python ]
from properties import *
t = """
BEGIN_PROVIDER [ $T, $X_pdmc_block_walk $D1 ]
&BEGIN_PROVIDER [ $T, $X_pdmc_block_walk_kahan $D2 ]
&BEGIN_PROVIDER [ $T, $X_2_pdmc_block_walk $D1 ]
&BEGIN_PROVIDER [ $T, $X_2_pdmc_block_walk_kahan $D2 ]
implicit none
BEGIN_DOC
! pdMC averages of $X. Computed in E_loc_pdmc_block_walk
END_DOC
$X_pdmc_block_walk = 0.d0
$X_pdmc_block_walk_kahan = 0.d0
$X_2_pdmc_block_walk = 0.d0
$X_2_pdmc_block_walk_kahan = 0.d0
END_PROVIDER
"""
for p in properties:
if p[1] != 'e_loc':
if p[2] == "":
D1 = ""
D2 = ", (3)"
else:
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_pdmc_block_walk ]
&BEGIN_PROVIDER [ double precision, E_loc_2_pdmc_block_walk ]
&BEGIN_PROVIDER [ double precision, E_loc_pdmc_block_walk_kahan, (3) ]
&BEGIN_PROVIDER [ double precision, E_loc_2_pdmc_block_walk_kahan, (3) ]
implicit none
include '../types.F'
BEGIN_DOC
! Properties averaged over the block using the PDMC method
END_DOC
real, allocatable :: elec_coord_tmp(:,:,:)
integer :: mod_align
double precision :: E_loc_save(walk_num_dmc_max)
double precision :: E_loc_save_tmp(walk_num_dmc_max)
double precision :: psi_value_save(walk_num)
double precision :: psi_value_save_tmp(walk_num)
double precision :: pdmc_weight(walk_num)
double precision, allocatable :: psi_grad_psi_inv_save(:,:,:)
double precision, allocatable :: psi_grad_psi_inv_save_tmp(:,:,:)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: psi_grad_psi_inv_save
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: psi_grad_psi_inv_save_tmp
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: E_loc_save
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: E_loc_save_tmp
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: psi_value_save
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: psi_value_save_tmp
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: pdmc_weight
allocate ( psi_grad_psi_inv_save(elec_num_8,3,walk_num) , &
psi_grad_psi_inv_save_tmp(elec_num_8,3,walk_num) , &
elec_coord_tmp(mod_align(elec_num+1),3,walk_num) )
psi_value_save = 0.d0
psi_value_save_tmp = 0.d0
pdmc_weight = 1.d0
! Initialization
if (vmc_algo /= t_Brownian) then
call abrt(irp_here,'PDMC should run with Brownian algorithm')
endif
integer :: k, i_walk, i_step
BEGIN_SHELL [ /usr/bin/python ]
from properties import *
t = """
if (calc_$X) then
!DIR$ VECTOR ALIGNED
$X_pdmc_block_walk = 0.d0
!DIR$ VECTOR ALIGNED
$X_pdmc_block_walk_kahan = 0.d0
!DIR$ VECTOR ALIGNED
$X_2_pdmc_block_walk = 0.d0
!DIR$ VECTOR ALIGNED
$X_2_pdmc_block_walk_kahan = 0.d0
endif
"""
for p in properties:
print t.replace("$X",p[1])
END_SHELL
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
real, external :: accep_rate
double precision :: delta, thr
thr = 2.d0/time_step_sq
logical :: first_loop
first_loop = .True.
if (walk_num > 1) then
call abrt(irp_here,'walk_num > 1')
endif
do while (loop)
! Every walker makes a step
do i_walk=1,walk_num
if (.not.first_loop) then
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
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
psi_value = psi_value_save(i_walk)
E_loc = E_loc_save(i_walk)
enddo
SOFT_TOUCH elec_coord psi_grad_psi_inv_x psi_grad_psi_inv_y psi_grad_psi_inv_z psi_value E_loc
else
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
psi_value_save(i_walk) = psi_value
E_loc_save(i_walk) = E_loc
endif
double precision :: p,q
real :: delta_x
logical :: accepted
call brownian_step(p,q,accepted,delta_x)
if ( psi_value * psi_value_save(i_walk) >= 0.d0 ) then
delta = ((E_loc+E_loc_save(i_walk))*0.5d0 - E_ref) * p
if ( delta > thr ) then
pdmc_weight(i_walk) = dexp(-dtime_step*thr)
else if ( delta < -thr ) then
pdmc_weight(i_walk) = dexp(dtime_step*thr)
else
pdmc_weight(i_walk) = dexp(-dtime_step*delta)
endif
elec_coord(elec_num+1,1) += p*time_step
elec_coord(elec_num+1,2) = E_loc
elec_coord(elec_num+1,3) = pdmc_weight(i_walk) * pop_weight_mult
do l=1,3
do i=1,elec_num+1
elec_coord_full(i,l,i_walk) = elec_coord(i,l)
enddo
enddo
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
psi_value_save(i_walk) = psi_value
E_loc_save(i_walk) = E_loc
if (calc_E_loc_zv) then
if (dabs(pdmc_weight(i_walk)*pop_weight_mult) > 1.d-6) then
E_loc_zv = (E_ref-E_loc)/(pdmc_weight(i_walk)*pop_weight_mult)
else
E_loc_zv = 0.d0
endif
TOUCH E_loc_zv
endif
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_pdmc_block_walk += $X * pop_weight_mult * pdmc_weight(i_walk)
! $X_2_pdmc_block_walk += $X_2 * pop_weight_mult * pdmc_weight(i_walk)
! see http://en.wikipedia.org/wiki/Kahan_summation_algorithm
$X_pdmc_block_walk_kahan($D2 3) = $X * pop_weight_mult * pdmc_weight(i_walk) - $X_pdmc_block_walk_kahan($D2 1)
$X_pdmc_block_walk_kahan($D2 2) = $X_pdmc_block_walk $D1 + $X_pdmc_block_walk_kahan($D2 3)
$X_pdmc_block_walk_kahan($D2 1) = ($X_pdmc_block_walk_kahan($D2 2) - $X_pdmc_block_walk $D1 ) &
- $X_pdmc_block_walk_kahan($D2 3)
$X_pdmc_block_walk $D1 = $X_pdmc_block_walk_kahan($D2 2)
$X_2_pdmc_block_walk_kahan($D2 3) = $X_2 * pop_weight_mult * pdmc_weight(i_walk) - $X_2_pdmc_block_walk_kahan($D2 1)
$X_2_pdmc_block_walk_kahan($D2 2) = $X_2_pdmc_block_walk $D1 + $X_2_pdmc_block_walk_kahan($D2 3)
$X_2_pdmc_block_walk_kahan($D2 1) = ($X_2_pdmc_block_walk_kahan($D2 2) - $X_2_pdmc_block_walk $D1 ) &
- $X_2_pdmc_block_walk_kahan($D2 3)
$X_2_pdmc_block_walk $D1 = $X_2_pdmc_block_walk_kahan($D2 2)
endif
"""
for p in properties:
if p[2] == "":
D1 = ""
D2 = ""
else:
D1 = "("+":"*(p[2].count(',')+1)+")"
D2 = ":"*(p[2].count(',')+1)+","
print t.replace("$X",p[1]).replace("$D1",D1).replace("$D2",D2)
END_SHELL
block_weight += pop_weight_mult * pdmc_weight(i_walk)
else
pdmc_weight(i_walk) = 0.d0
endif
enddo
! Move to the next projection step
if (pdmc_projection > 0) then
pdmc_projection_step = mod(pdmc_projection_step,pdmc_projection)+1
else
pdmc_projection_step = 1
endif
! Eventually, recompute the weight of the population
if (pdmc_projection_step == 1) then
pop_weight_mult = 1.d0
do k=1,pdmc_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(pdmc_projection_step)
! Compute the new weight of the population
double precision :: sum_weight
sum_weight = 0.d0
do k=1,walk_num
sum_weight += pdmc_weight(k)
enddo
pop_weight(pdmc_projection_step) = sum_weight/dble(walk_num)
! Update the running population weight
pop_weight_mult *= pop_weight(pdmc_projection_step)
call system_clock(cpu1, count_rate, count_max)
if (cpu1 < cpu0) then
cpu1 = cpu1+cpu0
endif
loop = dble(cpu1-cpu0) < 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
SOFT_TOUCH elec_coord_full pop_weight_mult
first_loop = .False.
enddo
double precision :: factor
factor = 1.d0/block_weight
SOFT_TOUCH block_weight
BEGIN_SHELL [ /usr/bin/python ]
from properties import *
t = """
if (calc_$X) then
$X_pdmc_block_walk *= factor
$X_2_pdmc_block_walk *= factor
endif
"""
for p in properties:
print t.replace("$X",p[1])
END_SHELL
deallocate ( elec_coord_tmp, psi_grad_psi_inv_save, psi_grad_psi_inv_save_tmp )
END_PROVIDER
BEGIN_PROVIDER [ double precision, pop_weight_mult ]
implicit none
BEGIN_DOC
! Population weight of PDMC
END_DOC
pop_weight_mult = pop_weight(pdmc_projection)
END_PROVIDER
BEGIN_PROVIDER [ integer, pdmc_projection ]
&BEGIN_PROVIDER [ integer, pdmc_projection_step ]
implicit none
BEGIN_DOC
! Number of projection steps for PDMC
END_DOC
real :: pdmc_projection_time
pdmc_projection_time = 1.
call get_simulation_srmc_projection_time(pdmc_projection_time)
pdmc_projection = int( pdmc_projection_time/time_step)
pdmc_projection_step = 0
END_PROVIDER
BEGIN_PROVIDER [ double precision, pop_weight, (0:pdmc_projection+1) ]
implicit none
BEGIN_DOC
! Population weight of PDMC
END_DOC
pop_weight(:) = 1.d0
END_PROVIDER

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

@ -15,7 +15,7 @@ subroutine reconfigure(ipos,w)
tmp = tmp + w(k)
enddo
dwalk_num = dble(walk_num)/tmp
integer :: kp, km
kp=0
km=0

1
src/SAMPLING/srmc_step.irp.f
View File

@ -169,6 +169,7 @@ END_SHELL
psi_value_save(i_walk) = psi_value
E_loc_save(i_walk) = E_loc
BEGIN_SHELL [ /usr/bin/python ]
from properties import *
t = """

4
src/simulation.irp.f
View File

@ -162,8 +162,10 @@ BEGIN_PROVIDER [ integer, qmc_method ]
qmc_method = t_SRMC
else if (method == types(t_FKMC)) then
qmc_method = t_FKMC
else if (method == types(t_PDMC)) then
qmc_method = t_PDMC
else
call abrt(irp_here, 'Method should be ( VMC | DMC | SRMC | FKMC )')
call abrt(irp_here, 'Method should be ( VMC | DMC | SRMC | FKMC | PDMC )')
endif
call cinfo(irp_here,'qmc_method',trim(method))

49
src/types.F
View File

@ -7,30 +7,41 @@
integer, parameter :: t_DMC = 8
integer, parameter :: t_SRMC = 9
integer, parameter :: t_FKMC = 10
integer, parameter :: t_PDMC = 11
integer, parameter :: t_Simple = 11
integer, parameter :: t_None = 12
integer, parameter :: t_Core = 14
integer, parameter :: t_Simple = 21
integer, parameter :: t_None = 22
integer, parameter :: t_Core = 24
integer, parameter :: t_Stopped = 0
integer, parameter :: t_Queued = 1
integer, parameter :: t_Running = 2
integer, parameter :: t_Stopping = 3
character*(32) :: types(15) = &
(/ ' ', &
' ', &
'Brownian ', &
'Langevin ', &
' ', &
' ', &
'VMC ', &
'DMC ', &
'SRMC ', &
'FKMC ', &
'Simple ', &
'None ', &
' ', &
'Core ', &
' '/)
character*(32) :: types(25) = &
(/ ' ', &
' ', &
'Brownian ', &
'Langevin ', &
' ', &
' ', &
'VMC ', &
'DMC ', &
'SRMC ', &
'FKMC ', &
'PDMC ', &
' ', &
' ', &
' ', &
' ', &
' ', &
' ', &
' ', &
' ', &
' ', &
'Simple ', &
'None ', &
' ', &
'Core ', &
' '/)