From a196526f8c130490626a6f5727bf03d6fcd2eb18 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 30 Mar 2016 21:33:38 +0200 Subject: [PATCH] Improved MRCC --- plugins/MRCC_CASSD/EZFIO.cfg | 4 +- plugins/MRCC_CASSD/mrcc_cassd.irp.f | 46 ++++++++++++++++--- plugins/MRCC_Utils/mrcc_general.irp.f | 12 +++-- plugins/MRCC_Utils/mrcc_utils.irp.f | 63 +++++++++++++++++---------- 4 files changed, 91 insertions(+), 34 deletions(-) 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..4b82bc5f 100644 --- a/plugins/MRCC_CASSD/mrcc_cassd.irp.f +++ b/plugins/MRCC_CASSD/mrcc_cassd.irp.f @@ -1,13 +1,47 @@ 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 + call diagonalize_ci_dressed(lambda) + 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 +52,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_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 c41a3431..f321cfa2 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -4,29 +4,44 @@ BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ] ! 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 + 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-8 ) 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 - if(dabs(ihpsi_current(k) * psi_non_ref_coef(i,k)) < 1d-5) then - i_pert_count +=1 - else - lambda_mrcc(k,i) = psi_non_ref_coef(i,k)/ihpsi_current(k) - endif - enddo - enddo + 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 + 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) + print*,'lambda min/max = ',maxval(dabs(lambda_mrcc)), minval(dabs(lambda_mrcc)) - 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 @@ -174,17 +189,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