scemama/qp_plugins_scemama
1
0
Fork 0
mirror of https://gitlab.com/scemama/qp_plugins_scemama.git synced 2024-11-07 14:43:41 +01:00
Code Issues Releases Wiki Activity
2ea84d4ffc
qp_plugins_scemama/devel/svdwf/buildpsi_diagSVDit_Anthony_v4.irp.f

777 lines
23 KiB
FortranFixed
Raw Normal View History

svd save
2021-07-28 17:19:18 +02:00
program buildpsi_diagSVDit_Anthony_v4
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
double precision :: h12
integer :: i, j, k, l
double precision :: norm_psi
double precision, allocatable :: Uref(:,:), Dref(:), Vtref(:,:), Aref(:,:), Vref(:,:)
double precision :: E0, tol_energy, ctmp, Ept2, E_prev
double precision, allocatable :: H_diag(:,:), Hkl(:,:), H0(:,:), H(:,:,:,:)
double precision, allocatable :: psi_postsvd(:,:)
integer :: n_TSVD, it_svd, it_svd_max
integer :: ii, jj, n_perturb
double precision :: norm, lambda, E_perturb
double precision, allocatable :: eigval0(:), eigvec0(:,:,:), H_tmp(:,:,:,:)
double precision, allocatable :: norm_row(:), norm_col(:), Utmp(:,:), Vtmp(:,:)
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 i_H_j(det1, det2, N_int, h12)
i_state = 1
! ---------------------------------------------------------------------------------------
! construct the initial CISD matrix
print *, ' ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~'
print *, ' CI matrix:', n_det_alpha_unique,'x',n_det_beta_unique
print *, ' N det :', N_det
print *, ' ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~'
norm_psi = 0.d0
do k = 1, N_det
norm_psi = norm_psi + psi_bilinear_matrix_values(k,i_state) &
* psi_bilinear_matrix_values(k,i_state)
enddo
print *, ' initial norm = ', norm_psi
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(max(n_det_beta_unique,n_det_alpha_unique)) )
allocate( Dref(min(n_det_beta_unique,n_det_alpha_unique)) )
allocate( Vref(n_det_beta_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 *, ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
do l = 1, n_det_beta_unique
do i = 1, n_det_beta_unique
Vref(i,l) = Vtref(l,i)
enddo
enddo
deallocate( Vtref )
! Truncated rank
!n_TSVD = 100
!call write_int(6,n_TSVD, 'Rank of psi')
!________________________________________________________________________________________________________
!
! loop over SVD iterations
!________________________________________________________________________________________________________
tol_energy = 1.d0
it_svd = 0
it_svd_max = 100
E_prev = 0.d0
!print *, ci_energy(1)
allocate(H_diag(n_det_alpha_unique,n_det_beta_unique))
allocate(psi_postsvd(n_det_alpha_unique,n_det_beta_unique))
do while( ( it_svd .lt. it_svd_max) .and. ( tol_energy .gt. 1d-6 ) )
it_svd = it_svd + 1
! Truncated rank
!n_TSVD = min(n_det_alpha_unique,n_det_beta_unique)
!do i = min(n_det_alpha_unique,n_det_beta_unique), 10, -1
! if( dabs(Dref(i)) .lt. 1d-2 ) then
! n_TSVD = n_TSVD - 1
! else
! exit
! endif
!enddo
!do i = 1, min(n_det_alpha_unique,n_det_beta_unique)
! print *, i, Dref(i)
!enddo
!call write_int(6,n_TSVD, 'Rank of psi')
n_TSVD = 30 !min(n_TSVD,100)
call write_int(6,n_TSVD, 'Rank of psi')
print *, '-- Compute H --'
allocate(H(n_TSVD,n_TSVD,n_det_alpha_unique,n_det_beta_unique))
call const_H_uv_modif(Uref, Vref, H, H_diag, n_TSVD)
!call const_H_uv(Uref, Vref, H, H_diag, n_TSVD)
! E0 = < psi_0 | H | psi_0 >
!norm = 0.d0
!do j = 1, n_TSVD
! norm = norm + Dref(j)*Dref(j)
!enddo
!Dref = Dref / dsqrt(norm)
!E0 = 0.d0
!do j = 1, n_TSVD
! do i = 1, n_TSVD
! E0 = E0 + Dref(i) * H(i,i,j,j) * Dref(j)
! enddo
!enddo
!E0 = E0 + nuclear_repulsion
!print *,' E0 bef diag =', E0
allocate( H_tmp(n_TSVD,n_TSVD,n_TSVD,n_TSVD) )
do l = 1, n_TSVD
do k = 1, n_TSVD
do j = 1, n_TSVD
do i = 1, n_TSVD
H_tmp(i,j,k,l) = H(i,j,k,l)
enddo
enddo
enddo
enddo
allocate( eigval0(n_TSVD**2) , eigvec0(n_TSVD,n_TSVD,n_TSVD**2) )
eigvec0 = 0.d0
!print *, ' --- Diag post-SVD --- '
call lapack_diag(eigval0, eigvec0, H_tmp, n_TSVD**2, n_TSVD**2)
!print *, 'eig =', eigval0(1) + nuclear_repulsion
deallocate(H_tmp, eigval0)
print *, ' --- first SVD --- '
Dref = 0.d0
call perform_newpostSVD(n_TSVD, eigvec0(1,1,1), size(eigvec0,1), Uref, Vref, Dref)
deallocate(eigvec0)
print *, ' --- Compute H --- '
call const_H_uv_modif(Uref, Vref, H, H_diag, n_TSVD)
!call const_H_uv(Uref, Vref, H, H_diag, n_TSVD)
! E0 = < psi_0 | H | psi_0 >
E0 = 0.d0
norm = 0.d0
do j = 1, n_det_beta_unique
do i = 1, n_TSVD
E0 = E0 + Dref(i) * H(i,i,j,j) * Dref(j)
enddo
norm = norm + Dref(j)*Dref(j)
enddo
E0 = E0 + nuclear_repulsion
print *,' E0 aft diag =', E0
! -----------------------------------------------------------------
!print *, ' --- Perturbation --- '
psi_postsvd = 0.d0
do i = 1, n_TSVD
psi_postsvd(i,i) = Dref(i)
enddo
lambda = 1.d0
Ept2 = 0.d0
do j = 1, n_TSVD
do i = n_TSVD+1, n_det_alpha_unique
ctmp = 0.d0
do k = 1, n_TSVD
ctmp = ctmp + H(k,k,i,j) * Dref(k)
enddo
psi_postsvd(i,j) = lambda * ctmp / (E0 - (H_diag(i,j)+nuclear_repulsion) )
Ept2 += ctmp * ctmp / (E0 - (H_diag(i,j)+nuclear_repulsion) )
enddo
enddo
do j = n_TSVD+1, n_det_beta_unique
do i = 1, n_TSVD
ctmp = 0.d0
do k = 1, n_TSVD
ctmp = ctmp + H(k,k,i,j) * Dref(k)
enddo
psi_postsvd(i,j) = lambda * ctmp / (E0 - (H_diag(i,j)+nuclear_repulsion) )
Ept2 += ctmp * ctmp / (E0 - (H_diag(i,j)+nuclear_repulsion) )
enddo
enddo
deallocate( H )
! -----------------------------------------------------------------
n_perturb = n_TSVD + n_TSVD
allocate( norm_row(n_det_alpha_unique) , norm_col(n_det_beta_unique) )
do i = 1, n_det_alpha_unique
norm_row(i) = 0.d0
do j = 1, n_det_beta_unique
norm_row(i) += psi_postsvd(i,j) * psi_postsvd(i,j)
enddo
enddo
do j = 1, n_det_beta_unique
norm_col(j) = 0.d0
do i = 1, n_det_alpha_unique
norm_col(j) += psi_postsvd(i,j) * psi_postsvd(i,j)
enddo
enddo
allocate( Utmp(n_det_alpha_unique,n_perturb) , Vtmp(n_det_beta_unique,n_perturb) )
do i = 1, n_perturb
ii = MAXLOC( norm_row , DIM=1 )
jj = MAXLOC( norm_col , DIM=1 )
if( (norm_row(ii).lt.1.d-12) .or. (norm_col(jj).lt.1.d-12) ) then
print *, ' !!!!!!! '
print *, ii, norm_row(ii)
print *, jj, norm_col(jj)
stop
endif
Utmp(:,i) = Uref(:,ii)
Vtmp(:,i) = Vref(:,jj)
norm_row(ii) = 0.d0
norm_col(jj) = 0.d0
enddo
deallocate( norm_row , norm_col )
print *, ' --- Compute H in n_perturb space --- '
allocate(H(n_perturb,n_perturb,n_det_alpha_unique,n_det_beta_unique))
call const_H_TSVD(Uref, Vref, H, n_perturb)
allocate( H_tmp(n_perturb,n_perturb,n_perturb,n_perturb) )
do l = 1, n_perturb
do k = 1, n_perturb
do j = 1, n_perturb
do i = 1, n_perturb
H_tmp(i,j,k,l) = H(i,j,k,l)
enddo
enddo
enddo
enddo
deallocate( H )
allocate( eigval0(n_perturb**2) , eigvec0(n_perturb,n_perturb,n_perturb**2) )
eigvec0 = 0.d0
call lapack_diag(eigval0, eigvec0, H_tmp, n_perturb**2, n_perturb**2)
E_perturb = eigval0(1) + nuclear_repulsion
print *, ' diag in n_perturb space: ', E_perturb
deallocate(H_tmp, eigval0)
print *, ' --- second SVD --- '
Uref = 0.d0
Vref = 0.d0
Dref = 0.d0
do l = 1, n_perturb
Uref(1:n_det_alpha_unique,l) = Utmp(1:n_det_alpha_unique,l)
Vref(1:n_det_beta_unique ,l) = Vtmp(1:n_det_beta_unique ,l)
enddo
deallocate( Utmp , Vtmp )
call perform_perturbSVD(n_TSVD, n_perturb, eigvec0(1,1,1), size(eigvec0,1), Uref, Vref, Dref)
deallocate(eigvec0)
tol_energy = dabs(E_prev - E0)
print '(I5, 2(2X,I5), 3(3X, F20.10))', it_svd, n_TSVD, n_perturb, E0, E_perturb, tol_energy
write(222,'(I5, 2(2X,I5), 3(3X, F20.10))') it_svd, n_TSVD, n_perturb, E0, E_perturb, tol_energy
E_prev = E0
!print *, ' --- SVD --- '
!call perform_newpostSVD(n_TSVD, psi_postsvd, size(psi_postsvd,1), Uref, Vref, Dref)
!call perform_newSVD(n_TSVD, psi_postsvd, size(psi_postsvd,1), Uref, Vref, Dref)
end do
end
subroutine perform_newpostSVD(n_TSVD, psi_postsvd, LDP, Uref, Vref, Dref)
integer, intent(in) :: n_TSVD, LDP
double precision, intent(in) :: psi_postsvd(LDP,n_TSVD)
double precision, intent(inout) :: Uref(n_det_alpha_unique,n_det_alpha_unique)
double precision, intent(inout) :: Vref(n_det_beta_unique ,n_det_beta_unique)
!double precision, intent(inout) :: Dref(max(n_det_beta_unique,n_det_alpha_unique))
double precision, intent(inout) :: Dref(min(n_det_beta_unique,n_det_alpha_unique))
integer :: mm, nn, i, j, l, na, nb
double precision :: err0, err_norm, err_tmp, norm_tmp
double precision, allocatable :: S_mat(:,:), SxVt(:,:)
double precision, allocatable :: U_svd(:,:), V_svd(:,:)
double precision, allocatable :: U_newsvd(:,:), V_newsvd(:,:), Vt_newsvd(:,:), D_newsvd(:), A_newsvd(:,:)
mm = n_det_alpha_unique
nn = n_det_beta_unique
allocate( U_svd(mm,n_TSVD) , V_svd(nn,n_TSVD) , S_mat(n_TSVD,n_TSVD) )
U_svd(1:mm,1:n_TSVD) = Uref(1:mm,1:n_TSVD)
V_svd(1:nn,1:n_TSVD) = Vref(1:nn,1:n_TSVD)
S_mat(1:n_TSVD,1:n_TSVD) = psi_postsvd(1:n_TSVD,1:n_TSVD)
! first compute S_mat x transpose(V_svd)
allocate( SxVt(n_TSVD,nn) )
call dgemm( 'N', 'T', n_TSVD, nn, n_TSVD, 1.d0 &
, S_mat , size(S_mat,1) &
, V_svd , size(V_svd,1) &
, 0.d0, SxVt, size(SxVt ,1) )
deallocate(S_mat)
! then compute U_svd x SxVt
allocate( A_newsvd(mm,nn) )
call dgemm( 'N', 'N', mm, nn, n_TSVD, 1.d0 &
, U_svd , size(U_svd ,1) &
, SxVt , size(SxVt ,1) &
, 0.d0, A_newsvd, size(A_newsvd,1) )
deallocate( SxVt )
! perform new SVD
!allocate( U_newsvd(mm,mm), Vt_newsvd(nn,nn), D_newsvd(max(mm,nn)) )
allocate( U_newsvd(mm,mm), Vt_newsvd(nn,nn), D_newsvd(min(mm,nn)) )
call svd_s( A_newsvd, size(A_newsvd,1), &
U_newsvd, size(U_newsvd,1), &
D_newsvd, &
Vt_newsvd, size(Vt_newsvd,1), &
mm, nn)
deallocate(A_newsvd)
allocate( V_newsvd(nn,nn) )
do l = 1, nn
do j = 1, nn
V_newsvd(j,l) = Vt_newsvd(l,j)
enddo
enddo
deallocate(Vt_newsvd)
!do l = 1, n_TSVD
! Dref(l) = D_newsvd(l)
! Uref(1:mm,l) = U_newsvd(1:mm,l)
! Vref(1:nn,l) = V_newsvd(1:nn,l)
!enddo
Dref(1:n_TSVD) = D_newsvd(1:n_TSVD)
Uref(1:mm,1:mm) = U_newsvd(1:mm,1:mm)
Vref(1:nn,1:nn) = V_newsvd(1:nn,1:nn)
deallocate(U_newsvd)
deallocate(V_newsvd)
deallocate(D_newsvd)
end subroutine perform_newpostSVD
subroutine perform_perturbSVD(n_TSVD, n_perturb, psi_postsvd, LDP, Uref, Vref, Dref)
integer, intent(in) :: n_TSVD, LDP, n_perturb
double precision, intent(in) :: psi_postsvd(LDP,n_TSVD)
double precision, intent(inout) :: Uref(n_det_alpha_unique,n_det_alpha_unique)
double precision, intent(inout) :: Vref(n_det_beta_unique ,n_det_beta_unique)
double precision, intent(inout) :: Dref(min(n_det_beta_unique,n_det_alpha_unique))
integer :: mm, nn, i, j, l, na, nb
double precision :: err0, err_norm, err_tmp, norm_tmp
double precision, allocatable :: S_mat(:,:), SxVt(:,:)
double precision, allocatable :: U_svd(:,:), V_svd(:,:)
double precision, allocatable :: U_newsvd(:,:), V_newsvd(:,:), Vt_newsvd(:,:), D_newsvd(:), A_newsvd(:,:)
mm = n_det_alpha_unique
nn = n_det_beta_unique
allocate( U_svd(mm,n_perturb) , V_svd(nn,n_perturb) , S_mat(n_perturb,n_perturb) )
U_svd(1:mm,1:n_perturb) = Uref(1:mm,1:n_perturb)
V_svd(1:nn,1:n_perturb) = Vref(1:nn,1:n_perturb)
S_mat(1:n_perturb,1:n_perturb) = psi_postsvd(1:n_perturb,1:n_perturb)
! first compute S_mat x transpose(V_svd)
allocate( SxVt(n_perturb,nn) )
call dgemm( 'N', 'T', n_perturb, nn, n_perturb, 1.d0 &
, S_mat , size(S_mat,1) &
, V_svd , size(V_svd,1) &
, 0.d0, SxVt, size(SxVt ,1) )
deallocate(S_mat)
! then compute U_svd x SxVt
allocate( A_newsvd(mm,nn) )
call dgemm( 'N', 'N', mm, nn, n_perturb, 1.d0 &
, U_svd , size(U_svd ,1) &
, SxVt , size(SxVt ,1) &
, 0.d0, A_newsvd, size(A_newsvd,1) )
deallocate( SxVt )
! perform new SVD
allocate( U_newsvd(mm,mm), Vt_newsvd(nn,nn), D_newsvd(min(mm,nn)) )
call svd_s( A_newsvd, size(A_newsvd,1), &
U_newsvd, size(U_newsvd,1), &
D_newsvd, &
Vt_newsvd, size(Vt_newsvd,1), &
mm, nn)
deallocate(A_newsvd)
allocate( V_newsvd(nn,nn) )
do l = 1, nn
do j = 1, nn
V_newsvd(j,l) = Vt_newsvd(l,j)
enddo
enddo
deallocate(Vt_newsvd)
Dref(1:n_TSVD) = D_newsvd(1:n_TSVD)
Uref(1:mm,1:mm) = U_newsvd(1:mm,1:mm)
Vref(1:nn,1:nn) = V_newsvd(1:nn,1:nn)
deallocate(U_newsvd)
deallocate(V_newsvd)
deallocate(D_newsvd)
end subroutine perform_perturbSVD
subroutine perform_newSVD(n_TSVD, psi_postsvd, LDP, Uref, Vref, Dref)
integer, intent(in) :: n_TSVD, LDP
double precision, intent(in) :: psi_postsvd(LDP,n_TSVD)
double precision, intent(inout) :: Uref(n_det_alpha_unique,n_det_alpha_unique)
double precision, intent(inout) :: Vref(n_det_beta_unique ,n_det_beta_unique)
double precision, intent(inout) :: Dref(min(n_det_beta_unique,n_det_alpha_unique))
integer :: mm, nn, i, j, l, na, nb
double precision :: err0, err_norm, err_tmp, norm_tmp
double precision, allocatable :: S_mat(:,:), SxVt(:,:)
double precision, allocatable :: U_svd(:,:), V_svd(:,:)
double precision, allocatable :: U_newsvd(:,:), V_newsvd(:,:), Vt_newsvd(:,:), D_newsvd(:), A_newsvd(:,:)
mm = n_det_alpha_unique
nn = n_det_beta_unique
allocate( U_svd(mm,mm) , V_svd(nn,nn) , S_mat(mm,nn) )
U_svd(1:mm,1:mm) = Uref(1:mm,1:mm)
V_svd(1:nn,1:nn) = Vref(1:nn,1:nn)
norm_tmp = 0.d0
do i = 1, nn
do j = 1, mm
S_mat(j,i) = psi_postsvd(j,i)
norm_tmp += psi_postsvd(j,i) * psi_postsvd(j,i)
enddo
enddo
norm_tmp = 1.d0 / dsqrt(norm_tmp)
do i = 1, nn
do j = 1, mm
S_mat(j,i) = S_mat(j,i) * norm_tmp
enddo
enddo
! first compute S_mat x transpose(V_svd)
allocate( SxVt(mm,nn) )
call dgemm( 'N', 'T', mm, nn, nn, 1.d0 &
, S_mat , size(S_mat,1) &
, V_svd , size(V_svd,1) &
, 0.d0, SxVt, size(SxVt ,1) )
deallocate(S_mat)
! then compute U_svd x SxVt
allocate( A_newsvd(mm,nn) )
call dgemm( 'N', 'N', mm, nn, mm, 1.d0 &
, U_svd , size(U_svd ,1) &
, SxVt , size(SxVt ,1) &
, 0.d0, A_newsvd, size(A_newsvd,1) )
deallocate( SxVt )
! perform new SVD
allocate( U_newsvd(mm,mm), Vt_newsvd(nn,nn), D_newsvd(min(mm,nn)) )
call svd_s( A_newsvd, size(A_newsvd,1), U_newsvd, size(U_newsvd,1), D_newsvd, Vt_newsvd, size(Vt_newsvd,1), mm, nn)
deallocate(A_newsvd)
allocate( V_newsvd(nn,nn) )
do l = 1, nn
do j = 1, nn
V_newsvd(j,l) = Vt_newsvd(l,j)
enddo
enddo
deallocate(Vt_newsvd)
!Dref(1:n_TSVD) = D_newsvd(1:n_TSVD)
Dref = D_newsvd
Uref(1:mm,1:mm) = U_newsvd(1:mm,1:mm)
Vref(1:nn,1:nn) = V_newsvd(1:nn,1:nn)
deallocate(U_newsvd)
deallocate(V_newsvd)
deallocate(D_newsvd)
return
end subroutine perform_newSVD
subroutine const_H_uv_modif(Uref, Vref, H, H_diag, n_TSVD)
USE OMP_LIB
implicit none
integer, intent(in) :: n_TSVD
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) :: H(n_TSVD,n_TSVD, n_det_alpha_unique, n_det_beta_unique)
double precision, intent(out) :: H_diag(n_det_alpha_unique,n_det_beta_unique)
integer(bit_kind) :: det1(N_int,2), det2(N_int,2)
integer :: i, j, k, l, degree, n, m, na, nb
double precision :: h12
double precision, allocatable :: H0(:,:,:,:)
double precision, allocatable :: H1(:,:,:,:)
double precision, allocatable :: tmp3(:,:,:)
double precision, allocatable :: tmp1(:,:), tmp0(:,:), tmp4(:,:)
double precision :: c_tmp
na = n_det_alpha_unique
nb = n_det_beta_unique
call wall_time(t0)
allocate( H0(na,nb,n_TSVD,n_TSVD) )
allocate( tmp3(na,nb,nb) )
H0 = 0.d0
tmp3 = 0.d0
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(i,j,k,l,m,n,det1,det2,degree,h12,c_tmp,tmp1,tmp0,tmp4) &
!$OMP SHARED(na,nb,psi_det_alpha_unique,psi_det_beta_unique, &
!$OMP N_int,tmp3,Uref,Vref,H_diag,H0,n_TSVD)
allocate(tmp1(na,na), tmp0(nb,nb), tmp4(nb,nb))
do i = 1, na
do m = 1, na
tmp1(m,i) = Uref(i,m)
enddo
enddo
!$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 > 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 > 2) cycle
call i_H_j(det1, det2, N_int, h12)
if (h12 == 0.d0) cycle
do m = 1, na
tmp3(m,j,l) = tmp3(m,j,l) + h12 * tmp1(m,i) * tmp1(m,k)
enddo
do n = 1, n_TSVD
c_tmp = h12 * Vref(j,n)
if (c_tmp == 0.d0) cycle
do m = 1, n_TSVD
H0(k,l,m,n) = H0(k,l,m,n) + c_tmp * tmp1(m,i)
enddo
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP DO
do m = 1, na
do l = 1, nb
do j = 1, nb
tmp4(j,l) = tmp3(m,j,l)
enddo
enddo
call DGEMM('N','N',nb,nb,nb,1.d0, &
tmp4, size(tmp4,1), &
Vref, size(Vref,1), &
0.d0, tmp0, size(tmp0,1))
do n = 1, nb
H_diag(m,n) = 0.d0
do j = 1, nb
H_diag(m,n) = H_diag(m,n) + tmp0(j,n) * Vref(j,n)
enddo
enddo
enddo
!$OMP END DO
deallocate(tmp1, tmp0)
deallocate(tmp4)
!$OMP END PARALLEL
deallocate(tmp3)
call wall_time(t1)
! H0(na,nb,n_TSVD,n_TSVD)
allocate( H1(nb,n_TSVD,n_TSVD,na) )
call DGEMM('T', 'N', nb * n_TSVD * n_TSVD, na, na, 1.d0, &
H0 , size(H0,1) , &
Uref, size(Uref,1), 0.d0, &
H1 , size(H1,1) * size(H1,2) * size(H1,3) )
deallocate( H0 )
! (l,p,q,r) -> (p,q,r,s)
call DGEMM('T','N', n_TSVD * n_TSVD * na, nb, nb, &
1.d0, H1, size(H1,1), Vref, size(Vref,1), 0.d0, H, size(H,1)*size(H,2)*size(H,3))
! do j=1,n_TSVD
! do i=1,n_TSVD
! print *, H_diag(i,j), H(i,j,i,j)
! enddo
! enddo
deallocate(H1)
call wall_time(t2)
print *, 't=', t1-t0, t2-t1
double precision :: t0, t1, t2
! stop
end subroutine const_H_uv_modif
subroutine const_H_TSVD(Uref, Vref, H, n_TSVD)
USE OMP_LIB
implicit none
integer, intent(in) :: n_TSVD
double precision, intent(in) :: Uref(n_det_alpha_unique,n_TSVD)
double precision, intent(in) :: Vref(n_det_beta_unique ,n_TSVD)
double precision, intent(out) :: H(n_TSVD,n_TSVD, n_det_alpha_unique, n_det_beta_unique)
integer(bit_kind) :: det1(N_int,2), det2(N_int,2)
integer :: i, j, k, l, degree, n, m, na, nb
double precision :: h12
double precision, allocatable :: H0(:,:,:,:)
double precision, allocatable :: H1(:,:,:,:)
double precision, allocatable :: tmp1(:,:)
double precision :: c_tmp
double precision :: t0, t1, t2
na = n_det_alpha_unique
nb = n_det_beta_unique
call wall_time(t0)
allocate( H0(na,nb,n_TSVD,n_TSVD) )
H0 = 0.d0
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(i,j,k,l,m,n,det1,det2,degree,h12,c_tmp,tmp1) &
!$OMP SHARED(na,nb,psi_det_alpha_unique,psi_det_beta_unique &
!$OMP ,N_int,Uref,Vref,H0,n_TSVD)
allocate(tmp1(na,na))
do i = 1, na
do m = 1, na
tmp1(m,i) = Uref(i,m)
enddo
enddo
!$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 > 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 > 2) cycle
call i_H_j(det1, det2, N_int, h12)
if (h12 == 0.d0) cycle
do n = 1, n_TSVD
c_tmp = h12 * Vref(j,n)
if (c_tmp == 0.d0) cycle
do m = 1, n_TSVD
H0(k,l,m,n) = H0(k,l,m,n) + c_tmp * tmp1(m,i)
enddo
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
deallocate(tmp1)
!$OMP END PARALLEL
call wall_time(t1)
allocate( H1(nb,n_TSVD,n_TSVD,na) )
call DGEMM('T','N', nb * n_TSVD * n_TSVD, na, na, &
1.d0, H0, size(H0,1), Uref, size(Uref,1), 0.d0, H1, size(H1,1)*size(H1,2)*size(H1,3))
deallocate( H0 )
! (l,p,q,r) -> (p,q,r,s)
call DGEMM('T','N', n_TSVD * n_TSVD * na, nb, nb, &
1.d0, H1, size(H1,1), Vref, size(Vref,1), 0.d0, H, size(H,1)*size(H,2)*size(H,3))
deallocate(H1)
call wall_time(t2)
print *, 't=', t1-t0, t2-t1
end subroutine const_H_TSVD
Reference in New Issue Copy Permalink