LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-23 03:53:29 +01:00
Code Releases Activity

Merge branch 'dev' of github.com:QuantumPackage/qp2 into dev
Some checks failed
continuous-integration/drone/push Build is failing
Details

Browse Source
This commit is contained in:
Anthony Scemama 2022-12-23 18:31:55 +01:00
commit 9ae324001e
10 changed files with 1848 additions and 589 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
config/bull.cfg
View File

@ -6,7 +6,7 @@
# --align=32 : Align all provided arrays on a 32-byte boundary
#
[COMMON]
FC : mpiifort -shared-libgcc -shared-intel -fpic
FC : mpiifort -fpic -xCORE-AVX2
LAPACK_LIB : -mkl=parallel
IRPF90 : irpf90
IRPF90_FLAGS : --ninja --align=32 -DMPI
@ -31,7 +31,7 @@ OPENMP : 1 ; Append OpenMP flags
# -ftz : Flushes denormal results to zero
#
[OPT]
FCFLAGS : -xCORE-AVX2 -O2 -ip -ftz -g -traceback -qopt-prefetch=5 -qopt-prefetch-issue-excl-hint -unroll-aggressive
FCFLAGS : -O2 -ip -ftz -g -traceback -qopt-prefetch=5 -qopt-prefetch-issue-excl-hint -unroll-aggressive
# Profiling flags
#################
#

4
src/csf/cfgCI_utils.c
View File

@ -253,9 +253,9 @@ void generateAllBFs(int64_t Isomo, int64_t MS, Tree *bftree, int *NBF, int *NSOM
buildTreeDriver(bftree, *NSOMO, MS, NBF);
}
void ortho_qr_csf(double *overlapMatrix, int lda, double *orthoMatrix, int rows, int cols);
//void ortho_qr_csf(double *overlapMatrix, int lda, double *orthoMatrix, int rows, int cols);
// QR to orthogonalize CSFs does not work
//void gramSchmidt_qp(double *overlapMatrix, int rows, int cols, double *orthoMatrix){
// int i,j;
// //for(j=0;j<cols;++j){

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

File diff suppressed because it is too large Load Diff

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

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

@ -38,6 +38,7 @@ subroutine obtain_connected_J_givenI(idxI, givenI, connectedI, idxs_connectedI,
integer :: holetype(mo_num)
integer :: end_index
integer :: Nsomo_I
integer :: listall(N_int*bit_kind_size), nelall
!
! 2 2 1 1 0 0 : 1 1 0 0 0 0
@ -65,9 +66,12 @@ subroutine obtain_connected_J_givenI(idxI, givenI, connectedI, idxs_connectedI,
! Since CFGs are sorted wrt to seniority
! we don't have to search the full CFG list
Isomo = givenI(1,1)
Idomo = givenI(1,2)
Nsomo_I = POPCNT(Isomo)
Nsomo_I = 0
do i=1,N_int
Isomo = givenI(i,1)
Idomo = givenI(i,2)
Nsomo_I += POPCNT(Isomo)
end do
end_index = min(N_configuration,cfg_seniority_index(min(Nsomo_I+6,elec_num))-1)
if(end_index .LT. 0) end_index= N_configuration
!end_index = N_configuration
@ -83,17 +87,24 @@ subroutine obtain_connected_J_givenI(idxI, givenI, connectedI, idxs_connectedI,
! idxs_connectedI(nconnectedI)=i
! cycle
!endif
Isomo = givenI(1,1)
Idomo = givenI(1,2)
Jsomo = psi_configuration(1,1,i)
Jdomo = psi_configuration(1,2,i)
ndiffSOMO = 0
ndiffDOMO = 0
nxordiffSOMODOMO = 0
do ii=1,N_int
Isomo = givenI(ii,1)
Idomo = givenI(ii,2)
Jsomo = psi_configuration(ii,1,i)
Jdomo = psi_configuration(ii,2,i)
diffSOMO = IEOR(Isomo,Jsomo)
ndiffSOMO = POPCNT(diffSOMO)
ndiffSOMO += POPCNT(diffSOMO)
diffDOMO = IEOR(Idomo,Jdomo)
xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO)
ndiffDOMO = POPCNT(diffDOMO)
nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO)
nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO
ndiffDOMO += POPCNT(diffDOMO)
nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO)
nxordiffSOMODOMO += POPCNT(diffSOMO) + POPCNT(diffDOMO)
end do
if((nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then
!-------
! MONO |
@ -144,25 +155,45 @@ subroutine obtain_connected_J_givenI(idxI, givenI, connectedI, idxs_connectedI,
! find out all pq holes possible
nholes = 0
! holes in SOMO
Isomo = psi_configuration(1,1,i)
Idomo = psi_configuration(1,2,i)
do iii = 1,n_act_orb
ii = list_act(iii)
if(POPCNT(IAND(Isomo,IBSET(0_8,ii-1))) .EQ. 1) then
!Isomo = psi_configuration(1,1,i)
!Idomo = psi_configuration(1,2,i)
!do iii = 1,n_act_orb
! ii = list_act(iii)
! if(POPCNT(IAND(Isomo,IBSET(0_8,ii-1))) .EQ. 1) then
! nholes += 1
! listholes(nholes) = ii
! holetype(nholes) = 1
! endif
!end do
call bitstring_to_list(psi_configuration(1,1,i),listall,nelall,N_int)
do iii=1,nelall
nholes += 1
listholes(nholes) = ii
listholes(nholes) = listall(iii)
holetype(nholes) = 1
endif
end do
! holes in DOMO
do iii = 1,n_act_orb
ii = list_act(iii)
if(POPCNT(IAND(Idomo,IBSET(0_8,ii-1))) .EQ. 1) then
!do iii = 1,n_act_orb
! ii = list_act(iii)
! if(POPCNT(IAND(Idomo,IBSET(0_8,ii-1))) .EQ. 1) then
! nholes += 1
! listholes(nholes) = ii
! holetype(nholes) = 2
! endif
!end do
call bitstring_to_list(psi_configuration(1,2,i),listall,nelall,N_int)
do iii=1,nelall
if(listall(iii) .gt. n_core_orb)then
nholes += 1
listholes(nholes) = ii
listholes(nholes) = listall(iii)
holetype(nholes) = 2
endif
end do
ntotalconnectedI += max(1,(psi_config_data(i,2)-psi_config_data(i,1)+1)*nholes)
endif
end do
@ -199,6 +230,8 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI
integer*8 :: Isomo
integer*8 :: Jdomo
integer*8 :: Jsomo
integer(bit_kind) :: Jcfg(N_int,2)
integer(bit_kind) :: Icfg(N_int,2)
integer*8 :: IJsomo
integer*8 :: diffSOMO
integer*8 :: diffDOMO
@ -209,132 +242,261 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI
integer :: iii,ii,i,j,k,l,p,q,nsomoJ,nsomoalpha,starti,endi,extyp,nholes
integer :: listholes(mo_num)
integer :: holetype(mo_num)
integer :: end_index
integer :: Nsomo_alpha
integer :: end_index, ishift
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
! Since CFGs are sorted wrt to seniority
! we don't have to search the full CFG list
Isomo = Ialpha(1,1)
Idomo = Ialpha(1,2)
Nsomo_alpha = POPCNT(Isomo)
!Isomo = Ialpha(1,1)
!Idomo = Ialpha(1,2)
!Nsomo_alpha = POPCNT(Isomo)
Icfg = Ialpha
Nsomo_alpha = 0
!print *," Ialpha="
do ii=1,N_int
Isomo = Ialpha(ii,1)
Idomo = Ialpha(ii,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
if(end_index .LT. 0 .OR. end_index .lt. idxI) end_index= N_configuration
end_index = N_configuration
p = 0
q = 0
if (N_int > 1) stop 'obtain_connected_i_foralpha : N_int > 1'
!if (N_int > 1) stop 'obtain_connected_i_foralpha : N_int > 1'
do i=idxI,end_index
Isomo = Ialpha(1,1)
Idomo = Ialpha(1,2)
Jsomo = psi_configuration(1,1,i)
Jdomo = psi_configuration(1,2,i)
! Check for Minimal alpha electrons (MS)
if(POPCNT(Isomo).lt.MS)then
if(Nsomo_alpha .lt. MS)then
cycle
endif
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)
!if(idxI.eq.1)then
! print *," \t idxI=",i," diffS=",ndiffSOMO," popJs=", POPCNT(Jsomo)," popIs=",POPCNT(Isomo)
!endif
ndiffSOMO += POPCNT(diffSOMO)
diffDOMO = IEOR(Idomo,Jdomo)
xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO)
ndiffDOMO = POPCNT(diffDOMO)
nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO)
nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO
ndiffDOMO += POPCNT(diffDOMO)
nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO)
nxordiffSOMODOMO += POPCNT(diffSOMO) + POPCNT(diffDOMO)
end do
!if(idxI.eq.218)then
! print *,"I=",idxI,"Nsomo_alpha=",Nsomo_alpha,"nxordiffSOMODOMO(4)=",nxordiffSOMODOMO, " ndiffSOMO(2)=",ndiffSOMO, " ndiffDOMO=",ndiffDOMO
!endif
!Jcfg = psi_configuration(:,:,i)
!print *,"nxordiffSOMODOMO(4)=",nxordiffSOMODOMO, " ndiffSOMO(2)=",ndiffSOMO
if((nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then
select case(ndiffDOMO)
case (0)
! SOMO -> VMO
!print *,"obt SOMO -> VMO"
extyp = 3
!if(N_int .eq. 1) then
! IJsomo = IEOR(Isomo, Jsomo)
! p = TRAILZ(IAND(Isomo,IJsomo)) + 1
! IJsomo = IBCLR(IJsomo,p-1)
! q = TRAILZ(IJsomo) + 1
! !print *," p=",p," q=",q
! !call get_single_excitation_cfg(Jcfg, Icfg, p, q, N_int)
!else
! Find p
do ii=1,N_int
Isomo = Ialpha(ii,1)
Jsomo = psi_configuration(ii,1,i)
IJsomo = IEOR(Isomo, Jsomo)
!IRP_IF WITHOUT_TRAILZ
! p = (popcnt(ieor( IAND(Isomo,IJsomo) , IAND(Isomo,IJsomo) -1))-1) + 1
!IRP_ELSE
p = TRAILZ(IAND(Isomo,IJsomo)) + 1
!IRP_ENDIF
IJsomo = IBCLR(IJsomo,p-1)
!IRP_IF WITHOUT_TRAILZ
! q = (popcnt(ieor(IJsomo,IJsomo-1))-1) + 1
!IRP_ELSE
q = TRAILZ(IJsomo) + 1
!IRP_ENDIF
if(popcnt(IAND(Isomo,IJsomo)) > 0)then
p = TRAILZ(IAND(Isomo,IJsomo)) + 1 + (ii-1) * bit_kind_size
EXIT
endif
end do
! Find q
do ii=1,N_int
Isomo = Ialpha(ii,1)
Jsomo = psi_configuration(ii,1,i)
IJsomo = IEOR(Isomo, Jsomo)
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-1) * bit_kind_size
EXIT
endif
enddo
!endif
!assert ( p == pp)
!assert ( q == qq)
!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"
extyp = 2
!IRP_IF WITHOUT_TRAILZ
! p = (popcnt(ieor( IEOR(Idomo,Jdomo),IEOR(Idomo,Jdomo) -1))-1) + 1
!IRP_ELSE
p = TRAILZ(IEOR(Idomo,Jdomo)) + 1
!IRP_ENDIF
Isomo = IEOR(Isomo, Jsomo)
Isomo = IBCLR(Isomo,p-1)
!IRP_IF WITHOUT_TRAILZ
! q = (popcnt(ieor(Isomo,Isomo-1))-1) + 1
!IRP_ELSE
q = TRAILZ(Isomo) + 1
!IRP_ENDIF
!if(N_int.eq.1)then
! p = TRAILZ(IEOR(Idomo,Jdomo)) + 1
! Isomo = IEOR(Isomo, Jsomo)
! Isomo = IBCLR(Isomo,p-1)
! q = TRAILZ(Isomo) + 1
!else
! Find p
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
p = TRAILZ(IEOR(Idomo,Jdomo)) + 1 + (ii-1) * bit_kind_size
EXIT
endif
end do
! Find q
do ii=1,N_int
Isomo = Ialpha(ii,1)
Jsomo = psi_configuration(ii,1,i)
IJsomo = IEOR(Isomo, Jsomo)
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-1) * bit_kind_size
EXIT
endif
end do
!endif
!assert ( p == pp)
!assert ( q == qq)
else
! SOMO -> SOMO
!print *,"obt SOMO -> SOMO"
extyp = 1
!IRP_IF WITHOUT_TRAILZ
! q = (popcnt(ieor( IEOR(Idomo,Jdomo), IEOR(Idomo,Jdomo)-1))-1) + 1
!IRP_ELSE
q = TRAILZ(IEOR(Idomo,Jdomo)) + 1
!IRP_ENDIF
Isomo = IEOR(Isomo, Jsomo)
Isomo = IBCLR(Isomo,q-1)
!IRP_IF WITHOUT_TRAILZ
! p = (popcnt(ieor(Isomo,Isomo-1))-1) + 1
!IRP_ELSE
p = TRAILZ(Isomo) + 1
!IRP_ENDIF
! Check for Minimal alpha electrons (MS)
!if(POPCNT(Isomo).lt.MS)then
! cycle
!if(N_int.eq.1)then
! q = TRAILZ(IEOR(Idomo,Jdomo)) + 1
! Isomo = IEOR(Isomo, Jsomo)
! Isomo = IBCLR(Isomo,q-1)
! p = TRAILZ(Isomo) + 1
! ! Check for Minimal alpha electrons (MS)
! !if(POPCNT(Isomo).lt.MS)then
! ! cycle
! !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-1) * bit_kind_size
EXIT
endif
enddo
! Find q
do ii=1,N_int
Isomo = Ialpha(ii,1)
Jsomo = psi_configuration(ii,1,i)
IJsomo = IEOR(Isomo, Jsomo)
iint = shiftr(q-1,bit_kind_shift) + 1
ipos = q-shiftl((iint-1),bit_kind_shift)
if(iint .eq. ii)then
IJsomo = IBCLR(IJsomo,ipos-1)
endif
!print *,"ii=",ii," Isomo=",Isomo
if(popcnt(IJsomo) > 0)then
p = TRAILZ(IJsomo) + 1 + (ii-1) * bit_kind_size
EXIT
endif
enddo
!endif
end if
!assert ( p == pp)
!assert ( q == qq)
endif
!print *," 2--- p=",p," q=",q
case (2)
! DOMO -> SOMO
!print *,"obt DOMO -> SOMO"
extyp = 4
!if(N_int.eq.1)then
! IJsomo = IEOR(Isomo, Jsomo)
! p = TRAILZ(IAND(Jsomo,IJsomo)) + 1
! IJsomo = IBCLR(IJsomo,p-1)
! q = TRAILZ(IJsomo) + 1
!else
! Find p
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)
IJsomo = IEOR(Isomo, Jsomo)
!IRP_IF WITHOUT_TRAILZ
! p = (popcnt(ieor( IAND(Jsomo,IJsomo), IAND(Jsomo,IJsomo)-1))-1) + 1
!IRP_ELSE
p = TRAILZ(IAND(Jsomo,IJsomo)) + 1
!IRP_ENDIF
IJsomo = IBCLR(IJsomo,p-1)
!IRP_IF WITHOUT_TRAILZ
! q = (popcnt(ieor( IJsomo , IJsomo -1))-1) + 1
!IRP_ELSE
q = TRAILZ(IJsomo) + 1
!IRP_ENDIF
if(popcnt(IAND(Jsomo,IJsomo)) > 0)then
p = TRAILZ(IAND(Jsomo,IJsomo)) + 1 + (ii-1) * bit_kind_size
EXIT
endif
enddo
! Find q
do ii=1,N_int
Isomo = Ialpha(ii,1)
Jsomo = psi_configuration(ii,1,i)
IJsomo = IEOR(Isomo, Jsomo)
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-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)
connectedI(k,2,nconnectedI) = psi_configuration(k,2,i)
do ii=1,N_int
connectedI(ii,1,nconnectedI) = psi_configuration(ii,1,i)
connectedI(ii,2,nconnectedI) = psi_configuration(ii,2,i)
enddo
idxs_connectedI(nconnectedI)=starti
excitationIds(1,nconnectedI)=p
@ -343,28 +505,51 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI
diagfactors(nconnectedI) = 1.0d0
else if((ndiffSOMO + ndiffDOMO) .EQ. 0) then
! find out all pq holes possible
!print *,"I = ",i
!print *,"I somo= ",psi_configuration(1,1,i), " domo=", psi_configuration(1,2,i)
!print *,"alp somo= ",Ialpha(1,1), " domo=", Ialpha(1,2)
nholes = 0
! holes in SOMO
Isomo = psi_configuration(1,1,i)
Idomo = psi_configuration(1,2,i)
do iii = 1,n_act_orb
ii = list_act(iii)
if(POPCNT(IAND(Isomo,IBSET(0_8,ii-1))) .EQ. 1) then
!Isomo = psi_configuration(1,1,i)
!Idomo = psi_configuration(1,2,i)
!do iii = 1,n_act_orb
! ii = list_act(iii)
! if(POPCNT(IAND(Isomo,IBSET(0_8,ii-1))) .EQ. 1) then
! nholes += 1
! listholes(nholes) = ii
! holetype(nholes) = 1
! endif
!end do
call bitstring_to_list(psi_configuration(1,1,i),listall,nelall,N_int)
do iii=1,nelall
nholes += 1
listholes(nholes) = ii
listholes(nholes) = listall(iii)
holetype(nholes) = 1
endif
end do
! holes in DOMO
do iii = 1,n_act_orb
ii = list_act(iii)
if(POPCNT(IAND(Idomo,IBSET(0_8,ii-1))) .EQ. 1) then
!do iii = 1,n_act_orb
! ii = list_act(iii)
! if(POPCNT(IAND(Idomo,IBSET(0_8,ii-1))) .EQ. 1) then
! nholes += 1
! listholes(nholes) = ii
! holetype(nholes) = 2
! endif
!end do
nelall=0
listall=0
call bitstring_to_list(psi_configuration(1,2,i),listall,nelall,N_int)
do iii=1,nelall
if(listall(iii) .gt. n_core_orb)then
nholes += 1
listholes(nholes) = ii
listholes(nholes) = listall(iii)
holetype(nholes) = 2
endif
end do
do k=1,nholes
p = listholes(k)
q = p
@ -372,6 +557,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
@ -382,6 +568,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
@ -390,8 +577,10 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI
excitationTypes(nconnectedI) = extyp
diagfactors(nconnectedI) = 2.0d0
endif
!print *,excitationIds(1,nconnectedI), excitationIds(2,nconnectedI)
enddo
endif
end do
!print *,"nconnectedExtradiag=",nconnectedExtradiag," nconnectedDiad=",nconnectedDiag
end subroutine obtain_connected_I_foralpha

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

@ -146,7 +146,6 @@
ncfgprev = cfg_seniority_index(i+2)
end do
!print *," ^^^^^ N_CSF = ",n_CSF," N_CFG=",N_configuration
END_PROVIDER
@ -832,7 +831,7 @@ subroutine calculate_preconditioner_cfg(diag_energies)
! the configurations in psi_configuration
! returns : diag_energies :
END_DOC
integer :: i,j,k,kk,l,p,q,noccp,noccq, ii, jj
integer :: i,j,k,kk,l,p,q,noccp,noccq, ii, jj, iii
real*8,intent(out) :: diag_energies(n_CSF)
integer :: nholes
integer :: nvmos
@ -859,6 +858,7 @@ subroutine calculate_preconditioner_cfg(diag_energies)
real*8 :: hpp
real*8 :: meCC
real*8 :: core_act_contrib
integer :: listall(N_int*bit_kind_size), nelall
!PROVIDE h_core_ri
PROVIDE core_fock_operator
@ -875,11 +875,11 @@ subroutine calculate_preconditioner_cfg(diag_energies)
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))
!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)
@ -888,48 +888,63 @@ subroutine calculate_preconditioner_cfg(diag_energies)
! find out all pq holes possible
nholes = 0
listholes = -1
! holes in SOMO
!do k = 1,mo_num
do kk = 1,n_act_orb
k = list_act(kk)
if(POPCNT(IAND(Isomo,IBSET(0_8,k-1))) .EQ. 1) then
!do kk = 1,n_act_orb
! k = list_act(kk)
! if(POPCNT(IAND(Isomo,IBSET(0_8,k-1))) .EQ. 1) then
! nholes += 1
! listholes(nholes) = k
! holetype(nholes) = 1
! endif
!enddo
call bitstring_to_list(psi_configuration(1,1,i),listall,nelall,N_int)
do iii=1,nelall
nholes += 1
listholes(nholes) = k
listholes(nholes) = listall(iii)
holetype(nholes) = 1
endif
enddo
end do
! 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
do kk = 1,n_act_orb
k = list_act(kk)
if(POPCNT(IAND(Idomo,IBSET(0_8,k-1))) .EQ. 1) then
!do kk = 1,n_act_orb
! k = list_act(kk)
! if(POPCNT(IAND(Idomo,IBSET(0_8,k-1))) .EQ. 1) then
! nholes += 1
! listholes(nholes) = k
! holetype(nholes) = 2
! endif
!enddo
call bitstring_to_list(psi_configuration(1,2,i),listall,nelall,N_int)
do iii=1,nelall
if(listall(iii) .gt. n_core_orb)then
nholes += 1
listholes(nholes) = k
listholes(nholes) = listall(iii)
holetype(nholes) = 2
endif
enddo
end do
! 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
do kk = 1,n_act_orb
k = list_act(kk)
!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
!!! 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
!!do kk = 1,n_act_orb
!! k = list_act(kk)
!! !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)
@ -1219,27 +1234,30 @@ 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(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
@ -1253,40 +1271,155 @@ 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))
tmpp += POPCNT(Isomo(ii))
end do
mask = ISHFT(1_bit_kind,ipos+1) - 1
Isomotmp = IAND(Isomo(iint),mask)
!pmodel = tmpp + POPCNT(mask) - POPCNT(XOR(Isomotmp,mask))
pmodel = tmpp + POPCNT(Isomotmp)
!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))
tmpq += POPCNT(Isomo(ii))
end do
mask = ISHFT(1_bit_kind,ipos+1) - 1
Isomotmp = IAND(Isomo(iint),mask)
!qmodel = tmpq + POPCNT(mask) - POPCNT(XOR(Isomotmp,mask))
qmodel = tmpq + POPCNT(Isomotmp)
!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))
tmpp += POPCNT(Jsomo(ii))
end do
mask = ISHFT(1_bit_kind,ipos+1) - 1
Jsomotmp = IAND(Jsomo(iint),mask)
!pmodel = tmpp + POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask))
pmodel = tmpp + POPCNT(Jsomotmp)
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))
tmpq += POPCNT(Jsomo(ii))
end do
mask = ISHFT(1_bit_kind,ipos+1) - 1
Jsomotmp = IAND(Jsomo(iint),mask)
!qmodel = tmpq + POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask))
qmodel = tmpq + POPCNT(Jsomotmp)
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))
tmpp += POPCNT(Isomo(ii))
end do
mask = ISHFT(1_bit_kind,ipos+1) - 1
Isomotmp = IAND(Isomo(iint),mask)
!pmodel = tmpp + POPCNT(mask) - POPCNT(XOR(Isomotmp,mask))
pmodel = tmpp + POPCNT(Isomotmp)
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))
tmpq += POPCNT(Jsomo(ii))
end do
mask = ISHFT(1_bit_kind,ipos+1) - 1
Jsomotmp = IAND(Jsomo(iint),mask)
!qmodel = tmpq + POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) + 1
qmodel = tmpq + POPCNT(Jsomotmp) + 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))
tmpp += POPCNT(Isomo(ii))
end do
mask = ISHFT(1_bit_kind,ipos+1) - 1
Isomotmp = IAND(Isomo(iint),mask)
!pmodel = tmpp + POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) + 1
pmodel = tmpp + POPCNT(Isomotmp) + 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))
tmpq += POPCNT(Jsomo(ii))
end do
mask = ISHFT(1_bit_kind,ipos+1) - 1
Jsomotmp = IAND(Jsomo(iint),mask)
!qmodel = tmpq + POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask))
qmodel = tmpq + POPCNT(Jsomotmp)
endif
case (4)
! DOMO -> SOMO
@ -1294,19 +1427,75 @@ 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))
tmpp += POPCNT(Jsomo(ii))
end do
mask = ISHFT(1_bit_kind,ipos+1) - 1
Jsomotmp = IAND(Jsomo(iint),mask)
!pmodel = tmpp + POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask))
pmodel = tmpp + POPCNT(Jsomotmp)
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))
tmpq += POPCNT(Isomo(ii))
end do
mask = ISHFT(1_bit_kind,ipos+1) - 1
Isomotmp = IAND(Isomo(iint),mask)
!qmodel = tmpq + POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) + 1
qmodel = tmpq + POPCNT(Isomotmp) + 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))
tmpp += POPCNT(Jsomo(ii))
end do
mask = ISHFT(1_bit_kind,ipos+1) - 1
Jsomotmp = IAND(Jsomo(iint),mask)
!pmodel = tmpp + POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) + 1
pmodel = tmpp + POPCNT(Jsomotmp) + 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))
tmpq += POPCNT(Isomo(ii))
end do
mask = ISHFT(1_bit_kind,ipos+1) - 1
Isomotmp = IAND(Isomo(iint),mask)
!qmodel = tmpq + POPCNT(mask) - POPCNT(XOR(Isomotmp,mask))
qmodel = tmpq + POPCNT(Isomotmp)
endif
case default
print *,"something is wrong in convertOrbIdsToModelSpaceIds"
@ -1364,8 +1553,13 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
integer :: rowsTKI
integer :: noccpp
integer :: istart_cfg, iend_cfg, num_threads_max
integer :: iint, jint, ipos, jpos, Nsomo_I, iii
integer :: nconnectedJ,nconnectedtotalmax,nconnectedmaxJ,maxnalphas,ntotJ
integer*8 :: MS, Isomo, Idomo, Jsomo, Jdomo, Ialpha, Ibeta
integer*8 :: MS,Ialpha, Ibeta
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 :: moi, moj, mok, mol, starti, endi, startj, endj, cnti, cntj, cntk
real*8 :: norm_coef_cfg, fac2eints
real*8 :: norm_coef_det
@ -1380,6 +1574,8 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
real*8,dimension(:),allocatable:: diag_energies
real*8 :: tmpvar, tmptot
real*8 :: core_act_contrib
integer :: listall(N_int*bit_kind_size), nelall
integer :: countelec
integer(omp_lock_kind), allocatable :: lock(:)
call omp_set_max_active_levels(1)
@ -1408,8 +1604,8 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
!nconnectedtotalmax = 1000
!nconnectedmaxJ = 1000
maxnalphas = elec_num*mo_num
Icfg(1,1) = psi_configuration(1,1,1)
Icfg(1,2) = psi_configuration(1,2,1)
Icfg(:,1) = psi_configuration(:,1,1)
Icfg(:,2) = psi_configuration(:,2,1)
allocate(listconnectedJ(N_INT,2,max(sze,10000)))
allocate(idslistconnectedJ(max(sze,10000)))
call obtain_connected_J_givenI(1, Icfg, listconnectedJ, idslistconnectedJ, nconnectedmaxJ, nconnectedtotalmax)
@ -1441,6 +1637,7 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
!$OMP shared(istart_cfg, iend_cfg, psi_configuration, mo_num, psi_config_data,&
!$OMP N_int, N_st, psi_out, psi_in, h_core_ri, core_energy, h_act_ri, AIJpqContainer,&
!$OMP pp, sze, NalphaIcfg_list,alphasIcfg_list, bit_tmp, &
!$OMP qq, iint, jint, ipos, jpos, nelall, listall, Nsomo_I, countelec,&
!$OMP AIJpqMatrixDimsList, diag_energies, n_CSF, lock, NBFmax,nconnectedtotalmax, nconnectedmaxJ,maxnalphas,&
!$OMP n_core_orb, n_act_orb, list_act, n, list_core, list_core_is_built,core_act_contrib, num_threads_max,&
!$OMP n_core_orb_is_built, mo_integrals_map, mo_integrals_map_is_built)
@ -1463,10 +1660,12 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
! else
! cycle
Icfg(1,1) = psi_configuration(1,1,i)
Icfg(1,2) = psi_configuration(1,2,i)
Isomo = Icfg(1,1)
Idomo = Icfg(1,2)
do ii=1,N_INT
Icfg(ii,1) = psi_configuration(ii,1,i)
Icfg(ii,2) = psi_configuration(ii,2,i)
Isomo(ii) = Icfg(ii,1)
Idomo(ii) = Icfg(ii,2)
enddo
NSOMOI = getNSOMO(Icfg)
! find out all pq holes possible
@ -1477,42 +1676,86 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
! list_core_inact
! bitmasks
!do k = 1,mo_num
do kk = 1,n_act_orb
k = list_act(kk)
if(POPCNT(IAND(Isomo,IBSET(0_8,k-1))) .EQ. 1) then
! do kk = 1,n_act_orb
! k = list_act(kk)
! 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 = 1,mo_num
! do kk = 1,n_act_orb
! k = list_act(kk)
! if(POPCNT(IAND(Idomo,IBSET(0_8,k-1))) .EQ. 1) then
! nholes += 1
! listholes(nholes) = k
! holetype(nholes) = 2
! endif
! enddo
! ! find vmos
! do kk = 1,n_act_orb
! k = list_act(kk)
! !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
! find out all pq holes possible
nholes = 0
call bitstring_to_list(Isomo,listall,nelall,N_int)
do iii=1,nelall
nholes += 1
listholes(nholes) = k
listholes(nholes) = listall(iii)
holetype(nholes) = 1
endif
enddo
! holes in DOMO
!do k = 1,mo_num
do kk = 1,n_act_orb
k = list_act(kk)
if(POPCNT(IAND(Idomo,IBSET(0_8,k-1))) .EQ. 1) then
end do
Nsomo_I = nelall
call bitstring_to_list(Idomo,listall,nelall,N_int)
do iii=1,nelall
if(listall(iii) .gt. n_core_orb)then
nholes += 1
listholes(nholes) = k
listholes(nholes) = listall(iii)
holetype(nholes) = 2
endif
enddo
end do
! find vmos
listvmos = -1
vmotype = -1
nvmos = 0
do kk = 1,n_act_orb
k = list_act(kk)
!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
! find vmos
! Take into account N_int
do ii = 1, n_act_orb
iii = list_act(ii)
iint = shiftr(iii-1,bit_kind_shift) + 1
ipos = iii-shiftl((iint-1),bit_kind_shift)-1
if(IAND(Idomo(iint),(IBSET(0_8,ipos))) .EQ. 0) then
if(IAND(Isomo(iint),(IBSET(0_8,ipos))) .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
listvmos(nvmos) = iii
vmotype(nvmos) = 1
else if(POPCNT(IAND(Isomo(iint),(IBSET(0_8,ipos)))) .EQ. 1) then
nvmos += 1
listvmos(nvmos) = iii
vmotype(nvmos) = 2
end if
enddo
end if
end do
! Icsf ids
@ -1531,16 +1774,31 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
extype = excitationTypes_single(j)
! Off diagonal terms
call convertOrbIdsToModelSpaceIds(Icfg, singlesI(1,1,j), p, q, extype, pmodel, qmodel)
Jsomo = singlesI(1,1,j)
Jdomo = singlesI(1,2,j)
do ii=1,N_INT
Jsomo(ii) = singlesI(1,1,j)
Jdomo(ii) = singlesI(1,2,j)
enddo
! Get actual p pos
pp = p
iint = shiftr(pp-1,bit_kind_shift) + 1
ipos = pp-shiftl((iint-1),bit_kind_shift)-1
! Get actual q pos
qq = q
jint = shiftr(qq-1,bit_kind_shift) + 1
jpos = qq-shiftl((jint-1),bit_kind_shift)-1
! 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
!if(POPCNT(IAND(Jsomo,IBSET(0_8,q-1))) .EQ. 1 .AND. POPCNT(IAND(Isomo,IBSET(0_8,q-1))) .EQ. 0) then
if(POPCNT(IAND(Jsomo(jint),IBSET(0_8,jpos))) .EQ. 1 .AND. POPCNT(IAND(Isomo(jint),IBSET(0_8,jpos))) .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
!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
if((POPCNT(IAND(Jdomo(jint),IBSET(0_8,jpos))) .EQ. 1 .AND. POPCNT(IAND(Idomo(jint),IBSET(0_8,jpos))) .EQ. 0) .AND.&
POPCNT(IAND(Isomo(jint),IBSET(0_8,jpos))) .EQ. 0) then
nholes += 1
listholes(nholes) = q
holetype(nholes) = 2
@ -1576,10 +1834,12 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
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
!if(POPCNT(IAND(Jsomo,IBSET(0_8,q-1))) .EQ. 1 .AND. POPCNT(IAND(Isomo,IBSET(0_8,q-1))) .EQ. 0) then
if(POPCNT(IAND(Jsomo(jint),IBSET(0_8,jpos))) .EQ. 1 .AND. POPCNT(IAND(Isomo(jint),IBSET(0_8,jpos))) .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
if((POPCNT(IAND(Jdomo(jint),IBSET(0_8,jpos))) .EQ. 1 .AND. POPCNT(IAND(Idomo(jint),IBSET(0_8,jpos))) .EQ. 0) .AND.&
POPCNT(IAND(Isomo(jint),IBSET(0_8,jpos))) .EQ. 0) then
nholes -= 1
endif
enddo
@ -1591,6 +1851,9 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
deallocate(excitationTypes_single)
!print *," singles part psi(1,1)=",psi_out(1,1)
!do i=1,n_CSF
! print *,"i=",i," psi(i)=",psi_out(1,i)
!enddo
allocate(listconnectedJ(N_INT,2,max(sze,10000)))
allocate(alphas_Icfg(N_INT,2,max(sze,10000)))
@ -1605,7 +1868,6 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
!!!====================!!!
!!! Double Excitations !!!
!!!====================!!!
! Loop over all selected configurations
!$OMP DO SCHEDULE(static)
do i = istart_cfg,iend_cfg
@ -1615,8 +1877,10 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
! else
! cycle
Icfg(1,1) = psi_configuration(1,1,i)
Icfg(1,2) = psi_configuration(1,2,i)
do ii=1,N_INT
Icfg(ii,1) = psi_configuration(ii,1,i)
Icfg(ii,2) = psi_configuration(ii,2,i)
enddo
starti = psi_config_data(i,1)
endi = psi_config_data(i,2)
@ -1627,14 +1891,15 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
Nalphas_Icfg = NalphaIcfg_list(i)
alphas_Icfg(1:n_int,1:2,1:Nalphas_Icfg) = alphasIcfg_list(1:n_int,1:2,i,1:Nalphas_Icfg)
if(Nalphas_Icfg .GT. maxnalphas) then
print *,"Nalpha > maxnalpha"
endif
!if(Nalphas_Icfg .GT. maxnalphas) then
! print *,"Nalpha > maxnalpha"
!endif
call obtain_connected_J_givenI(i, Icfg, listconnectedJ, idslistconnectedJ, nconnectedJ, ntotJ)
!call obtain_connected_J_givenI(i, Icfg, listconnectedJ, idslistconnectedJ, nconnectedJ, ntotJ)
! TODO : remove doubly excited for return
!print *,"I=",i," isomo=",psi_configuration(1,1,i)," idomo=",psi_configuration(1,2,i), " psiout=",psi_out(1,5)
!print *,"I=",i,"isomo=",psi_configuration(1,1,i),psi_configuration(2,1,i),POPCNT(psi_configuration(1,1,i)),POPCNT(psi_configuration(2,1,i)),&
!"idomo=",psi_configuration(1,2,i),psi_configuration(2,2,i),POPCNT(psi_configuration(1,2,i)),POPCNT(psi_configuration(2,2,i)), "Nalphas_Icfg=",Nalphas_Icfg
do k = 1,Nalphas_Icfg
! Now generate all singly excited with respect to a given alpha CFG
@ -1645,15 +1910,18 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
call obtain_connected_I_foralpha(i, alphas_Icfg(1,1,k), connectedI_alpha, idxs_connectedI_alpha, &
nconnectedI, excitationIds, excitationTypes, diagfactors)
!if(i .EQ. 218) then
! 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
! !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)), " NconnectedI=",nconnectedI
!endif
if(nconnectedI .EQ. 0) then
cycle
endif
!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))
!endif
! Here we do 2x the loop. One to count for the size of the matrix, then we compute.
totcolsTKI = 0
rowsTKI = -1
@ -1663,15 +1931,30 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
p = excitationIds(1,j)
q = excitationIds(2,j)
extype = excitationTypes(j)
!print *,"K=",k,"j=",j, "countelec=",countelec," p=",p," q=",q, " extype=",extype, "NSOMOalpha=",NSOMOalpha," NSOMOI=",NSOMOI, "alphas_Icfg(1,1,k)=",alphas_Icfg(1,1,k), &
!alphas_Icfg(2,1,k), " domo=",alphas_Icfg(1,2,k), alphas_Icfg(2,2,k), " connected somo=",connectedI_alpha(1,1,j), &
!connectedI_alpha(2,1,j), " domo=",connectedI_alpha(1,2,j), connectedI_alpha(2,2,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
rowsikpq = AIJpqMatrixDimsList(NSOMOalpha,extype,pmodel,qmodel,1)
colsikpq = AIJpqMatrixDimsList(NSOMOalpha,extype,pmodel,qmodel,2)
!if(i.eq.218)then
! print *,"j=",j," k=",k,"p=",p,"q=",q,"NSOMOalpha=",NSOMOalpha, "pmodel=",pmodel,"qmodel=",qmodel, "extype=",extype,&
! "conn somo=",connectedI_alpha(1,1,j),connectedI_alpha(2,1,j),&
! "conn domo=",connectedI_alpha(1,2,j),connectedI_alpha(2,2,j)
! do m=1,colsikpq
! print *,idxs_connectedI_alpha(j)+m-1
! enddo
!endif
!print *,"j=",j," Nsomo=",NSOMOalpha," rowsikpq=",rowsikpq," colsikpq=",colsikpq, " p=",pmodel," q=",qmodel, " extyp=",extype
totcolsTKI += colsikpq
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)
@ -1701,10 +1984,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
@ -1743,16 +2026,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
@ -1787,6 +2078,10 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
!$OMP END DO
!$OMP END PARALLEL
!print *," ----- "
!do i=1,sze
! print *,"i=",i," psi_out(i)=",psi_out(1,i)
!end do
call omp_set_max_active_levels(4)
deallocate(diag_energies)

0
src/dav_general_mat/test_dav.irp.f → src/dav_general_mat/test_dav.irp.f.example
View File

13
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), &
@ -325,7 +330,7 @@ subroutine davidson_diag_cfg_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
enddo
enddo
!tmpU =0.0d0
!tmpU(1,2)=1.0d0
!tmpU(1,1)=1.0d0
double precision :: irp_rdtsc
double precision :: ticks_0, ticks_1
integer*8 :: irp_imax
@ -348,9 +353,9 @@ subroutine davidson_diag_cfg_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
!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
! !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)

94
src/determinants/slater_rules.irp.f
View File

@ -83,7 +83,7 @@ subroutine get_excitation(det1,det2,exc,degree,phase,Nint)
! exc(1,1,1) = q
! exc(1,2,1) = p
! T^alpha_pq : exc(0,1,2) = 1
! T^beta_pq : exc(0,1,2) = 1
! exc(0,2,2) = 1
! exc(1,1,2) = q
! exc(1,2,2) = p
@ -434,6 +434,98 @@ subroutine get_single_excitation(det1,det2,exc,phase,Nint)
end
subroutine get_single_excitation_cfg(cfg1,cfg2,p,q,Nint)
use bitmasks
implicit none
BEGIN_DOC
! Returns the excitation operator between two singly excited configurations.
END_DOC
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: cfg1(Nint,2)
integer(bit_kind), intent(in) :: cfg2(Nint,2)
integer, intent(out) :: p, q
integer :: tz
integer :: l, ispin, idx_hole, idx_particle, ishift
integer :: nperm
integer :: i,j,k,m,n
integer :: high, low
integer :: a,b,c,d
integer(bit_kind) :: hole, particle, tmp
integer :: exc(0:2,2,2)
ASSERT (Nint > 0)
nperm = 0
p = 0
q = 0
exc(0,1,1) = 0
exc(0,2,1) = 0
exc(0,1,2) = 0
exc(0,2,2) = 0
do ispin = 1,2
ishift = 1-bit_kind_size
do l=1,Nint
ishift = ishift + bit_kind_size
if (cfg1(l,ispin) == cfg2(l,ispin)) then
cycle
endif
tmp = xor( cfg1(l,ispin), cfg2(l,ispin) )
particle = iand(tmp, cfg2(l,ispin))
hole = iand(tmp, cfg1(l,ispin))
if (particle /= 0_bit_kind) then
tz = trailz(particle)
exc(0,2,ispin) = 1
exc(1,2,ispin) = tz+ishift
!print *,"part ",tz+ishift, " ispin=",ispin
endif
if (hole /= 0_bit_kind) then
tz = trailz(hole)
exc(0,1,ispin) = 1
exc(1,1,ispin) = tz+ishift
!print *,"hole ",tz+ishift, " ispin=",ispin
endif
if ( iand(exc(0,1,ispin),exc(0,2,ispin)) /= 1) then ! exc(0,1,ispin)/=1 and exc(0,2,ispin) /= 1
cycle
endif
high = max(exc(1,1,ispin), exc(1,2,ispin))-1
low = min(exc(1,1,ispin), exc(1,2,ispin))
ASSERT (low >= 0)
ASSERT (high > 0)
k = shiftr(high,bit_kind_shift)+1
j = shiftr(low,bit_kind_shift)+1
m = iand(high,bit_kind_size-1)
n = iand(low,bit_kind_size-1)
if (j==k) then
nperm = nperm + popcnt(iand(cfg1(j,ispin), &
iand( shiftl(1_bit_kind,m)-1_bit_kind, &
not(shiftl(1_bit_kind,n))+1_bit_kind)) )
else
nperm = nperm + popcnt( &
iand(cfg1(j,ispin), &
iand(not(0_bit_kind), &
(not(shiftl(1_bit_kind,n)) + 1_bit_kind) ))) &
+ popcnt(iand(cfg1(k,ispin), &
(shiftl(1_bit_kind,m) - 1_bit_kind ) ))
do i=j+1,k-1
nperm = nperm + popcnt(cfg1(i,ispin))
end do
endif
! Set p and q
q = max(exc(1,1,1),exc(1,1,2))
p = max(exc(1,2,1),exc(1,2,2))
return
enddo
enddo
end
subroutine bitstring_to_list_ab( string, list, n_elements, Nint)
use bitmasks
implicit none

4
src/utils/linear_algebra.irp.f
View File

@ -1136,7 +1136,6 @@ subroutine ortho_svd(A,LDA,m,n)
end
! QR to orthonormalize CSFs does not work :-(
!subroutine ortho_qr_withB(A,LDA,B,m,n)
! implicit none
! BEGIN_DOC
@ -1223,7 +1222,7 @@ end
!
! !deallocate(WORK,TAU)
!end
!
!subroutine ortho_qr_csf(A, LDA, B, m, n) bind(C, name="ortho_qr_csf")
! use iso_c_binding
! integer(c_int32_t), value :: LDA
@ -1234,6 +1233,7 @@ end
! call ortho_qr_withB(A,LDA,B,m,n)
!end subroutine ortho_qr_csf
subroutine ortho_qr(A,LDA,m,n)
implicit none
BEGIN_DOC