diff --git a/src/csf/configuration_CI_sigma_helpers.irp.f b/src/csf/configuration_CI_sigma_helpers.irp.f index 581498c5..3794e8bb 100644 --- a/src/csf/configuration_CI_sigma_helpers.irp.f +++ b/src/csf/configuration_CI_sigma_helpers.irp.f @@ -671,11 +671,21 @@ use bitmasks Jdomo(jint) = IBSET(Jdomo(jint),jpos) endif + do ii=1, N_int + Jcfg(i,1) = Jsomo(i) + Jcfg(i,2) = Jdomo(i) + enddo + call bitstring_to_list(Jcfg,listall,nelall,N_int) Nsomo_J = nelall ! Check for Minimal alpha electrons (MS) if(Nsomo_J.lt.MS)then + if(vmotype(j) == 1)then + Jsomo(jint) = IBCLR(Jsomo(jint),jpos) + else if(vmotype(j) == 2)then + Jdomo(jint) = IBCLR(Jdomo(jint),jpos) + endif cycle endif @@ -715,12 +725,19 @@ use bitmasks endif end do - if(pqAlreadyGenQ) cycle + if(pqAlreadyGenQ) then + if(vmotype(j) == 1)then + Jsomo(jint) = IBCLR(Jsomo(jint),jpos) + else if(vmotype(j) == 2)then + Jdomo(jint) = IBCLR(Jdomo(jint),jpos) + endif + cycle + endif pqExistsQ = .FALSE. if(.NOT. pqExistsQ) then - tableUniqueAlphas(p,q) = .TRUE. + tableUniqueAlphas(pp,qq) = .TRUE. endif @@ -740,10 +757,13 @@ use bitmasks !print *,tableUniqueAlphas(:,:) ! prune list of alphas - Isomo = Icfg(1,1) - Idomo = Icfg(1,2) - Jsomo = Icfg(1,1) - Jdomo = Icfg(1,2) + do i=1, N_int + Isomo(i) = iand(act_bitmask(i,1),psi_configuration(i,1,idxI)) + Idomo(i) = iand(act_bitmask(i,1),psi_configuration(i,2,idxI)) + Jsomo(i) = Isomo(i) + Jdomo(i) = Idomo(i) + enddo + NalphaIcfg = 0 do i = 1, nholes pp = listholes(i) @@ -753,6 +773,7 @@ use bitmasks Jsomo(iint) = IBCLR(Jsomo(iint),ipos) else if(holetype(i) == 2)then Jdomo(iint) = IBCLR(Jdomo(iint),ipos) + Jsomo(iint) = IBSET(Jsomo(iint),ipos) endif do j = 1, nvmos @@ -763,14 +784,28 @@ use bitmasks Jsomo(jint) = IBSET(Jsomo(jint),jpos) else if(vmotype(j) == 2)then Jdomo(jint) = IBSET(Jdomo(jint),jpos) + Jsomo(jint) = IBCLR(Jsomo(jint),jpos) + endif + if(pp .EQ. qq) then + if(vmotype(j) == 1)then + Jsomo(jint) = IBCLR(Jsomo(jint),jpos) + else if(vmotype(j) == 2)then + Jdomo(jint) = IBCLR(Jdomo(jint),jpos) + endif + cycle endif - if(pp .EQ. qq) cycle if(tableUniqueAlphas(pp,qq)) then + do ii=1, N_int + Jcfg(i,1) = Jsomo(i) + Jcfg(i,2) = Jdomo(i) + enddo + call bitstring_to_list(Jcfg,listall,nelall,N_int) Nsomo_J = nelall if(Nsomo_J .ge. NSOMOMin) then + !print *," Idx = ",idxI, "p = ",pp, " q = ",qq," Jsomo=",Jsomo(1), " Jdomo=",IOR(Jdomo(1),ISHFT(1_8,n_core_orb)-1) NalphaIcfg += 1 alphasIcfg_list(:,1,idxI,NalphaIcfg) = Jcfg(:,1) if(n_core_orb .le. 63)then @@ -794,12 +829,14 @@ use bitmasks Jsomo(jint) = IBCLR(Jsomo(jint),jpos) else if(vmotype(j) == 2)then Jdomo(jint) = IBCLR(Jdomo(jint),jpos) + Jsomo(jint) = IBSET(Jsomo(jint),jpos) endif end do if(holetype(i) == 1)then Jsomo(iint) = IBSET(Jsomo(iint),ipos) else if(holetype(i) == 2)then Jdomo(iint) = IBSET(Jdomo(iint),ipos) + Jsomo(iint) = IBCLR(Jsomo(iint),ipos) endif end do @@ -863,6 +900,7 @@ END_PROVIDER integer,intent(in) :: idxI ! The id of the Ith CFG integer(bit_kind),intent(in) :: Icfg(N_int,2) + integer(bit_kind) :: Jcfg(N_int,2) integer,intent(out) :: NalphaIcfg integer(bit_kind),intent(out) :: alphasIcfg(N_int,2,*) logical,dimension(:,:),allocatable :: tableUniqueAlphas @@ -872,74 +910,84 @@ END_PROVIDER 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(bit_kind) :: Idomo(N_int), Idomop(N_int), Idomoq(N_int) + integer(bit_kind) :: Isomo(N_int), Isomop(N_int), Isomoq(N_int) + integer(bit_kind) :: Jdomo(N_int), Jdomop(N_int), Jdomoq(N_int) + integer(bit_kind) :: Jsomo(N_int), Jsomop(N_int), Jsomoq(N_int) + integer(bit_kind) :: diffDOMO(N_int), xordiffSOMODOMO(N_int), diffSOMO(N_int) integer :: ndiffSOMO integer :: ndiffDOMO integer :: nxordiffSOMODOMO integer :: ndiffAll integer :: i, ii integer :: j, jj + integer :: iii, iint, jint, ipos, jpos + integer :: i_s, i_d integer :: k, kk integer :: kstart integer :: kend - integer :: Nsomo_I - integer :: hole - integer :: p - integer :: q + integer :: Nsomo_I, Nsomo_J + integer :: hole, n_core_orb_64 + integer :: p, pp, p_s + integer :: q, qq, q_s integer :: countalphas logical :: pqAlreadyGenQ logical :: pqExistsQ logical :: ppExistsQ - Isomo = iand(act_bitmask(1,1),Icfg(1,1)) - Idomo = iand(act_bitmask(1,1),Icfg(1,2)) + integer :: listall(N_int*bit_kind_size), nelall + + do i=1, N_int + Isomo(i) = iand(act_bitmask(i,1),Icfg(i,1)) + Idomo(i) = iand(act_bitmask(i,1),Icfg(i,2)) + enddo + !print*,"Input cfg" !call debug_spindet(Isomo,1) !call debug_spindet(Idomo,1) ! find out all pq holes possible - nholes = 0 - ! holes in SOMO - do ii = 1,n_act_orb - i = list_act(ii) - 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 ii = 1,n_act_orb - i = list_act(ii) - if(POPCNT(IAND(Idomo,IBSET(0_8,i-1))) .EQ. 1) then - nholes += 1 - listholes(nholes) = i - holetype(nholes) = 2 - endif - end do + nholes = 0 + call bitstring_to_list(Isomo,listall,nelall,N_int) + + do iii=1,nelall + nholes += 1 + listholes(nholes) = listall(iii) + holetype(nholes) = 1 + end do + + Nsomo_I = nelall + + call bitstring_to_list(Idomo,listall,nelall,N_int) + + do iii=1,nelall + if(listall(iii) .gt. n_core_orb)then + nholes += 1 + listholes(nholes) = listall(iii) + holetype(nholes) = 2 + endif + end do + ! find vmos - listvmos = -1 - vmotype = -1 - nvmos = 0 - do ii = 1,n_act_orb - i = list_act(ii) - !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 + ! Take into account N_int + do ii = 1, n_act_orb + iii = list_act(ii) + iint = shiftr(iii-1,bit_kind_shift) + 1 + ipos = iii-shiftl((iint-1),bit_kind_shift)-1 + + if(IAND(Idomo(iint),(IBSET(0_8,ipos))) .EQ. 0) then + if(IAND(Isomo(iint),(IBSET(0_8,ipos))) .EQ. 0) then + nvmos += 1 + listvmos(nvmos) = iii + vmotype(nvmos) = 1 + else if(POPCNT(IAND(Isomo(iint),(IBSET(0_8,ipos)))) .EQ. 1) then + nvmos += 1 + listvmos(nvmos) = iii + vmotype(nvmos) = 2 + end if + end if + end do + !print *,"Nvmo=",nvmos !print *,listvmos @@ -948,10 +996,15 @@ END_PROVIDER allocate(tableUniqueAlphas(mo_num,mo_num)) tableUniqueAlphas = .FALSE. + ! Now find the allowed (p,q) excitations + do i=1, N_int + Isomo(i) = iand(act_bitmask(i,1),Icfg(i,1)) + Idomo(i) = iand(act_bitmask(i,1),Icfg(i,2)) + Jsomo(i) = Isomo(i) + Jdomo(i) = Idomo(i) + enddo + ! Now find the allowed (p,q) excitations - Isomo = iand(act_bitmask(1,1),Icfg(1,1)) - Idomo = iand(act_bitmask(1,1),Icfg(1,2)) - Nsomo_I = POPCNT(Isomo) if(Nsomo_I .EQ. 0) then kstart = 1 else @@ -971,41 +1024,40 @@ END_PROVIDER !enddo do i = 1,nholes - p = listholes(i) + pp = listholes(i) + iint = shiftr(pp-1,bit_kind_shift) + 1 + ipos = pp-shiftl((iint-1),bit_kind_shift)-1 + if(holetype(i) == 1)then + Jsomo(iint) = IBCLR(Jsomo(iint),ipos) + else if(holetype(i) == 2)then + Jdomo(iint) = IBCLR(Jdomo(iint),ipos) + endif + 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 - kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-2))) - kend = idxI-1 - 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) - kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-4))) - kend = idxI-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) - kstart = cfg_seniority_index(Nsomo_I) - kend = idxI-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) - kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-2))) - kend = idxI-1 - else - print*,"Something went wrong in obtain_associated_alphaI" + qq = listvmos(j) + jint = shiftr(qq-1,bit_kind_shift) + 1 + jpos = qq-shiftl((iint-1),bit_kind_shift)-1 + if(vmotype(j) == 1)then + Jsomo(jint) = IBSET(Jsomo(jint),jpos) + else if(vmotype(j) == 2)then + Jdomo(jint) = IBSET(Jdomo(jint),jpos) + endif + + do ii=1, N_int + Jcfg(i,1) = Jsomo(i) + Jcfg(i,2) = Jdomo(i) + enddo + + call bitstring_to_list(Jcfg,listall,nelall,N_int) + Nsomo_J = nelall + + if(pp .EQ. qq) then + if(vmotype(j) == 1)then + Jsomo(jint) = IBCLR(Jsomo(jint),jpos) + else if(vmotype(j) == 2)then + Jdomo(jint) = IBCLR(Jdomo(jint),jpos) + endif + cycle endif ! Again, we don't have to search from 1 @@ -1016,14 +1068,21 @@ END_PROVIDER pqAlreadyGenQ = .FALSE. ! First check if it can be generated before do k = kstart, kend - diffSOMO = IEOR(Jsomo,iand(act_bitmask(1,1),psi_configuration(1,1,k))) - ndiffSOMO = POPCNT(diffSOMO) + ndiffSOMO = 0 + ndiffDOMO = 0 + nxordiffSOMODOMO = 0 + do ii = 1, N_int + Jsomo = Jcfg(ii,1) + Jdomo = Jcfg(ii,2) + diffSOMO = IEOR(Jsomo,iand(reunion_of_act_virt_bitmask(ii,1),psi_configuration(ii,1,k))) + ndiffSOMO += POPCNT(diffSOMO(ii)) + diffDOMO = IEOR(Jdomo,iand(reunion_of_act_virt_bitmask(ii,2),psi_configuration(ii,2,k))) + xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) + ndiffDOMO += POPCNT(diffDOMO(ii)) + nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO(ii)) + nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO + end do if((ndiffSOMO .NE. 0) .AND. (ndiffSOMO .NE. 2)) cycle - diffDOMO = IEOR(Jdomo,iand(act_bitmask(1,1),psi_configuration(1,2,k))) - xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) - ndiffDOMO = POPCNT(diffDOMO) - nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO) - nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO !if(POPCNT(IEOR(diffSOMO,diffDOMO)) .LE. 1 .AND. ndiffDOMO .LT. 3) then if((ndiffSOMO+ndiffDOMO) .EQ. 0) then pqAlreadyGenQ = .TRUE. @@ -1033,19 +1092,20 @@ END_PROVIDER if((nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then pqAlreadyGenQ = .TRUE. !EXIT - !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 + if(pqAlreadyGenQ) then + if(vmotype(j) == 1)then + Jsomo(jint) = IBCLR(Jsomo(jint),jpos) + else if(vmotype(j) == 2)then + Jdomo(jint) = IBCLR(Jdomo(jint),jpos) + endif + cycle + endif pqExistsQ = .FALSE. ! now check if this exists in the selected list @@ -1066,53 +1126,67 @@ END_PROVIDER !call debug_spindet(Jsomo,1) !call debug_spindet(Jdomo,1) endif + if(vmotype(j) == 1)then + Jsomo(jint) = IBCLR(Jsomo(jint),jpos) + else if(vmotype(j) == 2)then + Jdomo(jint) = IBCLR(Jdomo(jint),jpos) + endif end do + if(holetype(i) == 1)then + Jsomo(iint) = IBSET(Jsomo(iint),ipos) + else if(holetype(i) == 2)then + Jdomo(iint) = IBSET(Jdomo(iint),ipos) + endif end do !print *,tableUniqueAlphas(:,:) ! prune list of alphas - Isomo = Icfg(1,1) - Idomo = Icfg(1,2) - Jsomo = Icfg(1,1) - Jdomo = Icfg(1,2) + do i=1, N_int + Isomo(i) = iand(act_bitmask(i,1),Icfg(i,1)) + Idomo(i) = iand(act_bitmask(i,1),Icfg(i,2)) + Jsomo(i) = Isomo(i) + Jdomo(i) = Idomo(i) + enddo + 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 + pp = listholes(i) + iint = shiftr(pp-1,bit_kind_shift) + 1 + ipos = pp-shiftl((iint-1),bit_kind_shift)-1 + if(holetype(i) == 1)then + Jsomo(iint) = IBCLR(Jsomo(iint),ipos) + else if(holetype(i) == 2)then + Jdomo(iint) = IBCLR(Jdomo(iint),ipos) + endif + do j = 1, nvmos + qq = listvmos(j) + jint = shiftr(qq-1,bit_kind_shift) + 1 + jpos = qq-shiftl((iint-1),bit_kind_shift)-1 + if(vmotype(j) == 1)then + Jsomo(jint) = IBSET(Jsomo(jint),jpos) + else if(vmotype(j) == 2)then + Jdomo(jint) = IBSET(Jdomo(jint),jpos) + endif + if(pp .EQ. qq) cycle + if(tableUniqueAlphas(pp,qq)) then ! SOMO NalphaIcfg += 1 - !print *,i,j,"|",NalphaIcfg - alphasIcfg(1,1,NalphaIcfg) = Jsomo - alphasIcfg(1,2,NalphaIcfg) = IOR(Jdomo,ISHFT(1_8,n_core_orb)-1) + alphasIcfg_list(:,1,idxI,NalphaIcfg) = Jcfg(:,1) + if(n_core_orb .le. 63)then + alphasIcfg_list(1,2,idxI,NalphaIcfg) = IOR(Jcfg(1,2),ISHFT(1_8,n_core_orb)-1) + else + n_core_orb_64 = n_core_orb + do ii=1,N_int + if(n_core_orb_64 .gt. 0)then + alphasIcfg_list(ii,2,idxI,NalphaIcfg) = IOR(Jcfg(ii,2),ISHFT(1_8,n_core_orb_64)-1) + else + alphasIcfg_list(ii,2,idxI,NalphaIcfg) = Jcfg(ii,2) + endif + n_core_orb_64 = ISHFT(n_core_orb_64,-6) + end do + endif !print *,"I = ",idxI, " Na=",NalphaIcfg," - ",Jsomo, IOR(Jdomo,ISHFT(1_8,n_core_orb)-1) endif end do @@ -1123,12 +1197,22 @@ END_PROVIDER Isomo = iand(act_bitmask(1,1),Icfg(1,1)) Idomo = iand(act_bitmask(1,1),Icfg(1,2)) do k = 1, idxI-1 - diffSOMO = IEOR(Isomo,iand(act_bitmask(1,1),psi_configuration(1,1,k))) - diffDOMO = IEOR(Idomo,iand(act_bitmask(1,1),psi_configuration(1,2,k))) - xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) - ndiffSOMO = POPCNT(diffSOMO) - ndiffDOMO = POPCNT(diffDOMO) - nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO) + do ii=1,N_int + diffSOMO = IEOR(Icfg(ii,1),iand(act_bitmask(ii,1),psi_configuration(ii,1,k))) + ndiffSOMO += POPCNT(diffSOMO(ii)) + end do + ! ndiffSOMO cannot be 0 (I /= k) + ! if ndiffSOMO /= 2 then it has to be greater than 2 and hense + ! this Icfg could not have been generated before. + if (ndiffSOMO /= 2) cycle + ndiffDOMO = 0 + nxordiffSOMODOMO = 0 + do ii=1,N_int + diffDOMO = IEOR(Icfg(ii,2),iand(act_bitmask(ii,1),psi_configuration(ii,2,k))) + xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) + ndiffDOMO += POPCNT(diffDOMO(ii)) + nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO(ii)) + end do if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then ppExistsQ = .TRUE. EXIT @@ -1141,8 +1225,8 @@ END_PROVIDER !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) + alphasIcfg_list(:,1,idxI,NalphaIcfg) = Icfg(:,1) + alphasIcfg_list(:,2,idxI,NalphaIcfg) = Icfg(:,2) endif end subroutine