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

Fixed bug in calculation of Nsomo_J.

This commit is contained in:
v1j4y 2022-12-16 16:02:00 +01:00
parent bc69ac42b1
commit 2e0c96fd6b
3 changed files with 42 additions and 37 deletions

View File

@ -591,6 +591,7 @@ use bitmasks
! find out all pq holes possible
nholes = 0
listholes=-1
call bitstring_to_list(Isomo,listall,nelall,N_int)
do iii=1,nelall
@ -638,12 +639,19 @@ use bitmasks
tableUniqueAlphas = .FALSE.
! Now find the allowed (p,q) excitations
do i=1, N_int
Isomo(i) = iand(reunion_of_act_virt_bitmask(i,1),psi_configuration(i,1,idxI))
Idomo(i) = iand(reunion_of_act_virt_bitmask(i,2),psi_configuration(i,2,idxI))
Jsomo(i) = Isomo(i)
Jdomo(i) = Idomo(i)
do ii=1, N_int
!Isomo(ii) = iand(reunion_of_act_virt_bitmask(i,1),psi_configuration(ii,1,idxI))
!Idomo(ii) = iand(reunion_of_act_virt_bitmask(i,2),psi_configuration(ii,2,idxI))
Isomo(ii) = psi_configuration(ii,1,idxI)
Idomo(ii) = psi_configuration(ii,2,idxI)
Jsomo(ii) = Isomo(ii)
Jdomo(ii) = Idomo(ii)
enddo
!print *,"I=",idxI
!print *,"Isomo=",Isomo(1)!, Isomo(2)
!print *,"Idomo=",Idomo(1)!, Idomo(2)
!print *,listholes
!print *,listvmos
if(Nsomo_I .EQ. 0) then
kstart = 1
@ -789,6 +797,7 @@ use bitmasks
do j = 1, nvmos
qq = listvmos(j)
if(pp.eq.qq) cycle
jint = shiftr(qq-1,bit_kind_shift) + 1
jpos = qq-shiftl((jint-1),bit_kind_shift)-1
if(vmotype(j) == 1)then
@ -797,31 +806,25 @@ use bitmasks
Jdomo(jint) = IBSET(Jdomo(jint),jpos)
Jsomo(jint) = IBCLR(Jsomo(jint),jpos)
endif
if(pp .EQ. qq) then
if(vmotype(j) == 1)then
Jsomo(jint) = IBCLR(Jsomo(jint),jpos)
else if(vmotype(j) == 2)then
Jdomo(jint) = IBCLR(Jdomo(jint),jpos)
Jsomo(jint) = IBSET(Jsomo(jint),jpos)
endif
cycle
endif
if(tableUniqueAlphas(pp,qq)) then
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,listall,nelall,N_int)
!Nsomo_J = nelall
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
if(n_core_orb .le. 64)then
alphasIcfg_list(1,2,idxI,NalphaIcfg) = IOR(Jcfg(1,2),ISHFT(1_8,n_core_orb)-1)
else
n_core_orb_64 = n_core_orb

View File

@ -1892,8 +1892,10 @@ 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),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),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
@ -1913,9 +1915,9 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
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
if(i.eq.1)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

View File

@ -329,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,1)=1.0d0
tmpU =0.0d0
tmpU(1,1)=1.0d0
double precision :: irp_rdtsc
double precision :: ticks_0, ticks_1
integer*8 :: irp_imax
@ -345,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