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/psiSVD_naiv1by1_v1.irp.f
Abdallah Ammar 2fd2fcac5d svd save
2021-07-28 17:19:18 +02:00

528 lines
16 KiB
Fortran

program psiSVD_naiv1by1_v1
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
double precision :: x, y, h12
double precision, allocatable :: Uref(:,:), Dref(:), Vtref(:,:), Aref(:,:), Vref(:,:)
integer :: rank_max
double precision :: E0, overlop, Ept2, Em, norm
double precision, allocatable :: H0(:,:)
double precision, allocatable :: eigvec0(:,:), eigval0(:), coeff_psi(:), coeff_tmp(:)
integer :: ii, jj, ia, ib, ja, jb
double precision, allocatable :: Hdiag(:), H0_1d(:), H0_tmp(:,:)
integer :: na_new, nb_new, ind_new, ind_gs
double precision :: ctmp, coeff_new
double precision, allocatable :: epsil(:), epsil_energ(:), check_ov(:)
double precision, allocatable :: Uezfio(:,:,:), Dezfio(:,:), Vezfio(:,:,:)
integer :: n_selected, n_toselect, n_tmp, na_max, nb_max
integer, allocatable :: numalpha_selected(:), numbeta_selected(:)
integer, allocatable :: numalpha_tmp(:), numbeta_tmp(:)
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.
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 )
deallocate( Aref )
! ---------------------------------------------------------------------------------------
! *** PARAMETERS *** !
na_max = n_det_alpha_unique
nb_max = n_det_beta_unique
! *** ***** *** !
print *, ' na_max = ', na_max
print *, ' nb_max = ', nb_max
! ---------------------------------------------------------------------------------------
! initial wavefunction: psi_0
n_selected = 1
allocate( numalpha_selected(n_selected) , numbeta_selected(n_selected) )
numalpha_selected(1) = 1
numbeta_selected (1) = 1
! ---------------------------------------------------------------------------------------
! ---------------------------------------------------------------------------------------
! construnc the initial basis to select phi_1 from the FSVD
n_toselect = min(na_max,nb_max) - n_selected
print *, ' toselect = ', n_toselect
print *, ' to trun = ', n_det_alpha_unique*n_det_beta_unique - na_max*nb_max
! ---------------------------------------------------------------------------------------
! ---------------------------------------------------------------------------------------
! read < u_k v_l | H | u_k v_l > for all vectors
allocate( Hdiag(n_toselect) , H0(n_selected,n_selected) )
n_tmp = n_det_alpha_unique * n_det_beta_unique - 1
open( unit=11, FILE="klHkl_v1.dat", ACTION="READ")
!open( unit=11, FILE="klHkl_v2.dat", ACTION="READ")
read(11,*) i, i, E0
H0(1,1) = E0
do i = 1, n_tmp
read(11,*) ia, ib, ctmp
if( ia .eq. ib ) then
ii = ia - 1
Hdiag(ii) = ctmp
!print *, ia, ib , Hdiag(ia-1)
endif
enddo
close(11)
! ---------------------------------------------------------------------------------------
E0 = E0 + nuclear_repulsion
Em = E0
print*, ' space dimen = ', n_selected
print*, ' ground state Em = ', Em
print*, ' ground state E0 = ', E0
na_new = 1
nb_new = 1
ind_new = 0
!________________________________________________________________________________________________________
!
! increase the size of psi0 iteratively
!________________________________________________________________________________________________________
! *** PARAMETERS *** !
rank_max = min( na_max , nb_max ) - 1
! *** ***** *** !
if( rank_max .gt. n_toselect ) then
print *, " rank_max should be less then n_toselect"
stop
endif
do while( n_selected .lt. rank_max )
call SYSTEM_CLOCK(COUNT=W_tbeg_it, COUNT_RATE=W_ir)
print*, ' '
print*, ' new iteration '
if( n_toselect .lt. 1 ) then
print*, ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
print*, ' no more vectors to construct a new basis '
print*, ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
stop
else
! ---------------------------------------------------------------------------------------
! select a new vector
na_new += 1
nb_new += 1
ind_new += 1
print *, ' best vector', na_new, nb_new
! < psi_old | H | delta_psi >
allocate( H0_1d(n_selected) )
call const_H0_1d(na_new, nb_new, na_max, nb_max, n_selected, Uref, Vref, numalpha_selected, numbeta_selected, H0_1d)
! ---------------------------------------------------------------------------------------
! ---------------------------------------------------------------------------------------
! new psi
allocate( numalpha_tmp(n_selected), numbeta_tmp(n_selected) )
allocate( H0_tmp(n_selected,n_selected) )
numalpha_tmp(:) = numalpha_selected(:)
numbeta_tmp (:) = numbeta_selected (:)
H0_tmp (:,:) = H0 (:,:)
deallocate( numalpha_selected, numbeta_selected, H0 )
n_tmp = n_selected
n_selected = n_selected + 1
allocate( numalpha_selected(n_selected) , numbeta_selected(n_selected) )
allocate( H0(n_selected,n_selected) )
H0(:,:) = 0.d0
do l = 1, n_tmp
numalpha_selected(l) = numalpha_tmp(l)
numbeta_selected (l) = numbeta_tmp (l)
enddo
H0(1:n_tmp,1:n_tmp) = H0_tmp(1:n_tmp,1:n_tmp)
deallocate( numalpha_tmp, numbeta_tmp, H0_tmp )
numalpha_selected(n_selected) = na_new
numbeta_selected (n_selected) = nb_new
H0(1:n_tmp,n_selected) = H0_1d(1:n_tmp)
H0(n_selected,1:n_tmp) = H0_1d(1:n_tmp)
deallocate( H0_1d )
H0(n_selected,n_selected) = Hdiag(ind_new)
! ---------------------------------------------------------------------------------------
! ---------------------------------------------------------------------------------------
! energy without diag
! < psi | psi >
norm = 0.d0
do j = 1, n_selected
ja = numalpha_selected(j)
jb = numbeta_selected (j)
if(ja.eq.jb) norm = norm + Dref(ja)*Dref(jb)
enddo
! < psi | H | psi >
Em = 0.d0
do j = 1, n_selected
ja = numalpha_selected(j)
jb = numbeta_selected (j)
if(ja.eq.jb) then
do i = 1, n_selected
ia = numalpha_selected(i)
ib = numbeta_selected (i)
if(ia.eq.ib) Em = Em + Dref(ja) * H0(j,i) * Dref(ia)
enddo
endif
enddo
! Em = < psi | H | psi > / < psi | psi >
Em = Em / norm + nuclear_repulsion
! ---------------------------------------------------------------------------------------
! ---------------------------------------------------------------------------------------
! energy with diag
allocate( eigvec0(n_selected,n_selected), eigval0(n_selected) )
call lapack_diag(eigval0, eigvec0, H0, n_selected, n_selected)
! get the postsvd ground state
allocate( check_ov(n_selected) )
do l = 1, n_selected
overlop = 0.d0
do i = 1, n_selected
ia = numalpha_selected(i)
ib = numbeta_selected (i)
if( ia .eq. ib ) overlop = overlop + eigvec0(i,l) * Dref(ia)
enddo
check_ov(l) = dabs(overlop)
enddo
ind_gs = MAXLOC( check_ov, DIM=1 )
overlop = check_ov(ind_gs)
E0 = eigval0(ind_gs)+nuclear_repulsion
coeff_psi(:) = eigvec0(:,ind_gs)
deallocate( check_ov, eigval0, eigvec0 )
! ---------------------------------------------------------------------------------------
print*, ' space dimen = ', n_selected
print*, ' E bef diag = ', Em
print*, ' E aft diag = ', E0
print*, ' overlop = ', overlop
print*, ' index = ', ind_gs
write(211, '( 3(I5,3X), 4(F15.8,3X) )') n_selected, na_new, nb_new, Em, E0, overlop
! ---------------------------------------------------------------------------------------
! remove selected pair | na_new nb_new >
n_toselect = n_toselect - 1
print*, ' rank to select = ', n_toselect
! ---------------------------------------------------------------------------------------
endif
call SYSTEM_CLOCK(COUNT=W_tend_it, COUNT_RATE=W_ir)
W_tot_time_it = real(W_tend_it-W_tbeg_it, kind=8) / real(W_ir, kind=8)
print*, " "
print*, " elapsed time (min) = ", W_tot_time_it/60.d0
end do
!________________________________________________________________________________________________________
!________________________________________________________________________________________________________
! ***************************************************************************************************
! save to ezfion
!allocate( Uezfio(n_det_alpha_unique,rank0,1), Dezfio(rank0,1), Vezfio(n_det_beta_unique,rank0,1) )
!do l = 1, rank0
! Dezfio(l,1) = coeff_psi(l)
! Uezfio(:,l,1) = U0(:,l)
! Vezfio(:,l,1) = V0(:,l)
!enddo
!call ezfio_set_spindeterminants_n_det(N_det)
!call ezfio_set_spindeterminants_n_states(N_states)
!call ezfio_set_spindeterminants_n_det_alpha(n_det_alpha_unique)
!call ezfio_set_spindeterminants_n_det_beta(n_det_beta_unique)
!call ezfio_set_spindeterminants_psi_coef_matrix_rows(psi_bilinear_matrix_rows)
!call ezfio_set_spindeterminants_psi_coef_matrix_columns(psi_bilinear_matrix_columns)
!call ezfio_set_spindeterminants_psi_coef_matrix_values(psi_bilinear_matrix_values)
!call ezfio_set_spindeterminants_n_svd_coefs(rank0)
!call ezfio_set_spindeterminants_psi_svd_alpha(Uezfio)
!call ezfio_set_spindeterminants_psi_svd_beta(Vezfio )
!call ezfio_set_spindeterminants_psi_svd_coefs(Dezfio)
!deallocate( Uezfio, Dezfio, Vezfio )
! ***************************************************************************************************
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 *, ' ___________________________________________________________________'
deallocate( Dref )
deallocate( Uref, Vref )
deallocate( coeff_psi )
deallocate( numalpha_selected, numbeta_selected )
deallocate( H0, Hdiag )
end
subroutine const_H0_1d(na_new, nb_new, na_max, nb_max, n_selected, Uref, Vref, numalpha_selected, numbeta_selected, H0_1d)
implicit none
integer, intent(in) :: na_new, nb_new, na_max, nb_max, n_selected
integer, intent(in) :: numalpha_selected(n_selected), numbeta_selected(n_selected)
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) :: H0_1d(n_selected)
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
double precision :: h12
double precision, allocatable :: Hmat_kl(:,:), tmp1(:,:), tmp2(:,:)
double precision, allocatable :: U1d(:), V1d(:)
double precision, allocatable :: Utmp(:,:), Vtmp(:,:)
double precision :: ti, tf
print *, ""
print *, " start const_H0_1d"
call wall_time(ti)
na = n_det_alpha_unique
nb = n_det_beta_unique
allocate( U1d(na) , V1d(nb) )
U1d(1:na) = Uref(1:na,na_new)
V1d(1:nb) = Vref(1:nb,nb_new)
allocate( tmp1(na,nb) )
tmp1 = 0.d0
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(i,j,k,l,h12,det1,det2,degree,tmp2) &
!$OMP SHARED(na,nb,psi_det_alpha_unique,psi_det_beta_unique, &
!$OMP N_int,U1d,V1d,tmp1)
allocate( tmp2(na,nb) )
tmp2 = 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
tmp2(i,j) += h12 * U1d(k) * V1d(l)
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP CRITICAL
do j = 1, nb
do i = 1, na
tmp1(i,j) += tmp2(i,j)
enddo
enddo
!$OMP END CRITICAL
deallocate( tmp2 )
!$OMP END PARALLEL
deallocate( U1d , V1d )
! tmp2(j,m) = sum_i tmp1(i,j) x Uref(i,m)
allocate( Utmp(na,na_max) )
Utmp(1:na,1:na_max) = Uref(1:na,1:na_max)
allocate( tmp2(nb,na_max) )
call DGEMM('T', 'N', nb, na_max, na, 1.d0, &
tmp1, size(tmp1,1), Utmp, size(Utmp,1), &
0.d0, tmp2, size(tmp2,1) )
deallocate( tmp1 )
deallocate( Utmp )
! Hmat_kl(m,n) = sum_j tmp2(j,m) x Vref(j,n)
allocate( Vtmp(nb,nb_max) )
Vtmp(1:nb,1:nb_max) = Vref(1:nb,1:nb_max)
allocate( Hmat_kl(na_max,nb_max) )
call DGEMM('T', 'N', na_max, nb_max, nb, 1.d0, &
tmp2, size(tmp2,1), Vtmp, size(Vtmp,1), &
0.d0, Hmat_kl, size(Hmat_kl,1) )
deallocate( tmp2 )
deallocate( Vtmp )
do m = 1, n_selected
ii = numalpha_selected(m)
jj = numbeta_selected (m)
H0_1d(m) = Hmat_kl(ii,jj)
enddo
deallocate( Hmat_kl )
call wall_time(tf)
print *, " end const_H0_1d after (min) ", (tf-ti)/60.
print *, ""
return
end subroutine const_H0_1d