Changed norm of MOs when the fitcusp is used

ezfio_config/qmc.config

@@ -33,7 +33,6 @@ nuclei
   nucl_label          character*(32) (nuclei_nucl_num)
   nucl_charge         real           (nuclei_nucl_num)
   nucl_coord          real           (nuclei_nucl_num,3)
-  nucl_fitcusp_radius real           (nuclei_nucl_num)
 
 spindeterminants
   n_det_alpha     integer
@@ -55,7 +54,7 @@ simulation
   equilibration       logical
   http_server         character*(128)
   do_jast             logical
-  do_nucl_fitcusp     logical
+  nucl_fitcusp_factor real           
   method              character*(32)
   block_time          integer
   sampling            character*(32)

ocaml/Default.ml

@@ -1,11 +1,17 @@
 open Core.Std;;
 
 
-let simulation_do_nucl_fitcusp  = lazy(
+let simulation_nucl_fitcusp_factor  = lazy(
-  if (not (Ezfio.has_pseudo_do_pseudo ())) then
+  let default = 
-    not (Ezfio.get_pseudo_do_pseudo ())
+    0.5
+  in
+  if (Ezfio.has_pseudo_do_pseudo ()) then
+    if (Ezfio.get_pseudo_do_pseudo ()) then
+      0.
     else
-    true
+      default
+  else
+    default
 )
 
 let electrons_elec_walk_num         = lazy ( 30          )
@@ -26,9 +32,9 @@ let reset_defaults () =
       "/jastrow/jast_type"             ;
       "/simulation/block_time"         ;
       "/simulation/ci_threshold"       ;
-      "/simulation/do_nucl_fitcusp"    ;
       "/simulation/method"             ;
       "/simulation/sampling"           ;
       "/simulation/stop_time"          ;
-      "/simulation/time_step"          ]
+      "/simulation/time_step"          ;
+      "/simulation/nucl_fitcusp_factor" ]
 

ocaml/Input.ml

@@ -66,81 +66,69 @@ end = struct
 
 end
 
-module Fitcusp : sig
 
-  type t = bool
+
+module Fitcusp_factor : sig
+
+  type t = float
   val doc   : string
   val read  : unit -> t
   val write : t -> unit
-  val to_bool : t -> bool
+  val to_float : t -> float
-  val of_bool : bool -> t
+  val of_float : float -> t
-  val to_int  : t -> int
-  val of_int  : int -> t
   val to_string : t -> string
   val of_string : string -> t
 
 end = struct
 
-  type t = bool
+  type t = float
 
-  let doc = "Correct wave function to verify electron-nucleus cusp condition"
+  let doc = "Correct wave function to verify electron-nucleus cusp condition. 
+Fit is done for r < r_c(f) where r_c(f) = (1s orbital radius) x f. Value of f"
 
-  let of_bool x = x 
+  let of_float x = 
+    if (x < 0.) then
+       failwith "Fitcusp_factor should be >= 0.";
+    if (x > 10.) then
+       failwith "Fitcusp_factor is too large.";
+    x 
 
-  let to_bool x = x 
+  let to_float x = x 
 
   let read () = 
-    let _ =
+    ignore @@
-      Lazy.force Qputils.ezfio_filename
+      Lazy.force Qputils.ezfio_filename ;
+    if (not (Ezfio.has_simulation_nucl_fitcusp_factor ())) then
+      begin
+        let factor =
+          Lazy.force Default.simulation_nucl_fitcusp_factor ;
         in
-    if (not (Ezfio.has_simulation_do_nucl_fitcusp ())) then
+        Ezfio.set_simulation_nucl_fitcusp_factor factor
-      Lazy.force Default.simulation_do_nucl_fitcusp
+      end ;
-      |> Ezfio.set_simulation_do_nucl_fitcusp ;
+    Ezfio.get_simulation_nucl_fitcusp_factor () 
-    Ezfio.get_simulation_do_nucl_fitcusp ()
+    |> of_float 
-    |> of_bool
 
 
   let write t =
     let _ =
       Lazy.force Qputils.ezfio_filename
     in
-    let () =
+    to_float t
-      match (Pseudo.read () |> Pseudo.to_bool, to_bool t) with
+    |> Ezfio.set_simulation_nucl_fitcusp_factor 
-      | (true, true) -> failwith "Pseudopotentials and Fitcusp are incompatible"
-      | _ -> ()
-    in
-    to_bool t
-    |> Ezfio.set_simulation_do_nucl_fitcusp
 
 
   let to_string t =
-    to_bool t
+    to_float t
-    |> Bool.to_string
+    |> Float.to_string
 
 
   let of_string t =
     try
-      String.lowercase t
+      Float.of_string t
-      |> Bool.of_string 
+      |> of_float
-      |> of_bool
     with
     | Invalid_argument msg -> failwith msg
 
-
-  let to_int t =
-    let t = 
-      to_bool t
-    in
-    if t then 1
-    else 0
-
-
-  let of_int = function
-    | 0 -> false
-    | 1 -> true
-    | _ -> failwith "Expected 0 or 1"
-
-
 end
 
 module Block_time : sig
@@ -855,8 +843,6 @@ let validate () =
      Time_step.read ()
   and jast_type = 
     Jastrow_type.read ()
-  and do_fitcusp =
-    Fitcusp.read ()
   and do_pseudo =
     Pseudo.read ()
   in
@@ -915,13 +901,23 @@ let validate () =
     | _ -> ()
   in
 
-  (* Fitcusp is not recommended with pseudo *)
+  (* Fitcusp is incompatible with pseudo *)
   let () =
-    match (Pseudo.to_bool do_pseudo, Fitcusp.to_bool do_fitcusp) with
+    let f = 
-    | (true, true) -> warn "Fitcusp is incompatible with Pseudopotentials"
+      Fitcusp_factor.read ()
+      |> Fitcusp_factor.to_float
+    in
+    match (Pseudo.to_bool do_pseudo, f > 0.) with
+    | (true, true) ->
+      begin
+        warn "Electron-nucleus cusp fitting is incompatible with Pseudopotentials.";
+        Fitcusp_factor.of_float 0.
+        |> Fitcusp_factor.write
+      end
     | _ -> ()
   in
 
+  
   (* Other Checks *)
   let () =
     let _ =

ocaml/Md5.ml

@@ -37,10 +37,9 @@ let files_to_track = [
   "mo_basis/mo_tot_num" ;
   "nuclei/nucl_charge.gz" ;
   "nuclei/nucl_coord.gz" ;
-  "nuclei/nucl_fitcusp_radius.gz" ;
   "nuclei/nucl_num" ;
   "simulation/ci_threshold" ;
-  "simulation/do_nucl_fitcusp" ;
+  "simulation/nucl_fitcusp_factor" ;
   "simulation/jast_a_up_dn" ;
   "simulation/jast_a_up_up" ;
   "simulation/jast_b_up_dn" ;

ocaml/Qmcchem_edit.ml

@@ -20,7 +20,7 @@ type field =
  | Walk_num
  | Walk_num_tot
  | Stop_time
- | Fitcusp
+ | Fitcusp_factor
  | Method
  | Sampling
  | Ref_energy
@@ -54,8 +54,8 @@ let get field =
    option_to_string Walk_num_tot.read Walk_num_tot.to_string Walk_num_tot.doc  
  | Stop_time    ->
    option_to_string Stop_time.read    Stop_time.to_string    Stop_time.doc  
- | Fitcusp      ->
+ | Fitcusp_factor ->
-   option_to_string Fitcusp.read      Fitcusp.to_string      Fitcusp.doc
+   option_to_string Fitcusp_factor.read      Fitcusp_factor.to_string      Fitcusp_factor.doc
  | Method       ->
    option_to_string Method.read       Method.to_string       Method.doc  
  | Sampling     ->
@@ -136,7 +136,7 @@ let run ~c ?f ?t ?l ?m ?e ?s ?ts ?w ?wt ?n ?j ?p ?input ezfio_filename =
   handle_option         Input.Method.(of_string, write) m;
   handle_option      Input.Stop_time.(of_int   , write) t;
   handle_option       Input.Sampling.(of_string, write) s;
-  handle_option      Input.Fitcusp.(of_int   , write) f;
+  handle_option Input.Fitcusp_factor.(of_float , write) f;
   handle_option      Input.Time_step.(of_float , write) ts;
   handle_option       Input.Walk_num.(of_int   , write) w;
   handle_option   Input.Walk_num_tot.(of_int   , write) wt;
@@ -155,7 +155,7 @@ let run ~c ?f ?t ?l ?m ?e ?s ?ts ?w ?wt ?n ?j ?p ?input ezfio_filename =
       Ref_energy           ;
       Walk_num             ;
       Walk_num_tot         ;
-      Fitcusp              ;
+      Fitcusp_factor       ;
       CI_threshold         ;
       Jastrow_type         ;
       Properties           ;
@@ -214,7 +214,7 @@ let run ~c ?f ?t ?l ?m ?e ?s ?ts ?w ?wt ?n ?j ?p ?input ezfio_filename =
          begin
             match f with
             |  Stop_time             ->  Stop_time.(of_string            s  |>  write)
-            |  Fitcusp               ->  Fitcusp.(of_string              s  |>  write)
+            |  Fitcusp_factor        ->  Fitcusp_factor.(of_string       s  |>  write)
             |  Block_time            ->  Block_time.(of_string           s  |>  write)
             |  Method                ->  Method.(of_string               s  |>  write)
             |  Ref_energy            ->  Ref_energy.(of_string           s  |>  write)
@@ -271,8 +271,8 @@ let spec =
   empty
   +> flag "c"  no_arg
      ~doc:(" Clear blocks")
-  +> flag "f"  (optional int)
+  +> flag "f"  (optional float)
-     ~doc:("0|1 "^Input.Fitcusp.doc)
+     ~doc:("float "^Input.Fitcusp_factor.doc)
   +> flag "t"  (optional int)
      ~doc:("seconds "^Input.Stop_time.doc)
   +> flag "l"  (optional int)

src/AO/ao_full.irp.f

@@ -73,6 +73,7 @@ BEGIN_PROVIDER [ logical, primitives_reduced ]
       PROVIDE  ao_power
       PROVIDE  ao_coef
       PROVIDE  ao_nucl
+      PROVIDE  mo_fitcusp_normalization_before
       do i=1,ao_num
         if (ao_oned_p(i) /= 0.) then
           l=ao_power(i,1)+ao_power(i,2)+ao_power(i,3)

src/ezfio_interface.irp.f

@@ -6,7 +6,6 @@ data = [ \
 ("nuclei_nucl_num"                 , "integer"       , ""                      ),
 ("nuclei_nucl_charge"              , "real"          , "(nucl_num)"            ),
 ("nuclei_nucl_coord"               , "real"          , "(nucl_num,3)"          ),
-("nuclei_nucl_fitcusp_radius"      , "real"          , "(nucl_num)"            ),
 ("mo_basis_mo_coef"                , "real"          , "(ao_num,mo_tot_num)" ),
 ("electrons_elec_fitcusp_radius"   , "real"          , ""                      ),
 ("electrons_elec_alpha_num"        , "integer"       , ""                      ),
@@ -38,9 +37,9 @@ data = [ \
 ("simulation_time_step"            , "real"          , ""                      ),
 ("simulation_srmc_projection_time" , "real"          , ""                      ),
 ("simulation_method"               , "character*(32)", ""                      ),
+("simulation_nucl_fitcusp_factor"  , "real"          , ""                     ),
 ("simulation_save_data"            , "logical"       , ""                      ),
 ("simulation_print_level"          , "integer"       , ""                      ),
-("simulation_do_nucl_fitcusp"      , "logical"       , ""                      ),
 ("simulation_sampling"             , "character*(32)", ""                      ),
 ("simulation_ci_threshold"         , "double precision"          , ""                      ),
 ("simulation_http_server"          , "character*(128)", ""                      ),

src/mo.irp.f

@@ -47,6 +47,7 @@ END_PROVIDER
   END_DOC
   mo_scale = 1.d0/(0.4d0*log(float(elec_num+1)))
   mo_norm = mo_scale*mo_scale
+
 END_PROVIDER
 
 
@@ -273,6 +274,15 @@ END_PROVIDER
     enddo
   endif
   
+  do i=1,mo_num
+    do j=1,elec_num
+        mo_value_transp(i,j)  *= mo_cusp_rescale(i)
+        mo_grad_transp_x(i,j) *= mo_cusp_rescale(i)
+        mo_grad_transp_y(i,j) *= mo_cusp_rescale(i)
+        mo_grad_transp_z(i,j) *= mo_cusp_rescale(i)
+        mo_lapl_transp(i,j)   *= mo_cusp_rescale(i)
+    enddo
+  enddo
 END_PROVIDER
 
 
@@ -401,6 +411,7 @@ BEGIN_PROVIDER [ double precision , mo_value_at_nucl, (mo_num_8,nucl_num) ]
   integer                        :: i, j, k, l
   real                           :: ao_value_at_nucl_no_S(ao_num)
   
+  PROVIDE mo_fitcusp_normalization_before
   do k=1,nucl_num
     point(1) = nucl_coord(k,1)
     point(2) = nucl_coord(k,2)
@@ -466,6 +477,99 @@ END_PROVIDER
   FREE ao_value_p ao_grad_p ao_lapl_p ao_axis_grad_p ao_oned_grad_p ao_oned_prim_grad_p ao_oned_lapl_p ao_axis_lapl_p ao_oned_prim_lapl_p ao_oned_p ao_oned_prim_p ao_axis_p ao_axis_power_p
   SOFT_TOUCH point
   
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, mo_fitcusp_normalization_before, (mo_tot_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Renormalization factor of MOs due to cusp fitting
+ END_DOC
+ include 'constants.F'
+ integer :: i,j,k,l
+ double precision :: dr, r, f, t
+ integer, save :: ifirst = 0
+
+  if (ifirst == 0) then
+    ifirst = 1
+    mo_fitcusp_normalization_before = 0.d0
+    do k=1,nucl_num
+      dr = nucl_fitcusp_radius(k)*1.d-2
+      point(1) = nucl_coord(k,1)
+      point(2) = nucl_coord(k,2)
+      point(3) = nucl_coord(k,3)-dr
+      do l=1,101
+        r = point(3) + dr
+        point(3) = r
+        TOUCH point
+        f = dfour_pi*r*r*dr
+        do i=1,mo_tot_num
+          mo_fitcusp_normalization_before(i) += f*mo_value_p(i)**2
+        enddo
+      enddo
+    enddo
+  endif
+
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, mo_fitcusp_normalization_after, (mo_tot_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Renormalization factor of MOs due to cusp fitting
+ END_DOC
+ include 'constants.F'
+ integer :: i,j,k,l
+ double precision :: dr, r, f, t, t2
+ integer, save :: ifirst = 0
+
+  PROVIDE primitives_reduced
+  if (ifirst == 0) then
+    ifirst = 1
+    mo_fitcusp_normalization_after = 0.d0
+    do k=1,nucl_num
+      dr = nucl_fitcusp_radius(k)*1.d-2
+      point(1) = nucl_coord(k,1)
+      point(2) = nucl_coord(k,2)
+      point(3) = nucl_coord(k,3)- dr
+      do l=1,101
+        point(3) = point(3)+ dr
+        TOUCH point nucl_fitcusp_param primitives_reduced mo_coef
+        r = point(3)
+        f = dfour_pi*r*r*dr
+        do i=1,mo_tot_num
+          t = 0.d0
+          do j=1,ao_num
+            if ( (ao_nucl(j) /= k).or.(ao_power(j,4) > 0) ) then
+              t = t + mo_coef(j,i) * ao_value_p(j) 
+            endif
+          enddo
+          t = t + nucl_fitcusp_param(1,i,k) +                        &
+              r * (nucl_fitcusp_param(2,i,k) +                       &
+              r * (nucl_fitcusp_param(3,i,k) +                       &
+              r *  nucl_fitcusp_param(4,i,k) ))
+          mo_fitcusp_normalization_after(i) += t*t*f
+        enddo
+      enddo
+    enddo
+  endif
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ real, mo_cusp_rescale, (mo_tot_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Rescaling coefficient to normalize MOs after applying fitcusp
+ END_DOC
+ integer :: i
+ if (do_nucl_fitcusp) then
+   do i=1,mo_tot_num
+!     mo_cusp_rescale(i) = dsqrt(mo_fitcusp_normalization_before(i) / mo_fitcusp_normalization_after(i))
+     mo_cusp_rescale(i) = 1.d0/dsqrt(1.d0 - mo_fitcusp_normalization_before(i) + mo_fitcusp_normalization_after(i))
+   enddo
+ else
+     mo_cusp_rescale = 1.d0
+ endif
+
 END_PROVIDER
 
 

src/nuclei.irp.f

@@ -125,17 +125,20 @@ BEGIN_PROVIDER [ real, nucl_fitcusp_radius, (nucl_num) ]
   BEGIN_DOC
 ! Distance threshold for the fit
   END_DOC
-  real                           :: def(nucl_num)
+  real                           :: def(nucl_num), factor
   integer                        :: k
+  real, parameter                :: a = 1.74891
+  real, parameter                :: b = 0.126057
   
-  if (.not.do_nucl_fitcusp) then
+  if (.not. do_nucl_fitcusp) then
     nucl_fitcusp_radius = 0.d0
     return
   endif
+
   do k=1,nucl_num
-    nucl_fitcusp_radius(k) = .5/nucl_charge(k)
+    nucl_fitcusp_radius(k) = nucl_fitcusp_factor/(a*nucl_charge(k)+b)
   enddo
-  call get_nuclei_nucl_fitcusp_radius(nucl_fitcusp_radius)
+
   ! Avoid dummy atoms
   do k=1,nucl_num
     if (nucl_charge(k) < 5.d-1) then

src/simulation.irp.f

@@ -250,17 +250,22 @@ BEGIN_PROVIDER [ character*(64), hostname]
 END_PROVIDER
 
 
-BEGIN_PROVIDER [ logical, do_nucl_fitcusp ]
+ BEGIN_PROVIDER [ real, nucl_fitcusp_factor ]
+&BEGIN_PROVIDER [ logical, do_nucl_fitcusp ]
  implicit none
  BEGIN_DOC
-   ! If true, do the fit of the electron-nucleus cusp
+ ! The electron-nucleus cusp fitting is done between 0 and r_c,
+ ! where r_c is chosen as nucl_fitcusp_factor * (radius_of_1s AO)
  END_DOC
-   do_nucl_fitcusp = .True.
+ nucl_fitcusp_factor = 0.
-   call get_simulation_do_nucl_fitcusp(do_nucl_fitcusp)
+ call get_simulation_nucl_fitcusp_factor(nucl_fitcusp_factor)
-   call linfo(irp_here,'do_nucl_fitcusp',do_nucl_fitcusp)
+ do_nucl_fitcusp = nucl_fitcusp_factor > 0.
+ call info(irp_here,'nucl_fitcusp_factor',nucl_fitcusp_factor)
+
 END_PROVIDER
 
 
+ 
 BEGIN_PROVIDER [ integer, vmc_algo ]
    implicit none
    include 'types.F'