From a196526f8c130490626a6f5727bf03d6fcd2eb18 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 30 Mar 2016 21:33:38 +0200 Subject: [PATCH 1/8] 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 From 71e7c51608be11bdd0d338e2c3455df11ae62a13 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 30 Mar 2016 21:34:53 +0200 Subject: [PATCH 2/8] Added mrcc_noiter --- plugins/MRCC_CASSD/mrcc_noiter.irp.f | 41 ++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) create mode 100644 plugins/MRCC_CASSD/mrcc_noiter.irp.f 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 + From 2cd7922b81b5f9845543767991c555f2025aa91b Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 31 Mar 2016 10:59:46 +0200 Subject: [PATCH 3/8] changed lambda mrcc --- plugins/MRCC_Utils/mrcc_utils.irp.f | 12 ++++++------ tests/bats/qp.bats | 6 +++--- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index 79e8dd32..62cc82cc 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -22,12 +22,12 @@ BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ] i_pert_count += 1 lambda_mrcc(k,i) = 0.d0 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 - 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 +! endif enddo enddo diff --git a/tests/bats/qp.bats b/tests/bats/qp.bats index 1be43412..17e2d859 100644 --- a/tests/bats/qp.bats +++ b/tests/bats/qp.bats @@ -53,7 +53,7 @@ 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 @@ -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.2288641014933E+02 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.9483705 } @test "FCI H2O VDZ pseudo" { From 650c67911f63e5e5eed6fa736cdb35f9127fbb98 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 31 Mar 2016 11:05:39 +0200 Subject: [PATCH 4/8] Repaired tests --- tests/bats/qp.bats | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/bats/qp.bats b/tests/bats/qp.bats index 17e2d859..b6bb38ca 100644 --- a/tests/bats/qp.bats +++ b/tests/bats/qp.bats @@ -56,7 +56,7 @@ function run_HF() { 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.2288641014933E+02 1.e-4 + eq $energy -76.2288641014933 1.e-4 } @@ -166,7 +166,7 @@ function run_all_1h_1p() { } @test "SCF H2O VDZ pseudo" { - run_HF h2o_pseudo.ezfio -16.9483705 + run_HF h2o_pseudo.ezfio -16.9483708496728 } @test "FCI H2O VDZ pseudo" { From af3038f9951b3f82bea4bc78c5e45b7ed697d00a Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 31 Mar 2016 11:08:24 +0200 Subject: [PATCH 5/8] Repaired tests --- tests/bats/qp.bats | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/bats/qp.bats b/tests/bats/qp.bats index b6bb38ca..0977801d 100644 --- a/tests/bats/qp.bats +++ b/tests/bats/qp.bats @@ -166,7 +166,7 @@ function run_all_1h_1p() { } @test "SCF H2O VDZ pseudo" { - run_HF h2o_pseudo.ezfio -16.9483708496728 + run_HF h2o_pseudo.ezfio -0.169483703903844E+02 } @test "FCI H2O VDZ pseudo" { From 39d4186d2844a3eda9c881bdd95ed8a76b3edba5 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 31 Mar 2016 11:51:25 +0200 Subject: [PATCH 6/8] Removed extra MRCC iteration --- plugins/MRCC_CASSD/mrcc_cassd.irp.f | 1 - 1 file changed, 1 deletion(-) diff --git a/plugins/MRCC_CASSD/mrcc_cassd.irp.f b/plugins/MRCC_CASSD/mrcc_cassd.irp.f index 4b82bc5f..c5ed175f 100644 --- a/plugins/MRCC_CASSD/mrcc_cassd.irp.f +++ b/plugins/MRCC_CASSD/mrcc_cassd.irp.f @@ -18,7 +18,6 @@ subroutine run 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 From 1ae3e9b5f472362347f55d4393e67aef99bcc760 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 31 Mar 2016 16:48:23 +0200 Subject: [PATCH 7/8] Official Lambda_MRCC --- plugins/MRCC_Utils/mrcc_utils.irp.f | 61 +++++++++++++---------------- plugins/Psiref_CAS/psi_ref.irp.f | 1 + 2 files changed, 28 insertions(+), 34 deletions(-) diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index 62cc82cc..80fdd4c7 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -3,50 +3,43 @@ BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ] BEGIN_DOC ! cm/ or perturbative 1/Delta_E(m) END_DOC - integer :: i,k + integer :: i,k double precision :: ihpsi_current(N_states) - integer :: i_pert_count - + 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) - 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 - 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 -! endif + + 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) = 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 - 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 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 From fcd0b39077c4c4a960df331209439a33a4f10d8c Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 31 Mar 2016 16:51:58 +0200 Subject: [PATCH 8/8] Updated tests --- tests/bats/qp.bats | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/bats/qp.bats b/tests/bats/qp.bats index 0977801d..6f1c916b 100644 --- a/tests/bats/qp.bats +++ b/tests/bats/qp.bats @@ -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.2288641014933 1.e-4 + 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.169483703903844E+02 + run_HF h2o_pseudo.ezfio -16.9483708495521 } @test "FCI H2O VDZ pseudo" {