1
0
mirror of https://gitlab.com/scemama/qp_plugins_scemama.git synced 2024-12-22 20:34:05 +01:00
qp_plugins_scemama/devel/svdwf/printSQ_ab_T_gd_v0.irp.f
2021-11-02 16:18:07 +01:00

403 lines
13 KiB
Fortran

program printSQ_ab_T_gd_v0
implicit none
BEGIN_DOC
! perturbative approach to build psi_postsvd
END_DOC
read_wf = .True.
TOUCH read_wf
PROVIDE N_int
call run()
end
subroutine run
USE OMP_LIB
USE bitmasks
implicit none
integer(bit_kind) :: det1(N_int,2), det2(N_int,2)
integer :: degree, i_state
double precision :: h12
integer :: i, j, k, l
double precision, allocatable :: T(:,:,:,:)
double precision, allocatable :: Uref(:,:), Dref(:), Vtref(:,:), Aref(:,:), Vref(:,:)
double precision :: ti, tf
integer :: nb_taches
!$OMP PARALLEL
nb_taches = OMP_GET_NUM_THREADS()
!$OMP END PARALLEL
call wall_time(ti)
i_state = 1
det1(:,1) = psi_det_alpha_unique(:,1)
det2(:,1) = psi_det_alpha_unique(:,1)
det1(:,2) = psi_det_beta_unique(:,1)
det2(:,2) = psi_det_beta_unique(:,1)
call get_excitation_degree_spin(det1(1,1),det2(1,1),degree,N_int)
call get_excitation_degree(det1,det2,degree,N_int)
call i_H_j(det1, det2, N_int, h12)
! ---------------------------------------------------------------------------------------
! ---------------------------------------------------------------------------------------
! ---------------------------------------------------------------------------------------
! construct the initial CI matrix
print *, ' ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~'
print *, ' CI matrix :', n_det_alpha_unique,'x',n_det_beta_unique
print *, ' N det :', N_det
print *, ' ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~'
allocate( Aref(n_det_alpha_unique,n_det_beta_unique) )
Aref(:,:) = 0.d0
do k = 1, N_det
i = psi_bilinear_matrix_rows(k)
j = psi_bilinear_matrix_columns(k)
Aref(i,j) = psi_bilinear_matrix_values(k,i_state)
enddo
allocate( Uref(n_det_alpha_unique,n_det_alpha_unique) )
allocate( Dref(min(n_det_alpha_unique,n_det_beta_unique)) )
allocate( Vtref(n_det_beta_unique,n_det_beta_unique) )
call svd_s(Aref, size(Aref,1), Uref, size(Uref,1), Dref, Vtref &
, size(Vtref,1), n_det_alpha_unique, n_det_beta_unique)
print *, " ---------- Full SVD is performed ----------------- "
allocate( Vref(n_det_beta_unique,n_det_beta_unique) )
do l = 1, n_det_beta_unique
do i = 1, n_det_beta_unique
Vref(i,l) = Vtref(l,i)
enddo
enddo
deallocate( Aref , Vtref , Dref )
! ---------------------------------------------------------------------------------------
! ---------------------------------------------------------------------------------------
! ---------------------------------------------------------------------------------------
! ---------------------------------------------------------------------------------------
allocate( T(n_det_alpha_unique,n_det_beta_unique,n_det_alpha_unique,n_det_beta_unique) )
call const_2b(T)
call const_ab_T_gd(Uref, Vref, T)
deallocate( T )
deallocate( Uref , Vref )
call wall_time(tf)
print *, ' ___________________________________________________________________'
print *, ' '
!print *, " Execution avec ", nb_taches, " threads"
print *, " Execution avec 1 threads"
print *, " total elapsed time (min) = ", (tf-ti)/60.d0
print *, ' ___________________________________________________________________'
end
!/-\! !/-\! !/-\! !/-\! !/-\! !/-\! !/-\! !/-\! !/-\!
!---! !---! !---! !---! !---! !---! !---! !---! !---!
! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !
!___! !___! !___! !___! !___! !___! !___! !___! !___!
!___! !___! !___! !___! !___! !___! !___! !___! !___!
!/-\! !/-\! !/-\! !/-\! !/-\! !/-\! !/-\! !/-\! !/-\!
!---! !---! !---! !---! !---! !---! !---! !---! !---!
! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !
!___! !___! !___! !___! !___! !___! !___! !___! !___!
!___! !___! !___! !___! !___! !___! !___! !___! !___!
!/-\! !/-\! !/-\! !/-\! !/-\! !/-\! !/-\! !/-\! !/-\!
!---! !---! !---! !---! !---! !---! !---! !---! !---!
! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !
!___! !___! !___! !___! !___! !___! !___! !___! !___!
!___! !___! !___! !___! !___! !___! !___! !___! !___!
!/-\! !/-\! !/-\! !/-\! !/-\! !/-\! !/-\! !/-\! !/-\!
!---! !---! !---! !---! !---! !---! !---! !---! !---!
! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !
!___! !___! !___! !___! !___! !___! !___! !___! !___!
!___! !___! !___! !___! !___! !___! !___! !___! !___!
subroutine const_2b(T)
USE bitmasks
implicit none
double precision, external :: get_two_e_integral
double precision, intent(out) :: T(n_det_alpha_unique,n_det_beta_unique,n_det_alpha_unique,n_det_beta_unique )
integer :: na, nb
integer :: i, k, j, l
integer(bit_kind) :: psi_ij(N_int,2), psi_kl(N_int,2)
double precision :: phase
integer :: degree, h1, h2, p1, p2, s1, s2, e1, e2
integer :: ii, jj
integer :: exc(0:2,2,2)
integer :: occ(N_int*bit_kind_size,2), n_occ_alpha
double precision :: two_body_fact
na = n_det_alpha_unique
nb = n_det_beta_unique
T(:,:,:,:) = 0.d0
! -----------------------------------------------------------------------------------------------------------------
do i = 1, na
psi_ij(1:N_int,1) = psi_det_alpha_unique(1:N_int,i)
do j = 1, nb
psi_ij(1:N_int,2) = psi_det_beta_unique(1:N_int,j)
call bitstring_to_list(psi_ij(1,1), occ(1,1), n_occ_alpha, N_int)
call bitstring_to_list(psi_ij(1,2), occ(1,2), n_occ_alpha, N_int)
do k = 1, na
psi_kl(1:N_int,1) = psi_det_alpha_unique(1:N_int,k)
do l = 1, nb
psi_kl(1:N_int,2) = psi_det_beta_unique(1:N_int,l)
call get_excitation_degree(psi_ij, psi_kl, degree, N_int)
two_body_fact = 0.d0
if(degree .eq. 2) then
call get_double_excitation(psi_ij, psi_kl, exc, phase, N_int)
call decode_exc(exc, degree, h1, p1, h2, p2, s1, s2)
select case(s1+s2)
case(2,4)
two_body_fact += phase * get_two_e_integral(h1, h2, p1, p2, mo_integrals_map)
two_body_fact -= phase * get_two_e_integral(h1, h2, p2, p1, mo_integrals_map)
case(3)
two_body_fact += 0.5d0 * phase * get_two_e_integral(h1, h2, p1, p2, mo_integrals_map)
two_body_fact += 0.5d0 * phase * get_two_e_integral(h2, h1, p2, p1, mo_integrals_map)
end select
else if(degree .eq. 1) then
call get_single_excitation(psi_ij, psi_kl, exc, phase, N_int)
call decode_exc(exc, degree, h1, p1, h2, p2, s1, s2)
select case(s1)
case(1)
do ii = 1, elec_alpha_num
p2 = occ(ii,1)
h2 = p2
two_body_fact += 0.5d0 * phase * get_two_e_integral(h1, h2, p1, p2, mo_integrals_map)
two_body_fact -= 0.5d0 * phase * get_two_e_integral(h1, h2, p2, p1, mo_integrals_map)
two_body_fact += 0.5d0 * phase * get_two_e_integral(h2, h1, p2, p1, mo_integrals_map)
two_body_fact -= 0.5d0 * phase * get_two_e_integral(h2, h1, p1, p2, mo_integrals_map)
enddo
do ii = 1, elec_beta_num
p2 = occ(ii,2)
h2 = p2
two_body_fact += 0.5d0 * phase * get_two_e_integral(h1, h2, p1, p2, mo_integrals_map)
two_body_fact += 0.5d0 * phase * get_two_e_integral(h2, h1, p2, p1, mo_integrals_map)
enddo
case(2)
do ii = 1, elec_alpha_num
p2 = occ(ii,1)
h2 = p2
two_body_fact += 0.5d0 * phase * get_two_e_integral(h1, h2, p1, p2, mo_integrals_map)
two_body_fact += 0.5d0 * phase * get_two_e_integral(h2, h1, p2, p1, mo_integrals_map)
enddo
do ii = 1, elec_beta_num
p2 = occ(ii,2)
h2 = p2
two_body_fact += 0.5d0 * phase * get_two_e_integral(h1, h2, p1, p2, mo_integrals_map)
two_body_fact -= 0.5d0 * phase * get_two_e_integral(h1, h2, p2, p1, mo_integrals_map)
two_body_fact += 0.5d0 * phase * get_two_e_integral(h2, h1, p2, p1, mo_integrals_map)
two_body_fact -= 0.5d0 * phase * get_two_e_integral(h2, h1, p1, p2, mo_integrals_map)
enddo
end select
else if(degree .eq. 0) then
do ii = 1, elec_alpha_num
e1 = occ(ii,1)
do jj = 1, elec_alpha_num
e2 = occ(jj,1)
two_body_fact += 0.5d0 * get_two_e_integral(e1, e2, e1, e2, mo_integrals_map)
two_body_fact -= 0.5d0 * get_two_e_integral(e1, e2, e2, e1, mo_integrals_map)
enddo
do jj = 1, elec_beta_num
e2 = occ(jj,2)
two_body_fact += 0.5d0 * get_two_e_integral(e1, e2, e1, e2, mo_integrals_map)
two_body_fact += 0.5d0 * get_two_e_integral(e2, e1, e2, e1, mo_integrals_map)
enddo
enddo
do ii = 1, elec_beta_num
e1 = occ(ii,2)
do jj = 1, elec_beta_num
e2 = occ(jj,2)
two_body_fact += 0.5d0 * get_two_e_integral(e1, e2, e1, e2, mo_integrals_map)
two_body_fact -= 0.5d0 * get_two_e_integral(e1, e2, e2, e1, mo_integrals_map)
enddo
enddo
end if
T(i,j,k,l) = two_body_fact
enddo
enddo
enddo
enddo
! -----------------------------------------------------------------------------------------------------------------
return
end subroutine const_2b
subroutine const_ab_T_gd(Uref, Vref, T)
implicit none
double precision, intent(in) :: Uref(n_det_alpha_unique,n_det_alpha_unique)
double precision, intent(in) :: Vref(n_det_beta_unique ,n_det_beta_unique)
double precision, intent(in) :: T(n_det_alpha_unique,n_det_beta_unique,n_det_alpha_unique,n_det_beta_unique )
integer :: na, nb
integer :: i, j, k, l, ii, jj, kk, ll
double precision :: h12, T_abgd, sum_i, sum_j, sum_k, sum_l
double precision, allocatable :: tmp1(:,:,:,:), tmp2(:,:,:,:)
double precision :: t1, t2
print *, ""
print *, " start const_ab_T_gd"
call wall_time(t1)
na = n_det_alpha_unique
nb = n_det_beta_unique
allocate( tmp1(na,nb,na,nb) )
tmp1(:,:,:,:) = 0.d0
do ii = 1, na
do l = 1, nb
do j = 1, nb
do k = 1, na
sum_i = 0.d0
do i = 1, na
sum_i = sum_i + T(i,j,k,l) * Uref(i,ii)
enddo
tmp1(ii,j,k,l) = sum_i
enddo
enddo
enddo
enddo
allocate( tmp2(na,nb,na,nb) )
tmp2(:,:,:,:) = 0.d0
do ii = 1, na
do jj = 1, nb
do k = 1, na
do l = 1, nb
sum_j = 0.d0
do j = 1, nb
sum_j = sum_j + tmp1(ii,j,k,l) * Vref(j,jj)
enddo
tmp2(ii,jj,k,l) = sum_j
enddo
enddo
enddo
enddo
deallocate(tmp1)
allocate( tmp1(na,nb,na,nb) )
tmp1(:,:,:,:) = 0.d0
do ii = 1, na
do jj = 1, nb
do kk = 1, na
do l = 1, nb
sum_k = 0.d0
do k = 1, na
sum_k = sum_k + tmp2(ii,jj,k,l) * Uref(k,kk)
enddo
tmp1(ii,jj,kk,l) = sum_k
enddo
enddo
enddo
enddo
deallocate(tmp2)
allocate( tmp2(na,nb,na,nb) )
tmp2(:,:,:,:) = 0.d0
do ii = 1, na
do jj = 1, nb
do kk = 1, na
do ll = 1, nb
sum_l = 0.d0
do l = 1, nb
sum_l = sum_l + tmp1(ii,jj,kk,l) * Vref(l,ll)
enddo
tmp2(ii,jj,kk,ll) = sum_l
enddo
enddo
enddo
enddo
deallocate(tmp1)
open(UNIT=11, FILE="ab_T_gd_v0.dat", ACTION="WRITE")
do ii = 1, na
do jj = 1, nb
do kk = 1, na
do ll = 1, nb
T_abgd = tmp2(ii,jj,kk,ll)
write(11, '((F15.8))') T_abgd
enddo
enddo
enddo
enddo
close(11)
deallocate(tmp2)
call wall_time(t2)
print *, " end const_ab_T_gd after (min) ", (t2-t1)/60.
print *, ""
return
end subroutine const_ab_T_gd