diff --git a/stable/utilities/extract_amplitudes.irp.f b/stable/utilities/extract_amplitudes.irp.f new file mode 100644 index 0000000..5c15266 --- /dev/null +++ b/stable/utilities/extract_amplitudes.irp.f @@ -0,0 +1,236 @@ +program extract_amplitudes + implicit none + BEGIN_DOC +! Print the T1 and T2 amplitudes into a file. + END_DOC + + read_wf = .True. + TOUCH read_wf + call run +end + +subroutine run + implicit none + BEGIN_DOC +! Compute t1 and T2 amplitudes + END_DOC + double precision, allocatable :: t1_a(:,:), t1_b(:,:) + double precision, allocatable :: t2_aa(:,:,:,:) + double precision, allocatable :: t2_ab(:,:,:,:) + double precision, allocatable :: t2_ba(:,:,:,:) + double precision, allocatable :: t2_bb(:,:,:,:) + double precision :: phase, norm + integer :: exc(0:2,2,2), degree + integer :: h1,h2,p1,p2,s1,s2, istate + integer :: i,j,k,l,a,b,c,d + + istate = 1 + + allocate ( & + t1_a(elec_alpha_num,mo_num), & + t1_b(elec_alpha_num,mo_num), & + t2_aa(elec_alpha_num,elec_alpha_num,mo_num,mo_num), & + t2_ab(elec_alpha_num,elec_alpha_num,mo_num,mo_num), & + t2_ba(elec_alpha_num,elec_alpha_num,mo_num,mo_num), & + t2_bb(elec_alpha_num,elec_alpha_num,mo_num,mo_num) ) + + t1_a = 0.d0 + t1_b = 0.d0 + t2_aa = 0.d0 + t2_ab = 0.d0 + t2_ba = 0.d0 + t2_bb = 0.d0 + + norm = 1.d0 / psi_coef_sorted(1,istate) + + do k=1,N_det + + call get_excitation_degree(hf_bitmask,psi_det(1,1,k),degree,N_int) + + select case (degree) + + case (1) + call get_excitation(hf_bitmask,psi_det(1,1,k),exc,degree,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + + if (h1 > elec_alpha_num) then + print *, irp_here, "h1 > elec_alpha_num" + endif + if (p1 <= elec_alpha_num) then + print *, irp_here, "p1 <= elec_alpha_num" + endif + + if (s1 == 1) then + t1_a(h1,p1) += phase * psi_coef(k,istate) * norm + else if (s1 == 2) then + t1_b(h1,p1) += phase * psi_coef(k,istate) * norm + else + print *, irp_here, ': Bug!', s1, s2 + endif + + + case (2) + call get_excitation(hf_bitmask,psi_det(1,1,k),exc,degree,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + + if (h1 > elec_alpha_num) then + print *, irp_here, "h1 > elec_alpha_num" + endif + if (p1 <= elec_alpha_num) then + print *, irp_here, "p1 <= elec_alpha_num" + endif + if (h2 > elec_alpha_num) then + print *, irp_here, "h2 > elec_alpha_num" + endif + if (p2 <= elec_alpha_num) then + print *, irp_here, "p2 <= elec_alpha_num" + endif + + if ( (s1 == 1).and.(s2 == 1) ) then + t2_aa(h1,h2,p1,p2) += phase * psi_coef(k,istate) * norm + else if ( (s1 == 1).and.(s2 == 2) ) then + t2_ab(h1,h2,p1,p2) += phase * psi_coef(k,istate) * norm + else if ( (s1 == 2).and.(s2 == 1) ) then + t2_ba(h1,h2,p1,p2) += phase * psi_coef(k,istate) * norm + else if ( (s1 == 2).and.(s2 == 2) ) then + t2_bb(h1,h2,p1,p2) += phase * psi_coef(k,istate) * norm + else + print *, irp_here, ': Bug!' + endif + + end select + + enddo + + 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 + 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_aa(i,j,a,b) += t1_a(i,b)*t1_a(j,a) +! t2_bb(i,j,a,b) += t1_b(i,b)*t1_b(j,a) + + t2_ab(i,j,a,b) -= t1_a(i,a)*t1_b(j,b) + t2_ba(i,j,a,b) -= t1_b(i,a)*t1_a(j,b) + + enddo + enddo + enddo + enddo + + + integer :: iunit, getunitandopen + + iunit = getunitandopen('t1','w') + + write(iunit,*) '# T1_a' + do i=1,elec_alpha_num + do a=elec_alpha_num+1, mo_num + if (t1_a(i,a) /= 0.d0) then + write(iunit,'(I4,X,I4,X,E20.10)') i, a, t1_a(i,a) + endif + enddo + enddo + + write(iunit,*) '# T1_b' + do i=1,elec_alpha_num + do a=elec_alpha_num+1, mo_num + if (t1_b(i,a) /= 0.d0) then + write(iunit,'(I4,X,I4,X,E20.10)') i, a, t1_b(i,a) + endif + enddo + enddo + + + iunit = getunitandopen('t2','w') + + write(iunit,*) '# T2_aa' + do i=1,elec_alpha_num + do j=1,elec_alpha_num + do a=elec_alpha_num+1, mo_num + do b=elec_alpha_num+1, mo_num + if (t2_aa(i,j,a,b) /= 0.d0) then + write(iunit,'(4(I4,X),E20.10)') i, j, a, b, t2_aa(i,j,a,b) + endif + enddo + enddo + enddo + enddo + + + write(iunit,*) '# T2_ab' + do i=1,elec_alpha_num + do j=1,elec_alpha_num + do a=elec_alpha_num+1, mo_num + do b=elec_alpha_num+1, mo_num + if (t2_ab(i,j,a,b) /= 0.d0) then + write(iunit,'(4(I4,X),E20.10)') i, j, a, b, t2_ab(i,j,a,b) + endif + enddo + enddo + enddo + enddo + + + write(iunit,*) '# T2_ba' + do i=1,elec_alpha_num + do j=1,elec_alpha_num + do a=elec_alpha_num+1, mo_num + do b=elec_alpha_num+1, mo_num + if (t2_ba(i,j,a,b) /= 0.d0) then + write(iunit,'(4(I4,X),E20.10)') i, j, a, b, t2_ba(i,j,a,b) + endif + enddo + enddo + enddo + enddo + + + write(iunit,*) '# T2_bb' + do i=1,elec_alpha_num + do j=1,elec_alpha_num + do a=elec_alpha_num+1, mo_num + do b=elec_alpha_num+1, mo_num + if (t2_bb(i,j,a,b) /= 0.d0) then + write(iunit,'(4(I4,X),E20.10)') i, j, a, b, t2_bb(i,j,a,b) + endif + enddo + enddo + enddo + enddo + + + double precision :: e_cor + double precision, external :: get_two_e_integral + + 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 + ( & + 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) + t2_ba(i,j,a,b) + & + t1_a(i,a) * t1_b(j,b) + & + t1_b(i,a) * t1_a(j,b) ) * & + get_two_e_integral(i,j,a,b,mo_integrals_map) + enddo + enddo + enddo + enddo + e_cor = e_cor + print '(A,F15.10)', 'HF energy: ', hf_energy + print '(A,F15.10)', 'corr energy: ', e_cor + print '(A,F15.10)', 'total energy: ', hf_energy + e_cor + + deallocate(t1_a,t1_b,t2_aa,t2_ab,t2_ba,t2_bb) + +end