From 4cb4d5e41696dd1eb8dd3b970328a6fd1fae0719 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 8 Jun 2017 22:15:42 +0200 Subject: [PATCH] Working on Slater dressing --- config/gfortran_debug.cfg | 2 +- plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f | 19 +- plugins/Hartree_Fock/diagonalize_fock.irp.f | 172 ++++++++++-------- plugins/Hartree_Fock_SlaterDressed/EZFIO.cfg | 14 ++ .../LinearSystem.irp.f | 53 +++--- .../SCF_dressed.irp.f | 69 +++++++ .../Hartree_Fock_SlaterDressed/at_nucl.irp.f | 2 +- .../integrals.irp.f | 101 ++++++++-- .../Hartree_Fock_SlaterDressed/slater.irp.f | 43 +++++ plugins/mrcepa0/dressing.irp.f | 27 +-- 10 files changed, 364 insertions(+), 138 deletions(-) create mode 100644 plugins/Hartree_Fock_SlaterDressed/EZFIO.cfg create mode 100644 plugins/Hartree_Fock_SlaterDressed/SCF_dressed.irp.f create mode 100644 plugins/Hartree_Fock_SlaterDressed/slater.irp.f diff --git a/config/gfortran_debug.cfg b/config/gfortran_debug.cfg index 214c3fbe..91a12345 100644 --- a/config/gfortran_debug.cfg +++ b/config/gfortran_debug.cfg @@ -13,7 +13,7 @@ FC : gfortran -g -ffree-line-length-none -I . LAPACK_LIB : -llapack -lblas IRPF90 : irpf90 -IRPF90_FLAGS : --ninja --align=32 +IRPF90_FLAGS : --ninja --align=32 --assert # Global options ################ diff --git a/plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f b/plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f index 032b7cd1..988528f3 100644 --- a/plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f +++ b/plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f @@ -216,8 +216,25 @@ END_DOC double precision, allocatable :: AF(:,:) allocate (AF(dim_DIIS+1,dim_DIIS+1)) double precision :: rcond, ferr, berr - integer :: iwork(dim_DIIS+1) + integer :: iwork(dim_DIIS+1), lwork + call dsysvx('N','U',dim_DIIS+1,1, & + B_matrix_DIIS,size(B_matrix_DIIS,1), & + AF, size(AF,1), & + ipiv, & + C_vector_DIIS,size(C_vector_DIIS,1), & + X_vector_DIIS,size(X_vector_DIIS,1), & + rcond, & + ferr, & + berr, & + scratch,-1, & + iwork, & + info & + ) + lwork = int(scratch(1,1)) + deallocate(scratch) + allocate(scratch(lwork,1)) + call dsysvx('N','U',dim_DIIS+1,1, & B_matrix_DIIS,size(B_matrix_DIIS,1), & AF, size(AF,1), & diff --git a/plugins/Hartree_Fock/diagonalize_fock.irp.f b/plugins/Hartree_Fock/diagonalize_fock.irp.f index 78d08aa3..b303b537 100644 --- a/plugins/Hartree_Fock/diagonalize_fock.irp.f +++ b/plugins/Hartree_Fock/diagonalize_fock.irp.f @@ -1,115 +1,129 @@ - BEGIN_PROVIDER [ double precision, diagonal_Fock_matrix_mo, (mo_tot_num) ] + BEGIN_PROVIDER [ double precision, diagonal_Fock_matrix_mo, (ao_num) ] &BEGIN_PROVIDER [ double precision, eigenvectors_Fock_matrix_mo, (ao_num_align,mo_tot_num) ] implicit none BEGIN_DOC ! Diagonal Fock matrix in the MO basis END_DOC - integer :: i,j, m + integer :: i,j integer :: liwork, lwork, n, info - integer, allocatable :: iwork(:), isuppz(:) - double precision, allocatable :: work(:), F(:,:), F2(:,:) - integer :: iorb,jorb + integer, allocatable :: iwork(:) + double precision, allocatable :: work(:), F(:,:), S(:,:) - allocate( F(mo_tot_num,mo_tot_num),F2(mo_tot_num,mo_tot_num), isuppz(2*mo_tot_num) ) + allocate( F(mo_tot_num,mo_tot_num) ) do j=1,mo_tot_num do i=1,mo_tot_num F(i,j) = Fock_matrix_mo(i,j) enddo enddo if(no_oa_or_av_opt)then + integer :: iorb,jorb do i = 1, n_act_orb iorb = list_act(i) - ASSERT (iorb > 0) - ASSERT (iorb <= mo_tot_num) do j = 1, n_inact_orb jorb = list_inact(j) - ASSERT (jorb > 0) - ASSERT (jorb <= mo_tot_num) F(iorb,jorb) = 0.d0 F(jorb,iorb) = 0.d0 enddo do j = 1, n_virt_orb jorb = list_virt(j) - ASSERT (jorb > 0) - ASSERT (jorb <= mo_tot_num) F(iorb,jorb) = 0.d0 F(jorb,iorb) = 0.d0 enddo do j = 1, n_core_orb jorb = list_core(j) - ASSERT (jorb > 0) - ASSERT (jorb <= mo_tot_num) F(iorb,jorb) = 0.d0 F(jorb,iorb) = 0.d0 enddo enddo - endif - - - - - ! Insert level shift here - do i = elec_beta_num+1, elec_alpha_num - F(i,i) += 0.5d0*level_shift - enddo - - do i = elec_alpha_num+1, mo_tot_num - F(i,i) += level_shift - enddo - - n = mo_tot_num - lwork = 1+6*n + 2*n*n - liwork = 10*n - - allocate(work(lwork)) - allocate(iwork(liwork) ) - - call dsyevr('V', 'A', 'U', mo_tot_num, F, size(F,1), & - -100.d0, 100.d0, 1, mo_tot_num, 0.d0, & - m, diagonal_Fock_matrix_mo, & - F2, size(F2,1), & - isuppz, work, lwork, iwork, liwork, info) - - if (info /= 0) then - print *, irp_here//' DSYEV failed : ', info - stop 1 - endif - - call dgemm('N','N',ao_num,mo_tot_num,mo_tot_num, 1.d0, & - mo_coef, size(mo_coef,1), F2, size(F2,1), & - 0.d0, eigenvectors_Fock_matrix_mo, size(eigenvectors_Fock_matrix_mo,1)) - deallocate(work, F2, F) - deallocate(iwork, isuppz) - - - - -END_PROVIDER + endif -BEGIN_PROVIDER [double precision, diagonal_Fock_matrix_mo_sum, (mo_tot_num)] - implicit none - BEGIN_DOC - ! diagonal element of the fock matrix calculated as the sum over all the interactions - ! with all the electrons in the RHF determinant - ! diagonal_Fock_matrix_mo_sum(i) = sum_{j=1, N_elec} 2 J_ij -K_ij - END_DOC - integer :: i,j - double precision :: accu - do j = 1,elec_alpha_num - accu = 0.d0 - do i = 1, elec_alpha_num - accu += 2.d0 * mo_bielec_integral_jj_from_ao(i,j) - mo_bielec_integral_jj_exchange_from_ao(i,j) - enddo - diagonal_Fock_matrix_mo_sum(j) = accu + mo_mono_elec_integral(j,j) - enddo - do j = elec_alpha_num+1,mo_tot_num - accu = 0.d0 - do i = 1, elec_alpha_num - accu += 2.d0 * mo_bielec_integral_jj_from_ao(i,j) - mo_bielec_integral_jj_exchange_from_ao(i,j) - enddo - diagonal_Fock_matrix_mo_sum(j) = accu + mo_mono_elec_integral(j,j) - enddo + + + + ! Insert level shift here + do i = elec_beta_num+1, elec_alpha_num + F(i,i) += 0.5d0*level_shift + enddo + + do i = elec_alpha_num+1, mo_tot_num + F(i,i) += level_shift + enddo + + n = mo_tot_num + lwork = 1+6*n + 2*n*n + liwork = 3 + 5*n + + allocate(work(lwork)) + allocate(iwork(liwork) ) + + lwork = -1 + liwork = -1 + + call dsyevd( 'V', 'U', mo_tot_num, F, & + size(F,1), diagonal_Fock_matrix_mo, & + work, lwork, iwork, liwork, info) + + if (info /= 0) then + print *, irp_here//' DSYEVD failed : ', info + stop 1 + endif + lwork = int(work(1)) + liwork = iwork(1) + deallocate(iwork) + deallocate(work) + + allocate(work(lwork)) + allocate(iwork(liwork) ) + call dsyevd( 'V', 'U', mo_tot_num, F, & + size(F,1), diagonal_Fock_matrix_mo, & + work, lwork, iwork, liwork, info) + deallocate(iwork) + + + if (info /= 0) then + call dsyev( 'V', 'L', mo_tot_num, F, & + size(F,1), diagonal_Fock_matrix_mo, & + work, lwork, info) + if (info /= 0) then + print *, irp_here//' DSYEV failed : ', info + stop 1 + endif + endif + + call dgemm('N','N',ao_num,mo_tot_num,mo_tot_num, 1.d0, & + mo_coef, size(mo_coef,1), F, size(F,1), & + 0.d0, eigenvectors_Fock_matrix_mo, size(eigenvectors_Fock_matrix_mo,1)) + deallocate(work, F) + + + +END_PROVIDER + +BEGIN_PROVIDER [double precision, diagonal_Fock_matrix_mo_sum, (mo_tot_num)] + implicit none + BEGIN_DOC + ! diagonal element of the fock matrix calculated as the sum over all the interactions + ! with all the electrons in the RHF determinant + ! diagonal_Fock_matrix_mo_sum(i) = sum_{j=1, N_elec} 2 J_ij -K_ij + END_DOC + integer :: i,j + double precision :: accu + do j = 1,elec_alpha_num + accu = 0.d0 + do i = 1, elec_alpha_num + accu += 2.d0 * mo_bielec_integral_jj_from_ao(i,j) - mo_bielec_integral_jj_exchange_from_ao(i,j) + enddo + diagonal_Fock_matrix_mo_sum(j) = accu + mo_mono_elec_integral(j,j) + enddo + do j = elec_alpha_num+1,mo_tot_num + accu = 0.d0 + do i = 1, elec_alpha_num + accu += 2.d0 * mo_bielec_integral_jj_from_ao(i,j) - mo_bielec_integral_jj_exchange_from_ao(i,j) + enddo + diagonal_Fock_matrix_mo_sum(j) = accu + mo_mono_elec_integral(j,j) + enddo + END_PROVIDER diff --git a/plugins/Hartree_Fock_SlaterDressed/EZFIO.cfg b/plugins/Hartree_Fock_SlaterDressed/EZFIO.cfg new file mode 100644 index 00000000..555066a2 --- /dev/null +++ b/plugins/Hartree_Fock_SlaterDressed/EZFIO.cfg @@ -0,0 +1,14 @@ +[slater_expo_ezfio] +type: double precision +doc: Exponents of the additional Slater functions +size: (nuclei.nucl_num) +interface: ezfio, provider + +[slater_coef_ezfio] +type: double precision +doc: Exponents of the additional Slater functions +size: (mo_basis.mo_tot_num,nuclei.nucl_num) +interface: ezfio, provider + + + diff --git a/plugins/Hartree_Fock_SlaterDressed/LinearSystem.irp.f b/plugins/Hartree_Fock_SlaterDressed/LinearSystem.irp.f index 579dee57..fd23d15f 100644 --- a/plugins/Hartree_Fock_SlaterDressed/LinearSystem.irp.f +++ b/plugins/Hartree_Fock_SlaterDressed/LinearSystem.irp.f @@ -6,21 +6,19 @@ BEGIN_PROVIDER [ double precision, cusp_A, (nucl_num, nucl_num) ] integer :: mu, A, B - do B=1,nucl_num - do A=1,nucl_num - cusp_A(A,B) = 0.d0 - if (A/=B) then - cusp_A(A,B) -= slater_value_at_nucl(A,B) - endif - do mu=1,ao_num - cusp_A(A,B) += slater_overlap(mu,B) * ao_value_at_nucl(mu,A) - enddo - enddo + cusp_A = 0.d0 + do A=1,nucl_num + cusp_A(A,A) = slater_expo(A)/nucl_charge(A) * slater_value_at_nucl(A,A) + do B=1,nucl_num + cusp_A(A,B) -= slater_value_at_nucl(B,A) + do mu=1,ao_num + cusp_A(A,B) += GauSlaOverlap_matrix(mu,B) * ao_value_at_nucl(mu,A) + enddo + enddo enddo - END_PROVIDER -BEGIN_PROVIDER [ double precision, cusp_C, (nucl_num, mo_tot_num) ] +BEGIN_PROVIDER [ double precision, cusp_B, (nucl_num, mo_tot_num) ] implicit none BEGIN_DOC ! Equations to solve : A.C = B @@ -30,20 +28,25 @@ BEGIN_PROVIDER [ double precision, cusp_C, (nucl_num, mo_tot_num) ] do i=1,mo_tot_num do A=1,nucl_num - cusp_C(A,i) = mo_value_at_nucl(i,A) + cusp_B(A,i) = mo_value_at_nucl(i,A) enddo enddo - - integer, allocatable :: ipiv(:) - allocate ( ipiv(nucl_num) ) - call dgegv(nucl_num, mo_tot_num, cusp_A, size(cusp_A,1), & - ipiv, cusp_C, size(cusp_C,1), info) - deallocate (ipiv) - - if (info /= 0) then - print *, 'Cusp : linear solve failed' - stop -1 - endif - +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, cusp_C, (nucl_num, mo_tot_num) ] + implicit none + BEGIN_DOC + ! Equations to solve : A.C = B + END_DOC + + double precision, allocatable :: AF(:,:) + integer :: info + allocate ( AF(nucl_num,nucl_num) ) + + call get_pseudo_inverse(cusp_A,nucl_num,nucl_num,AF,size(AF,1)) + call dgemm('N','N',nucl_num,mo_tot_num,nucl_num,1.d0, & + AF,size(AF,1), cusp_B, size(cusp_B,1), 0.d0, cusp_C, size(cusp_C,1)) + END_PROVIDER diff --git a/plugins/Hartree_Fock_SlaterDressed/SCF_dressed.irp.f b/plugins/Hartree_Fock_SlaterDressed/SCF_dressed.irp.f new file mode 100644 index 00000000..ee49405e --- /dev/null +++ b/plugins/Hartree_Fock_SlaterDressed/SCF_dressed.irp.f @@ -0,0 +1,69 @@ +program scf + BEGIN_DOC +! Produce `Hartree_Fock` MO orbital with Slater cusp dressing +! output: mo_basis.mo_tot_num mo_basis.mo_label mo_basis.ao_md5 mo_basis.mo_coef mo_basis.mo_occ +! output: hartree_fock.energy +! optional: mo_basis.mo_coef + END_DOC + call check_mos + call debug + call run +end + +subroutine check_mos + implicit none + BEGIN_DOC +! Create a MO guess if no MOs are present in the EZFIO directory + END_DOC + logical :: exists + PROVIDE ezfio_filename + call ezfio_has_mo_basis_mo_coef(exists) + if (.not.exists) then + print *, 'Please run SCF first' + stop + endif +end + +subroutine debug + implicit none + integer :: i + print *, 'A' + do i=1,nucl_num + print *, i, cusp_A(1:nucl_num, i) + enddo + print *, 'B' + do i=1,mo_tot_num + print *, i, cusp_B(1:nucl_num, i) + enddo + print *, 'C' + do i=1,mo_tot_num + print *, i, cusp_C(1:nucl_num, i) + enddo +end + +subroutine run + + BEGIN_DOC +! Run SCF calculation + END_DOC + + use bitmasks + implicit none + + double precision :: SCF_energy_before,SCF_energy_after,diag_H_mat_elem + double precision :: EHF + integer :: i_it, i, j, k + + EHF = HF_energy + + mo_label = "CuspDressed" + + call ezfio_set_Hartree_Fock_SlaterDressed_slater_coef_ezfio(cusp_B) +! Choose SCF algorithm + + +! call Roothaan_Hall_SCF + +end + + diff --git a/plugins/Hartree_Fock_SlaterDressed/at_nucl.irp.f b/plugins/Hartree_Fock_SlaterDressed/at_nucl.irp.f index 48179628..60c21fb1 100644 --- a/plugins/Hartree_Fock_SlaterDressed/at_nucl.irp.f +++ b/plugins/Hartree_Fock_SlaterDressed/at_nucl.irp.f @@ -58,7 +58,7 @@ BEGIN_PROVIDER [ double precision , slater_value_at_nucl, (nucl_num,nucl_num) ] expo = slater_expo(i)*slater_expo(i)*((x*x) + (y*y) + (z*z)) if (expo > 160.d0) cycle expo = dsqrt(expo) - slater_value_at_nucl(i,k) = dexp(-expo) + slater_value_at_nucl(i,k) = dexp(-expo) * slater_normalization(i) enddo enddo END_PROVIDER diff --git a/plugins/Hartree_Fock_SlaterDressed/integrals.irp.f b/plugins/Hartree_Fock_SlaterDressed/integrals.irp.f index d49bb45d..f6208e21 100644 --- a/plugins/Hartree_Fock_SlaterDressed/integrals.irp.f +++ b/plugins/Hartree_Fock_SlaterDressed/integrals.irp.f @@ -1,15 +1,17 @@ !***************************************************************************** -subroutine GauSlaOverlap(expGau,cGau,aGau,expSla,cSla) - -! Compute the overlap integral between a Gaussian function -! with arbitrary angular momemtum and a s-type Slater function - +subroutine GauSlaOverlap(expGau,cGau,aGau,expSla,cSla,result) implicit none + BEGIN_DOC + ! Compute the overlap integral between a Gaussian function + ! with arbitrary angular momemtum and a s-type Slater function + END_DOC + ! Input variables double precision,intent(in) :: expGau,expSla double precision,intent(in) :: cGau(3),cSla(3) integer,intent(in) :: aGau(3) + double precision,intent(out) :: result ! Final value of the integrals double precision :: ss,ps,ds @@ -82,13 +84,38 @@ subroutine GauSlaOverlap(expGau,cGau,aGau,expSla,cSla) dxzs = AxBx*AzBz*ds dyzs = AyBy*AzBz*ds -! Print result - write(*,'(A10,F16.10)') & - '(s|s) = ',ss - write(*,'(A10,F16.10,3X,A10,F16.10,3X,A10,F16.10)') & - '(px|s) = ',pxs,'(py|s) = ',pys,'(pz|s) = ',pzs - write(*,'(A10,F16.10,3X,A10,F16.10,3X,A10,F16.10,3X,A10,F16.10,3X,A10,F16.10,3X,A10,F16.10)') & - '(dx2|s) = ',dxxs,'(dy2|s) = ',dyys,'(dz2|s) = ',dzzs,'(dxy|s) = ',dxys,'(dxz|s) = ',dxzs,'(dyz|s) = ',dyzs + select case (sum(aGau)) + case (0) + result = ss + + case (1) + if (aGau(1) == 1) then + result = pxs + else if (aGau(2) == 1) then + result = pys + else if (aGau(3) == 1) then + result = pzs + endif + + case (2) + if (aGau(1) == 2) then + result = dxxs + else if (aGau(2) == 2) then + result = dyys + else if (aGau(3) == 2) then + result = dzzs + else if (aGau(1)+aGau(2) == 2) then + result = dxys + else if (aGau(1)+aGau(3) == 2) then + result = dxzs + else if (aGau(2)+aGau(3) == 2) then + result = dyzs + endif + + case default + stop 'GauSlaOverlap not implemented' + + end select end !***************************************************************************** @@ -97,11 +124,13 @@ end !***************************************************************************** subroutine GauSlaKinetic(expGau,cGau,aGau,expSla,cSla) -! Compute the kinetic energy integral between a Gaussian function -! with arbitrary angular momemtum and a s-type Slater function - implicit none + BEGIN_DOC + ! Compute the kinetic energy integral between a Gaussian function + ! with arbitrary angular momemtum and a s-type Slater function + END_DOC + ! Input variables double precision,intent(in) :: expGau,expSla double precision,intent(in) :: cGau(3),cSla(3) @@ -195,11 +224,13 @@ end !***************************************************************************** subroutine GauSlaNuclear(expGau,cGau,aGau,expSla,cSla,ZNuc,cNuc) -! Compute the nuclear attraction integral between a Gaussian function -! with arbitrary angular momemtum and a s-type Slater function - implicit none + BEGIN_DOC + ! Compute the nuclear attraction integral between a Gaussian function + ! with arbitrary angular momemtum and a s-type Slater function + END_DOC + ! Input variables double precision,intent(in) :: expGau,expSla double precision,intent(in) :: cGau(3),cSla(3) @@ -242,7 +273,8 @@ subroutine GauSlaNuclear(expGau,cGau,aGau,expSla,cSla,ZNuc,cNuc) end !***************************************************************************** double precision function BoysF0(t) - + implicit none + double precision, intent(in) :: t double precision :: pi pi = 4d0*atan(1d0) @@ -257,4 +289,35 @@ end !***************************************************************************** +BEGIN_PROVIDER [ double precision, GauSlaOverlap_matrix, (ao_num, nucl_num) ] + implicit none + BEGIN_DOC + ! overlap matrix + END_DOC + integer :: i,j,k + double precision :: cGau(3) + double precision :: cSla(3) + double precision :: expSla, res, expGau + integer :: aGau(3) + + do k=1,nucl_num + cSla(1:3) = nucl_coord_transp(1:3,k) + expSla = slater_expo(k) + + do i=1,ao_num + cGau(1:3) = nucl_coord_transp(1:3, ao_nucl(i)) + aGau(1:3) = ao_power(i,1:3) + GauSlaOverlap_matrix(i,k) = 0.d0 + + do j=1,ao_prim_num(i) + expGau = ao_expo_ordered_transp(j,i) + call GauSlaOverlap(expGau,cGau,aGau,expSla,cSla,res) + GauSlaOverlap_matrix(i,k) += ao_coef_normalized_ordered_transp(j,i) * res + enddo + + enddo + + enddo + +END_PROVIDER diff --git a/plugins/Hartree_Fock_SlaterDressed/slater.irp.f b/plugins/Hartree_Fock_SlaterDressed/slater.irp.f new file mode 100644 index 00000000..8a11b9b1 --- /dev/null +++ b/plugins/Hartree_Fock_SlaterDressed/slater.irp.f @@ -0,0 +1,43 @@ +BEGIN_PROVIDER [ double precision, slater_expo, (nucl_num) ] + implicit none + BEGIN_DOC + ! Exponents of the Slater functions + END_DOC + logical :: exists + call ezfio_has_Hartree_Fock_SlaterDressed_slater_expo_ezfio(exists) + if (exists) then + slater_expo(1:nucl_num) = slater_expo_ezfio(1:nucl_num) + else + slater_expo(1:nucl_num) = nucl_charge(1:nucl_num) + call ezfio_set_Hartree_Fock_SlaterDressed_slater_expo_ezfio(slater_expo) + endif +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, slater_coef, (nucl_num,mo_tot_num) ] + implicit none + BEGIN_DOC + ! Exponents of the Slater functions + END_DOC + logical :: exists + slater_coef = 0.d0 + call ezfio_has_Hartree_Fock_SlaterDressed_slater_coef_ezfio(exists) + if (exists) then + slater_coef = slater_coef_ezfio + else + call ezfio_set_Hartree_Fock_SlaterDressed_slater_coef_ezfio(slater_coef) + endif +END_PROVIDER + +BEGIN_PROVIDER [ double precision, slater_normalization, (nucl_num) ] + implicit none + BEGIN_DOC + ! Normalization of Slater functions : sqrt(expo^3/pi) + END_DOC + integer :: i + do i=1,nucl_num + slater_normalization(i) = dsqrt( slater_expo(i)**3/dacos(-1.d0) ) + enddo + +END_PROVIDER + diff --git a/plugins/mrcepa0/dressing.irp.f b/plugins/mrcepa0/dressing.irp.f index 7cfe4fa8..9524716b 100644 --- a/plugins/mrcepa0/dressing.irp.f +++ b/plugins/mrcepa0/dressing.irp.f @@ -252,25 +252,28 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen else if (perturbative_triples) then ! Linked - call get_delta_e_dyall_general_mp(psi_ref(1,1,i_I),tq(1,1,i_alpha),Delta_E_inv) - hka = hij_cache(idx_alpha(k_sd)) - do i_state=1,N_states - ASSERT (Delta_E_inv(i_state) < 0.d0) - dka(i_state) = hka / Delta_E_inv(i_state) - enddo + if (dabs(hka) > 1.d-12) then + call get_delta_e_dyall_general_mp(psi_ref(1,1,i_I),tq(1,1,i_alpha),Delta_E_inv) + + do i_state=1,N_states + ASSERT (Delta_E_inv(i_state) < 0.d0) + dka(i_state) = hka / Delta_E_inv(i_state) + enddo + endif endif if (perturbative_triples.and. (degree2 == 1) ) then - call get_delta_e_dyall_general_mp(psi_ref(1,1,i_I),tq(1,1,i_alpha),Delta_E_inv) call i_h_j(psi_ref(1,1,i_I),tmp_det,Nint,hka) hka = hij_cache(idx_alpha(k_sd)) - hka - - do i_state=1,N_states - ASSERT (Delta_E_inv(i_state) < 0.d0) - dka(i_state) = hka / Delta_E_inv(i_state) - enddo + if (dabs(hka) > 1.d-12) then + call get_delta_e_dyall_general_mp(psi_ref(1,1,i_I),tq(1,1,i_alpha),Delta_E_inv) + do i_state=1,N_states + ASSERT (Delta_E_inv(i_state) < 0.d0) + dka(i_state) = hka / Delta_E_inv(i_state) + enddo + endif endif