mirror of
https://gitlab.com/scemama/qp_plugins_scemama.git
synced 2024-06-02 03:15:25 +02:00
248 lines
7.5 KiB
Fortran
248 lines
7.5 KiB
Fortran
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_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
|
|
character*(32) :: buffer
|
|
|
|
call getenv('QP_STATE',buffer)
|
|
istate = 1
|
|
read(buffer,*,err=5,end=5) istate
|
|
5 continue
|
|
call write_int(6,istate,'State for amplitudes')
|
|
|
|
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_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_bb = 0.d0
|
|
|
|
integer :: i_ref
|
|
|
|
i_ref = maxloc( dabs(psi_coef(:,istate)), dim=1 )
|
|
|
|
norm = 1.d0 / psi_coef(i_ref,istate)
|
|
|
|
do k=1,N_det
|
|
|
|
call get_excitation_degree(psi_det(1,1,i_ref),psi_det(1,1,k),degree,N_int)
|
|
|
|
select case (degree)
|
|
|
|
case (1)
|
|
call get_excitation(psi_det(1,1,i_ref),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(psi_det(1,1,i_ref),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
|
|
t2_aa(h1,h2,p2,p1) -= phase * psi_coef(k,istate) * norm
|
|
t2_aa(h2,h1,p2,p1) += phase * psi_coef(k,istate) * norm
|
|
t2_aa(h2,h1,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 * 0.5d0
|
|
t2_ab(h2,h1,p2,p1) += phase * psi_coef(k,istate) * norm * 0.5d0
|
|
else if ( (s1 == 2).and.(s2 == 1) ) then
|
|
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,p2,p1) -= phase * psi_coef(k,istate) * norm
|
|
t2_bb(h2,h1,p2,p1) += phase * psi_coef(k,istate) * norm
|
|
t2_bb(h2,h1,p1,p2) -= phase * psi_coef(k,istate) * norm
|
|
else
|
|
print *, irp_here, ': Bug!'
|
|
stop -2
|
|
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
|
|
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)
|
|
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)
|
|
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
|
|
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_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
|
|
|
|
|
|
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
|
|
|
|
|
|
double precision :: E0, e_cor
|
|
double precision, external :: get_two_e_integral
|
|
|
|
e_cor = 0.d0
|
|
do b=elec_alpha_num+1, mo_num
|
|
do j=1,elec_alpha_num
|
|
! TODO : singles contributions
|
|
! do k=1,mo_num
|
|
! e_cor = e_cor + t1_a(j,b) * get_two_e_integral(k,j,k,b,mo_integrals_map)
|
|
! e_cor = e_cor + t1_b(j,b) * get_two_e_integral(k,j,k,b,mo_integrals_map)
|
|
! enddo
|
|
do a=elec_alpha_num+1, mo_num
|
|
do i=1,elec_alpha_num
|
|
e_cor = e_cor + 0.25d0*( &
|
|
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) - &
|
|
t1_a(i,b) * t1_a(j,a) - &
|
|
t1_b(i,b) * t1_b(j,a) &
|
|
) * ( &
|
|
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 + 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
|
|
enddo
|
|
enddo
|
|
e_cor = e_cor
|
|
|
|
double precision, external :: diag_h_mat_elem
|
|
E0 = diag_h_mat_elem(psi_det(1,1,i_ref),N_int) + nuclear_repulsion
|
|
print '(A,F15.10)', 'E0 : ', E0
|
|
print '(A,F15.10)', 'corr energy: ', e_cor
|
|
print '(A,F15.10)', 'total energy: ', E0 + e_cor
|
|
|
|
deallocate(t1_a,t1_b,t2_aa,t2_ab,t2_bb)
|
|
|
|
end
|