LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-24 12:33:30 +01:00
Code Releases Activity

Working on sigma_vector for CFG-CI.

Browse Source
This commit is contained in:
v1j4y 2021-03-01 18:28:45 +01:00
parent c1ef95cd60
commit c38d25fdba
3 changed files with 508 additions and 27 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
src/csf/configuration_CI_sigma_helpers.irp.f
View File

@ -1,4 +1,4 @@
subroutine obtain_associated_alphaI(idxI, Icfg, alphasIcfg, NalphaIcfg, factor_alphaI)
subroutine obtain_associated_alphaI(idxI, Icfg, alphasIcfg, NalphaIcfg)
implicit none
use bitmasks
BEGIN_DOC
@ -10,7 +10,6 @@
integer,intent(in) :: idxI ! The id of the Ith CFG
integer(bit_kind),intent(in) :: Icfg(N_int,2)
integer,intent(out) :: NalphaIcfg
real*8 ,intent(out) :: factor_alphaI(*)
integer(bit_kind),intent(out) :: alphasIcfg(N_int,2,*)
logical,dimension(:,:),allocatable :: tableUniqueAlphas
integer :: listholes(mo_num)
@ -293,19 +292,26 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod
! Type 3 - SOMO -> VMO
! Type 4 - DOMO -> SOMO
END_DOC
integer(bit_kind),intent(in) :: Ialpha(N_int,2)
integer(bit_kind),intent(in) :: Jcfg(N_int,2)
integer,intent(in) :: p,q
integer,intent(in) :: extype
integer,intent(out) :: pmodel,qmodel
integer*8 :: Isomo
integer*8 :: Idomo
integer*8 :: Jsomo
integer*8 :: Jdomo
integer*8 :: mask
integer*8 :: Isomotmp
integer*8 :: Jsomotmp
integer :: pos0,pos0prev
integer(bit_kind),intent(in) :: Ialpha(N_int,2)
integer(bit_kind),intent(in) :: Jcfg(N_int,2)
integer,intent(in) :: p,q
integer,intent(in) :: extype
integer,intent(out) :: pmodel,qmodel
!integer(bit_kind) :: Isomo(N_int)
!integer(bit_kind) :: Idomo(N_int)
!integer(bit_kind) :: Jsomo(N_int)
!integer(bit_kind) :: Jdomo(N_int)
integer*8 :: Isomo
integer*8 :: Idomo
integer*8 :: Jsomo
integer*8 :: Jdomo
integer*8 :: mask
integer :: iint, ipos
!integer(bit_kind) :: Isomotmp(N_int)
!integer(bit_kind) :: Jsomotmp(N_int)
integer*8 :: Isomotmp
integer*8 :: Jsomotmp
integer :: pos0,pos0prev
! TODO Flag (print) when model space indices is > 64
Isomo = Ialpha(1,1)
@ -317,15 +323,9 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod
qmodel = q
if(p .EQ. q) then
!print *,"input pq=",p,q,"extype=",extype
pmodel = 1
qmodel = 1
else
!print *,"input pq=",p,q,"extype=",extype
!call debug_spindet(Isomo,1)
!call debug_spindet(Idomo,1)
!call debug_spindet(Jsomo,1)
!call debug_spindet(Jdomo,1)
select case(extype)
case (1)
! SOMO -> SOMO

464
src/csf/sigma_vector.irp.f
View File

@ -708,3 +708,467 @@ end subroutine get_phase_qp_to_cfg
end do
end do
END_PROVIDER
subroutine calculate_preconditioner_cfg(diag_energies)
implicit none
use bitmasks
BEGIN_DOC
! Documentation for calculate_preconditioner
!
! Calculates the diagonal energies of
! the configurations in psi_configuration
! returns : diag_energies :
END_DOC
integer :: i,j,k,l,p,q,noccp,noccq, ii, jj
real*8,intent(out) :: diag_energies(n_CSF)
integer :: nholes
integer :: nvmos
integer :: listvmos(mo_num)
integer :: vmotype(mo_num) ! 1 -> VMO 2 -> SOMO
integer :: listholes(mo_num)
integer :: holetype(mo_num) ! 1-> SOMO 2->DOMO
integer*8 :: Idomo
integer*8 :: Isomo
integer*8 :: Jdomo
integer*8 :: Jsomo
integer*8 :: diffSOMO
integer*8 :: diffDOMO
integer :: NSOMOI
integer :: NSOMOJ
integer :: ndiffSOMO
integer :: ndiffDOMO
integer :: starti, endi, cnti, cntj, rows,cols
integer :: extype,pmodel,qmodel
integer(bit_kind) :: Icfg(N_INT,2)
integer(bit_kind) :: Jcfg(N_INT,2)
integer,external :: getNSOMO
real*8, external :: mo_two_e_integral
real*8 :: hpp
real*8 :: meCC
real*8 :: ecore
! initialize energies
diag_energies = 0.d0
! calculate core energy
!call get_core_energy(ecore)
!diag_energies = ecore
! calculate the core energy
!print *,"Core energy=",ref_bitmask_energy
do i=1,N_configuration
Isomo = psi_configuration(1,1,i)
Idomo = psi_configuration(1,2,i)
Icfg(1,1) = psi_configuration(1,1,i)
Icfg(1,2) = psi_configuration(1,2,i)
NSOMOI = getNSOMO(psi_configuration(:,:,i))
starti = psi_config_data(i,1)
endi = psi_config_data(i,2)
! find out all pq holes possible
nholes = 0
! holes in SOMO
!do k = n_core_orb+1,n_core_orb + n_act_orb
do k = 1,mo_num
if(POPCNT(IAND(Isomo,IBSET(0_8,k-1))) .EQ. 1) then
nholes += 1
listholes(nholes) = k
holetype(nholes) = 1
endif
enddo
! holes in DOMO
!do k = n_core_orb+1,n_core_orb + n_act_orb
!do k = 1+n_core_inact_orb,n_core_orb+n_core_inact_act_orb
do k = 1,mo_num
if(POPCNT(IAND(Idomo,IBSET(0_8,k-1))) .EQ. 1) then
nholes += 1
listholes(nholes) = k
holetype(nholes) = 2
endif
enddo
! find vmos
listvmos = -1
vmotype = -1
nvmos = 0
!do k = n_core_orb+1,n_core_orb + n_act_orb
do k = 1,mo_num
!print *,i,IBSET(0,i-1),POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))), POPCNT(IAND(Idomo,(IBSET(0_8,i-1))))
if(POPCNT(IAND(Isomo,(IBSET(0_8,k-1)))) .EQ. 0 .AND. POPCNT(IAND(Idomo,(IBSET(0_8,k-1)))) .EQ. 0) then
nvmos += 1
listvmos(nvmos) = k
vmotype(nvmos) = 0
else if(POPCNT(IAND(Isomo,(IBSET(0_8,k-1)))) .EQ. 1 .AND. POPCNT(IAND(Idomo,(IBSET(0_8,k-1)))) .EQ. 0 ) then
nvmos += 1
listvmos(nvmos) = k
vmotype(nvmos) = 1
end if
enddo
!print *,"I=",i
!call debug_spindet(psi_configuration(1,1,i),N_int)
!call debug_spindet(psi_configuration(1,2,i),N_int)
do k=1,nholes
p = listholes(k)
noccp = holetype(k)
! Calculate one-electron
! and two-electron coulomb terms
do l=1,nholes
q = listholes(l)
noccq = holetype(l)
!print *,"--------------- K=",p," L=",q
! one-electron term
if(p.EQ.q) then
hpp = noccq * h_core_ri(p,q)!mo_one_e_integrals(q,q)
else
hpp = 0.d0
endif
do j=starti,endi
! coulomb term
! (pp,qq) = <pq|pq>
if(p.EQ.q) then
diag_energies(j) += hpp !+ 0.5d0 * (noccp * noccq * mo_two_e_integral(p,q,p,q))
!print *,"hpp=",hpp,"diga= ",diag_energies(j)
! else
! diag_energies(j) += ! 0.5d0 * noccp * noccq * mo_two_e_integral(p,q,p,q)
! print *,"diga= ",diag_energies(j)
endif
enddo
enddo
enddo
enddo
end subroutine calculate_preconditioner_cfg
subroutine calculate_sigma_vector_cfg_nst(psi_out, psi_in, n_st, sze, istart, iend, ishift, istep)
implicit none
use bitmasks
BEGIN_DOC
! Documentation for sigma-vector calculation
!
! Calculates the result of the
! application of the hamiltonian to the
! wavefunction in CFG basis once
! TODO : Things prepare outside this routine
! 1. Touch the providers for
! a. ApqIJ containers
! b. DET to CSF transformation matrices
! 2. DET to CSF transcormation
! 2. CSF to DET back transcormation
! returns : psi_coef_out_det :
END_DOC
integer,intent(in) :: sze, istart,iend, istep, ishift, n_st
real*8,intent(in):: psi_in(sze,n_st)
real*8,intent(out):: psi_out(sze,n_st)
integer(bit_kind) :: Icfg(N_INT,2)
integer :: i,j,k,l,p,q,noccp,noccq, ii, jj, m, n, idxI, kk, nocck,orbk
integer(bit_kind) :: alphas_Icfg(N_INT,2,sze)
integer(bit_kind) :: singlesI(N_INT,2,sze)
integer(bit_kind) :: connectedI_alpha(N_INT,2,sze)
integer :: idxs_singlesI(sze)
integer :: idxs_connectedI_alpha(sze)
integer(bit_kind) :: psi_configuration_out(N_INT,2,sze)
real*8 :: psi_coef_out(n_CSF)
logical :: psi_coef_out_init(n_CSF)
integer :: excitationIds_single(2,sze)
integer :: excitationTypes_single(sze)
integer :: excitationIds(2,sze)
integer :: excitationTypes(sze)
real*8 :: diagfactors(sze)
integer :: nholes
integer :: nvmos
integer :: listvmos(mo_num)
integer :: vmotype(mo_num) ! 1 -> VMO 2 -> SOMO
integer :: listholes(mo_num)
integer :: holetype(mo_num) ! 1-> SOMO 2->DOMO
integer :: Nalphas_Icfg, nconnectedI, rowsikpq, colsikpq, nsinglesI
integer :: extype,NSOMOalpha,NSOMOI,NSOMOJ,pmodel,qmodel
integer :: getNSOMO
integer :: totcolsTKI
integer :: rowsTKI
integer :: noccpp
integer :: istart_cfg, iend_cfg
integer*8 :: MS, Isomo, Idomo, Jsomo, Jdomo, Ialpha, Ibeta
integer :: moi, moj, mok, mol, starti, endi, startj, endj, cnti, cntj, cntk
real*8 :: norm_coef_cfg, fac2eints
real*8 :: norm_coef_det
real*8 :: meCC1, meCC2, diagfac
real*8,dimension(:,:,:),allocatable :: TKI
real*8,dimension(:,:),allocatable :: GIJpqrs
real*8,dimension(:,:,:),allocatable :: TKIGIJ
real*8, external :: mo_two_e_integral
real*8, external :: get_two_e_integral
real*8 :: diag_energies(n_CSF)
call calculate_preconditioner_cfg(diag_energies)
MS = 0
norm_coef_cfg=0.d0
psi_out=0.d0
psi_coef_out_init = .False.
istart_cfg = psi_csf_to_config_data(istart)
iend_cfg = psi_csf_to_config_data(iend)
!!! Single Excitations !!!
do i=istart_cfg,iend_cfg
Icfg(1,1) = psi_configuration(1,1,i)
Icfg(1,2) = psi_configuration(1,2,i)
Isomo = Icfg(1,1)
Idomo = Icfg(1,2)
NSOMOI = getNSOMO(Icfg)
! find out all pq holes possible
nholes = 0
! holes in SOMO
! list_act
! list_core
! list_core_inact
! bitmasks
!do k = n_core_orb+1,n_core_orb + n_act_orb
do k = 1,mo_num
if(POPCNT(IAND(Isomo,IBSET(0_8,k-1))) .EQ. 1) then
nholes += 1
listholes(nholes) = k
holetype(nholes) = 1
endif
enddo
! holes in DOMO
!do k = n_core_orb+1,n_core_orb + n_act_orb
do k = 1,mo_num
if(POPCNT(IAND(Idomo,IBSET(0_8,k-1))) .EQ. 1) then
nholes += 1
listholes(nholes) = k
holetype(nholes) = 2
endif
enddo
! find vmos
listvmos = -1
vmotype = -1
nvmos = 0
!do k = n_core_orb+1,n_core_orb + n_act_orb
do k = 1,mo_num
!print *,i,IBSET(0,i-1),POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))), POPCNT(IAND(Idomo,(IBSET(0_8,i-1))))
if(POPCNT(IAND(Isomo,(IBSET(0_8,k-1)))) .EQ. 0 .AND. POPCNT(IAND(Idomo,(IBSET(0_8,k-1)))) .EQ. 0) then
nvmos += 1
listvmos(nvmos) = k
vmotype(nvmos) = 0
else if(POPCNT(IAND(Isomo,(IBSET(0_8,k-1)))) .EQ. 1 .AND. POPCNT(IAND(Idomo,(IBSET(0_8,k-1)))) .EQ. 0 ) then
nvmos += 1
listvmos(nvmos) = k
vmotype(nvmos) = 1
end if
enddo
! Icsf ids
starti = psi_config_data(i,1)
endi = psi_config_data(i,2)
NSOMOI = getNSOMO(Icfg)
call generate_all_singles_cfg_with_type(Icfg,singlesI,idxs_singlesI,excitationIds_single, &
excitationTypes_single,nsinglesI,N_int)
do j = 1,nsinglesI
idxI = idxs_singlesI(j)
NSOMOJ = getNSOMO(singlesI(:,:,j))
p = excitationIds_single(1,j)
q = excitationIds_single(2,j)
extype = excitationTypes_single(j)
! Off diagonal terms
call convertOrbIdsToModelSpaceIds(Icfg, singlesI(:,:,j), p, q, extype, pmodel, qmodel)
Jsomo = singlesI(1,1,j)
Jdomo = singlesI(1,2,j)
! Add the hole on J
if(POPCNT(IAND(Jsomo,IBSET(0_8,q-1))) .EQ. 1 .AND. POPCNT(IAND(Isomo,IBSET(0_8,q-1))) .EQ. 0) then
nholes += 1
listholes(nholes) = q
holetype(nholes) = 1
endif
if((POPCNT(IAND(Jdomo,IBSET(0_8,q-1))) .EQ. 1 .AND. POPCNT(IAND(Idomo,IBSET(0_8,q-1))) .EQ. 0) .AND. POPCNT(IAND(Isomo,IBSET(0_8,q-1))) .EQ. 0) then
nholes += 1
listholes(nholes) = q
holetype(nholes) = 2
endif
startj = psi_config_data(idxI,1)
endj = psi_config_data(idxI,2)
!!! One-electron contribution !!!
do kk = 1,n_st
cnti = 0
do ii = starti, endi
cnti += 1
cntj = 0
do jj = startj, endj
cntj += 1
meCC1 = AIJpqContainer(NSOMOI,NSOMOJ,extype,pmodel,qmodel,cnti,cntj)
psi_out(jj,kk) += meCC1 * psi_in(ii,kk) * h_core_ri(p,q)
psi_coef_out_init(jj) = .True.
enddo
enddo
enddo
! Undo setting in listholes
if(POPCNT(IAND(Jsomo,IBSET(0_8,q-1))) .EQ. 1 .AND. POPCNT(IAND(Isomo,IBSET(0_8,q-1))) .EQ. 0) then
nholes -= 1
endif
if((POPCNT(IAND(Jdomo,IBSET(0_8,q-1))) .EQ. 1 .AND. POPCNT(IAND(Idomo,IBSET(0_8,q-1))) .EQ. 0) .AND. POPCNT(IAND(Isomo,IBSET(0_8,q-1))) .EQ. 0) then
nholes -= 1
endif
enddo
enddo
!!! Double Excitations !!!
! Loop over all selected configurations
do i = istart_cfg,iend_cfg
Icfg(1,1) = psi_configuration(1,1,i)
Icfg(1,2) = psi_configuration(1,2,i)
starti = psi_config_data(i,1)
endi = psi_config_data(i,2)
! Returns all unique (checking the past) singly excited cfgs connected to I
call obtain_associated_alphaI(i, Icfg, alphas_Icfg, Nalphas_Icfg)
! TODO : remove doubly excited for return
! Here we do 2x the loop. One to count for the size of the matrix, then we compute.
do k = 1,Nalphas_Icfg
! Now generate all singly excited with respect to a given alpha CFG
call obtain_connected_I_foralpha(i,alphas_Icfg(1,1,k),connectedI_alpha,idxs_connectedI_alpha,nconnectedI,excitationIds,excitationTypes,diagfactors)
if(nconnectedI .EQ. 0) then
cycle
endif
totcolsTKI = 0
rowsTKI = -1
do j = 1,nconnectedI
NSOMOalpha = getNSOMO(alphas_Icfg(:,:,k))
NSOMOI = getNSOMO(connectedI_alpha(:,:,j))
p = excitationIds(1,j)
q = excitationIds(2,j)
extype = excitationTypes(j)
call convertOrbIdsToModelSpaceIds(alphas_Icfg(1,1,k), connectedI_alpha(1,1,j), p, q, extype, pmodel, qmodel)
! for E_pp E_rs and E_ppE_rr case
if(p.EQ.q) then
NSOMOalpha = NSOMOI
endif
rowsikpq = AIJpqMatrixDimsList(NSOMOalpha,NSOMOI,extype,pmodel,qmodel,1)
colsikpq = AIJpqMatrixDimsList(NSOMOalpha,NSOMOI,extype,pmodel,qmodel,2)
totcolsTKI += colsikpq
if(rowsTKI .LT. rowsikpq .AND. rowsTKI .NE. -1) then
print *,">",j,"Something is wrong in sigma-vector", rowsTKI, rowsikpq, "(p,q)=",pmodel,qmodel,"ex=",extype,"na=",NSOMOalpha," nI=",NSOMOI
!rowsTKI = rowsikpq
else
rowsTKI = rowsikpq
endif
enddo
allocate(TKI(rowsTKI,n_st,totcolsTKI)) ! coefficients of CSF
! Initialize the inegral container
! dims : (totcolsTKI, nconnectedI)
allocate(GIJpqrs(totcolsTKI,nconnectedI)) ! gpqrs
allocate(TKIGIJ(rowsTKI,n_st,nconnectedI)) ! gpqrs
totcolsTKI = 0
do j = 1,nconnectedI
NSOMOalpha = getNSOMO(alphas_Icfg(:,:,k))
NSOMOI = getNSOMO(connectedI_alpha(:,:,j))
p = excitationIds(1,j)
q = excitationIds(2,j)
extype = excitationTypes(j)
call convertOrbIdsToModelSpaceIds(alphas_Icfg(:,:,k), connectedI_alpha(:,:,j), p, q, extype, pmodel, qmodel)
rowsikpq = AIJpqMatrixDimsList(NSOMOalpha,NSOMOI,extype,pmodel,qmodel,1)
colsikpq = AIJpqMatrixDimsList(NSOMOalpha,NSOMOI,extype,pmodel,qmodel,2)
do kk = 1,n_st
do l = 1,rowsTKI
do m = 1,colsikpq
TKI(l,kk,totcolsTKI+m) = AIJpqContainer(NSOMOalpha,NSOMOI,extype,pmodel,qmodel,l,m) * psi_in(idxs_connectedI_alpha(j)+m-1,kk)
enddo
enddo
enddo
do m = 1,colsikpq
do l = 1,nconnectedI
! <ij|kl> = (ik|jl)
moi = excitationIds(1,j) ! p
mok = excitationIds(2,j) ! q
moj = excitationIds(2,l) ! s
mol = excitationIds(1,l) ! r
if(moi.EQ.mok .AND. moj.EQ.mol)then
diagfac = diagfactors(j)
diagfac *= diagfactors(l)
!print *,"integrals (",totcolsTKI+m,l,")",mok,moi,mol,moj, "|", diagfac
GIJpqrs(totcolsTKI+m,l) = diagfac*0.5d0*mo_two_e_integral(mok,mol,moi,moj) ! g(pq,sr) = <ps,qr>
else
diagfac = diagfactors(j)*diagfactors(l)
!print *,"integrals (",totcolsTKI+m,l,")",mok,moi,mol,moj, "|", diagfac
GIJpqrs(totcolsTKI+m,l) = diagfac*0.5d0*mo_two_e_integral(mok,mol,moi,moj) ! g(pq,sr) = <ps,qr>
!endif
endif
enddo
enddo
totcolsTKI += colsikpq
enddo
! Do big BLAS
! TODO TKI, size(TKI,1)*size(TKI,2)
call dgemm('N','N', rowsTKI*n_st, nconnectedI, totcolsTKI, 1.d0, &
TKI, size(TKI,1)*n_st, GIJpqrs, size(GIJpqrs,1), 0.d0, &
TKIGIJ , size(TKIGIJ,1)*n_st )
! Collect the result
totcolsTKI = 0
do j = 1,nconnectedI
NSOMOalpha = getNSOMO(alphas_Icfg(:,:,k))
NSOMOI = getNSOMO(connectedI_alpha(:,:,j))
p = excitationIds(1,j)
q = excitationIds(2,j)
extype = excitationTypes(j)
call convertOrbIdsToModelSpaceIds(alphas_Icfg(:,:,k), connectedI_alpha(:,:,j), p, q, extype, pmodel, qmodel)
rowsikpq = AIJpqMatrixDimsList(NSOMOalpha,NSOMOI,extype,pmodel,qmodel,1)
colsikpq = AIJpqMatrixDimsList(NSOMOalpha,NSOMOI,extype,pmodel,qmodel,2)
!print *,">j=",j,rowsikpq,colsikpq, ">>",totcolsTKI,",",idxs_connectedI_alpha(j)
do kk = 1,n_st
do m = 1,colsikpq
do l = 1,rowsTKI
psi_out(idxs_connectedI_alpha(j)+m-1,kk) += AIJpqContainer(NSOMOalpha,NSOMOI,extype,pmodel,qmodel,l,m) * TKIGIJ(l,kk,j)
psi_coef_out_init(idxs_connectedI_alpha(j)+m-1) = .True.
enddo
enddo
enddo
totcolsTKI += colsikpq
enddo
deallocate(TKI) ! coefficients of CSF
! Initialize the inegral container
! dims : (totcolsTKI, nconnectedI)
deallocate(GIJpqrs) ! gpqrs
deallocate(TKIGIJ) ! gpqrs
enddo ! loop over alphas
enddo ! loop over I
! Add the diagonal contribution
do i = 1,n_CSF
psi_out(i,1) += 1.0d0*diag_energies(i)*psi_in(i,1)
enddo
end subroutine calculate_sigma_vector_cfg_nst

29
src/davidson/diagonalization_hcsf_dressed.irp.f
View File

@ -88,7 +88,7 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
double precision, intent(out) :: energies(N_st_diag_in)
integer :: iter, N_st_diag
integer :: i,j,k,l,m
integer :: i,j,k,l,m,kk
logical, intent(inout) :: converged
double precision, external :: u_dot_v, u_dot_u
@ -285,7 +285,7 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
! Make random verctors eigenstates of S2
call convertWFfromDETtoCSF(N_st_diag,U,U_csf)
call convertWFfromCSFtoDET(N_st_diag,U_csf,U)
!call convertWFfromCSFtoDET(N_st_diag,U_csf,U)
do while (.not.converged)
itertot = itertot+1
@ -302,11 +302,28 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
! Compute |W_k> = \sum_i |i><i|H|u_k>
! -----------------------------------
call convertWFfromCSFtoDET(N_st_diag,U_csf(1,shift+1),U)
!call convertWFfromCSFtoDET(N_st_diag,U_csf(1,shift+1),U)
if ((sze > 100000).and.distributed_davidson) then
call H_u_0_nstates_zmq (W,U,N_st_diag,sze)
!call convertWFfromCSFtoDET(N_st_diag,U_csf(1,shift+1),U)
!call convertWFfromCSFtoDET(N_st_diag,W_csf(1,shift+1),W)
!call H_u_0_nstates_zmq (W,U,N_st_diag,sze)
!call convertWFfromDETtoCSF(N_st_diag,U,U_csf(1,shift+1))
!call convertWFfromDETtoCSF(N_st_diag,W,W_csf(1,shift+1))
!call calculate_sigma_vector_cfg_nst(W_csf(1,shift+1),U_csf(1,shift+1),N_st_diag,sze_csf,1,sze_csf,0,1)
do kk=1,N_st_diag
call calculate_sigma_vector_cfg_nst(W_csf(1,shift+kk),U_csf(1,shift+kk),1,sze_csf,1,sze_csf,0,1)
enddo
else
call H_u_0_nstates_openmp(W,U,N_st_diag,sze)
!call convertWFfromCSFtoDET(N_st_diag,U_csf(1,shift+1),U)
!call convertWFfromCSFtoDET(N_st_diag,W_csf(1,shift+1),W)
!call H_u_0_nstates_openmp(W,U,N_st_diag,sze)
!call convertWFfromDETtoCSF(N_st_diag,U,U_csf(1,shift+1))
!call convertWFfromDETtoCSF(N_st_diag,W,W_csf(1,shift+1))
!call calculate_sigma_vector_cfg_nst(W_csf(1,shift+1),U_csf(1,shift+1),N_st_diag,sze_csf,1,sze_csf,0,1)
do kk=1,N_st_diag
call calculate_sigma_vector_cfg_nst(W_csf(1,shift+kk),U_csf(1,shift+kk),1,sze_csf,1,sze_csf,0,1)
enddo
endif
else
! Already computed in update below
@ -350,7 +367,7 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
endif
endif
call convertWFfromDETtoCSF(N_st_diag,W,W_csf(1,shift+1))
!call convertWFfromDETtoCSF(N_st_diag,W,W_csf(1,shift+1))
! Compute h_kl = <u_k | W_l> = <u_k| H |u_l>
! -------------------------------------------