Working on normf

ocaml/Input_ao_basis.ml

@@ -13,6 +13,8 @@ module Ao_basis : sig
       ao_coef         : AO_coef.t array;
       ao_expo         : AO_expo.t array;
       ao_cartesian    : bool;
+      ao_normalized   : bool;
+      primitives_normalized : bool;
     } [@@deriving sexp]
   ;;
   val read : unit -> t option
@@ -34,6 +36,8 @@ end = struct
       ao_coef         : AO_coef.t array;
       ao_expo         : AO_expo.t array;
       ao_cartesian    : bool;
+      ao_normalized   : bool;
+      primitives_normalized : bool;
     } [@@deriving sexp]
   ;;
 
@@ -107,6 +111,24 @@ end = struct
     Ezfio.get_ao_basis_ao_cartesian ()
   ;;
 
+  let read_ao_normalized () =
+    if not (Ezfio.has_ao_basis_ao_normalized()) then
+       get_default "ao_normalized"
+       |> bool_of_string
+       |> Ezfio.set_ao_basis_ao_normalized
+    ;
+    Ezfio.get_ao_basis_ao_normalized ()
+  ;;
+
+  let read_primitives_normalized () =
+    if not (Ezfio.has_ao_basis_primitives_normalized()) then
+       get_default "primitives_normalized"
+       |> bool_of_string
+       |> Ezfio.set_ao_basis_primitives_normalized
+    ;
+    Ezfio.get_ao_basis_primitives_normalized ()
+  ;;
+
   let to_long_basis b =
     let ao_num = AO_number.to_int b.ao_num in
     let gto_array = Array.init (AO_number.to_int b.ao_num)
@@ -169,6 +191,8 @@ end = struct
          ao_coef         ;
          ao_expo         ;
          ao_cartesian    ;
+         ao_normalized   ;
+         primitives_normalized ; 
        } = b
      in
      write_md5 b ;
@@ -201,6 +225,8 @@ end = struct
      ~rank:2 ~dim:[| ao_num ; 3 |] ~data:ao_power) ;
 
      Ezfio.set_ao_basis_ao_cartesian(ao_cartesian);
+     Ezfio.set_ao_basis_ao_normalized(ao_normalized);
+     Ezfio.set_ao_basis_primitives_normalized(primitives_normalized);
      
      let ao_coef =
       Array.to_list ao_coef
@@ -233,6 +259,8 @@ end = struct
             ao_coef         = read_ao_coef () ;
             ao_expo         = read_ao_expo () ;
             ao_cartesian    = read_ao_cartesian () ;
+            ao_normalized   = read_ao_normalized () ;
+            primitives_normalized   = read_primitives_normalized () ;
           }
         in
         to_md5 result
@@ -340,7 +368,10 @@ end = struct
       in
       { ao_basis = name ;
         ao_num ; ao_prim_num ; ao_prim_num_max ; ao_nucl ;
-        ao_power ; ao_coef ; ao_expo ; ao_cartesian }
+        ao_power ; ao_coef ; ao_expo ; ao_cartesian ;
+        ao_normalized = bool_of_string @@ get_default "ao_normalized";
+        primitives_normalized = bool_of_string @@ get_default "primitives_normalized";
+        }
   ;;
 
   let reorder b = 
@@ -394,6 +425,14 @@ Cartesian coordinates (6d,10f,...) ::
 
   ao_cartesian = %s
 
+Use normalized primitive functions ::
+
+  primitives_normalized = %s
+
+Use normalized basis functions ::
+
+  ao_normalized = %s
+
 Basis set (read-only) ::
 
 %s
@@ -407,6 +446,8 @@ Basis set (read-only) ::
 
 "   (AO_basis_name.to_string b.ao_basis)
     (string_of_bool b.ao_cartesian)
+    (string_of_bool b.primitives_normalized)
+    (string_of_bool b.ao_normalized)
     (Basis.to_string short_basis
        |> String_ext.split ~on:'\n'
        |> List.map (fun x-> "  "^x)
@@ -434,16 +475,18 @@ Basis set (read-only) ::
 
   let to_string b =
     Printf.sprintf "
-ao_basis        = %s
+ao_basis                = %s
-ao_num          = %s
+ao_num                  = %s
-ao_prim_num     = %s
+ao_prim_num             = %s
-ao_prim_num_max = %s
+ao_prim_num_max         = %s
-ao_nucl         = %s
+ao_nucl                 = %s
-ao_power        = %s
+ao_power                = %s
-ao_coef         = %s
+ao_coef                 = %s
-ao_expo         = %s
+ao_expo                 = %s
-ao_cartesian    = %s
+ao_cartesian            = %s
-md5             = %s
+ao_normalized           = %s
+primitives_normalized   = %s
+md5                     = %s
 "
     (AO_basis_name.to_string b.ao_basis)
     (AO_number.to_string b.ao_num)
@@ -459,6 +502,8 @@ md5             = %s
     (b.ao_expo  |> Array.to_list |> List.map AO_expo.to_string
       |> String.concat ", ")
     (b.ao_cartesian |> string_of_bool)
+    (b.ao_normalized |> string_of_bool)
+    (b.primitives_normalized |> string_of_bool)
     (to_md5 b |> MD5.to_string )
 
   ;;

src/ao_basis/EZFIO.cfg

@@ -55,3 +55,15 @@ doc: If |true|, use |AOs| in Cartesian coordinates (6d,10f,...)
 interface: ezfio, provider
 default: false
 
+[ao_normalized]
+type: logical
+doc: Use normalized basis functions
+interface: ezfio, provider
+default: true
+
+[primitives_normalized]
+type: logical
+doc: Use normalized primitive functions
+interface: ezfio, provider
+default: true
+

src/ao_basis/aos.irp.f

@@ -20,25 +20,38 @@ END_PROVIDER
   C_A(2) = 0.d0
   C_A(3) = 0.d0
   ao_coef_normalized = 0.d0
+
   do i=1,ao_num
 
     powA(1) = ao_power(i,1)
     powA(2) = ao_power(i,2)
     powA(3) = ao_power(i,3)
 
-    do j=1,ao_prim_num(i)
+    ! Normalization of the primitives
-      call overlap_gaussian_xyz(C_A,C_A,ao_expo(i,j),ao_expo(i,j),powA,powA,overlap_x,overlap_y,overlap_z,norm,nz)
+    if (primitives_normalized) then
-      ao_coef_normalized(i,j) = ao_coef(i,j)/sqrt(norm)
+      do j=1,ao_prim_num(i)
-    enddo
+        call overlap_gaussian_xyz(C_A,C_A,ao_expo(i,j),ao_expo(i,j),powA,powA,overlap_x,overlap_y,overlap_z,norm,nz)
+        ao_coef_normalized(i,j) = ao_coef(i,j)/sqrt(norm)
+      enddo
+    else
+      do j=1,ao_prim_num(i)
+        ao_coef_normalized(i,j) = ao_coef(i,j)
+      enddo
+    endif
+
     ! Normalization of the contracted basis functions
-    norm = 0.d0
+    if (ao_normalized) then
-    do j=1,ao_prim_num(i)
+      norm = 0.d0
-     do k=1,ao_prim_num(i)
+      do j=1,ao_prim_num(i)
-      call overlap_gaussian_xyz(C_A,C_A,ao_expo(i,j),ao_expo(i,k),powA,powA,overlap_x,overlap_y,overlap_z,c,nz)
+        do k=1,ao_prim_num(i)
-      norm = norm+c*ao_coef_normalized(i,j)*ao_coef_normalized(i,k)
+          call overlap_gaussian_xyz(C_A,C_A,ao_expo(i,j),ao_expo(i,k),powA,powA,overlap_x,overlap_y,overlap_z,c,nz)
-     enddo
+          norm = norm+c*ao_coef_normalized(i,j)*ao_coef_normalized(i,k)
-    enddo
+        enddo
-    ao_coef_normalization_factor(i) = 1.d0/sqrt(norm)
+      enddo
+      ao_coef_normalization_factor(i) = 1.d0/sqrt(norm)
+    else
+      ao_coef_normalization_factor(i) = 1.d0
+    endif
   enddo
 
 END_PROVIDER

src/cipsi/selection.irp.f

@@ -22,14 +22,14 @@ END_PROVIDER
 subroutine update_pt2_and_variance_weights(pt2, variance, norm, N_st)
   implicit none
   BEGIN_DOC
-! Updates the rPT2- and Variance- matching weights.
+! Updates the PT2- and Variance- matching weights.
   END_DOC
   integer, intent(in)          :: N_st
   double precision, intent(in) :: pt2(N_st)
   double precision, intent(in) :: variance(N_st)
   double precision, intent(in) :: norm(N_st)
 
-  double precision :: avg, rpt2(N_st), element, dt, x
+  double precision :: avg, pt2_rpt2(N_st), element, dt, x
   integer          :: k
   integer, save    :: i_iter=0
   integer, parameter :: i_itermax = 1
@@ -49,12 +49,13 @@ subroutine update_pt2_and_variance_weights(pt2, variance, norm, N_st)
   dt = 1.0d0
 
   do k=1,N_st
-    rpt2(k) = pt2(k)/(1.d0 + norm(k))
+    ! PT2 + rPT2
+    pt2_rpt2(k) = pt2(k)* (1.d0 + 1.d0/(1.d0 + norm(k)))
   enddo
 
-  avg = sum(rpt2(1:N_st)) / dble(N_st) - 1.d-32 ! Avoid future division by zero
+  avg = sum(pt2_rpt2(1:N_st)) / dble(N_st) - 1.d-32 ! Avoid future division by zero
   do k=1,N_st
-    element = exp(dt*(rpt2(k)/avg -1.d0))
+    element = exp(dt*(pt2_rpt2(k)/avg -1.d0))
     element = min(2.0d0 , element)
     element = max(0.5d0 , element)
     memo_pt2(k,i_iter) = element
@@ -71,7 +72,7 @@ subroutine update_pt2_and_variance_weights(pt2, variance, norm, N_st)
   enddo
 
   threshold_davidson_pt2 = min(1.d-6, &
-     max(threshold_davidson, 1.e-1 * PT2_relative_error * minval(abs(rpt2(1:N_states)))) )
+     max(threshold_davidson, 1.e-1 * PT2_relative_error * minval(abs(pt2(1:N_states)))) )
 
   SOFT_TOUCH pt2_match_weight variance_match_weight threshold_davidson_pt2
 end