From 6778974ae566d1d1d7fdc1bfba69c28cef1398d8 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 10 Sep 2019 13:35:20 +0200 Subject: [PATCH] Guess amplitudes OK --- devel/cc/amplitude_guess.irp.f | 41 ++++++++++++++++++++-- devel/cc/{CC.irp.f => cc.irp.f} | 0 devel/cc/spatial_to_spin_ERI.irp.f | 42 ++++++++++++++--------- stable/utilities/extract_amplitudes.irp.f | 38 +++++++++++++------- 4 files changed, 90 insertions(+), 31 deletions(-) rename devel/cc/{CC.irp.f => cc.irp.f} (100%) diff --git a/devel/cc/amplitude_guess.irp.f b/devel/cc/amplitude_guess.irp.f index 986cf23..354efd5 100644 --- a/devel/cc/amplitude_guess.irp.f +++ b/devel/cc/amplitude_guess.irp.f @@ -19,6 +19,7 @@ BEGIN_PROVIDER [ double precision, t1_guess, (spin_occ_num,spin_vir_num) ] read(iunit,*,err=10) i, a, amplitude i = 2*i-1 a = 2*a-1 - spin_occ_num +print '(I4,I4,F20.10,X,F20.10)', i,a,amplitude,t1_guess(i,a) t1_guess(i,a) = amplitude enddo 10 continue @@ -26,10 +27,14 @@ BEGIN_PROVIDER [ double precision, t1_guess, (spin_occ_num,spin_vir_num) ] read(iunit,*,end=20) i, a, amplitude i = 2*i a = 2*a - spin_occ_num +print '(I4,I4,F20.10,X,F20.10)', i,a,amplitude,t1_guess(i,a) t1_guess(i,a) = amplitude enddo 20 continue close(iunit) + else if (cc_guess == 2) then + call random_number(t1_guess) + t1_guess *= 1.d-3 endif END_PROVIDER @@ -49,6 +54,7 @@ BEGIN_PROVIDER [ double precision, t2_guess, (spin_occ_num,spin_occ_num,spin_vir integer :: i, j, a, b double precision :: amplitude + t2_guess(:,:,:,:) = 0.d0 iunit = getunitandopen('t2','r') read(iunit,*) do @@ -57,28 +63,57 @@ BEGIN_PROVIDER [ double precision, t2_guess, (spin_occ_num,spin_occ_num,spin_vir j = 2*j-1 a = 2*a-1 - spin_occ_num b = 2*b-1 - spin_occ_num - t2_guess(i,j,a,b) = amplitude +print '(I4,I4,I4,I4,F20.10,X,F20.10)', i,j,a,b,amplitude,t2_guess(i,j,a,b) + t2_guess(i,j,a,b) = amplitude enddo 10 continue +print *, '' do read(iunit,*,err=20) i, j, a, b, amplitude i = 2*i j = 2*j a = 2*a - spin_occ_num b = 2*b - spin_occ_num - t2_guess(i,j,a,b) = amplitude +print '(I4,I4,I4,I4,F20.10,X,F20.10)', i,j,a,b,amplitude,t2_guess(i,j,a,b) + t2_guess(i,j,a,b) = amplitude enddo 20 continue +print *, '' do read(iunit,*,end=30) i, j, a, b, amplitude i = 2*i-1 j = 2*j a = 2*a-1 - spin_occ_num b = 2*b - spin_occ_num - t2_guess(i,j,a,b) = amplitude +print '(I4,I4,I4,I4,F20.10,X,F20.10)', i,j,a,b,amplitude,t2_guess(i,j,a,b) + t2_guess(i,j,a,b) = amplitude + t2_guess(i,j,b,a) = -amplitude + + i = i+1 + j = j-1 + a = a+1 + b = b-1 +print '(I4,I4,I4,I4,F20.10,X,F20.10)', i,j,a,b,amplitude,t2_guess(i,j,a,b) + t2_guess(i,j,a,b) = amplitude + t2_guess(i,j,b,a) = -amplitude enddo 30 continue close(iunit) + print *, 'Non-zero amplitudes:' + do i=1,spin_occ_num + do j=1,spin_occ_num + do a=1,spin_vir_num + do b=1,spin_vir_num + if (dabs(t2_guess(i,j,a,b)) > 1.d-16) then + print *, i,j,a,b,t2_guess(i,j,a,b) + endif + enddo + enddo + enddo + enddo + else if (cc_guess == 2) then + call random_number(t2_guess) + t2_guess *= 1.d-3 endif END_PROVIDER diff --git a/devel/cc/CC.irp.f b/devel/cc/cc.irp.f similarity index 100% rename from devel/cc/CC.irp.f rename to devel/cc/cc.irp.f diff --git a/devel/cc/spatial_to_spin_ERI.irp.f b/devel/cc/spatial_to_spin_ERI.irp.f index e09c0da..c8bb3af 100644 --- a/devel/cc/spatial_to_spin_ERI.irp.f +++ b/devel/cc/spatial_to_spin_ERI.irp.f @@ -21,24 +21,34 @@ BEGIN_PROVIDER [ double precision, dbERI, (spin_mo_num,spin_mo_num,spin_mo_num,s do r=1,spin_mo_num do q=1,spin_mo_num do p=1,spin_mo_num + pqrs = 0.d0 + pqsr = 0.d0 + if ( ( iand(p,1) == iand(r,1) ) .and. & + ( iand(q,1) == iand(s,1) ) ) then +! pqrs = Kronecker_delta(mod(p,2),mod(r,2)) & +! * Kronecker_delta(mod(q,2),mod(s,2)) & +! * get_two_e_integral( & + pqrs = get_two_e_integral( & + (p+1)/2, & + (q+1)/2, & + (r+1)/2, & + (s+1)/2, & + mo_two_e_integrals_in_map) + endif - pqrs = Kronecker_delta(mod(p,2),mod(r,2)) & - * Kronecker_delta(mod(q,2),mod(s,2)) & - * get_two_e_integral( & - (p+1)/2, & - (q+1)/2, & - (r+1)/2, & - (s+1)/2, & - mo_two_e_integrals_in_map) + if ( ( iand(p,1) == iand(s,1) ) .and. & + ( iand(q,1) == iand(r,1) ) ) then - pqsr = Kronecker_delta(mod(p,2),mod(s,2)) & - * Kronecker_delta(mod(q,2),mod(r,2)) & - * get_two_e_integral( & - (p+1)/2, & - (q+1)/2, & - (s+1)/2, & - (r+1)/2, & - mo_two_e_integrals_in_map) +! pqsr = Kronecker_delta(mod(p,2),mod(s,2)) & +! * Kronecker_delta(mod(q,2),mod(r,2)) & +! * get_two_e_integral( & + pqsr = get_two_e_integral( & + (p+1)/2, & + (q+1)/2, & + (s+1)/2, & + (r+1)/2, & + mo_two_e_integrals_in_map) + endif dbERI(p,q,r,s) = pqrs - pqsr enddo diff --git a/stable/utilities/extract_amplitudes.irp.f b/stable/utilities/extract_amplitudes.irp.f index a13680b..e025bb8 100644 --- a/stable/utilities/extract_amplitudes.irp.f +++ b/stable/utilities/extract_amplitudes.irp.f @@ -84,15 +84,24 @@ subroutine run endif if ( (s1 == 1).and.(s2 == 1) ) then - t2_aa(h1,h2,p1,p2) += phase * psi_coef(k,istate) * norm + t2_aa(h1,h2,p1,p2) += phase * psi_coef(k,istate) * norm * 0.5d0 + t2_aa(h1,h2,p2,p1) -= phase * psi_coef(k,istate) * norm * 0.5d0 + t2_aa(h2,h1,p2,p1) += phase * psi_coef(k,istate) * norm * 0.5d0 + t2_aa(h2,h1,p1,p2) -= phase * psi_coef(k,istate) * norm * 0.5d0 else if ( (s1 == 1).and.(s2 == 2) ) then - t2_ab(h1,h2,p1,p2) += phase * psi_coef(k,istate) * norm + t2_ab(h1,h2,p1,p2) += phase * psi_coef(k,istate) * norm * 0.5d0 + t2_ab(h2,h1,p2,p1) += phase * psi_coef(k,istate) * norm * 0.5d0 else if ( (s1 == 2).and.(s2 == 1) ) then - t2_ab(h2,h1,p2,p1) += phase * psi_coef(k,istate) * norm + print *, irp_here, ': Bug!' + stop -1 else if ( (s1 == 2).and.(s2 == 2) ) then - t2_bb(h1,h2,p1,p2) += phase * psi_coef(k,istate) * norm + t2_bb(h1,h2,p1,p2) += phase * psi_coef(k,istate) * norm * 0.5d0 + t2_bb(h1,h2,p2,p1) -= phase * psi_coef(k,istate) * norm * 0.5d0 + t2_bb(h2,h1,p2,p1) += phase * psi_coef(k,istate) * norm * 0.5d0 + t2_bb(h2,h1,p1,p2) -= phase * psi_coef(k,istate) * norm * 0.5d0 else print *, irp_here, ': Bug!' + stop -2 endif end select @@ -103,9 +112,14 @@ subroutine run do a=elec_alpha_num+1, mo_num do j=1,elec_alpha_num do i=1,elec_alpha_num - t2_aa(i,j,a,b) -= t1_a(i,a)*t1_a(j,b) - t2_bb(i,j,a,b) -= t1_b(i,a)*t1_b(j,b) - t2_ab(i,j,a,b) -= t1_a(i,a)*t1_b(j,b) + if ( (i == j).or.(a == b)) then + t2_aa(i,j,a,b) = 0.d0 + t2_bb(i,j,a,b) = 0.d0 + else + t2_aa(i,j,a,b) -= t1_a(i,a)*t1_a(j,b) + t2_bb(i,j,a,b) -= t1_b(i,a)*t1_b(j,b) + endif + t2_ab(i,j,a,b) -= 0.5d0* (t1_a(i,a)*t1_b(j,b) + t1_b(i,a)*t1_a(j,b)) enddo enddo enddo @@ -182,21 +196,21 @@ subroutine run double precision :: e_cor double precision, external :: get_two_e_integral - e_cor = 0.d0 + e_cor = 0.d0 do b=elec_alpha_num+1, mo_num do a=elec_alpha_num+1, mo_num do j=1,elec_alpha_num do i=1,elec_alpha_num ! The t1 terms are zero because HF - e_cor = e_cor + ( & + e_cor = e_cor + 0.5d0*( & t2_aa(i,j,a,b) + t2_bb(i,j,a,b) + & t1_a(i,a) * t1_a(j,b) + & t1_b(i,a) * t1_b(j,b) ) * ( & get_two_e_integral(i,j,a,b,mo_integrals_map) - & get_two_e_integral(i,j,b,a,mo_integrals_map) ) - e_cor = e_cor + ( & - t2_ab(i,j,a,b) + & - t1_a(i,a) * t1_b(j,b) ) * & + e_cor = e_cor + 1.0d0 * ( & + t2_ab(i,j,a,b) + & + t1_a(i,a) * t1_b(j,b) ) * & get_two_e_integral(i,j,a,b,mo_integrals_map) enddo enddo