1
0
mirror of https://gitlab.com/scemama/qp_plugins_scemama.git synced 2024-12-23 04:43:38 +01:00
qp_plugins_scemama/devel/svdwf/kl_H_kl_v2.irp.f
2021-11-02 16:18:07 +01:00

327 lines
8.4 KiB
Fortran

program kl_H_kl_v2
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
implicit none
integer(bit_kind) :: det1(N_int,2), det2(N_int,2)
integer :: degree, i_state
integer :: i, j, k, l, m, n
integer :: ii, ia, ib
double precision :: x, y, h12
double precision, allocatable :: Uref(:,:), Dref(:), Vtref(:,:), Aref(:,:), Vref(:,:)
double precision, allocatable :: Hdiag(:)
integer :: ibeg_alpha, ibeg_beta, iend_alpha, iend_beta
integer :: n_toselect, na_max, nb_max
integer, allocatable :: numalpha_toselect(:), numbeta_toselect(:)
double precision :: t_beg, t_end
integer(kind=8) :: W_tbeg, W_tend, W_tbeg_it, W_tend_it, W_ir
real(kind=8) :: W_tot_time, W_tot_time_it
integer :: nb_taches
!$OMP PARALLEL
nb_taches = OMP_GET_NUM_THREADS()
!$OMP END PARALLEL
call SYSTEM_CLOCK(COUNT=W_tbeg, COUNT_RATE=W_ir)
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 CISD matrix
print *, ' ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~'
print *, ' CISD 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
! ---------------------------------------------------------------------------------------
! ---------------------------------------------------------------------------------------
! perform a Full SVD
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 cpu_time(t_beg)
call svd_s(Aref, size(Aref,1), Uref, size(Uref,1), Dref, Vtref &
, size(Vtref,1), n_det_alpha_unique, n_det_beta_unique)
call cpu_time(t_end)
print *, " SVD is performed after (min)", (t_end-t_beg)/60.
deallocate( Aref , Dref )
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( Vtref )
ibeg_alpha = 1
iend_alpha = 100
!iend_alpha = n_det_alpha_unique
ibeg_beta = 1
iend_beta = 100
!iend_beta = n_det_beta_unique
na_max = iend_alpha - ibeg_alpha + 1
nb_max = iend_beta - ibeg_beta + 1
n_toselect = na_max * nb_max
print *, ' na_max = ', na_max
print *, ' nb_max = ', nb_max
print *, ' n_toselect = ', n_toselect
allocate( numalpha_toselect(n_toselect) , numbeta_toselect(n_toselect) )
k = 0
do i = ibeg_alpha, iend_alpha
do j = ibeg_beta, iend_beta
k = k + 1
numalpha_toselect(k) = i
numbeta_toselect (k) = j
enddo
enddo
if( k.ne.n_toselect ) then
print *, " error in numbering"
stop
endif
allocate( Hdiag(n_toselect) )
! get < u_k v_l | H | u_k v_l > for all vectors
call const_Hdiag(na_max, nb_max, n_toselect, Uref, Vref, numalpha_toselect, numbeta_toselect, Hdiag)
open( UNIT=11, FILE="klHkl_v2.dat", ACTION="WRITE")
do i = 1, n_toselect
write(11, '(2(I5,2X), 5X, E15.7)') numalpha_toselect(i), numbeta_toselect(i), Hdiag(i)
enddo
close(11)
deallocate( Uref, Vref )
deallocate( numalpha_toselect, numbeta_toselect, Hdiag )
call SYSTEM_CLOCK(COUNT=W_tend, COUNT_RATE=W_ir)
W_tot_time = real(W_tend - W_tbeg, kind=8) / real(W_ir, kind=8)
print *, ' ___________________________________________________________________'
print *, ' '
print *, " Execution avec ", nb_taches, " threads"
print *, " total elapsed time (min) = ", W_tot_time/60.d0
print *, ' ___________________________________________________________________'
end
subroutine const_Hdiag(na_max, nb_max, n_toselect, Uref, Vref, numalpha_toselect, numbeta_toselect, Hdiag)
implicit none
integer, intent(in) :: n_toselect, na_max, nb_max
integer, intent(in) :: numalpha_toselect(n_toselect), numbeta_toselect(n_toselect)
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(out) :: Hdiag(n_toselect)
integer(bit_kind) :: det1(N_int,2)
integer(bit_kind) :: det2(N_int,2)
integer :: degree, na, nb
integer :: i, j, k, l, ii, jj, m, n
double precision :: h12, xtmp
double precision, allocatable :: Hmat_diag(:,:), Vt(:,:), bl1_tmp(:,:,:)
double precision, allocatable :: Ut(:,:), tmp0(:,:,:) , Hmat_diag_tmp(:,:)
double precision :: t1, t2, t3, t4
print *, ""
print *, " start const_Hdiag"
call wall_time(t1)
na = n_det_alpha_unique
nb = n_det_beta_unique
allocate(Hmat_diag(na_max,nb_max))
Hmat_diag = 0.d0
allocate( bl1_tmp(na,na,nb_max) )
bl1_tmp = 0.d0
allocate( Vt(nb_max,nb) )
do i = 1, nb
do n = 1, nb_max
Vt(n,i) = Vref(i,n)
enddo
enddo
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(i,j,k,l,n,h12,det1,det2,degree,tmp0) &
!$OMP SHARED(na,nb,psi_det_alpha_unique,psi_det_beta_unique, &
!$OMP N_int,nb_max,Vt,bl1_tmp)
allocate( tmp0(na,na,nb_max) )
tmp0 = 0.d0
!$OMP DO
do l = 1, nb
det2(:,2) = psi_det_beta_unique(:,l)
do j = 1, nb
det1(:,2) = psi_det_beta_unique(:,j)
call get_excitation_degree_spin(det1(1,2),det2(1,2),degree,N_int)
if(degree .gt. 2) cycle
do k = 1, na
det2(:,1) = psi_det_alpha_unique(:,k)
do i = 1, na
det1(:,1) = psi_det_alpha_unique(:,i)
call get_excitation_degree(det1,det2,degree,N_int)
if(degree .gt. 2) cycle
call i_H_j(det1, det2, N_int, h12)
if( h12 .eq. 0.d0) cycle
do n = 1, nb_max
tmp0(i,k,n) += h12 * Vt(n,j) * Vt(n,l)
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP CRITICAL
do n = 1, nb_max
do k = 1, na
do i = 1, na
bl1_tmp(i,k,n) += tmp0(i,k,n)
enddo
enddo
enddo
!$OMP END CRITICAL
deallocate( tmp0 )
!$OMP END PARALLEL
deallocate(Vt)
call wall_time(t2)
print *, " end bl1_tmp after (min) ", (t2-t1)/60.
allocate( Ut(na,na_max) )
Ut(1:na,1:na_max) = Uref(1:na,1:na_max)
allocate( tmp0(na,nb_max,na_max) )
call DGEMM('T', 'N', na*nb_max, na_max, na, 1.d0, &
bl1_tmp, size(bl1_tmp,1), Ut, size(Ut,1), &
0.d0, tmp0, size(tmp0,1)*size(tmp0,2) )
deallocate( bl1_tmp )
call wall_time(t3)
print *, " end DGEMM after (min) ", (t3-t2)/60.
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(k,m,n,Hmat_diag_tmp) &
!$OMP SHARED(na,na_max,nb_max,Ut,tmp0,Hmat_diag)
allocate( Hmat_diag_tmp(na_max,nb_max) )
Hmat_diag_tmp = 0.d0
!$OMP DO
do n = 1, nb_max
do m = 1, na_max
do k = 1, na
Hmat_diag_tmp(m,n) += tmp0(k,n,m) * Ut(k,m)
enddo
enddo
enddo
!$OMP END DO
!$OMP CRITICAL
do n = 1, nb_max
do m = 1, na_max
Hmat_diag(m,n) += Hmat_diag_tmp(m,n)
enddo
enddo
!$OMP END CRITICAL
deallocate( Hmat_diag_tmp )
!$OMP END PARALLEL
deallocate( tmp0 , Ut )
Hdiag(:) = 0.d0
do m = 1, n_toselect
ii = numalpha_toselect(m)
jj = numbeta_toselect (m)
Hdiag(m) = Hmat_diag(ii,jj)
enddo
deallocate( Hmat_diag )
call wall_time(t4)
print *, " end const_Hdiag after (min) ", (t4-t3)/60.
print *, ""
print *, " total time (min) ", (t4-t1)/60.
print *, ""
return
end subroutine const_Hdiag