diff --git a/plugins/CAS_SD_ZMQ/cassd_zmq.irp.f b/plugins/CAS_SD_ZMQ/cassd_zmq.irp.f index 5b364400..0d079680 100644 --- a/plugins/CAS_SD_ZMQ/cassd_zmq.irp.f +++ b/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) ) - - print *, 'N_det = ', N_det - print *, 'N_states = ', N_states + 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) .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) diff --git a/plugins/Full_CI_ZMQ/fci_zmq.irp.f b/plugins/Full_CI_ZMQ/fci_zmq.irp.f index 585e40e3..e12033b4 100644 --- a/plugins/Full_CI_ZMQ/fci_zmq.irp.f +++ b/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 - CALL ezfio_get_hartree_fock_energy(hf_energy_ref) - END IF + 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 @@ -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) & + ) - - IF (correlation_energy_ratio_max .NE. 1.d0) THEN - 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 - + 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 - IF (correlation_energy_ratio_max .NE. 1.d0) THEN - print*, 'correlation_ratio = ', correlation_energy_var - ENDIF + print*, 'correlation_ratio = ', correlation_energy_ratio 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_selection(0, pt2)! pour non-stochastic + call ZMQ_pt2(pt2,relative_error) ! Stochastic PT2 + !call ZMQ_selection(0, pt2) ! Deterministic PT2 print *, 'Final step' print *, 'N_det = ', N_det print *, 'N_states = ', N_states diff --git a/plugins/Perturbation/EZFIO.cfg b/plugins/Perturbation/EZFIO.cfg index 19c6c52d..5a6ff5f7 100644 --- a/plugins/Perturbation/EZFIO.cfg +++ b/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) - 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. +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). (E_HF) is not required. interface: ezfio,provider,ocaml default: 1.00 @@ -35,5 +28,5 @@ default: 0.999 type: Threshold doc: Thresholds on selectors (fraction of the norm) for final PT2 calculation interface: ezfio,provider,ocaml -default: 1. +default: 1. diff --git a/plugins/Perturbation/var_pt2_ratio_provider.irp.f b/plugins/Perturbation/var_pt2_ratio_provider.irp.f new file mode 100644 index 00000000..43bcbf3e --- /dev/null +++ b/plugins/Perturbation/var_pt2_ratio_provider.irp.f @@ -0,0 +1,10 @@ +BEGIN_PROVIDER [ double precision, var_pt2_ratio ] + implicit none + BEGIN_DOC + ! The selection process stops when the energy ratio variational/(variational+PT2) + ! is equal to var_pt2_ratio + END_DOC + + var_pt2_ratio = correlation_energy_ratio_max +END_PROVIDER + diff --git a/src/Utils/util.irp.f b/src/Utils/util.irp.f index 4001e9df..d3b20a10 100644 --- a/src/Utils/util.irp.f +++ b/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, save :: memomax = 1 - - ASSERT (iand(n,1) /= 0) - if (n<=memomax) then - 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)) + integer :: k + double precision :: prod + prod=0.d0 + do k=2,n,2 + prod=prod+dlog(dfloat(k)) enddo + logfact2=prod + return end function