Added Rousset FKMC algorithm

ocaml/Default.ml

@@ -3,7 +3,7 @@ open Core.Std;;
 
 let simulation_nucl_fitcusp_factor  = lazy(
   let default = 
-    0.5
+    1.
   in
   if (Ezfio.has_pseudo_do_pseudo ()) then
     if (Ezfio.get_pseudo_do_pseudo ()) then
@@ -14,8 +14,8 @@ let simulation_nucl_fitcusp_factor  = lazy(
     default
 )
 
-let electrons_elec_walk_num         = lazy ( 30          )
+let electrons_elec_walk_num         = lazy ( 100         )
-let electrons_elec_walk_num_tot     = lazy ( 10000       )
+let electrons_elec_walk_num_tot     = lazy ( 1000        )
 let jastrow_jast_type               = lazy ( "None"      )
 let simulation_block_time           = lazy ( 30          )
 let simulation_ci_threshold         = lazy ( 1.e-8       )

ocaml/Input.ml

@@ -387,7 +387,7 @@ end
 
 module Method : sig
 
-  type t = VMC | DMC | SRMC
+  type t = VMC | DMC | SRMC | FKMC
   val doc : string
   val read  : unit -> t
   val write : t -> unit
@@ -396,21 +396,23 @@ module Method : sig
 
 end = struct
 
-  type t = VMC | DMC | SRMC
+  type t = VMC | DMC | SRMC | FKMC
 
-  let doc = "QMC Method : [ VMC | DMC | SRMC ]"
+  let doc = "QMC Method : [ VMC | DMC | SRMC | FKMC ]"
 
   let of_string = function
   | "VMC"  | "vmc"  -> VMC
   | "DMC"  | "dmc"  -> DMC
   | "SRMC" | "srmc" -> SRMC
-  | x -> failwith ("Method should be [ VMC | DMC | SRMC ], not "^x^".")
+  | "FKMC" | "fkmc" -> FKMC
+  | x -> failwith ("Method should be [ VMC | DMC | SRMC | FKMC ], not "^x^".")
 
 
   let to_string = function
   | VMC  -> "VMC"
   | DMC  -> "DMC"
   | SRMC -> "SRMC"
+  | FKMC -> "FKMC"
 
 
   let read () = 
@@ -851,10 +853,12 @@ let validate () =
   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) )
@@ -865,8 +869,9 @@ let validate () =
       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, _) ->
-        failwith "Lanvegin sampling is incompatible with DMC/SRMC"
+        failwith "Lanvegin sampling is incompatible with DMC"
   in
 
 
@@ -874,6 +879,7 @@ let validate () =
   let () =
     match (meth, Ref_energy.(read () |> to_float) ) with
     | (Method.SRMC,0.) 
+    | (Method.FKMC,0.) 
     | (Method.DMC,0.) -> failwith ("E_ref should not be zero in "^(Method.to_string meth) )
     | _          -> ()
   in
@@ -888,6 +894,7 @@ let validate () =
   let () =
     match (meth, Property.(calc E_loc)) with
     | (Method.SRMC, 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"
     | _ -> ()

src/SAMPLING/block.irp.f

@@ -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_FKMC) then
+  PROVIDE E_loc_fkmc_block_walk
+  if (calc_$X) then
+    $X_block_walk = $X_fkmc_block_walk
+    $X_2_block_walk = $X_2_fkmc_block_walk
+  endif
  endif
 
 END_PROVIDER

src/SAMPLING/dmc_step.irp.f

@@ -195,7 +195,7 @@ END_SHELL
   do k=1,walk_num_dmc
     sum_weight += dmc_weight(k)
   enddo
-  E0 = E_ref - log(real(walk_num_dmc)/real(walk_num)) * 0.1d0/dtime_step
+  E0 = E_ref - log(sum_weight/real(walk_num)) * 0.1d0 /dtime_step
 
 ! Branching
   integer                        :: ipos(walk_num_dmc_max), walk_num_dmc_new
@@ -214,20 +214,22 @@ END_SHELL
     ipos(k) = k
   enddo
 
+  walk_num_dmc_new = walk_num_dmc
   do k=1,walk_num_dmc
     r = qmc_ranf()
     if (dmc_weight(k) > 1.d0) then
       if ( 1.d0+r < dmc_weight(k) ) then
-        walk_num_dmc = walk_num_dmc+1
+        walk_num_dmc_new = walk_num_dmc_new+1
-        ipos(walk_num_dmc) = k
+        ipos(walk_num_dmc_new) = k
       endif
     else
       if ( r > dmc_weight(k) ) then
-        ipos(k) = ipos(walk_num_dmc)
+        ipos(k) = ipos(walk_num_dmc_new)
-        walk_num_dmc = walk_num_dmc-1
+        walk_num_dmc_new = walk_num_dmc_new-1
       endif
     endif
   enddo
+  walk_num_dmc = walk_num_dmc_new
 
   integer :: ipm
   do k=1,walk_num_dmc
@@ -328,7 +330,7 @@ BEGIN_PROVIDER [ integer, walk_num_dmc_max ]
  BEGIN_DOC
  ! Max number of walkers in DMC
  END_DOC
- walk_num_dmc_max = 3 * walk_num
+ walk_num_dmc_max = max(3 * walk_num, 30)
 END_PROVIDER
 
 

src/SAMPLING/fkmc_step.irp.f

@@ -0,0 +1,375 @@
+! Providers of *_fkmc_block_walk
+!==============================
+BEGIN_SHELL [ /usr/bin/python ]
+from properties import *
+
+t = """
+ BEGIN_PROVIDER [ $T, $X_fkmc_block_walk $D1 ]
+&BEGIN_PROVIDER [ $T, $X_fkmc_block_walk_kahan $D2 ]
+&BEGIN_PROVIDER [ $T, $X_2_fkmc_block_walk $D1 ]
+&BEGIN_PROVIDER [ $T, $X_2_fkmc_block_walk_kahan $D2 ]
+ implicit none
+ BEGIN_DOC  
+! fkMC averages of $X. Computed in E_loc_fkmc_block_walk
+ END_DOC
+ $X_fkmc_block_walk = 0.d0
+ $X_fkmc_block_walk_kahan = 0.d0
+ $X_2_fkmc_block_walk = 0.d0
+ $X_2_fkmc_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_fkmc_block_walk       ]
+&BEGIN_PROVIDER [ double precision, E_loc_2_fkmc_block_walk     ]
+&BEGIN_PROVIDER [ double precision, E_loc_fkmc_block_walk_kahan, (3) ]
+&BEGIN_PROVIDER [ double precision, E_loc_2_fkmc_block_walk_kahan, (3) ]
+ implicit none
+ include '../types.F'
+ BEGIN_DOC
+! Properties averaged over the block using the FKMC method
+ END_DOC
+
+  integer, parameter :: BIRTH=1, DEATH=2
+  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 :: fkmc_weight(walk_num)
+  double precision :: delta(walk_num)
+  double precision, allocatable :: psi_grad_psi_inv_save(:,:,:)
+  double precision, allocatable :: psi_grad_psi_inv_save_tmp(:,:,:)
+  double precision, allocatable :: fkmc_clock_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 :: fkmc_weight
+  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: fkmc_clock_tmp
+  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),             &
+       fkmc_clock_tmp(2,walk_num) )
+  psi_value_save = 0.d0
+  psi_value_save_tmp = 0.d0
+  fkmc_weight = 1.d0
+
+! Initialization
+ if (vmc_algo /= t_Brownian) then
+   call abrt(irp_here,'FKMC 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_fkmc_block_walk = 0.d0
+   !DIR$ VECTOR ALIGNED
+   $X_fkmc_block_walk_kahan = 0.d0
+   !DIR$ VECTOR ALIGNED
+   $X_2_fkmc_block_walk = 0.d0
+   !DIR$ VECTOR ALIGNED
+   $X_2_fkmc_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               :: thr
+
+ thr = 2.d0/time_step_sq
+
+ logical :: first_loop
+ first_loop = .True.
+
+ 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
+      E_loc_save(i_walk) = E_loc 
+      psi_value_save(i_walk) = psi_value
+    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(i_walk) = ((E_loc+E_loc_save(i_walk))*0.5d0 - E_ref) * p
+     if ( delta(i_walk) > thr ) then
+       delta(i_walk) = thr
+     else if ( delta(i_walk) < -thr ) then
+       delta(i_walk) = -thr
+     endif
+     fkmc_weight(i_walk) = dexp(-dtime_step*delta(i_walk))
+     elec_coord(elec_num+1,1) += p*time_step
+     elec_coord(elec_num+1,2)  = E_loc
+     elec_coord(elec_num+1,3)  = fkmc_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
+     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
+
+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_fkmc_block_walk    += $X * fkmc_weight(i_walk)
+   !  $X_2_fkmc_block_walk  += $X_2 * fkmc_weight(i_walk)
+   ! see http://en.wikipedia.org/wiki/Kahan_summation_algorithm
+   
+      $X_fkmc_block_walk_kahan($D2 3) = $X * fkmc_weight(i_walk) - $X_fkmc_block_walk_kahan($D2 1)
+      $X_fkmc_block_walk_kahan($D2 2) = $X_fkmc_block_walk $D1  + $X_fkmc_block_walk_kahan($D2 3)
+      $X_fkmc_block_walk_kahan($D2 1) = ($X_fkmc_block_walk_kahan($D2 2) - $X_fkmc_block_walk $D1 ) &
+          - $X_fkmc_block_walk_kahan($D2 3)
+      $X_fkmc_block_walk $D1  =  $X_fkmc_block_walk_kahan($D2 2) 
+   
+   
+      $X_2_fkmc_block_walk_kahan($D2 3) = $X_2 * fkmc_weight(i_walk) - $X_2_fkmc_block_walk_kahan($D2 1)
+      $X_2_fkmc_block_walk_kahan($D2 2) = $X_2_fkmc_block_walk $D1 + $X_2_fkmc_block_walk_kahan($D2 3)
+      $X_2_fkmc_block_walk_kahan($D2 1) = ($X_2_fkmc_block_walk_kahan($D2 2) - $X_2_fkmc_block_walk $D1 ) &
+          - $X_2_fkmc_block_walk_kahan($D2 3)
+      $X_2_fkmc_block_walk $D1 =  $X_2_fkmc_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 += fkmc_weight(i_walk)
+
+   else
+     fkmc_weight(i_walk) = 0.d0
+     delta(i_walk) = 1.d5
+   endif
+
+  enddo
+
+  ! Compute the new weight of the population
+  double precision :: sum_weight
+  sum_weight = 0.d0
+  do k=1,walk_num
+    sum_weight += fkmc_weight(k)
+  enddo
+
+  do k=1,walk_num
+    do l=1,3
+     do i=1,elec_num+1
+      elec_coord_tmp(i,l,k) = elec_coord_full(i,l,k)
+     enddo
+     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)
+    E_loc_save_tmp(k) = E_loc_save(k)
+    if (fkmc_weight(k) == 0.d0) then
+      fkmc_clock(DEATH,k) = -1.d0
+    endif
+    if ( delta(k) <= 0.d0 ) then
+      fkmc_clock_tmp(BIRTH,k) = fkmc_clock(BIRTH,k) +time_step * delta(k)
+      fkmc_clock_tmp(DEATH,k) = fkmc_clock(DEATH,k)
+    else
+      fkmc_clock_tmp(BIRTH,k) = fkmc_clock(BIRTH,k)
+      fkmc_clock_tmp(DEATH,k) = fkmc_clock(DEATH,k) -time_step * delta(k)
+    endif
+  enddo
+
+! Reconfiguration
+! ===============
+
+  ! Identify first which walkers will be killed to place branched walkers there
+  ! later
+  
+  double precision, external     :: qmc_ranf
+  integer                        :: ipm, m
+  integer                        :: killed(walk_num)
+
+  m=1
+  do k=1,walk_num
+    fkmc_clock(DEATH,k) = fkmc_clock_tmp(DEATH,k) 
+    if (fkmc_clock_tmp(DEATH,k) <= 0.d0) then
+      killed(m) = k
+      m += 1
+      fkmc_clock(DEATH,k) = -dlog(qmc_ranf())
+      fkmc_clock(BIRTH,k) = -dlog(qmc_ranf())
+      ipm = k
+      do while (ipm == k)
+        ipm = 1 + int (walk_num*qmc_ranf())
+      enddo
+      do l=1,3
+        do i=1,elec_num+1
+        elec_coord_full(i,l,k) = elec_coord_tmp(i,l,ipm)
+        enddo
+        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)
+      E_loc_save(k) = E_loc_save_tmp(ipm)
+    endif
+  enddo
+  killed(m) = 0
+
+  m=1
+  do k=1,walk_num
+    fkmc_clock(BIRTH,k) = fkmc_clock_tmp(BIRTH,k) 
+    if (fkmc_clock_tmp(BIRTH,k) <= 0.d0) then
+      fkmc_clock(BIRTH,k) = -dlog(qmc_ranf())
+      if (killed(m) == 0) then
+        ipm = k
+        do while (ipm == k)
+          ipm = 1 + int (walk_num*qmc_ranf())
+        enddo
+      else
+        ipm = killed(m)
+        m +=1
+      endif
+      fkmc_clock(BIRTH,ipm) = -dlog(qmc_ranf())
+      fkmc_clock(DEATH,ipm) = -dlog(qmc_ranf())
+      do l=1,3
+        do i=1,elec_num+1
+        elec_coord_full(i,l,ipm) = elec_coord_tmp(i,l,k)
+        enddo
+        do i=1,elec_num
+          psi_grad_psi_inv_save(i,l,ipm) = psi_grad_psi_inv_save_tmp(i,l,k)
+        enddo
+      enddo
+      psi_value_save(ipm) = psi_value_save_tmp(k)
+      E_loc_save(ipm) = E_loc_save_tmp(k)
+      
+    endif
+
+  enddo
+
+  
+
+  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
+
+! Update E_ref to take into account the weight of the population
+  E_ref -= dlog(sum_weight / dble(walk_num) ) / time_step 
+  SOFT_TOUCH elec_coord_full E_ref
+
+  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_fkmc_block_walk   *= factor
+   $X_2_fkmc_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, &
+   fkmc_clock_tmp )
+
+END_PROVIDER
+
+
+
+BEGIN_PROVIDER [ double precision, fkmc_clock, (2,walk_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Branching clocks for the FKMC algotithm. (1,:) is the birth clock and
+ ! (2,:) is the death clock.
+ END_DOC
+ integer :: i
+ double precision, external :: qmc_ranf
+ do i=1, walk_num
+   fkmc_clock(1,i) = -dlog(qmc_ranf())
+   fkmc_clock(2,i) = -dlog(qmc_ranf())
+ enddo
+
+END_PROVIDER
+

src/SAMPLING/srmc_step.irp.f

@@ -133,7 +133,9 @@ END_SHELL
           elec_coord(i,l) = elec_coord_full(i,l,i_walk)
         enddo
       enddo
-      SOFT_TOUCH elec_coord
+      TOUCH elec_coord
+      psi_value_save(i_walk) = psi_value
+      E_loc_save(i_walk) = E_loc
     endif
 
    double precision               :: p,q

src/mo.irp.f

@@ -47,7 +47,6 @@ END_PROVIDER
   END_DOC
   mo_scale = 1.d0/(0.4d0*log(float(elec_num+1)))
   mo_norm = mo_scale*mo_scale
-
 END_PROVIDER
 
 
@@ -536,7 +535,7 @@ BEGIN_PROVIDER [ double precision, mo_fitcusp_normalization_after, (mo_tot_num)
         TOUCH point nucl_fitcusp_param primitives_reduced mo_coef
         r = point(3)
         f = dfour_pi*r*r*dr
-        do i=1,mo_tot_num
+        do i=1,mo_num
           t = 0.d0
           do j=1,ao_num
             if ( (ao_nucl(j) /= k).or.(ao_power(j,4) > 0) ) then

src/simulation.irp.f

@@ -148,7 +148,7 @@ BEGIN_PROVIDER  [ integer, qmc_method ]
    implicit none
    include 'types.F'
    BEGIN_DOC
-   ! qmc_method : Calculation method. Can be t_VMC, t_DMC, t_SRMC
+   ! qmc_method : Calculation method. Can be t_VMC, t_DMC, t_SRMC, t_FKMC
    END_DOC
    character*(32)                 :: method
    method = types(t_VMC)
@@ -160,8 +160,10 @@ BEGIN_PROVIDER  [ integer, qmc_method ]
      qmc_method = t_DMC
    else if (method == types(t_SRMC)) then
      qmc_method = t_SRMC
+   else if (method == types(t_FKMC)) then
+     qmc_method = t_FKMC
    else
-     call abrt(irp_here, 'Method should be ( VMC | DMC | SRMC )')
+     call abrt(irp_here, 'Method should be ( VMC | DMC | SRMC | FKMC )')
    endif
    
    call cinfo(irp_here,'qmc_method',trim(method))
@@ -288,6 +290,9 @@ BEGIN_PROVIDER [ integer, vmc_algo ]
      if (qmc_method == t_SRMC) then
        stop 'Langevin incompatible with SRMC'
      endif
+     if (qmc_method == t_FKMC) then
+       stop 'Langevin incompatible with FKMC'
+     endif
    else if (Sampling == types(t_MTM)) then
      vmc_algo = t_MTM
    else

src/types.F

@@ -6,6 +6,7 @@
       integer, parameter :: t_VMC      = 7
       integer, parameter :: t_DMC      = 8
       integer, parameter :: t_SRMC     = 9
+      integer, parameter :: t_FKMC     = 10
 
       integer, parameter :: t_Simple   = 11
       integer, parameter :: t_None     = 12
@@ -26,7 +27,7 @@
           'VMC            ',   &
           'DMC            ',   &
           'SRMC           ',   &
-          '               ',   &
+          'FKMC           ',   &
           'Simple         ',   &
           'None           ',   &
           '               ',   &