From b7beac4557449fc556300d4ba35e3471f05470b4 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 23 Feb 2015 12:37:43 +0100 Subject: [PATCH] Added var_pt2_ratio program to stop at a constant ratio variational/(variational+PT2) --- data/ezfio_defaults | 1 + ocaml/Input_full_ci.ml | 26 ++++++++++++ src/Dets/diagonalize_CI.irp.f | 4 +- src/Full_CI/full_ci.ezfio_config | 1 + src/Full_CI/options.irp.f | 7 ++++ src/Full_CI/var_pt2_ratio.irp.f | 71 ++++++++++++++++++++++++++++++++ 6 files changed, 108 insertions(+), 2 deletions(-) create mode 100644 src/Full_CI/var_pt2_ratio.irp.f diff --git a/data/ezfio_defaults b/data/ezfio_defaults index 879ae78f..21163b55 100644 --- a/data/ezfio_defaults +++ b/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 diff --git a/ocaml/Input_full_ci.ml b/ocaml/Input_full_ci.ml index 0ee59f53..01548e88 100644 --- a/ocaml/Input_full_ci.ml +++ b/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 ;; diff --git a/src/Dets/diagonalize_CI.irp.f b/src/Dets/diagonalize_CI.irp.f index bd17f8bd..288f28c3 100644 --- a/src/Dets/diagonalize_CI.irp.f +++ b/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 diff --git a/src/Full_CI/full_ci.ezfio_config b/src/Full_CI/full_ci.ezfio_config index c9c4b3b6..65e9c393 100644 --- a/src/Full_CI/full_ci.ezfio_config +++ b/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 diff --git a/src/Full_CI/options.irp.f b/src/Full_CI/options.irp.f index 807f778a..553fbe9f 100644 --- a/src/Full_CI/options.irp.f +++ b/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 diff --git a/src/Full_CI/var_pt2_ratio.irp.f b/src/Full_CI/var_pt2_ratio.irp.f new file mode 100644 index 00000000..20395fa9 --- /dev/null +++ b/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