9
1
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-11-07 22:13:38 +01:00

Removed deletion of domos in get_phase.

This commit is contained in:
v1j4y 2021-03-01 12:41:42 +01:00
parent 4f5f612c5d
commit 40a3a30ea1
2 changed files with 302 additions and 18 deletions

View File

@ -1,3 +1,273 @@
subroutine obtain_associated_alphaI(idxI, Icfg, alphasIcfg, NalphaIcfg, factor_alphaI)
implicit none
use bitmasks
BEGIN_DOC
! Documentation for alphasI
! Returns the associated alpha's for
! the input configuration Icfg.
END_DOC
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)
integer :: holetype(mo_num) ! 1-> SOMO 2->DOMO
integer :: nholes
integer :: nvmos
integer :: listvmos(mo_num)
integer :: vmotype(mo_num) ! 1 -> VMO 2 -> SOMO
integer*8 :: Idomo
integer*8 :: Isomo
integer*8 :: Jdomo
integer*8 :: Jsomo
integer*8 :: diffSOMO
integer*8 :: diffDOMO
integer*8 :: xordiffSOMODOMO
integer :: ndiffSOMO
integer :: ndiffDOMO
integer :: nxordiffSOMODOMO
integer :: ndiffAll
integer :: i
integer :: j
integer :: k
integer :: hole
integer :: p
integer :: q
integer :: countalphas
logical :: pqAlreadyGenQ
logical :: pqExistsQ
logical :: ppExistsQ
Isomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,1))
Idomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,2))
!print*,"Input cfg"
!call debug_spindet(Isomo,1)
!call debug_spindet(Idomo,1)
!print*,n_act_orb, "monum=",mo_num," n_core=",n_core_orb
! find out all pq holes possible
nholes = 0
! holes in SOMO
do i = 1,mo_num
if(POPCNT(IAND(Isomo,IBSET(0_8,i-1))) .EQ. 1) then
nholes += 1
listholes(nholes) = i
holetype(nholes) = 1
endif
end do
! holes in DOMO
do i = 1,mo_num
if(POPCNT(IAND(Idomo,IBSET(0_8,i-1))) .EQ. 1) then
nholes += 1
listholes(nholes) = i
holetype(nholes) = 2
endif
end do
! find vmos
listvmos = -1
vmotype = -1
nvmos = 0
do i = 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,i-1)))) .EQ. 0 .AND. POPCNT(IAND(Idomo,(IBSET(0_8,i-1)))) .EQ. 0) then
nvmos += 1
listvmos(nvmos) = i
vmotype(nvmos) = 1
else if(POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))) .EQ. 1 .AND. POPCNT(IAND(Idomo,(IBSET(0_8,i-1)))) .EQ. 0 ) then
nvmos += 1
listvmos(nvmos) = i
vmotype(nvmos) = 2
end if
end do
!print *,"Nvmo=",nvmos
!print *,listvmos
!print *,vmotype
allocate(tableUniqueAlphas(mo_num,mo_num))
tableUniqueAlphas = .FALSE.
! Now find the allowed (p,q) excitations
Isomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,1))
Idomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,2))
!print *,"Isomo"
!call debug_spindet(Isomo,1)
!call debug_spindet(Idomo,1)
!print *,"Nholes=",nholes," Nvmos=",nvmos, " idxi=",idxI
!do i = 1,nholes
! print *,i,"->",listholes(i)
!enddo
!do i = 1,nvmos
! print *,i,"->",listvmos(i)
!enddo
! TODO cfg_seniority_index
do i = 1,nholes
p = listholes(i)
do j = 1,nvmos
q = listvmos(j)
if(p .EQ. q) cycle
if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 1) then
! SOMO -> VMO
Jsomo = IBCLR(Isomo,p-1)
Jsomo = IBSET(Jsomo,q-1)
Jdomo = Idomo
else if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 2) then
! SOMO -> SOMO
Jsomo = IBCLR(Isomo,p-1)
Jsomo = IBCLR(Jsomo,q-1)
Jdomo = IBSET(Idomo,q-1)
else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then
! DOMO -> VMO
Jsomo = IBSET(Isomo,p-1)
Jsomo = IBSET(Jsomo,q-1)
Jdomo = IBCLR(Idomo,p-1)
else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 2) then
! DOMO -> SOMO
Jsomo = IBSET(Isomo,p-1)
Jsomo = IBCLR(Jsomo,q-1)
Jdomo = IBCLR(Idomo,p-1)
Jdomo = IBSET(Jdomo,q-1)
else
print*,"Something went wrong in obtain_associated_alphaI"
endif
pqAlreadyGenQ = .FALSE.
! First check if it can be generated before
do k = 1, idxI-1
diffSOMO = IEOR(Jsomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k)))
diffDOMO = IEOR(Jdomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k)))
xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO)
ndiffSOMO = POPCNT(diffSOMO)
ndiffDOMO = POPCNT(diffDOMO)
nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO)
!if(POPCNT(IEOR(diffSOMO,diffDOMO)) .LE. 1 .AND. ndiffDOMO .LT. 3) then
if((ndiffSOMO+ndiffDOMO) .EQ. 0) then
pqAlreadyGenQ = .TRUE.
ppExistsQ = .TRUE.
EXIT
endif
if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then
pqAlreadyGenQ = .TRUE.
!ppExistsQ = .TRUE.
!print *,i,k,ndiffSOMO,ndiffDOMO
!call debug_spindet(Jsomo,1)
!call debug_spindet(Jdomo,1)
!call debug_spindet(iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k)),1)
!call debug_spindet(iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k)),1)
EXIT
endif
end do
!print *,"(,",p,",",q,")",pqAlreadyGenQ
if(pqAlreadyGenQ) cycle
pqExistsQ = .FALSE.
! now check if this exists in the selected list
!do k = idxI+1, N_configuration
! diffSOMO = IEOR(OR(reunion_of_act_virt_bitmask(1,1),Jsomo),psi_configuration(1,1,k))
! diffDOMO = IEOR(OR(reunion_of_act_virt_bitmask(1,1),Jdomo),psi_configuration(1,2,k))
! ndiffSOMO = POPCNT(diffSOMO)
! ndiffDOMO = POPCNT(diffDOMO)
! if((ndiffSOMO + ndiffDOMO) .EQ. 0) then
! pqExistsQ = .TRUE.
! EXIT
! endif
!end do
if(.NOT. pqExistsQ) then
tableUniqueAlphas(p,q) = .TRUE.
!print *,p,q
!call debug_spindet(Jsomo,1)
!call debug_spindet(Jdomo,1)
endif
end do
end do
!print *,tableUniqueAlphas(:,:)
! prune list of alphas
Isomo = Icfg(1,1)
Idomo = Icfg(1,2)
Jsomo = Icfg(1,1)
Jdomo = Icfg(1,2)
NalphaIcfg = 0
do i = 1, nholes
p = listholes(i)
do j = 1, nvmos
q = listvmos(j)
if(p .EQ. q) cycle
if(tableUniqueAlphas(p,q)) then
if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 1) then
! SOMO -> VMO
Jsomo = IBCLR(Isomo,p-1)
Jsomo = IBSET(Jsomo,q-1)
Jdomo = Idomo
else if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 2) then
! SOMO -> SOMO
Jsomo = IBCLR(Isomo,p-1)
Jsomo = IBCLR(Jsomo,q-1)
Jdomo = IBSET(Idomo,q-1)
else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then
! DOMO -> VMO
Jsomo = IBSET(Isomo,p-1)
Jsomo = IBSET(Jsomo,q-1)
Jdomo = IBCLR(Idomo,p-1)
else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 2) then
! DOMO -> SOMO
Jsomo = IBSET(Isomo,p-1)
Jsomo = IBCLR(Jsomo,q-1)
Jdomo = IBCLR(Idomo,p-1)
Jdomo = IBSET(Jdomo,q-1)
else
print*,"Something went wrong in obtain_associated_alphaI"
endif
! SOMO
NalphaIcfg += 1
!print *,p,q,"|",holetype(i),vmotype(j),NalphaIcfg
alphasIcfg(1,1,NalphaIcfg) = Jsomo
alphasIcfg(1,2,NalphaIcfg) = IOR(Jdomo,ISHFT(1_8,n_core_orb)-1)
endif
end do
end do
! Check if this Icfg has been previously generated as a mono
ppExistsQ = .False.
Isomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,1))
Idomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,2))
do k = 1, idxI-1
diffSOMO = IEOR(Isomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k)))
diffDOMO = IEOR(Idomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k)))
xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO)
ndiffSOMO = POPCNT(diffSOMO)
ndiffDOMO = POPCNT(diffDOMO)
nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO)
if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then
ppExistsQ = .TRUE.
EXIT
endif
end do
! Diagonal part (pp,qq)
if(nholes > 0 .AND. (.NOT. ppExistsQ))then
! SOMO
NalphaIcfg += 1
!print *,p,q,"|",holetype(i),vmotype(j),NalphaIcfg
!call debug_spindet(Idomo,1)
!call debug_spindet(Jdomo,1)
alphasIcfg(1,1,NalphaIcfg) = Icfg(1,1)
alphasIcfg(1,2,NalphaIcfg) = Icfg(1,2)
endif
end subroutine
function getNSOMO(Icfg) result(NSOMO) function getNSOMO(Icfg) result(NSOMO)
implicit none implicit none
integer(bit_kind),intent(in) :: Icfg(N_int,2) integer(bit_kind),intent(in) :: Icfg(N_int,2)
@ -47,6 +317,7 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod
qmodel = q qmodel = q
if(p .EQ. q) then if(p .EQ. q) then
!print *,"input pq=",p,q,"extype=",extype
pmodel = 1 pmodel = 1
qmodel = 1 qmodel = 1
else else
@ -80,26 +351,44 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod
! SOMO -> VMO ! SOMO -> VMO
!print *,"type -> SOMO -> VMO" !print *,"type -> SOMO -> VMO"
!Isomo = IEOR(Isomo,Jsomo) !Isomo = IEOR(Isomo,Jsomo)
mask = ISHFT(1_8,p) - 1 if(p.LT.q) then
Isomo = IAND(Isomo,mask) mask = ISHFT(1_8,p) - 1
pmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) Isomo = IAND(Isomo,mask)
mask = ISHFT(1_8,q) - 1 pmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask))
Jsomo = IAND(Jsomo,mask) mask = ISHFT(1_8,q) - 1
qmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) Jsomo = IAND(Jsomo,mask)
qmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,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))
endif
case (4) case (4)
! DOMO -> SOMO ! DOMO -> SOMO
! remove all domos except one at p ! remove all domos except one at p
!print *,"type -> DOMO -> SOMO" !print *,"type -> DOMO -> SOMO"
!Isomo = IEOR(Isomo,Jsomo) !Isomo = IEOR(Isomo,Jsomo)
mask = ISHFT(1_8,p) - 1 if(p.LT.q) then
Jsomo = IAND(Jsomo,mask) mask = ISHFT(1_8,p) - 1
pmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) Jsomo = IAND(Jsomo,mask)
mask = ISHFT(1_8,q) - 1 pmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask))
Isomo = IAND(Isomo,mask) mask = ISHFT(1_8,q) - 1
qmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) Isomo = IAND(Isomo,mask)
qmodel = POPCNT(mask) - POPCNT(XOR(Isomo,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))
endif
case default case default
print *,"something is wrong in convertOrbIdsToModelSpaceIds" print *,"something is wrong in convertOrbIdsToModelSpaceIds"
end select end select
endif endif
!print *,p,q,"model ids=",pmodel,qmodel !print *,p,q,"model ids=",pmodel,qmodel
end subroutine convertOrbIdsToModelSpaceIds end subroutine convertOrbIdsToModelSpaceIds

View File

@ -65,11 +65,6 @@ subroutine get_phase_qp_to_cfg(Ialpha, Ibeta, phaseout)
integer :: nbetas integer :: nbetas
integer :: count, k integer :: count, k
! Remove the DOMOs
mask2 = IAND(Ialpha,Ibeta)
deta = IEOR(Ialpha,mask2)
detb = IEOR(Ibeta ,mask2)
! Find how many alpha electrons there are in all the N_ints ! Find how many alpha electrons there are in all the N_ints
integer :: Na(N_int) integer :: Na(N_int)
do k=1,N_int do k=1,N_int