diff --git a/src/csf/configuration_CI_sigma_helpers.org b/src/csf/configuration_CI_sigma_helpers.org index 6677139a..88b27a7a 100644 --- a/src/csf/configuration_CI_sigma_helpers.org +++ b/src/csf/configuration_CI_sigma_helpers.org @@ -168,6 +168,12 @@ the input determinant \(|D_I\rangle\). !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) @@ -183,16 +189,16 @@ the input determinant \(|D_I\rangle\). 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 + !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. @@ -242,6 +248,7 @@ the input determinant \(|D_I\rangle\). 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 @@ -250,6 +257,33 @@ the input determinant \(|D_I\rangle\). 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 @@ -259,7 +293,7 @@ Next step is to obtain the connected CFGs \(|I\rangle\) that belong to the selec 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) +subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI, nconnectedI, excitationIds, excitationTypes, diagfactors) implicit none use bitmasks BEGIN_DOC @@ -284,6 +318,7 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, 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 @@ -291,34 +326,14 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI integer*8 :: IJsomo integer*8 :: diffSOMO integer*8 :: diffDOMO + integer*8 :: xordiffSOMODOMO integer :: ndiffSOMO integer :: ndiffDOMO integer :: nxordiffSOMODOMO - integer :: i,j,k,l,p,q,nsomoJ,nsomoalpha,starti,endi,extyp,nholes + integer :: ii,i,j,k,l,p,q,nsomoJ,nsomoalpha,starti,endi,extyp,nholes integer :: listholes(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 p = 0 @@ -335,11 +350,15 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI !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(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 @@ -565,17 +584,26 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod 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) - 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)) + 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 @@ -589,7 +617,7 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod print *,"something is wrong in convertOrbIdsToModelSpaceIds" end select endif - !print *,p,q,"model ids=",pmodel,qmodel + !print *,p,q,"model ids=",pmodel,qmodel end subroutine convertOrbIdsToModelSpaceIds #+end_src diff --git a/src/determinants/configuration_CI_sigma_helpers.org b/src/determinants/configuration_CI_sigma_helpers.org deleted file mode 100644 index 88b27a7a..00000000 --- a/src/determinants/configuration_CI_sigma_helpers.org +++ /dev/null @@ -1,657 +0,0 @@ -#+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 |