diff --git a/plugins/MRCC_CASSD/EZFIO.cfg b/plugins/MRCC_CASSD/EZFIO.cfg index e145c9e0..17ee7f29 100644 --- a/plugins/MRCC_CASSD/EZFIO.cfg +++ b/plugins/MRCC_CASSD/EZFIO.cfg @@ -7,11 +7,11 @@ interface: ezfio type: Threshold doc: Threshold on the convergence of the MRCC energy interface: ezfio,provider,ocaml -default: 1.e-7 +default: 1.e-5 [n_it_mrcc_max] type: Strictly_positive_int doc: Maximum number of MRCC iterations interface: ezfio,provider,ocaml -default: 20 +default: 10 diff --git a/plugins/MRCC_CASSD/mrcc_cassd.irp.f b/plugins/MRCC_CASSD/mrcc_cassd.irp.f index e784a167..c5ed175f 100644 --- a/plugins/MRCC_CASSD/mrcc_cassd.irp.f +++ b/plugins/MRCC_CASSD/mrcc_cassd.irp.f @@ -1,13 +1,46 @@ program mrcc implicit none - if (.not.read_wf) then - print *, 'read_wf has to be true.' - stop 1 - endif + read_wf = .True. + SOFT_TOUCH read_wf call print_cas_coefs - call run_mrcc + call set_generators_bitmasks_as_holes_and_particles + call run end +subroutine run + implicit none + + integer :: i + + double precision :: E_new, E_old, delta_e + integer :: iteration + double precision :: E_past(4), lambda + E_new = 0.d0 + delta_E = 1.d0 + iteration = 0 + lambda = 1.d0 + do while (delta_E > thresh_mrcc) + iteration += 1 + print *, '===========================' + print *, 'MRCC Iteration', iteration + print *, '===========================' + print *, '' + E_old = sum(ci_energy_dressed) + call write_double(6,ci_energy_dressed(1),"MRCC energy") + call diagonalize_ci_dressed(lambda) + E_new = sum(ci_energy_dressed) + delta_E = dabs(E_new - E_old) + call save_wavefunction + call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1)) + if (iteration > n_it_mrcc_max) then + exit + endif + enddo + call write_double(6,ci_energy_dressed(1),"Final MRCC energy") + +end + + subroutine print_cas_coefs implicit none @@ -18,7 +51,7 @@ subroutine print_cas_coefs print *, psi_cas_coef(i,:) call debug_det(psi_cas(1,1,i),N_int) enddo - call write_double(6,ci_energy(1),"Initial CI energy") + end diff --git a/plugins/MRCC_CASSD/mrcc_noiter.irp.f b/plugins/MRCC_CASSD/mrcc_noiter.irp.f new file mode 100644 index 00000000..949a0dab --- /dev/null +++ b/plugins/MRCC_CASSD/mrcc_noiter.irp.f @@ -0,0 +1,41 @@ +program mrcc_noiter + implicit none + read_wf = .True. + SOFT_TOUCH read_wf + call print_cas_coefs + call set_generators_bitmasks_as_holes_and_particles + call run +end + +subroutine run + implicit none + + double precision :: lambda + integer :: i,j + do j=1,N_states_diag + do i=1,N_det + psi_coef(i,j) = CI_eigenvectors_dressed(i,j) + enddo + enddo + SOFT_TOUCH psi_coef ci_energy_dressed + call write_double(6,ci_energy_dressed(1),"Final MRCC energy") + call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1)) + call save_wavefunction + +end + + +subroutine print_cas_coefs + implicit none + + integer :: i,j + print *, 'CAS' + print *, '===' + do i=1,N_det_cas + print *, psi_cas_coef(i,:) + call debug_det(psi_cas(1,1,i),N_int) + enddo + call write_double(6,ci_energy(1),"Initial CI energy") + +end + diff --git a/plugins/MRCC_Utils/mrcc_general.irp.f b/plugins/MRCC_Utils/mrcc_general.irp.f index 5d9acfc1..50343fdb 100644 --- a/plugins/MRCC_Utils/mrcc_general.irp.f +++ b/plugins/MRCC_Utils/mrcc_general.irp.f @@ -11,12 +11,13 @@ subroutine mrcc_iterations double precision :: E_new, E_old, delta_e integer :: iteration,i_oscillations - double precision :: E_past(4) + double precision :: E_past(4), lambda E_new = 0.d0 delta_E = 1.d0 iteration = 0 j = 1 i_oscillations = 0 + lambda = 1.d0 do while (delta_E > 1.d-7) iteration += 1 print *, '===========================' @@ -25,10 +26,15 @@ subroutine mrcc_iterations print *, '' E_old = sum(ci_energy_dressed) call write_double(6,ci_energy_dressed(1),"MRCC energy") - call diagonalize_ci_dressed + call diagonalize_ci_dressed(lambda) E_new = sum(ci_energy_dressed) delta_E = dabs(E_new - E_old) - +! if (E_new > E_old) then +! lambda = lambda * 0.7d0 +! else +! lambda = min(1.d0, lambda * 1.1d0) +! endif +! print *, 'energy lambda ', lambda E_past(j) = E_new j +=1 call save_wavefunction diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index 712daca1..80fdd4c7 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -1,89 +1,45 @@ -! BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ] -! implicit none -! BEGIN_DOC -! ! cm/ or perturbative 1/Delta_E(m) -! END_DOC -! integer :: i,k -! double precision :: ihpsi(N_states),ihpsi_current(N_states) -! integer :: i_pert_count -! -! i_pert_count = 0 -! lambda_mrcc = 0.d0 -! -! do i=1,N_det_non_ref -! call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref,size(psi_ref_coef,1), n_states, ihpsi_current) -! do k=1,N_states -! if (ihpsi_current(k) == 0.d0) then -! ihpsi_current(k) = 1.d-32 -! endif -! if(dabs(ihpsi_current(k) * psi_non_ref_coef(i,k)) < 1d-6) then -! i_pert_count +=1 -! else -! lambda_mrcc(k,i) = psi_non_ref_coef(i,k)/ihpsi_current(k) -! endif -! enddo -! enddo -! -! print*,'N_det_non_ref = ',N_det_non_ref -! print*,'Number of ignored determinants = ',i_pert_count -! print*,'psi_coef_ref_ratio = ',psi_ref_coef(2,1)/psi_ref_coef(1,1) -! END_PROVIDER - BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ] - implicit none - BEGIN_DOC - ! cm/ or perturbative 1/Delta_E(m) - END_DOC - integer :: i,k + implicit none + BEGIN_DOC + ! cm/ or perturbative 1/Delta_E(m) + END_DOC + integer :: i,k double precision :: ihpsi_current(N_states) integer :: i_pert_count double precision :: hii, lambda_pert - + i_pert_count = 0 lambda_mrcc = 0.d0 - - do i=1,N_det_non_ref - call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref, & - size(psi_ref_coef,1), N_states,ihpsi_current) - call i_H_j(psi_non_ref(1,1,i),psi_non_ref(1,1,i),N_int,hii) - do k=1,N_states - if (ihpsi_current(k) == 0.d0) then - ihpsi_current(k) = 1.d-32 - endif - lambda_mrcc(k,i) = psi_non_ref_coef(i,k)/ihpsi_current(k) - if ( dabs(psi_non_ref_coef(i,k)*ihpsi_current(k)) < 1.d-6 ) then - i_pert_count += 1 - lambda_mrcc(k,i) = 0.d0 -! lambda_pert = 1.d0 / (psi_ref_energy_diagonalized(k)-hii) -! if((ihpsi_current(k) * lambda_pert) < 0.5d0 * psi_non_ref_coef_restart(i,k) ) then -! lambda_mrcc(k,i) = 0.d0 -! endif - endif - double precision, parameter :: x = 2.d0 - if (lambda_mrcc(k,i) > x) then - lambda_mrcc(k,i) = x - else if (lambda_mrcc(k,i) < -x) then - lambda_mrcc(k,i) = -x + + do i=1,N_det_non_ref + call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref,& + size(psi_ref_coef,1), N_states,ihpsi_current) + call i_H_j(psi_non_ref(1,1,i),psi_non_ref(1,1,i),N_int,hii) + do k=1,N_states + if (ihpsi_current(k) == 0.d0) then + ihpsi_current(k) = 1.d-32 endif - enddo - enddo - + lambda_mrcc(k,i) = min(0.d0,psi_non_ref_coef(i,k)/ihpsi_current(k) ) + lambda_pert = 1.d0 / (psi_ref_energy_diagonalized(k)-hii) + if (lambda_pert / lambda_mrcc(k,i) < 0.5d0 ) then + i_pert_count += 1 + lambda_mrcc(k,i) = 0.d0 + endif + enddo + enddo + print*,'N_det_non_ref = ',N_det_non_ref - print*,'Number of ignored determinants = ',i_pert_count + print*,'Number of ignored determinants = ',i_pert_count print*,'psi_coef_ref_ratio = ',psi_ref_coef(2,1)/psi_ref_coef(1,1) - print*,'lambda min/max = ',maxval(dabs(lambda_mrcc)), minval(dabs(lambda_mrcc)) + print*,'lambda max = ',maxval(dabs(lambda_mrcc)) END_PROVIDER -!BEGIN_PROVIDER [ double precision, delta_ij_non_ref, (N_det_non_ref, N_det_non_ref,N_states) ] -!implicit none -!BEGIN_DOC -!! Dressing matrix in SD basis -!END_DOC -!delta_ij_non_ref = 0.d0 -!call H_apply_mrcc_simple(delta_ij_non_ref,N_det_non_ref) -!END_PROVIDER + + + + BEGIN_PROVIDER [ double precision, hij_mrcc, (N_det_non_ref,N_det_ref) ] implicit none @@ -218,17 +174,19 @@ BEGIN_PROVIDER [ double precision, CI_energy_dressed, (N_states_diag) ] END_PROVIDER -subroutine diagonalize_CI_dressed +subroutine diagonalize_CI_dressed(lambda) implicit none BEGIN_DOC ! Replace the coefficients of the CI states by the coefficients of the ! eigenstates of the CI matrix END_DOC + double precision, intent(in) :: lambda integer :: i,j do j=1,N_states_diag do i=1,N_det - psi_coef(i,j) = CI_eigenvectors_dressed(i,j) + psi_coef(i,j) = lambda * CI_eigenvectors_dressed(i,j) + (1.d0 - lambda) * psi_coef(i,j) enddo + call normalize(psi_coef(1,j), N_det) enddo SOFT_TOUCH psi_coef diff --git a/plugins/Psiref_CAS/psi_ref.irp.f b/plugins/Psiref_CAS/psi_ref.irp.f index f67f0587..d3dbda08 100644 --- a/plugins/Psiref_CAS/psi_ref.irp.f +++ b/plugins/Psiref_CAS/psi_ref.irp.f @@ -23,6 +23,7 @@ use bitmasks psi_ref_coef(i,k) = psi_cas_coef(i,k) enddo enddo + call normalize(psi_ref_coef, N_det_ref) END_PROVIDER diff --git a/tests/bats/qp.bats b/tests/bats/qp.bats index 1be43412..6f1c916b 100644 --- a/tests/bats/qp.bats +++ b/tests/bats/qp.bats @@ -53,10 +53,10 @@ function test_exe() { } function run_HF() { - thresh=1.e-8 + thresh=1.e-7 test_exe SCF || skip ezfio set_file $1 - ezfio set hartree_fock thresh_scf 1.e-10 + ezfio set hartree_fock thresh_scf 1.e-11 qp_run SCF $1 energy="$(ezfio get hartree_fock energy)" eq $energy $2 $thresh @@ -155,7 +155,7 @@ function run_all_1h_1p() { ezfio set determinants read_wf True qp_run mrcc_cassd $INPUT energy="$(ezfio get mrcc_cassd energy)" - eq $energy -76.2289109271715 1.E-3 + eq $energy -76.2286679037553 1.e-4 } @@ -166,7 +166,7 @@ function run_all_1h_1p() { } @test "SCF H2O VDZ pseudo" { - run_HF h2o_pseudo.ezfio -0.169483703904991E+02 + run_HF h2o_pseudo.ezfio -16.9483708495521 } @test "FCI H2O VDZ pseudo" {