Cleaned correlation_energy_ratio_max

plugins/CAS_SD_ZMQ/cassd_zmq.irp.f

@@ -9,6 +9,15 @@ program fci_zmq
   
   allocate (pt2(N_states))
   
+  double precision               :: hf_energy_ref
+  logical                        :: has
+  call ezfio_has_hartree_fock_energy(has)
+  if (has) then
+    call ezfio_get_hartree_fock_energy(hf_energy_ref)
+  else
+    hf_energy_ref = ref_bitmask_energy
+  endif
+
   pt2 = 1.d0
   threshold_davidson_in = threshold_davidson
   threshold_davidson = threshold_davidson_in * 100.d0
@@ -41,10 +50,23 @@ program fci_zmq
   print*,'Beginning the selection ...'
   E_CI_before(1:N_states) = CI_energy(1:N_states)
   
-  do while ( (N_det < N_det_max) .and. (maxval(abs(pt2(1:N_states))) > pt2_max) )
+  double precision :: correlation_energy_ratio
-    
+  correlation_energy_ratio = E_CI_before(1) - hf_energy_ref
-    print *,  'N_det          = ', N_det
+  correlation_energy_ratio = correlation_energy_ratio / (correlation_energy_ratio + pt2(1))
-    print *,  'N_states       = ', N_states
+
+  do while (                                                         &
+        (N_det < N_det_max) .and.                                    &
+        (maxval(abs(pt2(1:N_states))) > pt2_max) .and.               &
+        (correlation_energy_ratio < correlation_energy_ratio_max)    &
+        )
+
+    correlation_energy_ratio = E_CI_before(1) - hf_energy_ref
+    correlation_energy_ratio = correlation_energy_ratio / (correlation_energy_ratio + pt2(1))
+ 
+    print *,  'N_det             = ', N_det
+    print *,  'N_states          = ', N_states
+    print*,   'correlation_ratio = ', correlation_energy_ratio
+
     do k=1, N_states
       print*,'State ',k
       print *,  'PT2            = ', pt2(k)

plugins/Full_CI_ZMQ/fci_zmq.irp.f

@@ -10,11 +10,14 @@ program fci_zmq
   
   allocate (pt2(N_states))
 
-  IF (correlation_energy_ratio_max .NE. 1.d0) THEN
+  double precision               :: hf_energy_ref
- 
+  logical                        :: has
-    DOUBLE PRECISION :: hf_energy_ref
+  call ezfio_has_hartree_fock_energy(has)
-    CALL ezfio_get_hartree_fock_energy(hf_energy_ref)
+  if (has) then
-  END IF
+    call ezfio_get_hartree_fock_energy(hf_energy_ref)
+  else
+    hf_energy_ref = ref_bitmask_energy
+  endif
 
   pt2 = 1.d0
   threshold_davidson_in = threshold_davidson
@@ -47,28 +50,22 @@ program fci_zmq
   E_CI_before(1:N_states) = CI_energy(1:N_states)
   n_det_before = 0
   
+  double precision :: correlation_energy_ratio
+  correlation_energy_ratio = E_CI_before(1) - hf_energy_ref
+  correlation_energy_ratio = correlation_energy_ratio / (correlation_energy_ratio + pt2(1))
 
-  do while ( (N_det < N_det_max) .and. (maxval(abs(pt2(1:N_states))) > pt2_max) )
+  do while (                                                         &
+        (N_det < N_det_max) .and.                                    &
+        (maxval(abs(pt2(1:N_states))) > pt2_max) .and.               &
+        (correlation_energy_ratio < correlation_energy_ratio_max)    &
+        )
 
-
+    correlation_energy_ratio = E_CI_before(1) - hf_energy_ref
-    IF (correlation_energy_ratio_max .NE. 1.d0) THEN
+    correlation_energy_ratio = correlation_energy_ratio / (correlation_energy_ratio + pt2(1))
-      DOUBLE PRECISION :: correlation_energy_var, correlation_energy_var_ratio
-
-      correlation_energy_var = MAXVAL(E_CI_before) - hf_energy_ref
-      correlation_energy_var_ratio = correlation_energy_var / (correlation_energy_var + MAXVAL(pt2(:)))
-
-      IF (correlation_energy_ratio_max < correlation_energy_var_ratio) THEN
-        EXIT
-      ENDIF
-    
-    ENDIF
- 
  
     print *,  'N_det             = ', N_det
     print *,  'N_states          = ', N_states
-    IF (correlation_energy_ratio_max .NE. 1.d0) THEN
+    print*,   'correlation_ratio = ', correlation_energy_ratio
-    print*,   'correlation_ratio = ', correlation_energy_var 
-    ENDIF
 
     do k=1, N_states
       print*,'State ',k
@@ -76,6 +73,7 @@ program fci_zmq
       print *,  'E               = ', CI_energy(k)
       print *,  'E(before)+PT2   = ', E_CI_before(k)+pt2(k)
     enddo
+
     print *,  '-----'
     if(N_states.gt.1)then
       print*,'Variational Energy difference'
@@ -130,8 +128,8 @@ program fci_zmq
     double precision :: relative_error
     relative_error=1.d-3
     pt2 = 0.d0
-    call ZMQ_pt2(pt2,relative_error)
+    call ZMQ_pt2(pt2,relative_error) ! Stochastic PT2
-    !call ZMQ_selection(0, pt2)! pour non-stochastic
+    !call ZMQ_selection(0, pt2)      ! Deterministic PT2
     print *,  'Final step'
     print *,  'N_det    = ', N_det
     print *,  'N_states = ', N_states

plugins/Perturbation/EZFIO.cfg

@@ -11,17 +11,10 @@ doc: The selection process stops when the largest PT2 (for all the state) is low
 interface: ezfio,provider,ocaml
 default: 0.0001
 
-[var_pt2_ratio]
-type: Normalized_float
-doc: The selection process stops when the energy ratio variational/(variational+PT2)
-     is equal to var_pt2_ratio. (Obsolete. Need to be removed)
-interface: ezfio,provider,ocaml
-default: 0.75
-
 [correlation_energy_ratio_max]
 type: Normalized_float
-doc: The selection process stops at a fixed correlation ratio (usefull for getting same accuracy between molecules)
+doc: The selection process stops at a fixed correlation ratio (useful for getting same accuracy between molecules)
-     Defined as (E_CI-E_HF)/ (E_CI+PT2 - E_HF)  If Ratio eq 1 it will not be used. (E_HF) is not required. 
+     Defined as (E_CI-E_HF)/ (E_CI+PT2 - E_HF). (E_HF) is not required. 
 interface: ezfio,provider,ocaml
 default: 1.00
 

src/Utils/util.irp.f

@@ -234,30 +234,14 @@ double precision function dble_logfact(n) result(logfact2)
   ! n!!
   END_DOC
   integer                        :: n
-  double precision, save         :: memo(1:100)
+  integer :: k
-  integer, save                  :: memomax = 1
+  double precision :: prod
-  
+  prod=0.d0
-  ASSERT (iand(n,1) /= 0)
+  do k=2,n,2
-  if (n<=memomax) then
+   prod=prod+dlog(dfloat(k))
-    if (n<3) then
-      logfact2 = 0.d0
-    else
-      logfact2 = memo(n)
-    endif
-    return
-  endif
-  
-  integer                        :: i
-  memo(1) = 0.d0
-  do i=memomax+2,min(n,99),2
-    memo(i) = memo(i-2)+ dlog(dble(i))
-  enddo
-  memomax = min(n,99)
-  logfact2 = memo(memomax)
-  
-  do i=101,n,2
-    logfact2 += dlog(dble(i))
   enddo
+  logfact2=prod
+  return
   
 end function