Added var_pt2_ratio program to stop at a constant ratio variational/(variational+PT2)

data/ezfio_defaults

@@ -28,6 +28,7 @@ full_ci
   n_det_max_fci              10000
   pt2_max                    1.e-4
   do_pt2_end                 True
+  var_pt2_ratio              0.75
 
 cassd  
   n_det_max_cassd            10000

ocaml/Input_full_ci.ml

@@ -7,6 +7,7 @@ module Full_ci : sig
     { n_det_max_fci      : Det_number_max.t;
       pt2_max            : PT2_energy.t;
       do_pt2_end         : bool;
+      var_pt2_ratio      : Normalized_float.t;
     } with sexp
   ;;
   val read  : unit -> t option
@@ -19,6 +20,7 @@ end = struct
     { n_det_max_fci      : Det_number_max.t;
       pt2_max            : PT2_energy.t;
       do_pt2_end         : bool;
+      var_pt2_ratio      : Normalized_float.t;
     } with sexp
   ;;
 
@@ -39,6 +41,21 @@ end = struct
     |> Ezfio.set_full_ci_n_det_max_fci
   ;;
 
+  let read_var_pt2_ratio () =
+    if not (Ezfio.has_full_ci_var_pt2_ratio ()) then
+       get_default "var_pt2_ratio"
+       |> Float.of_string
+       |> Ezfio.set_full_ci_var_pt2_ratio
+    ;
+    Ezfio.get_full_ci_var_pt2_ratio ()
+    |> Normalized_float.of_float
+  ;;
+
+  let write_var_pt2_ratio ratio =
+    Normalized_float.to_float ratio
+    |> Ezfio.set_full_ci_var_pt2_ratio
+  ;;
+
   let read_pt2_max () = 
     if not (Ezfio.has_full_ci_pt2_max ()) then
        get_default "pt2_max"
@@ -73,6 +90,7 @@ end = struct
     { n_det_max_fci = read_n_det_max_fci ();
       pt2_max       = read_pt2_max ();
       do_pt2_end    = read_do_pt2_end ();
+      var_pt2_ratio = read_var_pt2_ratio ();
     }
   ;;
 
@@ -80,10 +98,12 @@ end = struct
   let write { n_det_max_fci ;
               pt2_max       ;
               do_pt2_end    ;
+              var_pt2_ratio ;
             } = 
      write_n_det_max_fci n_det_max_fci;
      write_pt2_max       pt2_max;
      write_do_pt2_end    do_pt2_end;
+     write_var_pt2_ratio var_pt2_ratio;
   ;;
 
   let to_string b =
@@ -91,10 +111,12 @@ end = struct
 n_det_max_fci = %s
 pt2_max = %s
 do_pt2_end = %s
+var_pt2_ratio = %s
 "
         (Det_number_max.to_string b.n_det_max_fci)
         (PT2_energy.to_string b.pt2_max)
         (Bool.to_string b.do_pt2_end)
+        (Normalized_float.to_string b.var_pt2_ratio)
     ;;
 
   let to_rst b =
@@ -111,10 +133,14 @@ Compute E(PT2) at the end ::
 
   do_pt2_end = %s
 
+Target energy ratio variational/(variational+PT2) ::
+
+  var_pt2_ratio = %s
 "
         (Det_number_max.to_string b.n_det_max_fci)
         (PT2_energy.to_string b.pt2_max)
         (Bool.to_string b.do_pt2_end)
+        (Normalized_float.to_string b.var_pt2_ratio)
     |> Rst_string.of_string
     ;;
 

src/Dets/diagonalize_CI.irp.f

@@ -69,10 +69,10 @@ END_PROVIDER
     i_state = 0
     do j=1,N_det
       call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2)
-      print *, 'j = ',j,s2, expected_s2
+!      print *, 'j = ',j,s2, expected_s2
       if(dabs(s2-expected_s2).le.0.3d0)then
        i_state += 1
-       print *,  'i_state = ',i_state
+!       print *,  'i_state = ',i_state
        do i=1,N_det
          CI_eigenvectors(i,i_state) = eigenvectors(i,j)
        enddo

src/Full_CI/full_ci.ezfio_config

@@ -4,3 +4,4 @@ full_ci
   do_pt2_end     logical
   energy         double precision
   energy_pt2     double precision
+  var_pt2_ratio  double precision

src/Full_CI/options.irp.f

@@ -21,5 +21,12 @@ T.set_doc       ( """The selection process stops when the largest PT2 (for all t
 is lower than pt2_max in absolute value""" )
 T.set_ezfio_name( "pt2_max" )
 print T
+
+T.set_type      ( "double precision" )
+T.set_name      ( "var_pt2_ratio" )
+T.set_doc       ( """The selection process stops when the energy ratio variational/(variational+PT2)
+is equal to var_pt2_ratio""" )
+T.set_ezfio_name( "var_pt2_ratio" )
+print T
 END_SHELL
 

src/Full_CI/var_pt2_ratio.irp.f

@@ -0,0 +1,71 @@
+program var_pt2_ratio_run
+  implicit none
+  integer                        :: i,k
+
+  
+  double precision, allocatable  :: pt2(:), norm_pert(:), H_pert_diag(:)
+  integer                        :: N_st, degree
+  N_st = N_states
+  allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st))
+  character*(64)                 :: perturbation
+  
+  double precision, allocatable  :: psi_det_save(:,:,:), psi_coef_save(:,:)
+
+  double precision :: E_fci, E_var, ratio, E_ref
+  integer :: Nmin, Nmax
+
+  pt2 = 1.d0
+  diag_algorithm = "Lapack"
+
+  ratio = 0.d0
+  Nmin=1
+  do while (ratio < var_pt2_ratio)
+    call H_apply_FCI(pt2, norm_pert, H_pert_diag,  N_st)
+    psi_det = psi_det_sorted
+    psi_coef = psi_coef_sorted
+    soft_touch N_det psi_det psi_coef
+    call diagonalize_CI
+    ratio = (CI_energy(1) - HF_energy) / (CI_energy(1)+pt2(1) - HF_energy)
+  enddo
+
+  threshold_selectors = 1.d0
+  threshold_generators = 0.999d0
+  call diagonalize_CI
+  call H_apply_FCI_PT2(pt2, norm_pert, H_pert_diag,  N_st)
+  E_ref = CI_energy(1) + pt2(1)
+  threshold_selectors = 0.999d0
+  threshold_generators = 0.99d0
+
+  Nmax=N_det
+  Nmin=1
+  do while (Nmax-Nmin > 1)
+    ratio = (CI_energy(1) - HF_energy) / (E_ref - HF_energy)
+
+    if (ratio < var_pt2_ratio) then
+      Nmin = N_det
+      Nmax = max(Nmax,Nmin+10)
+      ! Select new determinants
+      call H_apply_FCI(pt2, norm_pert, H_pert_diag,  N_st)
+    else
+      Nmax = N_det
+      N_det = Nmin + (Nmax-Nmin)/2
+    endif
+
+    psi_det = psi_det_sorted
+    psi_coef = psi_coef_sorted
+    soft_touch N_det psi_det psi_coef
+    call diagonalize_CI
+    call save_wavefunction
+    print *,  'Det min, Det max: ', Nmin, Nmax
+    print *,  'Ratio           : ', ratio, '  ~  ', var_pt2_ratio
+    print *,  'HF_energy = ', HF_energy
+    print *,  'Est FCI   = ', E_ref
+    print *,  'N_det     = ', N_det
+    print *,  'E         = ', CI_energy(1)
+    call ezfio_set_full_ci_energy(CI_energy)
+    if (abort_all) then
+      exit
+    endif
+  enddo
+  deallocate(pt2,norm_pert)
+end