From f2bbbfc2998f050d08c27a73b1b6cb909c87e573 Mon Sep 17 00:00:00 2001 From: v1j4y Date: Mon, 1 Mar 2021 12:28:45 +0100 Subject: [PATCH] Added my configuration_CI_sigma_helpers.org --- src/csf/configuration_CI_sigma_helpers.org | 103 ++- .../configuration_CI_sigma_helpers.org | 657 ++++++++++++++++++ 2 files changed, 742 insertions(+), 18 deletions(-) create mode 100644 src/determinants/configuration_CI_sigma_helpers.org diff --git a/src/csf/configuration_CI_sigma_helpers.org b/src/csf/configuration_CI_sigma_helpers.org index ede42b97..6677139a 100644 --- a/src/csf/configuration_CI_sigma_helpers.org +++ b/src/csf/configuration_CI_sigma_helpers.org @@ -43,8 +43,10 @@ the input determinant \(|D_I\rangle\). integer*8 :: Jsomo integer*8 :: diffSOMO integer*8 :: diffDOMO + integer*8 :: xordiffSOMODOMO integer :: ndiffSOMO integer :: ndiffDOMO + integer :: nxordiffSOMODOMO integer :: ndiffAll integer :: i integer :: j @@ -55,6 +57,7 @@ the input determinant \(|D_I\rangle\). 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" @@ -66,7 +69,7 @@ the input determinant \(|D_I\rangle\). ! find out all pq holes possible nholes = 0 ! 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 nholes += 1 listholes(nholes) = i @@ -74,7 +77,7 @@ the input determinant \(|D_I\rangle\). endif end do ! 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 nholes += 1 listholes(nholes) = i @@ -86,7 +89,7 @@ the input determinant \(|D_I\rangle\). listvmos = -1 vmotype = -1 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)))) 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 @@ -126,7 +129,7 @@ the input determinant \(|D_I\rangle\). p = listholes(i) do j = 1,nvmos q = listvmos(j) - if(p == q) cycle + if(p .EQ. q) cycle if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 1) then ! SOMO -> VMO Jsomo = IBCLR(Isomo,p-1) @@ -158,9 +161,12 @@ the input determinant \(|D_I\rangle\). 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) - 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. !print *,i,k,ndiffSOMO,ndiffDOMO !call debug_spindet(Jsomo,1) @@ -171,11 +177,13 @@ the input determinant \(|D_I\rangle\). endif end do + !print *,"(,",p,",",q,")",pqAlreadyGenQ + if(pqAlreadyGenQ) cycle pqExistsQ = .FALSE. ! 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)) diffDOMO = IEOR(OR(reunion_of_act_virt_bitmask(1,1),Jdomo),psi_configuration(1,2,k)) ndiffSOMO = POPCNT(diffSOMO) @@ -236,8 +244,6 @@ the input determinant \(|D_I\rangle\). NalphaIcfg += 1 !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,2,NalphaIcfg) = IOR(Jdomo,ISHFT(1_8,n_core_orb)-1) endif @@ -287,6 +293,7 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI integer*8 :: diffDOMO integer :: ndiffSOMO integer :: ndiffDOMO + integer :: nxordiffSOMODOMO integer :: i,j,k,l,p,q,nsomoJ,nsomoalpha,starti,endi,extyp,nholes integer :: listholes(mo_num) integer :: holetype(mo_num) @@ -316,7 +323,7 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI p = 0 q = 0 - do i=idxI+1,N_configuration + do i=idxI,N_configuration Isomo = Ialpha(1,1) Idomo = Ialpha(1,2) Jsomo = psi_configuration(1,1,i) @@ -330,8 +337,7 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI diffDOMO = IEOR(Idomo,Jdomo) ndiffSOMO = POPCNT(diffSOMO) ndiffDOMO = POPCNT(diffDOMO) - if((ndiffSOMO + ndiffDOMO) .EQ. 0) cycle - !print *,"-I--i=",i,diffSOMO,diffDOMO!Isomo,Jsomo,ndiffSOMO,ndiffDOMO + nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO) !print *,POPCNT(IEOR(diffSOMO,diffDOMO)), ndiffDOMO if(POPCNT(IEOR(diffSOMO,diffDOMO)) .LE. 1 .AND. ndiffDOMO .LT. 3) then !call debug_spindet(Isomo,1) @@ -390,7 +396,58 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI excitationIds(1,nconnectedI)=p excitationIds(2,nconnectedI)=q 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 end do @@ -468,6 +525,7 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod qmodel = q if(p .EQ. q) then + !print *,"input pq=",p,q,"extype=",extype pmodel = 1 qmodel = 1 else @@ -501,12 +559,12 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod ! SOMO -> VMO !print *,"type -> SOMO -> VMO" !Isomo = IEOR(Isomo,Jsomo) - 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)) + 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) case (4) ! DOMO -> SOMO ! remove all domos except one at p @@ -518,6 +576,15 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod mask = ISHFT(1_8,q) - 1 Isomo = IAND(Isomo,mask) qmodel = POPCNT(mask) - POPCNT(XOR(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 print *,"something is wrong in convertOrbIdsToModelSpaceIds" end select diff --git a/src/determinants/configuration_CI_sigma_helpers.org b/src/determinants/configuration_CI_sigma_helpers.org new file mode 100644 index 00000000..88b27a7a --- /dev/null +++ b/src/determinants/configuration_CI_sigma_helpers.org @@ -0,0 +1,657 @@ +#+title: Configuration Sigma Vector Helpers +#+author: Vijay Gopal Chilkuri +#+email: vijay.gopal.c@gmail.com + +* Generate the singly excited configurations on-the-fly + +The algorithm is based on the work by Garniron et. al. (see thesis Chap 5). +The basic idea is to generate \(|\alpha\rangle\)'s on-the-fly. + +The algorithm is based on the idea of splitting the list of \(|\alpha\rangle\)'s +into blocks associated with a selected determinant \(|D_I\rangle\). + +** Create a function to return a list of alphas + +Here we create a list of \(|\alpha\rangle\)'s associated with +the input determinant \(|D_I\rangle\). + +#+begin_src f90 :main no :tangle configuration_CI_sigma_helpers.irp.f + 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 +#+end_src + +** Given an \(\alpha\) CFG, return all the \(|I\rangle\) CFGs + +Next step is to obtain the connected CFGs \(|I\rangle\) that belong to the selected space +given a RI configuration \(|\alpha\rangle\). + +#+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, diagfactors) + implicit none + use bitmasks + BEGIN_DOC + ! Documentation for obtain_connected_I_foralpha + ! This function returns all those selected configurations + ! which are connected to the input configuration + ! Ialpha by a single excitation. + ! + ! The type of excitations are ordered as follows: + ! Type 1 - SOMO -> SOMO + ! Type 2 - DOMO -> VMO + ! Type 3 - SOMO -> VMO + ! Type 4 - DOMO -> SOMO + ! + ! Order of operators + ! \alpha> = a^\dag_p a_q |I> = E_pq |I> + END_DOC + integer ,intent(in) :: idxI + integer(bit_kind),intent(in) :: Ialpha(N_int,2) + integer(bit_kind),intent(out) :: connectedI(N_int,2,*) + integer ,intent(out) :: idxs_connectedI(*) + integer,intent(out) :: nconnectedI + integer,intent(out) :: excitationIds(2,*) + integer,intent(out) :: excitationTypes(*) + real*8 ,intent(out) :: diagfactors(*) + integer*8 :: Idomo + integer*8 :: Isomo + integer*8 :: Jdomo + integer*8 :: Jsomo + integer*8 :: IJsomo + integer*8 :: diffSOMO + integer*8 :: diffDOMO + integer*8 :: xordiffSOMODOMO + integer :: ndiffSOMO + integer :: ndiffDOMO + integer :: nxordiffSOMODOMO + integer :: ii,i,j,k,l,p,q,nsomoJ,nsomoalpha,starti,endi,extyp,nholes + integer :: listholes(mo_num) + integer :: holetype(mo_num) + + nconnectedI = 0 + + p = 0 + q = 0 + do i=idxI,N_configuration + Isomo = Ialpha(1,1) + Idomo = Ialpha(1,2) + Jsomo = psi_configuration(1,1,i) + Jdomo = psi_configuration(1,2,i) + !call debug_spindet(Isomo,1) + !call debug_spindet(Idomo,1) + !print *,"-J--i=",i,Idomo,Jdomo,">",N_configuration + !call debug_spindet(Jsomo,1) + !call debug_spindet(Jdomo,1) + diffSOMO = IEOR(Isomo,Jsomo) + diffDOMO = IEOR(Idomo,Jdomo) + xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) + ndiffSOMO = POPCNT(diffSOMO) + ndiffDOMO = POPCNT(diffDOMO) + nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO) + !print *,"-I--i=",i,ndiffSOMO,ndiffDOMO,nxordiffSOMODOMO!Isomo,Jsomo,ndiffSOMO,ndiffDOMO + !if((ndiffSOMO + ndiffDOMO) .EQ. 0) cycle + !print *,POPCNT(IEOR(diffSOMO,diffDOMO)), ndiffDOMO + !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(Idomo,1) + !print *,"-J--i=",i,Idomo,Jdomo,">",N_configuration + !call debug_spindet(Jsomo,1) + !call debug_spindet(Jdomo,1) + select case(ndiffDOMO) + case (0) + ! SOMO -> VMO + !print *,"obt SOMO -> VMO" + extyp = 3 + IJsomo = IEOR(Isomo, Jsomo) + p = TRAILZ(IAND(Isomo,IJsomo)) + 1 + IJsomo = IBCLR(IJsomo,p-1) + q = TRAILZ(IJsomo) + 1 + 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 + p = TRAILZ(IEOR(Idomo,Jdomo)) + 1 + Isomo = IEOR(Isomo, Jsomo) + Isomo = IBCLR(Isomo,p-1) + q = TRAILZ(Isomo) + 1 + else + ! SOMO -> SOMO + !print *,"obt SOMO -> SOMO" + extyp = 1 + q = TRAILZ(IEOR(Idomo,Jdomo)) + 1 + Isomo = IEOR(Isomo, Jsomo) + Isomo = IBCLR(Isomo,q-1) + p = TRAILZ(Isomo) + 1 + end if + case (2) + ! DOMO -> SOMO + !print *,"obt DOMO -> SOMO" + extyp = 4 + IJsomo = IEOR(Isomo, Jsomo) + p = TRAILZ(IAND(Jsomo,IJsomo)) + 1 + IJsomo = IBCLR(IJsomo,p-1) + q = TRAILZ(IJsomo) + 1 + case default + print *,"something went wront in get connectedI" + end select + 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 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 + end do + +end subroutine obtain_connected_I_foralpha +#+end_src + +#+begin_src fortran + print *,TRAILZ(8) + print *,IBCLR(8,TRAILZ(9)) + print *,TRAILZ(IBCLR(8,TRAILZ(9))) + +#+end_src + +#+RESULTS: +| 3 | +| 8 | +| 3 | + +** Function to get the NSOMOs (seniority) + +#+begin_src f90 :main no :tangle configuration_CI_sigma_helpers.irp.f + function getNSOMO(Icfg) result(NSOMO) + implicit none + integer(bit_kind),intent(in) :: Icfg(N_int,2) + integer :: NSOMO + integer :: i + NSOMO = 0 + do i = 1,N_int + NSOMO += POPCNT(Icfg(i,1)) + enddo + end function getNSOMO +#+end_src + +** Function to convert p,q to model space ids + +This function converts the real orbital ids \(i,j\) to model +space ids \(p,q\) which depend only on the number of somos. + +#+begin_src f90 :main no :tangle configuration_CI_sigma_helpers.irp.f +subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmodel) + implicit none + BEGIN_DOC + ! This function converts the orbital ids + ! in real space to those used in model space + ! in order to identify the matrices required + ! for the calculation of MEs. + ! + ! The type of excitations are ordered as follows: + ! Type 1 - SOMO -> SOMO + ! Type 2 - DOMO -> VMO + ! Type 3 - SOMO -> VMO + ! Type 4 - DOMO -> SOMO + END_DOC + integer(bit_kind),intent(in) :: Ialpha(N_int,2) + integer(bit_kind),intent(in) :: Jcfg(N_int,2) + integer,intent(in) :: p,q + integer,intent(in) :: extype + integer,intent(out) :: pmodel,qmodel + integer*8 :: Isomo + integer*8 :: Idomo + integer*8 :: Jsomo + integer*8 :: Jdomo + integer*8 :: mask + integer*8 :: Isomotmp + integer*8 :: Jsomotmp + integer :: pos0,pos0prev + + ! 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) + pos0prev = 0 + pmodel = p + qmodel = q + + if(p .EQ. q) then + !print *,"input pq=",p,q,"extype=",extype + pmodel = 1 + qmodel = 1 + 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) + case (1) + ! 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)) + 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)) + 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 + 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) + ! DOMO -> SOMO + ! remove all domos except one at p + !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 + 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 + print *,"something is wrong in convertOrbIdsToModelSpaceIds" + end select + endif + !print *,p,q,"model ids=",pmodel,qmodel +end subroutine convertOrbIdsToModelSpaceIds +#+end_src + +#+begin_src fortran + integer :: i + integer :: count + integer :: mask + integer :: isomo + count = 0 + mask = ISHFT(1_8,5)-1 + print *,mask + print *,POPCNT(mask) + isomo = 144 + isomo = IAND(isomo,mask) + print *,isomo + print *,XOR(isomo,mask) + print *,POPCNT(mask) - POPCNT(XOR(isomo,mask)) + +#+end_src + +#+RESULTS: +| 31 | +| 5 | +| 16 | +| 15 | +| 1 | + +#+begin_src fortran + print *,IBSET(0_8,4)-1 + print *,POPCNT(IBSET(0_8,4)-1) - POPCNT(IAND(716,IBSET(0_8,4)-1)) + print *,POPCNT(IBSET(0_8,8)-1) - POPCNT(IAND(716,IBSET(0_8,8)-1)) +#+end_src + +#+RESULTS: +| 15 | +| 2 | +| 4 |