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

Merge branch 'csf' of github.com:QuantumPackage/qp2 into csf

Browse Source
This commit is contained in:
Anthony Scemama 2021-03-01 18:33:35 +01:00
commit 46e425941c
6 changed files with 1168 additions and 242 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -1,3 +1,272 @@
subroutine obtain_associated_alphaI(idxI, Icfg, alphasIcfg, NalphaIcfg)
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
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)
@ -23,19 +292,26 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod
! Type 3 - SOMO -> VMO ! Type 3 - SOMO -> VMO
! Type 4 - DOMO -> SOMO ! Type 4 - DOMO -> SOMO
END_DOC END_DOC
integer(bit_kind),intent(in) :: Ialpha(N_int,2) integer(bit_kind),intent(in) :: Ialpha(N_int,2)
integer(bit_kind),intent(in) :: Jcfg(N_int,2) integer(bit_kind),intent(in) :: Jcfg(N_int,2)
integer,intent(in) :: p,q integer,intent(in) :: p,q
integer,intent(in) :: extype integer,intent(in) :: extype
integer,intent(out) :: pmodel,qmodel integer,intent(out) :: pmodel,qmodel
integer*8 :: Isomo !integer(bit_kind) :: Isomo(N_int)
integer*8 :: Idomo !integer(bit_kind) :: Idomo(N_int)
integer*8 :: Jsomo !integer(bit_kind) :: Jsomo(N_int)
integer*8 :: Jdomo !integer(bit_kind) :: Jdomo(N_int)
integer*8 :: mask integer*8 :: Isomo
integer*8 :: Isomotmp integer*8 :: Idomo
integer*8 :: Jsomotmp integer*8 :: Jsomo
integer :: pos0,pos0prev integer*8 :: Jdomo
integer*8 :: mask
integer :: iint, ipos
!integer(bit_kind) :: Isomotmp(N_int)
!integer(bit_kind) :: Jsomotmp(N_int)
integer*8 :: Isomotmp
integer*8 :: Jsomotmp
integer :: pos0,pos0prev
! TODO Flag (print) when model space indices is > 64 ! TODO Flag (print) when model space indices is > 64
Isomo = Ialpha(1,1) Isomo = Ialpha(1,1)
@ -50,11 +326,6 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod
pmodel = 1 pmodel = 1
qmodel = 1 qmodel = 1
else else
!print *,"input pq=",p,q,"extype=",extype
!call debug_spindet(Isomo,1)
!call debug_spindet(Idomo,1)
!call debug_spindet(Jsomo,1)
!call debug_spindet(Jdomo,1)
select case(extype) select case(extype)
case (1) case (1)
! SOMO -> SOMO ! SOMO -> SOMO
@ -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

211
src/csf/configuration_CI_sigma_helpers.org
View File

@ -43,8 +43,10 @@ the input determinant \(|D_I\rangle\).
integer*8 :: Jsomo integer*8 :: Jsomo
integer*8 :: diffSOMO integer*8 :: diffSOMO
integer*8 :: diffDOMO integer*8 :: diffDOMO
integer*8 :: xordiffSOMODOMO
integer :: ndiffSOMO integer :: ndiffSOMO
integer :: ndiffDOMO integer :: ndiffDOMO
integer :: nxordiffSOMODOMO
integer :: ndiffAll integer :: ndiffAll
integer :: i integer :: i
integer :: j integer :: j
@ -55,6 +57,7 @@ the input determinant \(|D_I\rangle\).
integer :: countalphas integer :: countalphas
logical :: pqAlreadyGenQ logical :: pqAlreadyGenQ
logical :: pqExistsQ logical :: pqExistsQ
logical :: ppExistsQ
Isomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,1)) Isomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,1))
Idomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,2)) Idomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,2))
!print*,"Input cfg" !print*,"Input cfg"
@ -66,7 +69,7 @@ the input determinant \(|D_I\rangle\).
! find out all pq holes possible ! find out all pq holes possible
nholes = 0 nholes = 0
! holes in SOMO ! holes in SOMO
do i = n_core_orb+1,n_core_orb + n_act_orb do i = 1,mo_num
if(POPCNT(IAND(Isomo,IBSET(0_8,i-1))) .EQ. 1) then if(POPCNT(IAND(Isomo,IBSET(0_8,i-1))) .EQ. 1) then
nholes += 1 nholes += 1
listholes(nholes) = i listholes(nholes) = i
@ -74,7 +77,7 @@ the input determinant \(|D_I\rangle\).
endif endif
end do end do
! holes in DOMO ! holes in DOMO
do i = n_core_orb+1,n_core_orb + n_act_orb do i = 1,mo_num
if(POPCNT(IAND(Idomo,IBSET(0_8,i-1))) .EQ. 1) then if(POPCNT(IAND(Idomo,IBSET(0_8,i-1))) .EQ. 1) then
nholes += 1 nholes += 1
listholes(nholes) = i listholes(nholes) = i
@ -86,7 +89,7 @@ the input determinant \(|D_I\rangle\).
listvmos = -1 listvmos = -1
vmotype = -1 vmotype = -1
nvmos = 0 nvmos = 0
do i = n_core_orb+1,n_core_orb + n_act_orb 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)))) !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 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 nvmos += 1
@ -126,7 +129,7 @@ the input determinant \(|D_I\rangle\).
p = listholes(i) p = listholes(i)
do j = 1,nvmos do j = 1,nvmos
q = listvmos(j) q = listvmos(j)
if(p == q) cycle if(p .EQ. q) cycle
if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 1) then if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 1) then
! SOMO -> VMO ! SOMO -> VMO
Jsomo = IBCLR(Isomo,p-1) Jsomo = IBCLR(Isomo,p-1)
@ -158,10 +161,19 @@ the input determinant \(|D_I\rangle\).
do k = 1, idxI-1 do k = 1, idxI-1
diffSOMO = IEOR(Jsomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k))) 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))) diffDOMO = IEOR(Jdomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k)))
xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO)
ndiffSOMO = POPCNT(diffSOMO) ndiffSOMO = POPCNT(diffSOMO)
ndiffDOMO = POPCNT(diffDOMO) ndiffDOMO = POPCNT(diffDOMO)
if(POPCNT(IEOR(diffSOMO,diffDOMO)) .LE. 1 .AND. ndiffDOMO .LT. 3) then nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO)
!if(POPCNT(IEOR(diffSOMO,diffDOMO)) .LE. 1 .AND. ndiffDOMO .LT. 3) then
if((ndiffSOMO+ndiffDOMO) .EQ. 0) then
pqAlreadyGenQ = .TRUE. 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 !print *,i,k,ndiffSOMO,ndiffDOMO
!call debug_spindet(Jsomo,1) !call debug_spindet(Jsomo,1)
!call debug_spindet(Jdomo,1) !call debug_spindet(Jdomo,1)
@ -171,20 +183,22 @@ the input determinant \(|D_I\rangle\).
endif endif
end do end do
!print *,"(,",p,",",q,")",pqAlreadyGenQ
if(pqAlreadyGenQ) cycle if(pqAlreadyGenQ) cycle
pqExistsQ = .FALSE. pqExistsQ = .FALSE.
! now check if this exists in the selected list ! now check if this exists in the selected list
do k = idxI, N_configuration !do k = idxI+1, N_configuration
diffSOMO = IEOR(OR(reunion_of_act_virt_bitmask(1,1),Jsomo),psi_configuration(1,1,k)) ! 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)) ! diffDOMO = IEOR(OR(reunion_of_act_virt_bitmask(1,1),Jdomo),psi_configuration(1,2,k))
ndiffSOMO = POPCNT(diffSOMO) ! ndiffSOMO = POPCNT(diffSOMO)
ndiffDOMO = POPCNT(diffDOMO) ! ndiffDOMO = POPCNT(diffDOMO)
if((ndiffSOMO + ndiffDOMO) .EQ. 0) then ! if((ndiffSOMO + ndiffDOMO) .EQ. 0) then
pqExistsQ = .TRUE. ! pqExistsQ = .TRUE.
EXIT ! EXIT
endif ! endif
end do !end do
if(.NOT. pqExistsQ) then if(.NOT. pqExistsQ) then
tableUniqueAlphas(p,q) = .TRUE. tableUniqueAlphas(p,q) = .TRUE.
@ -234,16 +248,42 @@ the input determinant \(|D_I\rangle\).
print*,"Something went wrong in obtain_associated_alphaI" print*,"Something went wrong in obtain_associated_alphaI"
endif endif
! SOMO
NalphaIcfg += 1 NalphaIcfg += 1
!print *,p,q,"|",holetype(i),vmotype(j),NalphaIcfg !print *,p,q,"|",holetype(i),vmotype(j),NalphaIcfg
!call debug_spindet(Idomo,1)
!call debug_spindet(Jdomo,1)
alphasIcfg(1,1,NalphaIcfg) = Jsomo alphasIcfg(1,1,NalphaIcfg) = Jsomo
alphasIcfg(1,2,NalphaIcfg) = IOR(Jdomo,ISHFT(1_8,n_core_orb)-1) alphasIcfg(1,2,NalphaIcfg) = IOR(Jdomo,ISHFT(1_8,n_core_orb)-1)
endif endif
end do end do
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 end subroutine
#+end_src #+end_src
@ -253,7 +293,7 @@ Next step is to obtain the connected CFGs \(|I\rangle\) that belong to the selec
given a RI configuration \(|\alpha\rangle\). given a RI configuration \(|\alpha\rangle\).
#+begin_src f90 :main no :tangle ../cfgCI/obtain_I_foralpha.irp.f #+begin_src f90 :main no :tangle ../cfgCI/obtain_I_foralpha.irp.f
subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI, nconnectedI, excitationIds, excitationTypes) subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI, nconnectedI, excitationIds, excitationTypes, diagfactors)
implicit none implicit none
use bitmasks use bitmasks
BEGIN_DOC BEGIN_DOC
@ -278,6 +318,7 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI
integer,intent(out) :: nconnectedI integer,intent(out) :: nconnectedI
integer,intent(out) :: excitationIds(2,*) integer,intent(out) :: excitationIds(2,*)
integer,intent(out) :: excitationTypes(*) integer,intent(out) :: excitationTypes(*)
real*8 ,intent(out) :: diagfactors(*)
integer*8 :: Idomo integer*8 :: Idomo
integer*8 :: Isomo integer*8 :: Isomo
integer*8 :: Jdomo integer*8 :: Jdomo
@ -285,38 +326,19 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI
integer*8 :: IJsomo integer*8 :: IJsomo
integer*8 :: diffSOMO integer*8 :: diffSOMO
integer*8 :: diffDOMO integer*8 :: diffDOMO
integer*8 :: xordiffSOMODOMO
integer :: ndiffSOMO integer :: ndiffSOMO
integer :: ndiffDOMO integer :: ndiffDOMO
integer :: i,j,k,l,p,q,nsomoJ,nsomoalpha,starti,endi,extyp,nholes integer :: nxordiffSOMODOMO
integer :: ii,i,j,k,l,p,q,nsomoJ,nsomoalpha,starti,endi,extyp,nholes
integer :: listholes(mo_num) integer :: listholes(mo_num)
integer :: holetype(mo_num) integer :: holetype(mo_num)
! find out all pq holes possible
nholes = 0
! holes in SOMO
Isomo = psi_configuration(1,1,idxI)
Idomo = psi_configuration(1,2,idxI)
do i = n_core_orb+1,n_core_orb + n_act_orb
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 = n_core_orb+1,n_core_orb + n_act_orb
if(POPCNT(IAND(Idomo,IBSET(0_8,i-1))) .EQ. 1) then
nholes += 1
listholes(nholes) = i
holetype(nholes) = 2
endif
end do
nconnectedI = 0 nconnectedI = 0
p = 0 p = 0
q = 0 q = 0
do i=idxI+1,N_configuration do i=idxI,N_configuration
Isomo = Ialpha(1,1) Isomo = Ialpha(1,1)
Idomo = Ialpha(1,2) Idomo = Ialpha(1,2)
Jsomo = psi_configuration(1,1,i) Jsomo = psi_configuration(1,1,i)
@ -328,12 +350,15 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI
!call debug_spindet(Jdomo,1) !call debug_spindet(Jdomo,1)
diffSOMO = IEOR(Isomo,Jsomo) diffSOMO = IEOR(Isomo,Jsomo)
diffDOMO = IEOR(Idomo,Jdomo) diffDOMO = IEOR(Idomo,Jdomo)
xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO)
ndiffSOMO = POPCNT(diffSOMO) ndiffSOMO = POPCNT(diffSOMO)
ndiffDOMO = POPCNT(diffDOMO) ndiffDOMO = POPCNT(diffDOMO)
if((ndiffSOMO + ndiffDOMO) .EQ. 0) cycle nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO)
!print *,"-I--i=",i,diffSOMO,diffDOMO!Isomo,Jsomo,ndiffSOMO,ndiffDOMO !print *,"-I--i=",i,ndiffSOMO,ndiffDOMO,nxordiffSOMODOMO!Isomo,Jsomo,ndiffSOMO,ndiffDOMO
!if((ndiffSOMO + ndiffDOMO) .EQ. 0) cycle
!print *,POPCNT(IEOR(diffSOMO,diffDOMO)), ndiffDOMO !print *,POPCNT(IEOR(diffSOMO,diffDOMO)), ndiffDOMO
if(POPCNT(IEOR(diffSOMO,diffDOMO)) .LE. 1 .AND. ndiffDOMO .LT. 3) then !if(POPCNT(IEOR(diffSOMO,diffDOMO)) .LE. 1 .AND. ndiffDOMO .LT. 3) then
if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then
!call debug_spindet(Isomo,1) !call debug_spindet(Isomo,1)
!call debug_spindet(Idomo,1) !call debug_spindet(Idomo,1)
!print *,"-J--i=",i,Idomo,Jdomo,">",N_configuration !print *,"-J--i=",i,Idomo,Jdomo,">",N_configuration
@ -390,7 +415,58 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI
excitationIds(1,nconnectedI)=p excitationIds(1,nconnectedI)=p
excitationIds(2,nconnectedI)=q excitationIds(2,nconnectedI)=q
excitationTypes(nconnectedI) = extyp excitationTypes(nconnectedI) = extyp
print *,"------ > output p,q in obt=",p,q diagfactors(nconnectedI) = 1.0d0
!print *,"------ > output p,q in obt=",p,q
else if((ndiffSOMO + ndiffDOMO) .EQ. 0) then
! find out all pq holes possible
nholes = 0
! holes in SOMO
Isomo = psi_configuration(1,1,i)
Idomo = psi_configuration(1,2,i)
do ii = 1,mo_num
if(POPCNT(IAND(Isomo,IBSET(0_8,ii-1))) .EQ. 1) then
nholes += 1
listholes(nholes) = ii
holetype(nholes) = 1
endif
end do
! holes in DOMO
do ii = 1,mo_num
if(POPCNT(IAND(Idomo,IBSET(0_8,ii-1))) .EQ. 1) then
nholes += 1
listholes(nholes) = ii
holetype(nholes) = 2
endif
end do
do k=1,nholes
p = listholes(k)
q = p
extyp = 1
if(holetype(k) .EQ. 1) then
starti = psi_config_data(i,1)
endi = psi_config_data(i,2)
nconnectedI += 1
connectedI(:,:,nconnectedI) = psi_configuration(:,:,i)
idxs_connectedI(nconnectedI)=starti
excitationIds(1,nconnectedI)=p
excitationIds(2,nconnectedI)=q
excitationTypes(nconnectedI) = extyp
diagfactors(nconnectedI) = 1.0d0
!print *,"------ > output p,q in obt=",p,q
else
starti = psi_config_data(i,1)
endi = psi_config_data(i,2)
nconnectedI += 1
connectedI(:,:,nconnectedI) = psi_configuration(:,:,i)
idxs_connectedI(nconnectedI)=starti
excitationIds(1,nconnectedI)=p
excitationIds(2,nconnectedI)=q
excitationTypes(nconnectedI) = extyp
diagfactors(nconnectedI) = 2.0d0
!print *,"------ > output p,q in obt=",p,q
endif
enddo
endif endif
end do end do
@ -468,6 +544,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
@ -501,28 +578,46 @@ 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
#+end_src #+end_src

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

@ -10,7 +10,7 @@ subroutine convertWFfromDETtoCSF(N_st,psi_coef_det_in, psi_coef_cfg_out)
integer, intent(in) :: N_st integer, intent(in) :: N_st
double precision, intent(in) :: psi_coef_det_in(N_det,N_st) double precision, intent(in) :: psi_coef_det_in(N_det,N_st)
double precision, intent(out) :: psi_coef_cfg_out(n_CSF,N_st) double precision, intent(out) :: psi_coef_cfg_out(n_CSF,N_st)
integer*8 :: Isomo, Idomo, mask integer*8 :: Isomo, Idomo
integer(bit_kind) :: Ialpha(N_int) ,Ibeta(N_int) integer(bit_kind) :: Ialpha(N_int) ,Ibeta(N_int)
integer :: rows, cols, i, j, k integer :: rows, cols, i, j, k
integer :: startdet, enddet integer :: startdet, enddet
@ -65,7 +65,7 @@ subroutine convertWFfromDETtoCSF(N_st,psi_coef_det_in, psi_coef_cfg_out)
countcsf += bfIcfg countcsf += bfIcfg
enddo enddo
end end subroutine convertWFfromDETtoCSF
subroutine convertWFfromCSFtoDET(N_st,psi_coef_cfg_in, psi_coef_det) subroutine convertWFfromCSFtoDET(N_st,psi_coef_cfg_in, psi_coef_det)
@ -88,11 +88,11 @@ subroutine convertWFfromCSFtoDET(N_st,psi_coef_cfg_in, psi_coef_det)
integer :: ndetI integer :: ndetI
integer :: getNSOMO integer :: getNSOMO
double precision,allocatable :: tempBuffer(:,:) double precision,allocatable :: tempBuffer(:,:)
double precision,allocatable :: tempCoeff (:,:) double precision,allocatable :: tempCoeff(:,:)
double precision :: phasedet double precision :: phasedet
integer :: idx integer :: idx
countcsf = 0 countcsf = 0
do i = 1,N_configuration do i = 1,N_configuration
startdet = psi_configuration_to_psi_det(1,i) startdet = psi_configuration_to_psi_det(1,i)
@ -138,11 +138,4 @@ subroutine convertWFfromCSFtoDET(N_st,psi_coef_cfg_in, psi_coef_det)
deallocate(tempBuffer) deallocate(tempBuffer)
enddo enddo
end end subroutine convertCSFtoDET

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

@ -1,9 +1,9 @@
BEGIN_PROVIDER [ integer, NSOMOMax] BEGIN_PROVIDER [ integer, NSOMOMax]
&BEGIN_PROVIDER [ integer, NCSFMax] &BEGIN_PROVIDER [ integer, NCSFMax]
&BEGIN_PROVIDER [ integer*8, NMO] &BEGIN_PROVIDER [ integer*8, NMO]
&BEGIN_PROVIDER [ integer, NBFMax] &BEGIN_PROVIDER [ integer, NBFMax]
&BEGIN_PROVIDER [ integer, n_CSF] &BEGIN_PROVIDER [ integer, n_CSF]
&BEGIN_PROVIDER [ integer, maxDetDimPerBF] &BEGIN_PROVIDER [ integer, maxDetDimPerBF]
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Documentation for NSOMOMax ! Documentation for NSOMOMax
@ -22,29 +22,38 @@
integer NSOMO integer NSOMO
integer dimcsfpercfg integer dimcsfpercfg
integer detDimperBF integer detDimperBF
real*8 :: coeff real*8 :: coeff
integer MS integer MS
integer ncfgpersomo integer ncfgpersomo
detDimperBF = 0 detDimperBF = 0
MS = elec_alpha_num-elec_beta_num MS = elec_alpha_num-elec_beta_num
!print *,"NSOMOMax=",NSOMOMax, cfg_seniority_index(0)
! number of cfgs = number of dets for 0 somos ! number of cfgs = number of dets for 0 somos
n_CSF = cfg_seniority_index(0)-1 n_CSF = cfg_seniority_index(0)-1
ncfgprev = cfg_seniority_index(0) ncfgprev = cfg_seniority_index(0)
do i = 0-iand(MS,1)+2, NSOMOMax,2 do i = 0-iand(MS,1)+2, NSOMOMax,2
if(cfg_seniority_index(i) .EQ. -1)then if(cfg_seniority_index(i) .EQ. -1)then
ncfgpersomo = N_configuration + 1 ncfgpersomo = N_configuration + 1
else else
ncfgpersomo = cfg_seniority_index(i) ncfgpersomo = cfg_seniority_index(i)
endif endif
ncfg = ncfgpersomo - ncfgprev ncfg = ncfgpersomo - ncfgprev
!detDimperBF = max(1,nint((binom(i,(i+1)/2)))) !detDimperBF = max(1,nint((binom(i,(i+1)/2))))
dimcsfpercfg = max(1,nint((binom(i-2,(i-2+1)/2)-binom(i-2,((i-2+1)/2)+1)))) dimcsfpercfg = max(1,nint((binom(i-2,(i-2+1)/2)-binom(i-2,((i-2+1)/2)+1))))
n_CSF += ncfg * dimcsfpercfg n_CSF += ncfg * dimcsfpercfg
!if(cfg_seniority_index(i+2) == -1) EXIT !print *,i,">(",ncfg,ncfgprev,ncfgpersomo,")",",",detDimperBF,">",dimcsfpercfg, " | dimbas= ", n_CSF
!if(detDimperBF > maxDetDimPerBF) maxDetDimPerBF = detDimperBF !if(cfg_seniority_index(i+2) == -1) EXIT
ncfgprev = cfg_seniority_index(i) !if(detDimperBF > maxDetDimPerBF) maxDetDimPerBF = detDimperBF
ncfgprev = cfg_seniority_index(i)
enddo enddo
END_PROVIDER if(NSOMOMax .EQ. elec_num)then
ncfgpersomo = N_configuration + 1
ncfg = ncfgpersomo - ncfgprev
dimcsfpercfg = max(1,nint((binom(i-2,(i-2+1)/2)-binom(i-2,((i-2+1)/2)+1))))
n_CSF += ncfg * dimcsfpercfg
!print *,i,">(",ncfg,ncfgprev,ncfgpersomo,")",",",detDimperBF,">",dimcsfpercfg, " | dimbas= ", n_CSF
endif
END_PROVIDER
subroutine get_phase_qp_to_cfg(Ialpha, Ibeta, phaseout) subroutine get_phase_qp_to_cfg(Ialpha, Ibeta, phaseout)
use bitmasks use bitmasks
@ -61,14 +70,13 @@ subroutine get_phase_qp_to_cfg(Ialpha, Ibeta, phaseout)
integer(bit_kind),intent(in) :: Ialpha(N_int) integer(bit_kind),intent(in) :: Ialpha(N_int)
integer(bit_kind),intent(in) :: Ibeta(N_int) integer(bit_kind),intent(in) :: Ibeta(N_int)
real*8,intent(out) :: phaseout real*8,intent(out) :: phaseout
integer(bit_kind) :: mask, mask2(N_int), deta(N_int), detb(N_int) integer(bit_kind) :: mask, deta(N_int), detb(N_int)
integer :: nbetas integer :: nbetas
integer :: count, k integer :: count, k
! Remove the DOMOs ! Initliaze deta and detb
mask2 = IAND(Ialpha,Ibeta) deta = Ialpha
deta = IEOR(Ialpha,mask2) detb = Ibeta
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)
@ -86,7 +94,7 @@ subroutine get_phase_qp_to_cfg(Ialpha, Ibeta, phaseout)
detb(k) = ibclr(detb(k),ipos) detb(k) = ibclr(detb(k),ipos)
! Create a mask will all MOs higher than the beta electron ! Create a mask will all MOs higher than the beta electron
mask = not(shiftl(1_bit_kind,ipos+1) - 1_bit_kind) mask = not(shiftl(1_bit_kind,ipos + 1) - 1_bit_kind)
! Apply the mask to the alpha string to count how many electrons to cross ! Apply the mask to the alpha string to count how many electrons to cross
nperm = popcnt( iand(mask, deta(k)) ) nperm = popcnt( iand(mask, deta(k)) )
@ -104,7 +112,111 @@ end subroutine get_phase_qp_to_cfg
BEGIN_PROVIDER [ integer, AIJpqMatrixDimsList, (0:NSOMOMax,0:NSOMOMax,4,NSOMOMax,NSOMOMax,2)] BEGIN_PROVIDER [ real*8, DetToCSFTransformationMatrix, (0:NSOMOMax,NBFMax,maxDetDimPerBF)]
&BEGIN_PROVIDER [ real*8, psi_coef_config, (n_CSF,1)]
&BEGIN_PROVIDER [ integer, psi_config_data, (N_configuration,2)]
&BEGIN_PROVIDER [ integer, psi_csf_to_config_data, (n_CSF)]
use cfunctions
implicit none
BEGIN_DOC
! Documentation for DetToCSFTransformationMatrix
! Provides the matrix of transformatons for the
! conversion between determinant to CSF basis (in BFs)
END_DOC
integer*8 :: Isomo, Idomo
integer(bit_kind) :: Ialpha(N_int),Ibeta(N_int)
integer :: rows, cols, i, j, k
integer :: startdet, enddet
integer*8 MS
integer ndetI
integer :: getNSOMO
real*8,dimension(:,:),allocatable :: tempBuffer
real*8,dimension(:),allocatable :: tempCoeff
real*8 :: norm_det1, phasedet
norm_det1 = 0.d0
MS = elec_alpha_num - elec_beta_num
print *,"Maxbfdim=",NBFMax
print *,"Maxdetdim=",maxDetDimPerBF
print *,"n_CSF=",n_CSF
print *,"N_configurations=",N_configuration
print *,"n_core_orb=",n_core_orb
! initialization
psi_coef_config = 0.d0
DetToCSFTransformationMatrix(0,:,:) = 1.d0
do i = 2-iand(elec_alpha_num-elec_beta_num,1), NSOMOMax,2
Isomo = IBSET(0_8, i) - 1_8
! rows = Ncsfs
! cols = Ndets
bfIcfg = max(1,nint((binom(i,(i+1)/2)-binom(i,((i+1)/2)+1))))
ndetI = max(1,nint((binom(i,(i+1)/2))))
allocate(tempBuffer(bfIcfg,ndetI))
call getCSFtoDETTransformationMatrix(Isomo, MS, NBFMax, maxDetDimPerBF, tempBuffer)
DetToCSFTransformationMatrix(i,:bfIcfg,:ndetI) = tempBuffer
deallocate(tempBuffer)
enddo
integer s, bfIcfg
integer countcsf
countcsf = 0
integer countdet
countdet = 0
integer istate
istate = 1
phasedet = 1.0d0
do i = 1,N_configuration
startdet = psi_configuration_to_psi_det(1,i)
enddet = psi_configuration_to_psi_det(2,i)
ndetI = enddet-startdet+1
allocate(tempCoeff(ndetI))
countdet = 1
do j = startdet, enddet
Ialpha = psi_det(:,1,psi_configuration_to_psi_det_data(j))
Ibeta = psi_det(:,2,psi_configuration_to_psi_det_data(j))
!call debug_spindet(Ialpha,1,1)
!call debug_spindet(Ibeta ,1,1)
call get_phase_qp_to_cfg(Ialpha, Ibeta, phasedet)
!print *,">>",Ialpha,Ibeta,phasedet
tempCoeff(countdet) = psi_coef(psi_configuration_to_psi_det_data(j), istate)*phasedet
!tempCoeff(countdet) = psi_coef(psi_configuration_to_psi_det_data(j), istate)
norm_det1 += tempCoeff(countdet)*tempCoeff(countdet)
countdet += 1
enddo
!print *,"dimcoef=",bfIcfg,norm_det1
!call printMatrix(tempCoeff,ndetI,1)
s = 0
do k=1,N_int
if (psi_configuration(k,1,i) == 0_bit_kind) cycle
s = s + popcnt(psi_configuration(k,1,i))
enddo
bfIcfg = max(1,nint((binom(s,(s+1)/2)-binom(s,((s+1)/2)+1))))
! perhaps blocking with CFGs of same seniority
! can be more efficient
allocate(tempBuffer(bfIcfg,ndetI))
tempBuffer = DetToCSFTransformationMatrix(s,:bfIcfg,:ndetI)
!print *,"csftodetdim=",bfIcfg,ndetI
!call printMatrix(tempBuffer,bfIcfg,ndetI)
call dgemm('N','N', bfIcfg, 1, ndetI, 1.d0, tempBuffer, size(tempBuffer,1), tempCoeff, size(tempCoeff,1), 0.d0, psi_coef_config(countcsf+1,1), size(psi_coef_config,1))
!call dgemv('N', NBFMax, maxDetDimPerBF, 1.d0, tempBuffer, size(tempBuffer,1), tempCoeff, 1, 0.d0, psi_coef_config(countcsf), 1)
!call printMatrix(psi_coef_config(countcsf+1,1),bfIcfg,1)
deallocate(tempCoeff)
deallocate(tempBuffer)
psi_config_data(i,1) = countcsf + 1
countcsf += bfIcfg
psi_config_data(i,2) = countcsf
psi_csf_to_config_data(countcsf) = i
enddo
print *,"Norm det=",norm_det1, size(psi_coef_config,1), " Dim csf=", countcsf
END_PROVIDER
BEGIN_PROVIDER [ integer, AIJpqMatrixDimsList, (0:NSOMOMax,0:NSOMOMax,4,NSOMOMax+1,NSOMOMax+1,2)]
&BEGIN_PROVIDER [ integer, rowsmax] &BEGIN_PROVIDER [ integer, rowsmax]
&BEGIN_PROVIDER [ integer, colsmax] &BEGIN_PROVIDER [ integer, colsmax]
use cfunctions use cfunctions
@ -127,9 +239,13 @@ end subroutine get_phase_qp_to_cfg
nsomomin = elec_alpha_num-elec_beta_num nsomomin = elec_alpha_num-elec_beta_num
rowsmax = 0 rowsmax = 0
colsmax = 0 colsmax = 0
print *,"NSOMOMax = ",NSOMOMax
!allocate(AIJpqMatrixDimsList(NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,2)) !allocate(AIJpqMatrixDimsList(NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,2))
! Type ! Type
! 1. SOMO -> SOMO ! 1. SOMO -> SOMO
!print *,"Doing SOMO->SOMO"
AIJpqMatrixDimsList(0,0,1,1,1,1) = 1
AIJpqMatrixDimsList(0,0,1,1,1,2) = 1
do i = 2-iand(nsomomin,1), NSOMOMax, 2 do i = 2-iand(nsomomin,1), NSOMOMax, 2
Isomo = ISHFT(1_8,i)-1 Isomo = ISHFT(1_8,i)-1
do j = i-2,i-2, 2 do j = i-2,i-2, 2
@ -157,6 +273,7 @@ end subroutine get_phase_qp_to_cfg
MS, & MS, &
rows, & rows, &
cols) cols)
!print *, "SOMO->SOMO \t",i,j,k,l,">",Isomo,Jsomo,">",rows, cols
if(rowsmax .LT. rows) then if(rowsmax .LT. rows) then
rowsmax = rows rowsmax = rows
end if end if
@ -172,6 +289,9 @@ end subroutine get_phase_qp_to_cfg
end do end do
! Type ! Type
! 2. DOMO -> VMO ! 2. DOMO -> VMO
!print *,"Doing DOMO->VMO"
AIJpqMatrixDimsList(0,0,2,1,1,1) = 1
AIJpqMatrixDimsList(0,0,2,1,1,2) = 1
do i = 0+iand(nsomomin,1), NSOMOMax, 2 do i = 0+iand(nsomomin,1), NSOMOMax, 2
Isomo = ISHFT(1_8,i)-1 Isomo = ISHFT(1_8,i)-1
tmpsomo = ISHFT(1_8,i+2)-1 tmpsomo = ISHFT(1_8,i+2)-1
@ -205,6 +325,7 @@ end subroutine get_phase_qp_to_cfg
MS, & MS, &
rows, & rows, &
cols) cols)
!print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols
if(rowsmax .LT. rows) then if(rowsmax .LT. rows) then
rowsmax = rows rowsmax = rows
end if end if
@ -221,6 +342,8 @@ end subroutine get_phase_qp_to_cfg
! Type ! Type
! 3. SOMO -> VMO ! 3. SOMO -> VMO
!print *,"Doing SOMO->VMO" !print *,"Doing SOMO->VMO"
AIJpqMatrixDimsList(0,0,3,1,1,1) = 1
AIJpqMatrixDimsList(0,0,3,1,1,2) = 1
do i = 2-iand(nsomomin,1), NSOMOMax, 2 do i = 2-iand(nsomomin,1), NSOMOMax, 2
Isomo = ISHFT(1_8,i)-1 Isomo = ISHFT(1_8,i)-1
do j = i,i, 2 do j = i,i, 2
@ -228,8 +351,8 @@ end subroutine get_phase_qp_to_cfg
if(j .GT. NSOMOMax .OR. j .LE. 0) then if(j .GT. NSOMOMax .OR. j .LE. 0) then
cycle cycle
end if end if
do k = 1,i do k = 1,i+1
do l = 1,i do l = 1,i+1
if(k .NE. l) then if(k .NE. l) then
Isomo = ISHFT(1_8,i+1)-1 Isomo = ISHFT(1_8,i+1)-1
Isomo = IBCLR(Isomo,l-1) Isomo = IBCLR(Isomo,l-1)
@ -244,6 +367,7 @@ end subroutine get_phase_qp_to_cfg
MS, & MS, &
rows, & rows, &
cols) cols)
!print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols
if(rowsmax .LT. rows) then if(rowsmax .LT. rows) then
rowsmax = rows rowsmax = rows
end if end if
@ -258,18 +382,20 @@ end subroutine get_phase_qp_to_cfg
end do end do
end do end do
! Type ! Type
! 4. DOMO -> VMO ! 4. DOMO -> SOMO
!print *,"Doing DOMO->SOMO" !print *,"Doing DOMO->SOMO"
AIJpqMatrixDimsList(0,0,4,1,1,1) = 1
AIJpqMatrixDimsList(0,0,4,1,1,2) = 1
do i = 2-iand(nsomomin,1), NSOMOMax, 2 do i = 2-iand(nsomomin,1), NSOMOMax, 2
do j = i,i, 2 do j = i,i, 2
if(j .GT. NSOMOMax .OR. j .LE. 0) then if(j .GT. NSOMOMax .OR. j .LE. 0) then
cycle cycle
end if end if
do k = 1,i do k = 1,i+1
do l = 1,i do l = 1,i+1
if(k .NE. l) then if(k .NE. l) then
Isomo = ISHFT(1_8,i+1)-1 Isomo = ISHFT(1_8,i+1)-1
Isomo = IBCLR(Isomo,k+1-1) Isomo = IBCLR(Isomo,k-1)
Jsomo = ISHFT(1_8,j+1)-1 Jsomo = ISHFT(1_8,j+1)-1
Jsomo = IBCLR(Jsomo,l-1) Jsomo = IBCLR(Jsomo,l-1)
else else
@ -281,6 +407,7 @@ end subroutine get_phase_qp_to_cfg
MS, & MS, &
rows, & rows, &
cols) cols)
!print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols
if(rowsmax .LT. rows) then if(rowsmax .LT. rows) then
rowsmax = rows rowsmax = rows
end if end if
@ -294,9 +421,10 @@ end subroutine get_phase_qp_to_cfg
end do end do
end do end do
end do end do
print *,"Rowsmax=",rowsmax," Colsmax=",colsmax
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ real*8, AIJpqContainer, (0:NSOMOMax,0:NSOMOMax,4,NSOMOMax,NSOMOMax,NBFMax,NBFMax)] BEGIN_PROVIDER [ real*8, AIJpqContainer, (0:NSOMOMax,0:NSOMOMax,4,NSOMOMax+1,NSOMOMax+1,NBFMax,NBFMax)]
use cfunctions use cfunctions
implicit none implicit none
BEGIN_DOC BEGIN_DOC
@ -332,13 +460,18 @@ end subroutine get_phase_qp_to_cfg
integer maxdim integer maxdim
!maxdim = max(rowsmax,colsmax) !maxdim = max(rowsmax,colsmax)
! allocate matrix ! allocate matrix
!print *,"rowsmax =",rowsmax," colsmax=",colsmax
!print *,"NSOMOMax = ",NSOMOMax
!allocate(AIJpqMatrixDimsList(NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,2)) !allocate(AIJpqMatrixDimsList(NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,2))
! Type ! Type
! 1. SOMO -> SOMO ! 1. SOMO -> SOMO
!print *,"Doing SOMO -> SOMO"
AIJpqContainer(0,0,1,1,1,1,1) = 1.0d0
do i = 2, NSOMOMax, 2 do i = 2, NSOMOMax, 2
Isomo = ISHFT(1_8,i)-1 Isomo = ISHFT(1_8,i)-1
do j = i-2,i-2, 2 do j = i-2,i-2, 2
if(j .GT. NSOMOMax .OR. j .LT. 0) cycle if(j .GT. NSOMOMax .OR. j .LT. 0) cycle
!print *,"i,j=",i,j
do k = 1,i do k = 1,i
do l = 1,i do l = 1,i
@ -355,6 +488,10 @@ end subroutine get_phase_qp_to_cfg
nsomoj = i nsomoj = i
endif endif
!print *,"k,l=",k,l
!call debug_spindet(Jsomo,1)
!call debug_spindet(Isomo,1)
AIJpqContainer(nsomoi,nsomoj,1,k,l,:,:) = 0.0d0 AIJpqContainer(nsomoi,nsomoj,1,k,l,:,:) = 0.0d0
call getApqIJMatrixDims(Isomo, & call getApqIJMatrixDims(Isomo, &
Jsomo, & Jsomo, &
@ -376,6 +513,8 @@ end subroutine get_phase_qp_to_cfg
meMatrix, & meMatrix, &
rows, & rows, &
cols) cols)
!print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols,">",rowsmax,colsmax
!call printMatrix(meMatrix,rows,cols)
! i -> j ! i -> j
do ri = 1,rows do ri = 1,rows
do ci = 1,cols do ci = 1,cols
@ -389,6 +528,8 @@ end subroutine get_phase_qp_to_cfg
end do end do
! Type ! Type
! 2. DOMO -> VMO ! 2. DOMO -> VMO
!print *,"Doing DOMO -> VMO"
AIJpqContainer(0,0,2,1,1,1,1) = 1.0d0
do i = 0, NSOMOMax, 2 do i = 0, NSOMOMax, 2
Isomo = ISHFT(1_8,i)-1 Isomo = ISHFT(1_8,i)-1
tmpsomo = ISHFT(1_8,i+2)-1 tmpsomo = ISHFT(1_8,i+2)-1
@ -414,6 +555,10 @@ end subroutine get_phase_qp_to_cfg
nsomoj = j nsomoj = j
endif endif
!print *,"k,l=",k,l
!call debug_spindet(Jsomo,1)
!call debug_spindet(Isomo,1)
AIJpqContainer(nsomoi,nsomoj,2,k,l,:,:) = 0.0d0 AIJpqContainer(nsomoi,nsomoj,2,k,l,:,:) = 0.0d0
call getApqIJMatrixDims(Isomo, & call getApqIJMatrixDims(Isomo, &
Jsomo, & Jsomo, &
@ -435,6 +580,8 @@ end subroutine get_phase_qp_to_cfg
meMatrix, & meMatrix, &
rows, & rows, &
cols) cols)
!print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols,">",rowsmax,colsmax
!call printMatrix(meMatrix,rows,cols)
! i -> j ! i -> j
do ri = 1,rows do ri = 1,rows
do ci = 1,cols do ci = 1,cols
@ -448,13 +595,14 @@ end subroutine get_phase_qp_to_cfg
end do end do
! Type ! Type
! 3. SOMO -> VMO ! 3. SOMO -> VMO
!print *,"Doing SOMO -> VMO"
do i = 2, NSOMOMax, 2 do i = 2, NSOMOMax, 2
Isomo = ISHFT(1_8,i)-1 Isomo = ISHFT(1_8,i)-1
do j = i,i, 2 do j = i,i, 2
Jsomo = ISHFT(1_8,j)-1 Jsomo = ISHFT(1_8,j)-1
if(j .GT. NSOMOMax .OR. j .LE. 0) cycle if(j .GT. NSOMOMax .OR. j .LE. 0) cycle
do k = 1,i do k = 1,i+1
do l = 1,i do l = 1,i+1
if(k .NE. l) then if(k .NE. l) then
Isomo = ISHFT(1_8,i+1)-1 Isomo = ISHFT(1_8,i+1)-1
Isomo = IBCLR(Isomo,l-1) Isomo = IBCLR(Isomo,l-1)
@ -465,6 +613,10 @@ end subroutine get_phase_qp_to_cfg
Jsomo = ISHFT(1_8,j)-1 Jsomo = ISHFT(1_8,j)-1
endif endif
!print *,"k,l=",k,l
!call debug_spindet(Jsomo,1)
!call debug_spindet(Isomo,1)
AIJpqContainer(i,j,3,k,l,:,:) = 0.0d0 AIJpqContainer(i,j,3,k,l,:,:) = 0.0d0
call getApqIJMatrixDims(Isomo, & call getApqIJMatrixDims(Isomo, &
Jsomo, & Jsomo, &
@ -486,6 +638,8 @@ end subroutine get_phase_qp_to_cfg
meMatrix, & meMatrix, &
rows, & rows, &
cols) cols)
!call printMatrix(meMatrix,rows,cols)
!print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols,">",rowsmax,colsmax
! i -> j ! i -> j
do ri = 1,rows do ri = 1,rows
do ci = 1,cols do ci = 1,cols
@ -499,18 +653,20 @@ end subroutine get_phase_qp_to_cfg
end do end do
! Type ! Type
! 4. DOMO -> SOMO ! 4. DOMO -> SOMO
!print *,"Doing DOMO -> SOMO"
AIJpqContainer(0,0,4,1,1,1,1) = 1.0d0
do i = 2, NSOMOMax, 2 do i = 2, NSOMOMax, 2
Isomo = ISHFT(1_8,i)-1 Isomo = ISHFT(1_8,i)-1
do j = i,i, 2 do j = i,i, 2
Jsomo = ISHFT(1_8,i)-1 Jsomo = ISHFT(1_8,i)-1
if(j .GT. NSOMOMax .OR. j .LE. 0) cycle if(j .GT. NSOMOMax .OR. j .LE. 0) cycle
do k = 1,i do k = 1,i+1
do l = 1,i do l = 1,i+1
if(k .NE. l) then if(k .NE. l) then
Isomo = ISHFT(1_8,i+1)-1 Isomo = ISHFT(1_8,i+1)-1
Isomo = IBCLR(Isomo,k-1) Isomo = IBCLR(Isomo,k-1)
Jsomo = ISHFT(1_8,j+1)-1 Jsomo = ISHFT(1_8,j+1)-1
Jsomo = IBCLR(Jsomo,l+1-1) Jsomo = IBCLR(Jsomo,l-1)
else else
Isomo = ISHFT(1_8,i)-1 Isomo = ISHFT(1_8,i)-1
Jsomo = ISHFT(1_8,j)-1 Jsomo = ISHFT(1_8,j)-1
@ -538,6 +694,8 @@ end subroutine get_phase_qp_to_cfg
meMatrix, & meMatrix, &
rows, & rows, &
cols) cols)
!call printMatrix(meMatrix,rows,cols)
!print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols,">",rowsmax,colsmax
! i -> j ! i -> j
do ri = 1,rows do ri = 1,rows
do ci = 1,cols do ci = 1,cols
@ -551,93 +709,466 @@ end subroutine get_phase_qp_to_cfg
end do end do
END_PROVIDER END_PROVIDER
subroutine calculate_preconditioner_cfg(diag_energies)
!!!!!!
BEGIN_PROVIDER [ real*8, DetToCSFTransformationMatrix, (0:NSOMOMax,NBFMax,maxDetDimPerBF)]
&BEGIN_PROVIDER [ real*8, psi_coef_config, (n_CSF)]
&BEGIN_PROVIDER [ integer, psi_config_data, (N_configuration,2)]
use cfunctions
use bitmasks
implicit none implicit none
use bitmasks
BEGIN_DOC BEGIN_DOC
! Documentation for DetToCSFTransformationMatrix ! Documentation for calculate_preconditioner
! Provides the matrix of transformatons for the !
! conversion between determinant to CSF basis (in BFs) ! Calculates the diagonal energies of
! the configurations in psi_configuration
! returns : diag_energies :
END_DOC END_DOC
integer(bit_kind) :: mask(N_int), Ialpha(N_int),Ibeta(N_int) integer :: i,j,k,l,p,q,noccp,noccq, ii, jj
integer :: rows, cols, i, j, k real*8,intent(out) :: diag_energies(n_CSF)
integer :: startdet, enddet integer :: nholes
integer*8 MS, Isomo, Idomo integer :: nvmos
integer ndetI integer :: listvmos(mo_num)
integer :: vmotype(mo_num) ! 1 -> VMO 2 -> SOMO
integer :: listholes(mo_num)
integer :: holetype(mo_num) ! 1-> SOMO 2->DOMO
integer*8 :: Idomo
integer*8 :: Isomo
integer*8 :: Jdomo
integer*8 :: Jsomo
integer*8 :: diffSOMO
integer*8 :: diffDOMO
integer :: NSOMOI
integer :: NSOMOJ
integer :: ndiffSOMO
integer :: ndiffDOMO
integer :: starti, endi, cnti, cntj, rows,cols
integer :: extype,pmodel,qmodel
integer(bit_kind) :: Icfg(N_INT,2)
integer(bit_kind) :: Jcfg(N_INT,2)
integer,external :: getNSOMO
real*8, external :: mo_two_e_integral
real*8 :: hpp
real*8 :: meCC
real*8 :: ecore
! initialize energies
diag_energies = 0.d0
! calculate core energy
!call get_core_energy(ecore)
!diag_energies = ecore
! calculate the core energy
!print *,"Core energy=",ref_bitmask_energy
do i=1,N_configuration
Isomo = psi_configuration(1,1,i)
Idomo = psi_configuration(1,2,i)
Icfg(1,1) = psi_configuration(1,1,i)
Icfg(1,2) = psi_configuration(1,2,i)
NSOMOI = getNSOMO(psi_configuration(:,:,i))
starti = psi_config_data(i,1)
endi = psi_config_data(i,2)
! find out all pq holes possible
nholes = 0
! holes in SOMO
!do k = n_core_orb+1,n_core_orb + n_act_orb
do k = 1,mo_num
if(POPCNT(IAND(Isomo,IBSET(0_8,k-1))) .EQ. 1) then
nholes += 1
listholes(nholes) = k
holetype(nholes) = 1
endif
enddo
! holes in DOMO
!do k = n_core_orb+1,n_core_orb + n_act_orb
!do k = 1+n_core_inact_orb,n_core_orb+n_core_inact_act_orb
do k = 1,mo_num
if(POPCNT(IAND(Idomo,IBSET(0_8,k-1))) .EQ. 1) then
nholes += 1
listholes(nholes) = k
holetype(nholes) = 2
endif
enddo
! find vmos
listvmos = -1
vmotype = -1
nvmos = 0
!do k = n_core_orb+1,n_core_orb + n_act_orb
do k = 1,mo_num
!print *,i,IBSET(0,i-1),POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))), POPCNT(IAND(Idomo,(IBSET(0_8,i-1))))
if(POPCNT(IAND(Isomo,(IBSET(0_8,k-1)))) .EQ. 0 .AND. POPCNT(IAND(Idomo,(IBSET(0_8,k-1)))) .EQ. 0) then
nvmos += 1
listvmos(nvmos) = k
vmotype(nvmos) = 0
else if(POPCNT(IAND(Isomo,(IBSET(0_8,k-1)))) .EQ. 1 .AND. POPCNT(IAND(Idomo,(IBSET(0_8,k-1)))) .EQ. 0 ) then
nvmos += 1
listvmos(nvmos) = k
vmotype(nvmos) = 1
end if
enddo
!print *,"I=",i
!call debug_spindet(psi_configuration(1,1,i),N_int)
!call debug_spindet(psi_configuration(1,2,i),N_int)
do k=1,nholes
p = listholes(k)
noccp = holetype(k)
! Calculate one-electron
! and two-electron coulomb terms
do l=1,nholes
q = listholes(l)
noccq = holetype(l)
!print *,"--------------- K=",p," L=",q
! one-electron term
if(p.EQ.q) then
hpp = noccq * h_core_ri(p,q)!mo_one_e_integrals(q,q)
else
hpp = 0.d0
endif
do j=starti,endi
! coulomb term
! (pp,qq) = <pq|pq>
if(p.EQ.q) then
diag_energies(j) += hpp !+ 0.5d0 * (noccp * noccq * mo_two_e_integral(p,q,p,q))
!print *,"hpp=",hpp,"diga= ",diag_energies(j)
! else
! diag_energies(j) += ! 0.5d0 * noccp * noccq * mo_two_e_integral(p,q,p,q)
! print *,"diga= ",diag_energies(j)
endif
enddo
enddo
enddo
enddo
end subroutine calculate_preconditioner_cfg
subroutine calculate_sigma_vector_cfg_nst(psi_out, psi_in, n_st, sze, istart, iend, ishift, istep)
implicit none
use bitmasks
BEGIN_DOC
! Documentation for sigma-vector calculation
!
! Calculates the result of the
! application of the hamiltonian to the
! wavefunction in CFG basis once
! TODO : Things prepare outside this routine
! 1. Touch the providers for
! a. ApqIJ containers
! b. DET to CSF transformation matrices
! 2. DET to CSF transcormation
! 2. CSF to DET back transcormation
! returns : psi_coef_out_det :
END_DOC
integer,intent(in) :: sze, istart,iend, istep, ishift, n_st
real*8,intent(in):: psi_in(sze,n_st)
real*8,intent(out):: psi_out(sze,n_st)
integer(bit_kind) :: Icfg(N_INT,2)
integer :: i,j,k,l,p,q,noccp,noccq, ii, jj, m, n, idxI, kk, nocck,orbk
integer(bit_kind) :: alphas_Icfg(N_INT,2,sze)
integer(bit_kind) :: singlesI(N_INT,2,sze)
integer(bit_kind) :: connectedI_alpha(N_INT,2,sze)
integer :: idxs_singlesI(sze)
integer :: idxs_connectedI_alpha(sze)
integer(bit_kind) :: psi_configuration_out(N_INT,2,sze)
real*8 :: psi_coef_out(n_CSF)
logical :: psi_coef_out_init(n_CSF)
integer :: excitationIds_single(2,sze)
integer :: excitationTypes_single(sze)
integer :: excitationIds(2,sze)
integer :: excitationTypes(sze)
real*8 :: diagfactors(sze)
integer :: nholes
integer :: nvmos
integer :: listvmos(mo_num)
integer :: vmotype(mo_num) ! 1 -> VMO 2 -> SOMO
integer :: listholes(mo_num)
integer :: holetype(mo_num) ! 1-> SOMO 2->DOMO
integer :: Nalphas_Icfg, nconnectedI, rowsikpq, colsikpq, nsinglesI
integer :: extype,NSOMOalpha,NSOMOI,NSOMOJ,pmodel,qmodel
integer :: getNSOMO integer :: getNSOMO
real*8,dimension(:,:),allocatable :: tempBuffer integer :: totcolsTKI
real*8,dimension(:),allocatable :: tempCoeff integer :: rowsTKI
real*8 :: norm_det1, phasedet integer :: noccpp
norm_det1 = 0.d0 integer :: istart_cfg, iend_cfg
MS = elec_alpha_num - elec_beta_num integer*8 :: MS, Isomo, Idomo, Jsomo, Jdomo, Ialpha, Ibeta
! initialization integer :: moi, moj, mok, mol, starti, endi, startj, endj, cnti, cntj, cntk
psi_coef_config = 0.d0 real*8 :: norm_coef_cfg, fac2eints
DetToCSFTransformationMatrix(0,:,:) = 1.d0 real*8 :: norm_coef_det
do i = 2-iand(elec_alpha_num-elec_beta_num,1), NSOMOMax,2 real*8 :: meCC1, meCC2, diagfac
Isomo = IBSET(0_8, i) - 1_8 real*8,dimension(:,:,:),allocatable :: TKI
! rows = Ncsfs real*8,dimension(:,:),allocatable :: GIJpqrs
! cols = Ndets real*8,dimension(:,:,:),allocatable :: TKIGIJ
bfIcfg = max(1,nint((binom(i,(i+1)/2)-binom(i,((i+1)/2)+1)))) real*8, external :: mo_two_e_integral
ndetI = max(1,nint((binom(i,(i+1)/2)))) real*8, external :: get_two_e_integral
real*8 :: diag_energies(n_CSF)
call calculate_preconditioner_cfg(diag_energies)
allocate(tempBuffer(bfIcfg,ndetI)) MS = 0
call getCSFtoDETTransformationMatrix(Isomo, MS, NBFMax, maxDetDimPerBF, tempBuffer) norm_coef_cfg=0.d0
DetToCSFTransformationMatrix(i,1:bfIcfg,1:ndetI) = tempBuffer(1:bfIcfg,1:ndetI)
deallocate(tempBuffer) psi_out=0.d0
psi_coef_out_init = .False.
istart_cfg = psi_csf_to_config_data(istart)
iend_cfg = psi_csf_to_config_data(iend)
!!! Single Excitations !!!
do i=istart_cfg,iend_cfg
Icfg(1,1) = psi_configuration(1,1,i)
Icfg(1,2) = psi_configuration(1,2,i)
Isomo = Icfg(1,1)
Idomo = Icfg(1,2)
NSOMOI = getNSOMO(Icfg)
! find out all pq holes possible
nholes = 0
! holes in SOMO
! list_act
! list_core
! list_core_inact
! bitmasks
!do k = n_core_orb+1,n_core_orb + n_act_orb
do k = 1,mo_num
if(POPCNT(IAND(Isomo,IBSET(0_8,k-1))) .EQ. 1) then
nholes += 1
listholes(nholes) = k
holetype(nholes) = 1
endif
enddo
! holes in DOMO
!do k = n_core_orb+1,n_core_orb + n_act_orb
do k = 1,mo_num
if(POPCNT(IAND(Idomo,IBSET(0_8,k-1))) .EQ. 1) then
nholes += 1
listholes(nholes) = k
holetype(nholes) = 2
endif
enddo
! find vmos
listvmos = -1
vmotype = -1
nvmos = 0
!do k = n_core_orb+1,n_core_orb + n_act_orb
do k = 1,mo_num
!print *,i,IBSET(0,i-1),POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))), POPCNT(IAND(Idomo,(IBSET(0_8,i-1))))
if(POPCNT(IAND(Isomo,(IBSET(0_8,k-1)))) .EQ. 0 .AND. POPCNT(IAND(Idomo,(IBSET(0_8,k-1)))) .EQ. 0) then
nvmos += 1
listvmos(nvmos) = k
vmotype(nvmos) = 0
else if(POPCNT(IAND(Isomo,(IBSET(0_8,k-1)))) .EQ. 1 .AND. POPCNT(IAND(Idomo,(IBSET(0_8,k-1)))) .EQ. 0 ) then
nvmos += 1
listvmos(nvmos) = k
vmotype(nvmos) = 1
end if
enddo
! Icsf ids
starti = psi_config_data(i,1)
endi = psi_config_data(i,2)
NSOMOI = getNSOMO(Icfg)
call generate_all_singles_cfg_with_type(Icfg,singlesI,idxs_singlesI,excitationIds_single, &
excitationTypes_single,nsinglesI,N_int)
do j = 1,nsinglesI
idxI = idxs_singlesI(j)
NSOMOJ = getNSOMO(singlesI(:,:,j))
p = excitationIds_single(1,j)
q = excitationIds_single(2,j)
extype = excitationTypes_single(j)
! Off diagonal terms
call convertOrbIdsToModelSpaceIds(Icfg, singlesI(:,:,j), p, q, extype, pmodel, qmodel)
Jsomo = singlesI(1,1,j)
Jdomo = singlesI(1,2,j)
! Add the hole on J
if(POPCNT(IAND(Jsomo,IBSET(0_8,q-1))) .EQ. 1 .AND. POPCNT(IAND(Isomo,IBSET(0_8,q-1))) .EQ. 0) then
nholes += 1
listholes(nholes) = q
holetype(nholes) = 1
endif
if((POPCNT(IAND(Jdomo,IBSET(0_8,q-1))) .EQ. 1 .AND. POPCNT(IAND(Idomo,IBSET(0_8,q-1))) .EQ. 0) .AND. POPCNT(IAND(Isomo,IBSET(0_8,q-1))) .EQ. 0) then
nholes += 1
listholes(nholes) = q
holetype(nholes) = 2
endif
startj = psi_config_data(idxI,1)
endj = psi_config_data(idxI,2)
!!! One-electron contribution !!!
do kk = 1,n_st
cnti = 0
do ii = starti, endi
cnti += 1
cntj = 0
do jj = startj, endj
cntj += 1
meCC1 = AIJpqContainer(NSOMOI,NSOMOJ,extype,pmodel,qmodel,cnti,cntj)
psi_out(jj,kk) += meCC1 * psi_in(ii,kk) * h_core_ri(p,q)
psi_coef_out_init(jj) = .True.
enddo
enddo
enddo
! Undo setting in listholes
if(POPCNT(IAND(Jsomo,IBSET(0_8,q-1))) .EQ. 1 .AND. POPCNT(IAND(Isomo,IBSET(0_8,q-1))) .EQ. 0) then
nholes -= 1
endif
if((POPCNT(IAND(Jdomo,IBSET(0_8,q-1))) .EQ. 1 .AND. POPCNT(IAND(Idomo,IBSET(0_8,q-1))) .EQ. 0) .AND. POPCNT(IAND(Isomo,IBSET(0_8,q-1))) .EQ. 0) then
nholes -= 1
endif
enddo
enddo enddo
integer s, bfIcfg !!! Double Excitations !!!
integer countcsf
countcsf = 0
integer countdet
countdet = 0
integer idx
integer istate
istate = 1
phasedet = 1.0d0
do i = 1,N_configuration
startdet = psi_configuration_to_psi_det(1,i)
enddet = psi_configuration_to_psi_det(2,i)
ndetI = enddet-startdet+1
allocate(tempCoeff(ndetI)) ! Loop over all selected configurations
countdet = 1 do i = istart_cfg,iend_cfg
do j = startdet, enddet
idx = psi_configuration_to_psi_det_data(j)
Ialpha(:) = psi_det(:,1,idx)
Ibeta(:) = psi_det(:,2,idx)
call get_phase_qp_to_cfg(Ialpha, Ibeta, phasedet)
tempCoeff(countdet) = psi_coef(idx, istate)*phasedet
norm_det1 += tempCoeff(countdet)*tempCoeff(countdet)
countdet += 1
enddo
s = 0 Icfg(1,1) = psi_configuration(1,1,i)
do k=1,N_int Icfg(1,2) = psi_configuration(1,2,i)
if (psi_configuration(k,1,i) == 0_bit_kind) cycle starti = psi_config_data(i,1)
s = s + popcnt(psi_configuration(k,1,i)) endi = psi_config_data(i,2)
enddo
bfIcfg = max(1,nint((binom(s,(s+1)/2)-binom(s,((s+1)/2)+1))))
! perhaps blocking with CFGs of same seniority ! Returns all unique (checking the past) singly excited cfgs connected to I
! can be more efficient call obtain_associated_alphaI(i, Icfg, alphas_Icfg, Nalphas_Icfg)
allocate(tempBuffer(bfIcfg,ndetI)) ! TODO : remove doubly excited for return
tempBuffer = DetToCSFTransformationMatrix(s,:bfIcfg,:ndetI) ! Here we do 2x the loop. One to count for the size of the matrix, then we compute.
do k = 1,Nalphas_Icfg
! Now generate all singly excited with respect to a given alpha CFG
call obtain_connected_I_foralpha(i,alphas_Icfg(1,1,k),connectedI_alpha,idxs_connectedI_alpha,nconnectedI,excitationIds,excitationTypes,diagfactors)
call dgemm('N','N', bfIcfg, 1, ndetI, 1.d0, tempBuffer, size(tempBuffer,1), tempCoeff, size(tempCoeff,1), 0.d0, psi_coef_config(countcsf+1), size(psi_coef_config,1)) if(nconnectedI .EQ. 0) then
!call dgemv('N', NBFMax, maxDetDimPerBF, 1.d0, tempBuffer, size(tempBuffer,1), tempCoeff, 1, 0.d0, psi_coef_config(countcsf), 1) cycle
endif
totcolsTKI = 0
rowsTKI = -1
do j = 1,nconnectedI
NSOMOalpha = getNSOMO(alphas_Icfg(:,:,k))
NSOMOI = getNSOMO(connectedI_alpha(:,:,j))
p = excitationIds(1,j)
q = excitationIds(2,j)
extype = excitationTypes(j)
call convertOrbIdsToModelSpaceIds(alphas_Icfg(1,1,k), connectedI_alpha(1,1,j), p, q, extype, pmodel, qmodel)
! for E_pp E_rs and E_ppE_rr case
if(p.EQ.q) then
NSOMOalpha = NSOMOI
endif
rowsikpq = AIJpqMatrixDimsList(NSOMOalpha,NSOMOI,extype,pmodel,qmodel,1)
colsikpq = AIJpqMatrixDimsList(NSOMOalpha,NSOMOI,extype,pmodel,qmodel,2)
totcolsTKI += colsikpq
if(rowsTKI .LT. rowsikpq .AND. rowsTKI .NE. -1) then
print *,">",j,"Something is wrong in sigma-vector", rowsTKI, rowsikpq, "(p,q)=",pmodel,qmodel,"ex=",extype,"na=",NSOMOalpha," nI=",NSOMOI
!rowsTKI = rowsikpq
else
rowsTKI = rowsikpq
endif
enddo
deallocate(tempCoeff) allocate(TKI(rowsTKI,n_st,totcolsTKI)) ! coefficients of CSF
deallocate(tempBuffer) ! Initialize the inegral container
psi_config_data(i,1) = countcsf + 1 ! dims : (totcolsTKI, nconnectedI)
countcsf += bfIcfg allocate(GIJpqrs(totcolsTKI,nconnectedI)) ! gpqrs
psi_config_data(i,2) = countcsf allocate(TKIGIJ(rowsTKI,n_st,nconnectedI)) ! gpqrs
totcolsTKI = 0
do j = 1,nconnectedI
NSOMOalpha = getNSOMO(alphas_Icfg(:,:,k))
NSOMOI = getNSOMO(connectedI_alpha(:,:,j))
p = excitationIds(1,j)
q = excitationIds(2,j)
extype = excitationTypes(j)
call convertOrbIdsToModelSpaceIds(alphas_Icfg(:,:,k), connectedI_alpha(:,:,j), p, q, extype, pmodel, qmodel)
rowsikpq = AIJpqMatrixDimsList(NSOMOalpha,NSOMOI,extype,pmodel,qmodel,1)
colsikpq = AIJpqMatrixDimsList(NSOMOalpha,NSOMOI,extype,pmodel,qmodel,2)
do kk = 1,n_st
do l = 1,rowsTKI
do m = 1,colsikpq
TKI(l,kk,totcolsTKI+m) = AIJpqContainer(NSOMOalpha,NSOMOI,extype,pmodel,qmodel,l,m) * psi_in(idxs_connectedI_alpha(j)+m-1,kk)
enddo
enddo
enddo
do m = 1,colsikpq
do l = 1,nconnectedI
! <ij|kl> = (ik|jl)
moi = excitationIds(1,j) ! p
mok = excitationIds(2,j) ! q
moj = excitationIds(2,l) ! s
mol = excitationIds(1,l) ! r
if(moi.EQ.mok .AND. moj.EQ.mol)then
diagfac = diagfactors(j)
diagfac *= diagfactors(l)
!print *,"integrals (",totcolsTKI+m,l,")",mok,moi,mol,moj, "|", diagfac
GIJpqrs(totcolsTKI+m,l) = diagfac*0.5d0*mo_two_e_integral(mok,mol,moi,moj) ! g(pq,sr) = <ps,qr>
else
diagfac = diagfactors(j)*diagfactors(l)
!print *,"integrals (",totcolsTKI+m,l,")",mok,moi,mol,moj, "|", diagfac
GIJpqrs(totcolsTKI+m,l) = diagfac*0.5d0*mo_two_e_integral(mok,mol,moi,moj) ! g(pq,sr) = <ps,qr>
!endif
endif
enddo
enddo
totcolsTKI += colsikpq
enddo
! Do big BLAS
! TODO TKI, size(TKI,1)*size(TKI,2)
call dgemm('N','N', rowsTKI*n_st, nconnectedI, totcolsTKI, 1.d0, &
TKI, size(TKI,1)*n_st, GIJpqrs, size(GIJpqrs,1), 0.d0, &
TKIGIJ , size(TKIGIJ,1)*n_st )
! Collect the result
totcolsTKI = 0
do j = 1,nconnectedI
NSOMOalpha = getNSOMO(alphas_Icfg(:,:,k))
NSOMOI = getNSOMO(connectedI_alpha(:,:,j))
p = excitationIds(1,j)
q = excitationIds(2,j)
extype = excitationTypes(j)
call convertOrbIdsToModelSpaceIds(alphas_Icfg(:,:,k), connectedI_alpha(:,:,j), p, q, extype, pmodel, qmodel)
rowsikpq = AIJpqMatrixDimsList(NSOMOalpha,NSOMOI,extype,pmodel,qmodel,1)
colsikpq = AIJpqMatrixDimsList(NSOMOalpha,NSOMOI,extype,pmodel,qmodel,2)
!print *,">j=",j,rowsikpq,colsikpq, ">>",totcolsTKI,",",idxs_connectedI_alpha(j)
do kk = 1,n_st
do m = 1,colsikpq
do l = 1,rowsTKI
psi_out(idxs_connectedI_alpha(j)+m-1,kk) += AIJpqContainer(NSOMOalpha,NSOMOI,extype,pmodel,qmodel,l,m) * TKIGIJ(l,kk,j)
psi_coef_out_init(idxs_connectedI_alpha(j)+m-1) = .True.
enddo
enddo
enddo
totcolsTKI += colsikpq
enddo
deallocate(TKI) ! coefficients of CSF
! Initialize the inegral container
! dims : (totcolsTKI, nconnectedI)
deallocate(GIJpqrs) ! gpqrs
deallocate(TKIGIJ) ! gpqrs
enddo ! loop over alphas
enddo ! loop over I
! Add the diagonal contribution
do i = 1,n_CSF
psi_out(i,1) += 1.0d0*diag_energies(i)*psi_in(i,1)
enddo enddo
END_PROVIDER
end subroutine calculate_sigma_vector_cfg_nst

1
src/csf/tree_utils.h
View File

@ -23,6 +23,7 @@ struct bin_tree {
}; };
#include "/opt/intel/oneapi/mkl/2021.1.1/include/mkl_cblas.h" #include "/opt/intel/oneapi/mkl/2021.1.1/include/mkl_cblas.h"
//#include "cblas.h"
#define MAX_SOMO 32 #define MAX_SOMO 32

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

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