mirror of
https://gitlab.com/scemama/qp_plugins_scemama.git
synced 2024-12-31 08:35:46 +01:00
Extract CC amplitudes
This commit is contained in:
parent
dd65b2677b
commit
ab8ff0c2ba
236
stable/utilities/extract_amplitudes.irp.f
Normal file
236
stable/utilities/extract_amplitudes.irp.f
Normal file
@ -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
|
||||
Loading…
Reference in New Issue
Block a user