PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-07-22 10:47:38 +02:00
Code Releases Activity

Fixed some bugs. Diagonal energy is OK. Some bugs still present.

Browse Source
This commit is contained in:
v1j4y 2022-12-15 17:15:30 +01:00
parent f291078945
commit 5622b9790d
5 changed files with 329 additions and 116 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -524,7 +524,7 @@ use bitmasks
!!!
!!!END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), alphasIcfg_list , (N_int,2,N_configuration,mo_num*(mo_num))]
BEGIN_PROVIDER [ integer(bit_kind), alphasIcfg_list , (N_int,2,N_configuration,mo_num*12)]
&BEGIN_PROVIDER [ integer, NalphaIcfg_list, (N_configuration) ]
implicit none
!use bitmasks
@ -549,7 +549,8 @@ use bitmasks
integer(bit_kind) :: Isomo(N_int), Isomop(N_int), Isomoq(N_int)
integer(bit_kind) :: Jdomo(N_int), Jdomop(N_int), Jdomoq(N_int)
integer(bit_kind) :: Jsomo(N_int), Jsomop(N_int), Jsomoq(N_int)
integer(bit_kind) :: diffDOMO(N_int), xordiffSOMODOMO(N_int), diffSOMO(N_int)
!integer(bit_kind) :: diffDOMO(N_int), xordiffSOMODOMO(N_int), diffSOMO(N_int)
integer(bit_kind) :: diffDOMO, xordiffSOMODOMO, diffSOMO
integer :: ndiffSOMO
integer :: ndiffDOMO
integer :: nxordiffSOMODOMO
@ -674,13 +675,15 @@ use bitmasks
Jsomo(jint) = IBCLR(Jsomo(jint),jpos)
endif
Nsomo_J=0
do ii=1, N_int
Jcfg(ii,1) = Jsomo(ii)
Jcfg(ii,2) = Jdomo(ii)
Nsomo_J += POPCNT(Jsomo(ii))
enddo
call bitstring_to_list(Jcfg,listall,nelall,N_int)
Nsomo_J = nelall
!call bitstring_to_list(Jcfg(1,1),listall,nelall,N_int)
!Nsomo_J = nelall
! Check for Minimal alpha electrons (MS)
if(Nsomo_J.lt.MS)then
@ -705,15 +708,13 @@ use bitmasks
ndiffDOMO = 0
nxordiffSOMODOMO = 0
do ii = 1, N_int
Jsomo = Jcfg(ii,1)
Jdomo = Jcfg(ii,2)
diffSOMO = IEOR(Jsomo,iand(reunion_of_act_virt_bitmask(ii,1),psi_configuration(ii,1,k)))
ndiffSOMO += POPCNT(diffSOMO(ii))
diffDOMO = IEOR(Jdomo,iand(reunion_of_act_virt_bitmask(ii,2),psi_configuration(ii,2,k)))
diffSOMO = IEOR(Jcfg(ii,1),iand(reunion_of_act_virt_bitmask(ii,1),psi_configuration(ii,1,k)))
ndiffSOMO += POPCNT(diffSOMO)
diffDOMO = IEOR(Jcfg(ii,2),iand(reunion_of_act_virt_bitmask(ii,2),psi_configuration(ii,2,k)))
xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO)
ndiffDOMO += POPCNT(diffDOMO(ii))
nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO(ii))
nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO
ndiffDOMO += POPCNT(diffDOMO)
nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO)
nxordiffSOMODOMO += POPCNT(diffSOMO) + POPCNT(diffDOMO)
end do
if((ndiffSOMO .ne. 0) .and. (ndiffSOMO .ne. 2)) cycle
@ -818,6 +819,7 @@ use bitmasks
if(Nsomo_J .ge. NSOMOMin) then
!print *," Idx = ",idxI, "p = ",pp, " q = ",qq," Jsomo=",Jsomo(1), " Jdomo=",IOR(Jdomo(1),ISHFT(1_8,n_core_orb)-1)
NalphaIcfg += 1
!print *," Idx = ",idxI, " Nalpha=",NalphaIcfg
alphasIcfg_list(:,1,idxI,NalphaIcfg) = Jcfg(:,1)
if(n_core_orb .le. 63)then
alphasIcfg_list(1,2,idxI,NalphaIcfg) = IOR(Jcfg(1,2),ISHFT(1_8,n_core_orb)-1)
@ -864,7 +866,7 @@ use bitmasks
ndiffSOMO = 0
do ii=1,N_int
diffSOMO = IEOR(Icfg(ii,1),iand(act_bitmask(ii,1),psi_configuration(ii,1,k)))
ndiffSOMO += POPCNT(diffSOMO(ii))
ndiffSOMO += POPCNT(diffSOMO)
end do
! ndiffSOMO cannot be 0 (I /= k)
! if ndiffSOMO /= 2 then it has to be greater than 2 and hense
@ -875,8 +877,8 @@ use bitmasks
do ii=1,N_int
diffDOMO = IEOR(Icfg(ii,2),iand(act_bitmask(ii,1),psi_configuration(ii,2,k)))
xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO)
ndiffDOMO += POPCNT(diffDOMO(ii))
nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO(ii))
ndiffDOMO += POPCNT(diffDOMO)
nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO)
end do
if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4)) then
ppExistsQ = .TRUE.
@ -926,7 +928,8 @@ END_PROVIDER
integer(bit_kind) :: Isomo(N_int), Isomop(N_int), Isomoq(N_int)
integer(bit_kind) :: Jdomo(N_int), Jdomop(N_int), Jdomoq(N_int)
integer(bit_kind) :: Jsomo(N_int), Jsomop(N_int), Jsomoq(N_int)
integer(bit_kind) :: diffDOMO(N_int), xordiffSOMODOMO(N_int), diffSOMO(N_int)
!integer(bit_kind) :: diffDOMO(N_int), xordiffSOMODOMO(N_int), diffSOMO(N_int)
integer(bit_kind) :: diffDOMO, xordiffSOMODOMO, diffSOMO
integer :: ndiffSOMO
integer :: ndiffDOMO
integer :: nxordiffSOMODOMO
@ -1084,15 +1087,13 @@ END_PROVIDER
ndiffDOMO = 0
nxordiffSOMODOMO = 0
do ii = 1, N_int
Jsomo = Jcfg(ii,1)
Jdomo = Jcfg(ii,2)
diffSOMO = IEOR(Jsomo,iand(reunion_of_act_virt_bitmask(ii,1),psi_configuration(ii,1,k)))
ndiffSOMO += POPCNT(diffSOMO(ii))
diffDOMO = IEOR(Jdomo,iand(reunion_of_act_virt_bitmask(ii,2),psi_configuration(ii,2,k)))
diffSOMO = IEOR(Jcfg(ii,1),iand(reunion_of_act_virt_bitmask(ii,1),psi_configuration(ii,1,k)))
ndiffSOMO += POPCNT(diffSOMO)
diffDOMO = IEOR(Jcfg(ii,2),iand(reunion_of_act_virt_bitmask(ii,2),psi_configuration(ii,2,k)))
xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO)
ndiffDOMO += POPCNT(diffDOMO(ii))
nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO(ii))
nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO
ndiffDOMO += POPCNT(diffDOMO)
nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO)
nxordiffSOMODOMO += POPCNT(diffSOMO) + POPCNT(diffDOMO)
end do
if((ndiffSOMO .NE. 0) .AND. (ndiffSOMO .NE. 2)) cycle
!if(POPCNT(IEOR(diffSOMO,diffDOMO)) .LE. 1 .AND. ndiffDOMO .LT. 3) then
@ -1211,7 +1212,7 @@ END_PROVIDER
do k = 1, idxI-1
do ii=1,N_int
diffSOMO = IEOR(Icfg(ii,1),iand(act_bitmask(ii,1),psi_configuration(ii,1,k)))
ndiffSOMO += POPCNT(diffSOMO(ii))
ndiffSOMO += POPCNT(diffSOMO)
end do
! ndiffSOMO cannot be 0 (I /= k)
! if ndiffSOMO /= 2 then it has to be greater than 2 and hense
@ -1222,8 +1223,8 @@ END_PROVIDER
do ii=1,N_int
diffDOMO = IEOR(Icfg(ii,2),iand(act_bitmask(ii,1),psi_configuration(ii,2,k)))
xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO)
ndiffDOMO += POPCNT(diffDOMO(ii))
nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO(ii))
ndiffDOMO += POPCNT(diffDOMO)
nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO)
end do
if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then
ppExistsQ = .TRUE.

1
src/csf/conversion.irp.f
View File

@ -114,6 +114,7 @@ subroutine convertWFfromCSFtoDET(N_st,psi_coef_cfg_in, psi_coef_det)
integer :: idx
integer MS
MS = elec_alpha_num-elec_beta_num
print *,"size=",size(tmp_psi_coef_det,1)," ",size(tmp_psi_coef_det,2)
countcsf = 0

72
src/csf/obtain_I_foralpha.irp.f
View File

@ -102,7 +102,7 @@ subroutine obtain_connected_J_givenI(idxI, givenI, connectedI, idxs_connectedI,
xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO)
ndiffDOMO += POPCNT(diffDOMO)
nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO)
nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO
nxordiffSOMODOMO += POPCNT(diffSOMO) + POPCNT(diffDOMO)
end do
if((nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then
@ -243,13 +243,16 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI
integer :: listholes(mo_num)
integer :: holetype(mo_num)
integer :: end_index, ishift
integer :: Nsomo_alpha, pp,qq, nperm
integer :: Nsomo_alpha, pp,qq, nperm, iint, ipos
integer*8 :: MS
integer :: exc(0:2,2,2), tz, m, n, high, low
integer :: listall(N_int*bit_kind_size), nelall
integer :: nconnectedExtradiag, nconnectedDiag
integer(bit_kind) :: hole, particle, tmp
MS = elec_alpha_num-elec_beta_num
nconnectedExtradiag=0
nconnectedDiag=0
nconnectedI = 0
end_index = N_configuration
@ -260,10 +263,12 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI
!Nsomo_alpha = POPCNT(Isomo)
Icfg = Ialpha
Nsomo_alpha = 0
!print *," Ialpha="
do i=1,N_int
Isomo = Ialpha(i,1)
Idomo = Ialpha(i,2)
Nsomo_alpha += POPCNT(Isomo)
!print *,Isomo, Idomo, "Nsomo=",Nsomo_alpha
end do
end_index = min(N_configuration,cfg_seniority_index(min(Nsomo_alpha+4,elec_num))-1)
if(end_index .LT. 0) end_index= N_configuration
@ -293,20 +298,25 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI
ndiffSOMO = 0
ndiffDOMO = 0
nxordiffSOMODOMO = 0
nsomoJ=0
nsomoalpha=0
do ii=1,N_int
Isomo = Ialpha(ii,1)
Idomo = Ialpha(ii,2)
Jsomo = psi_configuration(ii,1,i)
Jdomo = psi_configuration(ii,2,i)
nsomoJ += POPCNT(Jsomo)
nsomoalpha += POPCNT(Isomo)
diffSOMO = IEOR(Isomo,Jsomo)
ndiffSOMO += POPCNT(diffSOMO)
diffDOMO = IEOR(Idomo,Jdomo)
xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO)
ndiffDOMO += POPCNT(diffDOMO)
nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO)
nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO
nxordiffSOMODOMO += POPCNT(diffSOMO) + POPCNT(diffDOMO)
end do
Jcfg = psi_configuration(:,:,i)
!Jcfg = psi_configuration(:,:,i)
!print *,"nxordiffSOMODOMO(4)=",nxordiffSOMODOMO, " ndiffSOMO(2)=",ndiffSOMO
if((nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then
select case(ndiffDOMO)
@ -328,7 +338,7 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI
Jsomo = psi_configuration(ii,1,i)
IJsomo = IEOR(Isomo, Jsomo)
if(popcnt(IAND(Isomo,IJsomo)) > 0)then
p = TRAILZ(IAND(Isomo,IJsomo)) + 1 + ii * bit_kind_size
p = TRAILZ(IAND(Isomo,IJsomo)) + 1 + (ii-1) * bit_kind_size
EXIT
endif
end do
@ -337,22 +347,24 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI
Isomo = Ialpha(ii,1)
Jsomo = psi_configuration(ii,1,i)
IJsomo = IEOR(Isomo, Jsomo)
IJsomo = IBCLR(IJsomo,p-1)
iint = shiftr(p-1,bit_kind_shift) + 1
ipos = p-shiftl((iint-1),bit_kind_shift)
if(iint .eq. ii)then
IJsomo = IBCLR(IJsomo,ipos-1)
endif
if(popcnt(IJsomo) > 0)then
q = TRAILZ(IJsomo) + 1 + ii * bit_kind_size
q = TRAILZ(IJsomo) + 1 + (ii-1) * bit_kind_size
EXIT
endif
enddo
endif
!assert ( p == pp)
!assert ( q == qq)
!print *," --- p=",p," q=",q
!print *," 1--- p=",p," q=",q
case (1)
! DOMO -> VMO
! or
! SOMO -> SOMO
nsomoJ = POPCNT(Jsomo)
nsomoalpha = POPCNT(Isomo)
if(nsomoJ .GT. nsomoalpha) then
! DOMO -> VMO
!print *,"obt DOMO -> VMO"
@ -371,7 +383,7 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI
Idomo = Ialpha(ii,2)
Jdomo = psi_configuration(ii,2,i)
if(popcnt(IEOR(Idomo,Jdomo)) > 0)then
p = TRAILZ(IEOR(Idomo,Jdomo)) + 1 + ii * bit_kind_size
p = TRAILZ(IEOR(Idomo,Jdomo)) + 1 + (ii-1) * bit_kind_size
EXIT
endif
end do
@ -380,9 +392,13 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI
Isomo = Ialpha(ii,1)
Jsomo = psi_configuration(ii,1,i)
Isomo = IEOR(Isomo, Jsomo)
Isomo = IBCLR(Isomo,p-1)
iint = shiftr(p-1,bit_kind_shift) + 1
ipos = p-shiftl((iint-1),bit_kind_shift)
if(iint .eq. ii)then
Isomo = IBCLR(Isomo,ipos-1)
endif
if(popcnt(Isomo) > 0)then
q = TRAILZ(Isomo) + 1 + ii * bit_kind_size
q = TRAILZ(Isomo) + 1 + (ii-1) * bit_kind_size
EXIT
endif
end do
@ -404,13 +420,16 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI
!endif
else
! Find p
!print *,"Ialpha somo=",Ialpha(1,1), Ialpha(2,1)," Ialpha domo=",Ialpha(1,2), Ialpha(2,2)
!print *,"J somo=",psi_configuration(1,1,i), psi_configuration(2,1,i)," J domo=",psi_configuration(1,2,i),&
!psi_configuration(2,2,i)
do ii=1,N_int
Isomo = Ialpha(ii,1)
Jsomo = psi_configuration(ii,1,i)
Idomo = Ialpha(ii,2)
Jdomo = psi_configuration(ii,2,i)
if(popcnt(IEOR(Idomo,Jdomo)) > 0)then
q = TRAILZ(IEOR(Idomo,Jdomo)) + 1 + ii * bit_kind_size
q = TRAILZ(IEOR(Idomo,Jdomo)) + 1 + (ii-1) * bit_kind_size
EXIT
endif
enddo
@ -419,9 +438,14 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI
Isomo = Ialpha(ii,1)
Jsomo = psi_configuration(ii,1,i)
Isomo = IEOR(Isomo, Jsomo)
Isomo = IBCLR(Isomo,q-1)
iint = shiftr(q-1,bit_kind_shift) + 1
ipos = q-shiftl((iint-1),bit_kind_shift)
if(iint .eq. ii)then
Isomo = IBCLR(Isomo,ipos-1)
endif
!print *,"ii=",ii," Isomo=",Isomo
if(popcnt(Isomo) > 0)then
p = TRAILZ(Isomo) + 1 + ii * bit_kind_size
p = TRAILZ(Isomo) + 1 + (ii-1) * bit_kind_size
EXIT
endif
enddo
@ -429,6 +453,7 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI
!assert ( p == pp)
!assert ( q == qq)
endif
!print *," 2--- p=",p," q=",q
case (2)
! DOMO -> SOMO
!print *,"obt DOMO -> SOMO"
@ -447,7 +472,7 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI
Jdomo = psi_configuration(ii,2,i)
IJsomo = IEOR(Isomo, Jsomo)
if(popcnt(IAND(Jsomo,IJsomo)) > 0)then
p = TRAILZ(IAND(Jsomo,IJsomo)) + 1 + ii * bit_kind_size
p = TRAILZ(IAND(Jsomo,IJsomo)) + 1 + (ii-1) * bit_kind_size
EXIT
endif
enddo
@ -456,20 +481,26 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI
Isomo = Ialpha(ii,1)
Jsomo = psi_configuration(ii,1,i)
IJsomo = IEOR(Isomo, Jsomo)
IJsomo = IBCLR(IJsomo,p-1)
iint = shiftr(p-1,bit_kind_shift) + 1
ipos = p-shiftl((iint-1),bit_kind_shift)
if(iint .eq. ii)then
IJsomo = IBCLR(IJsomo,ipos-1)
endif
if(popcnt(IJsomo) > 0)then
q = TRAILZ(IJsomo) + 1 + ii * bit_kind_size
q = TRAILZ(IJsomo) + 1 + (ii-1) * bit_kind_size
EXIT
endif
enddo
endif
!assert ( p == pp)
!assert ( q == qq)
!print *," 3--- p=",p," q=",q
case default
print *,"something went wront in get connectedI"
end select
starti = psi_config_data(i,1)
endi = psi_config_data(i,2)
nconnectedExtradiag+=1
nconnectedI += 1
do k=1,N_int
connectedI(k,1,nconnectedI) = psi_configuration(k,1,i)
@ -534,6 +565,7 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI
if(holetype(k) .EQ. 1) then
starti = psi_config_data(i,1)
endi = psi_config_data(i,2)
nconnectedDiag+=1
nconnectedI += 1
connectedI(:,:,nconnectedI) = psi_configuration(:,:,i)
idxs_connectedI(nconnectedI)=starti
@ -544,6 +576,7 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI
else
starti = psi_config_data(i,1)
endi = psi_config_data(i,2)
nconnectedDiag+=1
nconnectedI += 1
connectedI(:,:,nconnectedI) = psi_configuration(:,:,i)
idxs_connectedI(nconnectedI)=starti
@ -556,5 +589,6 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI
enddo
endif
end do
!print *,"nconnectedExtradiag=",nconnectedExtradiag," nconnectedDiad=",nconnectedDiag
end subroutine obtain_connected_I_foralpha

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

@ -885,7 +885,7 @@ subroutine calculate_preconditioner_cfg(diag_energies)
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))
!NSOMOI = getNSOMO(psi_configuration(:,:,i))
starti = psi_config_data(i,1)
endi = psi_config_data(i,2)
@ -1239,27 +1239,34 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod
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(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 :: iint, ipos, ii
!integer(bit_kind) :: Isomotmp(N_int)
!integer(bit_kind) :: Jsomotmp(N_int)
integer*8 :: Isomotmp
integer*8 :: Jsomotmp
integer :: pos0,pos0prev
integer :: tmpp, tmpq
! TODO Flag (print) when model space indices is > 64
Isomo = Ialpha(1,1)
Idomo = Ialpha(1,2)
Jsomo = Jcfg(1,1)
Jdomo = Jcfg(1,2)
do ii=1,N_int
!Isomo = Ialpha(ii,1)
!Idomo = Ialpha(ii,2)
!Jsomo = Jcfg(ii,1)
!Jdomo = Jcfg(ii,2)
Isomo(ii) = Ialpha(ii,1)
Idomo(ii) = Ialpha(ii,2)
Jsomo(ii) = Jcfg(ii,1)
Jdomo(ii) = Jcfg(ii,2)
end do
pos0prev = 0
pmodel = p
qmodel = q
@ -1273,40 +1280,139 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod
! SOMO -> SOMO
! remove all domos
!print *,"type -> SOMO -> SOMO"
mask = ISHFT(1_8,p) - 1
Isomotmp = IAND(Isomo,mask)
pmodel = POPCNT(mask) - POPCNT(XOR(Isomotmp,mask))
mask = ISHFT(1_8,q) - 1
Isomotmp = IAND(Isomo,mask)
qmodel = POPCNT(mask) - POPCNT(XOR(Isomotmp,mask))
!mask = ISHFT(1_8,p) - 1
!Isomotmp = IAND(Isomo,mask)
!pmodel = POPCNT(mask) - POPCNT(XOR(Isomotmp,mask))
!mask = ISHFT(1_8,q) - 1
!Isomotmp = IAND(Isomo,mask)
!qmodel = POPCNT(mask) - POPCNT(XOR(Isomotmp,mask))
iint = shiftr(p-1,bit_kind_shift) + 1
ipos = p-shiftl((iint-1),bit_kind_shift)-1
tmpp = 0
!print *,"iint=",iint, " p=",p
do ii=1,iint-1
mask = ISHFT(1_bit_kind,-1)-1_bit_kind
Isomotmp = IAND(Isomo(ii),mask)
tmpp += POPCNT(mask) - POPCNT(XOR(Isomotmp,mask))
end do
mask = ISHFT(1_bit_kind,ipos+1) - 1
Isomotmp = IAND(Isomo(iint),mask)
pmodel = tmpp + POPCNT(mask) - POPCNT(XOR(Isomotmp,mask))
!print *,"iint=",iint, " ipos=",ipos,"pmodel=",pmodel, XOR(Isomotmp,mask),Isomo(iint)
iint = shiftr(q-1,bit_kind_shift) + 1
ipos = q-shiftl((iint-1),bit_kind_shift)-1
tmpq = 0
do ii=1,iint-1
mask = ISHFT(1_bit_kind,-1)-1_bit_kind
Isomotmp = IAND(Isomo(ii),mask)
tmpq += POPCNT(mask) - POPCNT(XOR(Isomotmp,mask))
end do
mask = ISHFT(1_bit_kind,ipos+1) - 1
Isomotmp = IAND(Isomo(iint),mask)
qmodel = tmpq + POPCNT(mask) - POPCNT(XOR(Isomotmp,mask))
!print *,"iint=",iint, " ipos=",ipos,"qmodel=",qmodel
case (2)
! DOMO -> VMO
! remove all domos except one at p
!print *,"type -> DOMO -> VMO"
mask = ISHFT(1_8,p) - 1
Jsomotmp = IAND(Jsomo,mask)
pmodel = POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask))
mask = ISHFT(1_8,q) - 1
Jsomotmp = IAND(Jsomo,mask)
qmodel = POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask))
!mask = ISHFT(1_8,p) - 1
!Jsomotmp = IAND(Jsomo,mask)
!pmodel = POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask))
!mask = ISHFT(1_8,q) - 1
!Jsomotmp = IAND(Jsomo,mask)
!qmodel = POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask))
iint = shiftr(p-1,bit_kind_shift) + 1
ipos = p-shiftl((iint-1),bit_kind_shift)-1
tmpp = 0
do ii=1,iint-1
mask = ISHFT(1_bit_kind,-1)-1_bit_kind
Jsomotmp = IAND(Jsomo(ii),mask)
tmpp += POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask))
end do
mask = ISHFT(1_bit_kind,ipos+1) - 1
Jsomotmp = IAND(Jsomo(iint),mask)
pmodel = tmpp + POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask))
iint = shiftr(q-1,bit_kind_shift) + 1
ipos = q-shiftl((iint-1),bit_kind_shift)-1
tmpq = 0
do ii=1,iint-1
mask = ISHFT(1_bit_kind,-1)-1_bit_kind
Jsomotmp = IAND(Jsomo(ii),mask)
tmpq += POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask))
end do
mask = ISHFT(1_bit_kind,ipos+1) - 1
Jsomotmp = IAND(Jsomo(iint),mask)
qmodel = tmpq + POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask))
case (3)
! SOMO -> VMO
!print *,"type -> SOMO -> VMO"
!Isomo = IEOR(Isomo,Jsomo)
if(p.LT.q) then
mask = ISHFT(1_8,p) - 1
Isomo = IAND(Isomo,mask)
pmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask))
mask = ISHFT(1_8,q) - 1
Jsomo = IAND(Jsomo,mask)
qmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) + 1
!mask = ISHFT(1_8,p) - 1
!Isomo = IAND(Isomo,mask)
!pmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask))
!mask = ISHFT(1_8,q) - 1
!Jsomo = IAND(Jsomo,mask)
!qmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) + 1
iint = shiftr(p-1,bit_kind_shift) + 1
ipos = p-shiftl((iint-1),bit_kind_shift)-1
tmpp = 0
do ii=1,iint-1
mask = ISHFT(1_bit_kind,-1)-1_bit_kind
Isomotmp = IAND(Isomo(ii),mask)
tmpp += POPCNT(mask) - POPCNT(XOR(Isomotmp,mask))
end do
mask = ISHFT(1_8,ipos+1) - 1
Isomotmp = IAND(Isomo(iint),mask)
pmodel = tmpp + POPCNT(mask) - POPCNT(XOR(Isomotmp,mask))
iint = shiftr(q-1,bit_kind_shift) + 1
ipos = q-shiftl((iint-1),bit_kind_shift)-1
tmpq = 0
do ii=1,iint-1
mask = ISHFT(1_bit_kind,-1)-1_bit_kind
Jsomotmp = IAND(Jsomo(ii),mask)
tmpq += POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask))
end do
mask = ISHFT(1_bit_kind,ipos+1) - 1
Jsomotmp = IAND(Jsomo(iint),mask)
qmodel = tmpq + POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) + 1
else
mask = ISHFT(1_8,p) - 1
Isomo = IAND(Isomo,mask)
pmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) + 1
mask = ISHFT(1_8,q) - 1
Jsomo = IAND(Jsomo,mask)
qmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask))
!mask = ISHFT(1_8,p) - 1
!Isomo = IAND(Isomo,mask)
!pmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) + 1
!mask = ISHFT(1_8,q) - 1
!Jsomo = IAND(Jsomo,mask)
!qmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask))
iint = shiftr(p-1,bit_kind_shift) + 1
ipos = p-shiftl((iint-1),bit_kind_shift)-1
tmpp = 0
do ii=1,iint-1
mask = ISHFT(1_bit_kind,-1)-1_bit_kind
Isomotmp = IAND(Isomo(ii),mask)
tmpp += POPCNT(mask) - POPCNT(XOR(Isomotmp,mask))
end do
mask = ISHFT(1_bit_kind,ipos+1) - 1
Isomotmp = IAND(Isomo(iint),mask)
pmodel = tmpp + POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) + 1
iint = shiftr(q-1,bit_kind_shift) + 1
ipos = q-shiftl((iint-1),bit_kind_shift)-1
tmpq = 0
do ii=1,iint-1
mask = ISHFT(1_bit_kind,-1)-1_bit_kind
Jsomotmp = IAND(Jsomo(ii),mask)
tmpq += POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask))
end do
mask = ISHFT(1_bit_kind,ipos+1) - 1
Jsomotmp = IAND(Jsomo(iint),mask)
qmodel = tmpq + POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask))
endif
case (4)
! DOMO -> SOMO
@ -1314,19 +1420,67 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod
!print *,"type -> DOMO -> SOMO"
!Isomo = IEOR(Isomo,Jsomo)
if(p.LT.q) then
mask = ISHFT(1_8,p) - 1
Jsomo = IAND(Jsomo,mask)
pmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask))
mask = ISHFT(1_8,q) - 1
Isomo = IAND(Isomo,mask)
qmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) + 1
!mask = ISHFT(1_8,p) - 1
!Jsomo = IAND(Jsomo,mask)
!pmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask))
!mask = ISHFT(1_8,q) - 1
!Isomo = IAND(Isomo,mask)
!qmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) + 1
iint = shiftr(p-1,bit_kind_shift) + 1
ipos = p-shiftl((iint-1),bit_kind_shift)-1
tmpp = 0
do ii=1,iint-1
mask = ISHFT(1_bit_kind,-1)-1_bit_kind
Jsomotmp = IAND(Jsomo(ii),mask)
tmpp += POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask))
end do
mask = ISHFT(1_bit_kind,ipos+1) - 1
Jsomotmp = IAND(Jsomo(iint),mask)
pmodel = tmpp + POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask))
iint = shiftr(q-1,bit_kind_shift) + 1
ipos = q-shiftl((iint-1),bit_kind_shift)-1
tmpq = 0
do ii=1,iint-1
mask = ISHFT(1_bit_kind,-1)-1_bit_kind
Isomotmp = IAND(Isomo(ii),mask)
tmpq += POPCNT(mask) - POPCNT(XOR(Isomotmp,mask))
end do
mask = ISHFT(1_bit_kind,ipos+1) - 1
Isomotmp = IAND(Isomo(iint),mask)
qmodel = tmpq + POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) + 1
else
mask = ISHFT(1_8,p) - 1
Jsomo = IAND(Jsomo,mask)
pmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) + 1
mask = ISHFT(1_8,q) - 1
Isomo = IAND(Isomo,mask)
qmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask))
!mask = ISHFT(1_8,p) - 1
!Jsomo = IAND(Jsomo,mask)
!pmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) + 1
!mask = ISHFT(1_8,q) - 1
!Isomo = IAND(Isomo,mask)
!qmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask))
iint = shiftr(p-1,bit_kind_shift) + 1
ipos = p-shiftl((iint-1),bit_kind_shift)-1
tmpp = 0
do ii=1,iint-1
mask = ISHFT(1_bit_kind,-1)-1_bit_kind
Jsomotmp = IAND(Jsomo(ii),mask)
tmpp += POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask))
end do
mask = ISHFT(1_bit_kind,ipos+1) - 1
Jsomotmp = IAND(Jsomo(iint),mask)
pmodel = tmpp + POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) + 1
iint = shiftr(q-1,bit_kind_shift) + 1
ipos = q-shiftl((iint-1),bit_kind_shift)-1
tmpq = 0
do ii=1,iint-1
mask = ISHFT(1_bit_kind,-1)-1_bit_kind
Isomotmp = IAND(Isomo(ii),mask)
tmpq += POPCNT(mask) - POPCNT(XOR(Isomotmp,mask))
end do
mask = ISHFT(1_bit_kind,ipos+1) - 1
Isomotmp = IAND(Isomo(iint),mask)
qmodel = tmpq + POPCNT(mask) - POPCNT(XOR(Isomotmp,mask))
endif
case default
print *,"something is wrong in convertOrbIdsToModelSpaceIds"
@ -1415,7 +1569,7 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
allocate(diag_energies(n_CSF))
call calculate_preconditioner_cfg(diag_energies)
!print *," diag energy =",diag_energies(1)
print *," diag energy =",diag_energies(1)
MS = 0
norm_coef_cfg=0.d0
@ -1590,6 +1744,7 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
call omp_set_lock(lock(jj))
do kk = 1,n_st
psi_out(kk,jj) = psi_out(kk,jj) + meCC1 * psi_in(kk,ii)
print *,"jj=",jj,'psi_out(kk)=',psi_out(kk,jj)
enddo
call omp_unset_lock(lock(jj))
enddo
@ -1666,7 +1821,8 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
nconnectedI, excitationIds, excitationTypes, diagfactors)
!if(i .EQ. 1) then
! print *,'k=',k,' kcfgSOMO=',alphas_Icfg(1,1,k),' ',POPCNT(alphas_Icfg(1,1,k)),' kcfgDOMO=',alphas_Icfg(1,2,k),' ',POPCNT(alphas_Icfg(1,2,k))
! print *,'k=',k,' kcfgSOMO=',alphas_Icfg(1,1,k),alphas_Icfg(2,1,k),' ',POPCNT(alphas_Icfg(1,1,k)),' &
! kcfgDOMO=',alphas_Icfg(1,2,k),alphas_Icfg(2,2,k),' ',POPCNT(alphas_Icfg(1,2,k)), " NconnectedI=",nconnectedI
!endif
@ -1682,9 +1838,13 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
NSOMOI = getNSOMO(connectedI_alpha(:,:,j))
p = excitationIds(1,j)
q = excitationIds(2,j)
!print *,"j=",j, " p=",p," q=",q
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(k.eq.722)then
! print *,"j=",j," k=",k,"p=",p,"q=",q,"NSOMOalpha=",NSOMOalpha, "pmodel=",pmodel,"qmodel=",qmodel, "extype=",extype
!endif
rowsikpq = AIJpqMatrixDimsList(NSOMOalpha,extype,pmodel,qmodel,1)
colsikpq = AIJpqMatrixDimsList(NSOMOalpha,extype,pmodel,qmodel,2)
!print *,"j=",j," Nsomo=",NSOMOalpha," rowsikpq=",rowsikpq," colsikpq=",colsikpq, " p=",pmodel," q=",qmodel, " extyp=",extype
@ -1692,6 +1852,10 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
rowsTKI = rowsikpq
enddo
!if(i.eq.1)then
! print *,"n_st=",n_st,"rowsTKI=",rowsTKI, " nconnectedI=",nconnectedI, &
! "totcolsTKI=",totcolsTKI
!endif
allocate(TKI(n_st,rowsTKI,totcolsTKI)) ! coefficients of CSF
! Initialize the integral container
! dims : (totcolsTKI, nconnectedI)
@ -1721,10 +1885,10 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
TKI(kk,l,totcolsTKI+m) = AIJpqContainer(l,m,pmodel,qmodel,extype,NSOMOalpha) &
* psi_in(kk,idxs_connectedI_alpha(j)+m-1)
enddo
!if(i.eq.1) then
! print *,AIJpqContainer(l,m,pmodel,qmodel,extype,NSOMOalpha)
!endif
enddo
!if(i.eq.1) then
! print *,"j=",j,"psi_in=",psi_in(1,idxs_connectedI_alpha(j)+m-1)
!endif
enddo
diagfactors_0 = diagfactors(j)*0.5d0
@ -1763,16 +1927,24 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
rowsTKI = rowsikpq
CCmattmp = 0.d0
!if(i.eq.1)then
! print *,"\t n_st=",n_st," colsikpq=",colsikpq," rowsTKI=",rowsTKI,&
! " | ",size(TKIGIJ,1),size(AIJpqContainer,1),size(CCmattmp,1)
!endif
call dgemm('N','N', n_st, colsikpq, rowsTKI, 1.d0, &
TKIGIJ(1,1,j), size(TKIGIJ,1), &
AIJpqContainer(1,1,pmodel,qmodel,extype,NSOMOalpha), &
size(AIJpqContainer,1), 0.d0, &
CCmattmp, size(CCmattmp,1) )
!print *,"j=",j,"colsikpq=",colsikpq, "sizeTIG=",size(TKIGIJ,1),"sizeaijpq=",size(AIJpqContainer,1)
do m = 1,colsikpq
call omp_set_lock(lock(idxs_connectedI_alpha(j)+m-1))
do kk = 1,n_st
psi_out(kk,idxs_connectedI_alpha(j)+m-1) += CCmattmp(kk,m)
!if(dabs(CCmattmp(kk,m)).gt.1e-10)then
! print *, CCmattmp(kk,m), " | ",idxs_connectedI_alpha(j)+m-1
!end if
enddo
call omp_unset_lock(lock(idxs_connectedI_alpha(j)+m-1))
enddo

35
src/davidson/diagonalization_hcfg.irp.f
View File

@ -112,6 +112,8 @@ subroutine davidson_diag_cfg_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
double precision, allocatable :: U(:,:), U_csf(:,:), overlap(:,:)
double precision, allocatable :: tmpU(:,:), tmpW(:,:)
double precision, pointer :: W(:,:), W_csf(:,:)
double precision, pointer :: W2(:,:), W_csf2(:,:)
double precision, allocatable :: U2(:,:), U_csf2(:,:)
logical :: disk_based
double precision :: energy_shift(N_st_diag_in*davidson_sze_max)
@ -234,12 +236,15 @@ subroutine davidson_diag_cfg_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
call c_f_pointer(ptr_w, W_csf, (/sze_csf,N_st_diag*itermax/))
else
allocate(W(sze,N_st_diag),W_csf(sze_csf,N_st_diag*itermax))
allocate(W2(sze,N_st_diag),W_csf2(sze_csf,N_st_diag*itermax))
endif
allocate( &
! Large
U(sze,N_st_diag), &
U2(sze,N_st_diag), &
U_csf(sze_csf,N_st_diag*itermax), &
U_csf2(sze_csf,N_st_diag*itermax), &
! Small
h(N_st_diag*itermax,N_st_diag*itermax), &
@ -324,8 +329,8 @@ subroutine davidson_diag_cfg_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
tmpU(kk,ii) = U_csf(ii,shift+kk)
enddo
enddo
!tmpU =0.0d0
!tmpU(1,2)=1.0d0
tmpU =0.0d0
tmpU(1,1)=1.0d0
double precision :: irp_rdtsc
double precision :: ticks_0, ticks_1
integer*8 :: irp_imax
@ -340,19 +345,19 @@ subroutine davidson_diag_cfg_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
enddo
enddo
!U_csf = 0.0d0
!U_csf(1,1) = 1.0d0
!u_in = 0.0d0
!call convertWFfromCSFtoDET(N_st_diag,tmpU,U2)
!call H_u_0_nstates_openmp(u_in,U2,N_st_diag,sze)
!call convertWFfromDETtoCSF(N_st_diag,u_in(1,1),W_csf2(1,1))
!do i=1,sze_csf
! print *,"I=",i," qp=",W_csf2(i,1)," my=",W_csf(i,1)," diff=",dabs(W_csf2(i,1))-dabs(W_csf(i,1))
! if(dabs(dabs(W_csf2(i,1))-dabs(W_csf(i,1))) .gt. 1.0e-10)then
! print *,"somo=",psi_configuration(1,1,i)," domo=",psi_configuration(1,2,i)," diff=",dabs(W_csf2(i,1))-dabs(W_csf(i,1))
! endif
!end do
!stop
U_csf = 0.0d0
U_csf(1,1) = 1.0d0
u_in = 0.0d0
call convertWFfromCSFtoDET(N_st_diag,tmpU,U2)
call H_u_0_nstates_openmp(u_in,U2,N_st_diag,sze)
call convertWFfromDETtoCSF(N_st_diag,u_in(1,1),W_csf2(1,1))
do i=1,sze_csf
print *,"I=",i," qp=",W_csf2(i,1)," my=",W_csf(i,1)," diff=",dabs(W_csf2(i,1))-dabs(W_csf(i,1))
!if(dabs(dabs(W_csf2(i,1))-dabs(W_csf(i,1))) .gt. 1.0e-10)then
! print *,"somo=",psi_configuration(1,1,i)," domo=",psi_configuration(1,2,i)," diff=",dabs(W_csf2(i,1))-dabs(W_csf(i,1))
!endif
end do
stop
deallocate(tmpW)
deallocate(tmpU)
endif