diff --git a/src/csf/configuration_CI_sigma_helpers.irp.f b/src/csf/configuration_CI_sigma_helpers.irp.f index f73362eb..eb532b03 100644 --- a/src/csf/configuration_CI_sigma_helpers.irp.f +++ b/src/csf/configuration_CI_sigma_helpers.irp.f @@ -1,3 +1,273 @@ + subroutine obtain_associated_alphaI(idxI, Icfg, alphasIcfg, NalphaIcfg, factor_alphaI) + implicit none + use bitmasks + BEGIN_DOC + ! Documentation for alphasI + ! Returns the associated alpha's for + ! the input configuration Icfg. + END_DOC + + integer,intent(in) :: idxI ! The id of the Ith CFG + integer(bit_kind),intent(in) :: Icfg(N_int,2) + integer,intent(out) :: NalphaIcfg + real*8 ,intent(out) :: factor_alphaI(*) + integer(bit_kind),intent(out) :: alphasIcfg(N_int,2,*) + logical,dimension(:,:),allocatable :: tableUniqueAlphas + integer :: listholes(mo_num) + integer :: holetype(mo_num) ! 1-> SOMO 2->DOMO + integer :: nholes + integer :: nvmos + integer :: listvmos(mo_num) + integer :: vmotype(mo_num) ! 1 -> VMO 2 -> SOMO + integer*8 :: Idomo + integer*8 :: Isomo + integer*8 :: Jdomo + integer*8 :: Jsomo + integer*8 :: diffSOMO + integer*8 :: diffDOMO + integer*8 :: xordiffSOMODOMO + integer :: ndiffSOMO + integer :: ndiffDOMO + integer :: nxordiffSOMODOMO + integer :: ndiffAll + integer :: i + integer :: j + integer :: k + integer :: hole + integer :: p + integer :: q + integer :: countalphas + logical :: pqAlreadyGenQ + logical :: pqExistsQ + logical :: ppExistsQ + Isomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,1)) + Idomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,2)) + !print*,"Input cfg" + !call debug_spindet(Isomo,1) + !call debug_spindet(Idomo,1) + + !print*,n_act_orb, "monum=",mo_num," n_core=",n_core_orb + + ! find out all pq holes possible + nholes = 0 + ! holes in SOMO + do i = 1,mo_num + if(POPCNT(IAND(Isomo,IBSET(0_8,i-1))) .EQ. 1) then + nholes += 1 + listholes(nholes) = i + holetype(nholes) = 1 + endif + end do + ! holes in DOMO + do i = 1,mo_num + if(POPCNT(IAND(Idomo,IBSET(0_8,i-1))) .EQ. 1) then + nholes += 1 + listholes(nholes) = i + holetype(nholes) = 2 + endif + end do + + ! find vmos + listvmos = -1 + vmotype = -1 + nvmos = 0 + do i = 1,mo_num + !print *,i,IBSET(0,i-1),POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))), POPCNT(IAND(Idomo,(IBSET(0_8,i-1)))) + if(POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))) .EQ. 0 .AND. POPCNT(IAND(Idomo,(IBSET(0_8,i-1)))) .EQ. 0) then + nvmos += 1 + listvmos(nvmos) = i + vmotype(nvmos) = 1 + else if(POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))) .EQ. 1 .AND. POPCNT(IAND(Idomo,(IBSET(0_8,i-1)))) .EQ. 0 ) then + nvmos += 1 + listvmos(nvmos) = i + vmotype(nvmos) = 2 + end if + end do + + !print *,"Nvmo=",nvmos + !print *,listvmos + !print *,vmotype + + allocate(tableUniqueAlphas(mo_num,mo_num)) + tableUniqueAlphas = .FALSE. + + ! Now find the allowed (p,q) excitations + Isomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,1)) + Idomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,2)) + !print *,"Isomo" + !call debug_spindet(Isomo,1) + !call debug_spindet(Idomo,1) + + !print *,"Nholes=",nholes," Nvmos=",nvmos, " idxi=",idxI + !do i = 1,nholes + ! print *,i,"->",listholes(i) + !enddo + !do i = 1,nvmos + ! print *,i,"->",listvmos(i) + !enddo + + ! TODO cfg_seniority_index + do i = 1,nholes + p = listholes(i) + do j = 1,nvmos + q = listvmos(j) + if(p .EQ. q) cycle + if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 1) then + ! SOMO -> VMO + Jsomo = IBCLR(Isomo,p-1) + Jsomo = IBSET(Jsomo,q-1) + Jdomo = Idomo + else if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 2) then + ! SOMO -> SOMO + Jsomo = IBCLR(Isomo,p-1) + Jsomo = IBCLR(Jsomo,q-1) + Jdomo = IBSET(Idomo,q-1) + else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then + ! DOMO -> VMO + Jsomo = IBSET(Isomo,p-1) + Jsomo = IBSET(Jsomo,q-1) + Jdomo = IBCLR(Idomo,p-1) + else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 2) then + ! DOMO -> SOMO + Jsomo = IBSET(Isomo,p-1) + Jsomo = IBCLR(Jsomo,q-1) + Jdomo = IBCLR(Idomo,p-1) + Jdomo = IBSET(Jdomo,q-1) + else + print*,"Something went wrong in obtain_associated_alphaI" + endif + + + pqAlreadyGenQ = .FALSE. + ! First check if it can be generated before + do k = 1, idxI-1 + diffSOMO = IEOR(Jsomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k))) + diffDOMO = IEOR(Jdomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k))) + xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) + ndiffSOMO = POPCNT(diffSOMO) + ndiffDOMO = POPCNT(diffDOMO) + nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO) + !if(POPCNT(IEOR(diffSOMO,diffDOMO)) .LE. 1 .AND. ndiffDOMO .LT. 3) then + if((ndiffSOMO+ndiffDOMO) .EQ. 0) then + pqAlreadyGenQ = .TRUE. + ppExistsQ = .TRUE. + EXIT + endif + if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then + pqAlreadyGenQ = .TRUE. + !ppExistsQ = .TRUE. + !print *,i,k,ndiffSOMO,ndiffDOMO + !call debug_spindet(Jsomo,1) + !call debug_spindet(Jdomo,1) + !call debug_spindet(iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k)),1) + !call debug_spindet(iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k)),1) + EXIT + endif + end do + + !print *,"(,",p,",",q,")",pqAlreadyGenQ + + if(pqAlreadyGenQ) cycle + + pqExistsQ = .FALSE. + ! now check if this exists in the selected list + !do k = idxI+1, N_configuration + ! diffSOMO = IEOR(OR(reunion_of_act_virt_bitmask(1,1),Jsomo),psi_configuration(1,1,k)) + ! diffDOMO = IEOR(OR(reunion_of_act_virt_bitmask(1,1),Jdomo),psi_configuration(1,2,k)) + ! ndiffSOMO = POPCNT(diffSOMO) + ! ndiffDOMO = POPCNT(diffDOMO) + ! if((ndiffSOMO + ndiffDOMO) .EQ. 0) then + ! pqExistsQ = .TRUE. + ! EXIT + ! endif + !end do + + if(.NOT. pqExistsQ) then + tableUniqueAlphas(p,q) = .TRUE. + !print *,p,q + !call debug_spindet(Jsomo,1) + !call debug_spindet(Jdomo,1) + endif + end do + end do + + !print *,tableUniqueAlphas(:,:) + + ! prune list of alphas + Isomo = Icfg(1,1) + Idomo = Icfg(1,2) + Jsomo = Icfg(1,1) + Jdomo = Icfg(1,2) + NalphaIcfg = 0 + do i = 1, nholes + p = listholes(i) + do j = 1, nvmos + q = listvmos(j) + if(p .EQ. q) cycle + if(tableUniqueAlphas(p,q)) then + if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 1) then + ! SOMO -> VMO + Jsomo = IBCLR(Isomo,p-1) + Jsomo = IBSET(Jsomo,q-1) + Jdomo = Idomo + else if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 2) then + ! SOMO -> SOMO + Jsomo = IBCLR(Isomo,p-1) + Jsomo = IBCLR(Jsomo,q-1) + Jdomo = IBSET(Idomo,q-1) + else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then + ! DOMO -> VMO + Jsomo = IBSET(Isomo,p-1) + Jsomo = IBSET(Jsomo,q-1) + Jdomo = IBCLR(Idomo,p-1) + else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 2) then + ! DOMO -> SOMO + Jsomo = IBSET(Isomo,p-1) + Jsomo = IBCLR(Jsomo,q-1) + Jdomo = IBCLR(Idomo,p-1) + Jdomo = IBSET(Jdomo,q-1) + else + print*,"Something went wrong in obtain_associated_alphaI" + endif + + ! SOMO + NalphaIcfg += 1 + !print *,p,q,"|",holetype(i),vmotype(j),NalphaIcfg + alphasIcfg(1,1,NalphaIcfg) = Jsomo + alphasIcfg(1,2,NalphaIcfg) = IOR(Jdomo,ISHFT(1_8,n_core_orb)-1) + endif + end do + end do + + ! Check if this Icfg has been previously generated as a mono + ppExistsQ = .False. + Isomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,1)) + Idomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,2)) + do k = 1, idxI-1 + diffSOMO = IEOR(Isomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k))) + diffDOMO = IEOR(Idomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k))) + xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) + ndiffSOMO = POPCNT(diffSOMO) + ndiffDOMO = POPCNT(diffDOMO) + nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO) + if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then + ppExistsQ = .TRUE. + EXIT + endif + end do + ! Diagonal part (pp,qq) + if(nholes > 0 .AND. (.NOT. ppExistsQ))then + ! SOMO + NalphaIcfg += 1 + !print *,p,q,"|",holetype(i),vmotype(j),NalphaIcfg + !call debug_spindet(Idomo,1) + !call debug_spindet(Jdomo,1) + alphasIcfg(1,1,NalphaIcfg) = Icfg(1,1) + alphasIcfg(1,2,NalphaIcfg) = Icfg(1,2) + endif + + end subroutine + function getNSOMO(Icfg) result(NSOMO) implicit none integer(bit_kind),intent(in) :: Icfg(N_int,2) @@ -47,6 +317,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 @@ -80,26 +351,44 @@ 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) + 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 + 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 + !print *,p,q,"model ids=",pmodel,qmodel end subroutine convertOrbIdsToModelSpaceIds diff --git a/src/csf/sigma_vector.irp.f b/src/csf/sigma_vector.irp.f index 9c1ac56a..e14f11bb 100644 --- a/src/csf/sigma_vector.irp.f +++ b/src/csf/sigma_vector.irp.f @@ -65,11 +65,6 @@ subroutine get_phase_qp_to_cfg(Ialpha, Ibeta, phaseout) integer :: nbetas integer :: count, k - ! Remove the DOMOs - mask2 = IAND(Ialpha,Ibeta) - deta = IEOR(Ialpha,mask2) - detb = IEOR(Ibeta ,mask2) - ! Find how many alpha electrons there are in all the N_ints integer :: Na(N_int) do k=1,N_int