From 194b1f750c0203db6a3b50cca5d36a626ccc8c44 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 17 Mar 2016 15:41:47 +0100 Subject: [PATCH] Changed norm of MOs when the fitcusp is used --- ezfio_config/qmc.config | 3 +- ocaml/Default.ml | 22 +++++--- ocaml/Input.ml | 96 +++++++++++++++++------------------ ocaml/Md5.ml | 3 +- ocaml/Qmcchem_edit.ml | 40 +++++++-------- src/AO/ao_full.irp.f | 1 + src/ezfio_interface.irp.f | 3 +- src/mo.irp.f | 104 ++++++++++++++++++++++++++++++++++++++ src/nuclei.irp.f | 11 ++-- src/simulation.irp.f | 21 +++++--- 10 files changed, 208 insertions(+), 96 deletions(-) diff --git a/ezfio_config/qmc.config b/ezfio_config/qmc.config index f1392ce..29aaac9 100644 --- a/ezfio_config/qmc.config +++ b/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) diff --git a/ocaml/Default.ml b/ocaml/Default.ml index 37904a3..65b542d 100644 --- a/ocaml/Default.ml +++ b/ocaml/Default.ml @@ -1,23 +1,29 @@ open Core.Std;; -let simulation_do_nucl_fitcusp = lazy( - if (not (Ezfio.has_pseudo_do_pseudo ())) then - not (Ezfio.get_pseudo_do_pseudo ()) +let simulation_nucl_fitcusp_factor = lazy( + let default = + 0.5 + in + if (Ezfio.has_pseudo_do_pseudo ()) then + if (Ezfio.get_pseudo_do_pseudo ()) then + 0. + else + default else - true + default ) let electrons_elec_walk_num = lazy ( 30 ) let electrons_elec_walk_num_tot = lazy ( 10000 ) let jastrow_jast_type = lazy ( "None" ) let simulation_block_time = lazy ( 30 ) -let simulation_ci_threshold = lazy ( 1.e-8 ) +let simulation_ci_threshold = lazy ( 1.e-8 ) let simulation_method = lazy ( "VMC" ) let simulation_sampling = lazy ( "Langevin" ) let simulation_stop_time = lazy ( 3600 ) let simulation_time_step = lazy ( 0.15 ) -let simulation_srmc_projection_time = lazy ( 1. ) +let simulation_srmc_projection_time = lazy ( 1. ) let reset_defaults () = List.iter ~f:(fun x -> Sys.remove ( (Lazy.force Qputils.ezfio_filename) ^ x)) @@ -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" ] diff --git a/ocaml/Input.ml b/ocaml/Input.ml index a61444b..ad9c770 100644 --- a/ocaml/Input.ml +++ b/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 of_bool : bool -> t - val to_int : t -> int - val of_int : int -> t + val to_float : t -> float + val of_float : float -> 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 _ = - Lazy.force Qputils.ezfio_filename - in - if (not (Ezfio.has_simulation_do_nucl_fitcusp ())) then - Lazy.force Default.simulation_do_nucl_fitcusp - |> Ezfio.set_simulation_do_nucl_fitcusp ; - Ezfio.get_simulation_do_nucl_fitcusp () - |> of_bool + ignore @@ + 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 + Ezfio.set_simulation_nucl_fitcusp_factor factor + end ; + Ezfio.get_simulation_nucl_fitcusp_factor () + |> of_float let write t = let _ = Lazy.force Qputils.ezfio_filename in - let () = - match (Pseudo.read () |> Pseudo.to_bool, to_bool t) with - | (true, true) -> failwith "Pseudopotentials and Fitcusp are incompatible" - | _ -> () - in - to_bool t - |> Ezfio.set_simulation_do_nucl_fitcusp + to_float t + |> Ezfio.set_simulation_nucl_fitcusp_factor let to_string t = - to_bool t - |> Bool.to_string + to_float t + |> Float.to_string let of_string t = try - String.lowercase t - |> Bool.of_string - |> of_bool + Float.of_string t + |> of_float 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 - | (true, true) -> warn "Fitcusp is incompatible with Pseudopotentials" + let f = + 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 _ = diff --git a/ocaml/Md5.ml b/ocaml/Md5.ml index 8abd5ca..48d0b9b 100644 --- a/ocaml/Md5.ml +++ b/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" ; diff --git a/ocaml/Qmcchem_edit.ml b/ocaml/Qmcchem_edit.ml index 9b75263..2f6a6e7 100644 --- a/ocaml/Qmcchem_edit.ml +++ b/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 -> - option_to_string Fitcusp.read Fitcusp.to_string Fitcusp.doc + | Fitcusp_factor -> + 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 -> @@ -130,19 +130,19 @@ let run ~c ?f ?t ?l ?m ?e ?s ?ts ?w ?wt ?n ?j ?p ?input ezfio_filename = in (); in - handle_option Input.Ref_energy.(of_float , write) e; - handle_option Input.Jastrow_type.(of_string, write) j; - handle_option Input.Block_time.(of_int , write) l; - 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.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; - handle_option Input.CI_threshold.(of_float , write) n; - handle_option Input.SRMC_projection_time.(of_float , write) p; - + handle_option Input.Ref_energy.(of_float , write) e; + handle_option Input.Jastrow_type.(of_string, write) j; + handle_option Input.Block_time.(of_int , write) l; + 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_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; + handle_option Input.CI_threshold.(of_float , write) n; + handle_option Input.SRMC_projection_time.(of_float , write) p; + let fields = [ @@ -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) - ~doc:("0|1 "^Input.Fitcusp.doc) + +> flag "f" (optional float) + ~doc:("float "^Input.Fitcusp_factor.doc) +> flag "t" (optional int) ~doc:("seconds "^Input.Stop_time.doc) +> flag "l" (optional int) diff --git a/src/AO/ao_full.irp.f b/src/AO/ao_full.irp.f index 3ff16ba..c59abb7 100644 --- a/src/AO/ao_full.irp.f +++ b/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) diff --git a/src/ezfio_interface.irp.f b/src/ezfio_interface.irp.f index 32e4f6d..08c9cfe 100644 --- a/src/ezfio_interface.irp.f +++ b/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)", "" ), diff --git a/src/mo.irp.f b/src/mo.irp.f index 233137a..2dfa9c6 100644 --- a/src/mo.irp.f +++ b/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 diff --git a/src/nuclei.irp.f b/src/nuclei.irp.f index 39d84a5..cdcbc8b 100644 --- a/src/nuclei.irp.f +++ b/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 diff --git a/src/simulation.irp.f b/src/simulation.irp.f index dfdb084..3608e08 100644 --- a/src/simulation.irp.f +++ b/src/simulation.irp.f @@ -250,16 +250,21 @@ BEGIN_PROVIDER [ character*(64), hostname] END_PROVIDER -BEGIN_PROVIDER [ logical, do_nucl_fitcusp ] - implicit none - BEGIN_DOC - ! If true, do the fit of the electron-nucleus cusp - END_DOC - do_nucl_fitcusp = .True. - call get_simulation_do_nucl_fitcusp(do_nucl_fitcusp) - call linfo(irp_here,'do_nucl_fitcusp',do_nucl_fitcusp) + BEGIN_PROVIDER [ real, nucl_fitcusp_factor ] +&BEGIN_PROVIDER [ logical, do_nucl_fitcusp ] + implicit none + BEGIN_DOC + ! 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 + nucl_fitcusp_factor = 0. + call get_simulation_nucl_fitcusp_factor(nucl_fitcusp_factor) + 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