diff --git a/src/csf/cfgCI_utils.c b/src/csf/cfgCI_utils.c index 3510db37..bad6434f 100644 --- a/src/csf/cfgCI_utils.c +++ b/src/csf/cfgCI_utils.c @@ -253,9 +253,9 @@ void generateAllBFs(int64_t Isomo, int64_t MS, Tree *bftree, int *NBF, int *NSOM buildTreeDriver(bftree, *NSOMO, MS, NBF); } -void ortho_qr_csf(double *overlapMatrix, int lda, double *orthoMatrix, int rows, int cols); +//void ortho_qr_csf(double *overlapMatrix, int lda, double *orthoMatrix, int rows, int cols); + -// QR to orthogonalize CSFs does not work //void gramSchmidt_qp(double *overlapMatrix, int rows, int cols, double *orthoMatrix){ // int i,j; // //for(j=0;j SOMO 2->DOMO +!!! integer :: nholes +!!! integer :: nvmos +!!! integer :: listvmos(mo_num) +!!! integer :: vmotype(mo_num) ! 1 -> VMO 2 -> SOMO +!!! integer*8 :: Idomo, Idomop, Idomoq +!!! integer*8 :: Isomo, Isomop, Isomoq +!!! integer*8 :: Jdomo, Jdomop, Jdomoq +!!! integer*8 :: Jsomo, Jsomop, Jsomoq +!!! integer*8 :: diffSOMO +!!! integer*8 :: diffDOMO +!!! integer*8 :: xordiffSOMODOMO +!!! integer :: ndiffSOMO +!!! integer :: ndiffDOMO +!!! integer :: nxordiffSOMODOMO +!!! integer :: ndiffAll +!!! integer :: i,ii,iii +!!! integer :: j,jj, i_s, i_d +!!! integer :: k,kk +!!! integer :: kstart +!!! integer :: kend +!!! 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 +!!! integer*8 :: MS +!!! integer :: listall(N_int*bit_kind_size), nelall +!!! +!!! double precision :: t0, t1 +!!! call wall_time(t0) +!!! +!!! MS = elec_alpha_num-elec_beta_num +!!! +!!! allocate(tableUniqueAlphas(mo_num,mo_num)) +!!! NalphaIcfg_list = 0 +!!! +!!! do idxI = 1, N_configuration +!!! +!!! Icfg = psi_configuration(:,:,idxI) +!!! Jcfg = psi_configuration(:,:,idxI) +!!! !print *," Jcfg somo=",Jcfg(1,1), " ", Jcfg(2,1) +!!! !print *," Jcfg domo=",Jcfg(1,2), " ", Jcfg(2,2) +!!! +!!! Isomo = iand(act_bitmask(1,1),Icfg(1,1)) +!!! Idomo = iand(act_bitmask(1,1),Icfg(1,2)) +!!! +!!! ! 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 +!!! call bitstring_to_list(psi_configuration(1,1,idxI),listall,nelall,N_int) +!!! +!!! !print *,'list somo' +!!! do iii=1,nelall +!!! nholes += 1 +!!! listholes(nholes) = listall(iii) +!!! !print *,listall(iii) +!!! holetype(nholes) = 1 +!!! end do +!!! +!!! Nsomo_I = nelall +!!! +!!! ! 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 +!!! +!!! !do iii=1,N_int +!!! ! print *,' iii=',iii, psi_configuration(iii,2,idxI), ' idxI=',idxI +!!! !end do +!!! call bitstring_to_list(psi_configuration(1,2,idxI),listall,nelall,N_int) +!!! +!!! !print *,'list domo ncore=',n_core_orb, ' nelall=',nelall +!!! do iii=1,nelall +!!! if(listall(iii) .gt. n_core_orb)then +!!! nholes += 1 +!!! listholes(nholes) = listall(iii) +!!! !print *,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) +!!! ! if(IAND(Idomo,(IBSET(0_8,i-1))) .EQ. 0) then +!!! ! if(IAND(Isomo,(IBSET(0_8,i-1))) .EQ. 0) then +!!! ! nvmos += 1 +!!! ! listvmos(nvmos) = i +!!! ! print *,'1 i=',i +!!! ! vmotype(nvmos) = 1 +!!! ! else if(POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))) .EQ. 1) then +!!! ! nvmos += 1 +!!! ! listvmos(nvmos) = i +!!! ! print *,'2 i=',i +!!! ! vmotype(nvmos) = 2 +!!! ! end if +!!! ! end if +!!! !end do +!!! !print *,'-----------' +!!! +!!! ! Take into account N_int +!!! do ii = 1, n_act_orb +!!! iii = list_act(ii) +!!! i_s = (1+((iii-1)/63)) +!!! i = iii - ( i_s -1 )*63 +!!! Isomo = iand(act_bitmask(i_s,1),Icfg(i_s,1)) +!!! Idomo = iand(act_bitmask(i_s,1),Icfg(i_s,2)) +!!! +!!! if(IAND(Idomo,(IBSET(0_8,i-1))) .EQ. 0) then +!!! if(IAND(Isomo,(IBSET(0_8,i-1))) .EQ. 0) then +!!! nvmos += 1 +!!! listvmos(nvmos) = iii +!!! vmotype(nvmos) = 1 +!!! else if(POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))) .EQ. 1) then +!!! nvmos += 1 +!!! listvmos(nvmos) = iii +!!! vmotype(nvmos) = 2 +!!! end if +!!! end if +!!! end do +!!! +!!! tableUniqueAlphas = .FALSE. +!!! +!!! ! 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 +!!! kstart = cfg_seniority_index(max(NSOMOMin,Nsomo_I-2)) +!!! endif +!!! kend = idxI-1 +!!! +!!! do i = 1,nholes +!!! pp = listholes(i) +!!! p_s = (1+((pp-1)/63)) +!!! p = pp - (p_s - 1)*63 +!!! !print *,' pp=',pp, ' p_s=',p_s, ' p=',p +!!! do j = 1,nvmos +!!! qq = listvmos(j) +!!! q_s = (1+((qq-1)/63)) +!!! q = qq - (q_s - 1)*63 +!!! !print *,' qq=',qq, ' q_s=',q_s, ' q=',q +!!! Isomop = iand(act_bitmask(i_s,1),Icfg(p_s,1)) +!!! Idomop = iand(act_bitmask(i_s,1),Icfg(p_s,2)) +!!! Isomop = iand(act_bitmask(i_s,1),Icfg(q_s,1)) +!!! Idomop = iand(act_bitmask(i_s,1),Icfg(q_s,2)) +!!! if(p .EQ. q) cycle +!!! if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 1) then +!!! ! SOMO -> VMO +!!! !print *,'SOMO -> VMO' +!!! if (p_s .eq. q_s) then +!!! Jsomop = IBCLR(Isomop,p-1) +!!! Jsomop = IBSET(Jsomop,q-1) +!!! Jsomoq = Jsomop +!!! else +!!! Jsomop = IBCLR(Isomop,p-1) +!!! Jsomoq = IBSET(Isomoq,q-1) +!!! endif +!!! +!!! ! Domo remains the same +!!! Jdomop = Idomop +!!! Jdomoq = Idomoq +!!! +!!! 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 +!!! !print *,'SOMO -> SOMO' +!!! if(p_s .eq. q_s) then +!!! Jsomop = IBCLR(Isomop,p-1) +!!! Jsomop = IBCLR(Jsomop,q-1) +!!! Jsomoq = Jsomop +!!! else +!!! Jsomop = IBCLR(Isomop,p-1) +!!! Jsomoq = IBCLR(Isomoq,q-1) +!!! endif +!!! +!!! Jdomoq = IBSET(Idomoq,q-1) +!!! +!!! ! Check for Minimal alpha electrons (MS) +!!! if(POPCNT(Jsomoq).ge.MS)then +!!! kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-4))) +!!! kend = idxI-1 +!!! else +!!! cycle +!!! endif +!!! else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then +!!! ! DOMO -> VMO +!!! !print *,'DOMO -> VMO', Isomop, p, q, Jsomop +!!! if(p_s .eq. q_s) then +!!! Jsomop = IBSET(Isomop,p-1) +!!! Jsomop = IBSET(Jsomop,q-1) +!!! Jsomoq = Jsomop +!!! else +!!! Jsomop = IBSET(Isomop,p-1) +!!! Jsomoq = IBSET(Jsomoq,q-1) +!!! endif +!!! !print *, 'Jsomop=', Jsomop +!!! +!!! Jdomop = IBCLR(Idomop,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 +!!! !print *,'DOMO -> SOMO' +!!! if(p_s .eq. q_s) then +!!! Jsomop = IBSET(Isomop,p-1) +!!! Jsomop = IBCLR(Jsomop,q-1) +!!! Jsomoq = Jsomop +!!! +!!! Jdomop = IBCLR(Idomop,p-1) +!!! Jdomop = IBSET(Jdomop,q-1) +!!! Jdomoq = Jdomop +!!! else +!!! Jsomop = IBSET(Isomop,p-1) +!!! Jsomoq = IBCLR(Jsomoq,q-1) +!!! +!!! Jdomop = IBCLR(Idomop,p-1) +!!! Jdomoq = IBSET(Jdomoq,q-1) +!!! endif +!!! +!!! kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-2))) +!!! kend = idxI-1 +!!! else +!!! print*,"Something went wrong in obtain_associated_alphaI" +!!! endif +!!! +!!! ! Save it to Jcfg +!!! !print *,i,j,"0| nalpha=",NalphaIcfg, " somo=",Jcfg(1,1),Jcfg(2,1) +!!! Jcfg(p_s,1) = Jsomop +!!! Jcfg(q_s,1) = Jsomoq +!!! Jcfg(p_s,2) = Jdomop +!!! Jcfg(q_s,2) = Jdomoq +!!! !print *,'p_s=',p_s,' q_s=', q_s +!!! !print *,'Jsomop=',Jsomop, ' Jsomoq=', Jsomoq, ' Jdomop=', Jdomop, ' Jdomoq=', Jdomo +!!! !print *,i,j,"1| nalpha=",NalphaIcfg, " somo=",Jcfg(1,1),Jcfg(2,1) +!!! call bitstring_to_list(Jcfg(1,1),listall,nelall,N_int) +!!! Nsomo_J = nelall +!!! +!!! ! Check for Minimal alpha electrons (MS) +!!! if(Nsomo_J.lt.MS)then +!!! cycle +!!! endif +!!! +!!! ! Again, we don't have to search from 1 +!!! ! we just use seniority to find the +!!! ! first index with NSOMO - 2 to NSOMO + 2 +!!! ! this is what is done in kstart, kend +!!! +!!! pqAlreadyGenQ = .FALSE. +!!! ! First check if it can be generated before +!!! do k = kstart, kend +!!! !diffSOMO = IEOR(Jsomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k))) +!!! !ndiffSOMO = POPCNT(diffSOMO) +!!! !if((ndiffSOMO .NE. 0) .AND. (ndiffSOMO .NE. 2)) cycle +!!! !diffDOMO = IEOR(Jdomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k))) +!!! !xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) +!!! !ndiffDOMO = POPCNT(diffDOMO) +!!! !nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO) +!!! !nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO +!!! +!!! 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) +!!! diffDOMO = IEOR(Jdomo,iand(reunion_of_act_virt_bitmask(ii,2),psi_configuration(ii,2,k))) +!!! xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) +!!! ndiffDOMO += POPCNT(diffDOMO) +!!! nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO) +!!! nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO +!!! end do +!!! +!!! if((ndiffSOMO .ne. 0) .and. (ndiffSOMO .ne. 2)) cycle +!!! +!!! if((ndiffSOMO+ndiffDOMO) .EQ. 0) then +!!! pqAlreadyGenQ = .TRUE. +!!! ppExistsQ = .TRUE. +!!! EXIT +!!! endif +!!! if((nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then +!!! pqAlreadyGenQ = .TRUE. +!!! EXIT +!!! endif +!!! end do +!!! +!!! if(pqAlreadyGenQ) cycle +!!! +!!! pqExistsQ = .FALSE. +!!! +!!! if(.NOT. pqExistsQ) then +!!! tableUniqueAlphas(p,q) = .TRUE. +!!! 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) +!!! pp = listholes(i) +!!! p_s = (1+((pp-1)/63)) +!!! p = pp - (p_s - 1)*63 +!!! do j = 1, nvmos +!!! !q = listvmos(j) +!!! qq = listvmos(j) +!!! q_s = (1+((qq-1)/63)) +!!! q = qq - (q_s - 1)*63 +!!! Isomop = iand(act_bitmask(i_s,1),Icfg(p_s,1)) +!!! Idomop = iand(act_bitmask(i_s,1),Icfg(p_s,2)) +!!! Isomoq = iand(act_bitmask(i_s,1),Icfg(q_s,1)) +!!! Idomoq = iand(act_bitmask(i_s,1),Icfg(q_s,2)) +!!! 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 +!!! if (p_s .eq. q_s) then +!!! Jsomop = IBCLR(Isomop,p-1) +!!! Jsomop = IBSET(Jsomop,q-1) +!!! Jsomoq = Jsomop +!!! else +!!! Jsomop = IBCLR(Isomop,p-1) +!!! Jsomoq = IBSET(Isomoq,q-1) +!!! endif +!!! +!!! ! Domo remains the same +!!! Jdomop = Idomop +!!! Jdomoq = Idomoq +!!! +!!! 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) +!!! +!!! if(p_s .eq. q_s) then +!!! Jsomop = IBCLR(Isomop,p-1) +!!! Jsomop = IBCLR(Jsomop,q-1) +!!! Jsomoq = Jsomop +!!! else +!!! Jsomop = IBCLR(Isomop,p-1) +!!! Jsomoq = IBCLR(Isomoq,q-1) +!!! endif +!!! +!!! Jdomoq = IBSET(Idomoq,q-1) +!!! +!!! if(POPCNT(Jsomoq).ge.MS)then +!!! kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-4))) +!!! kend = idxI-1 +!!! else +!!! cycle +!!! endif +!!! 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) +!!! +!!! if(p_s .eq. q_s) then +!!! Jsomop = IBSET(Isomop,p-1) +!!! Jsomop = IBSET(Jsomop,q-1) +!!! Jsomoq = Jsomop +!!! else +!!! Jsomop = IBSET(Isomop,p-1) +!!! Jsomoq = IBSET(Jsomoq,q-1) +!!! endif +!!! +!!! Jdomop = IBCLR(Idomop,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) +!!! if(p_s .eq. q_s) then +!!! Jsomop = IBSET(Isomop,p-1) +!!! Jsomop = IBCLR(Jsomop,q-1) +!!! Jsomoq = Jsomop +!!! +!!! Jdomop = IBCLR(Idomop,p-1) +!!! Jdomop = IBSET(Jdomop,q-1) +!!! Jdomoq = Jdomop +!!! else +!!! Jsomop = IBSET(Isomop,p-1) +!!! Jsomoq = IBCLR(Jsomoq,q-1) +!!! +!!! Jdomop = IBCLR(Idomop,p-1) +!!! Jdomoq = IBSET(Jdomoq,q-1) +!!! endif +!!! +!!! else +!!! print*,"Something went wrong in obtain_associated_alphaI" +!!! endif +!!! +!!! ! Save it to Jcfg +!!! Jcfg(p_s,1) = Jsomop +!!! Jcfg(q_s,1) = Jsomoq +!!! Jcfg(p_s,2) = Jdomop +!!! Jcfg(q_s,2) = Jdomoq +!!! +!!! ! SOMO +!!! !print *,i,j,"|",NalphaIcfg, Jsomo, IOR(Jdomo,ISHFT(1_8,n_core_orb)-1) +!!! if(POPCNT(Jsomo) .ge. NSOMOMin) then +!!! NalphaIcfg += 1 +!!! alphasIcfg_list(:,1,idxI,NalphaIcfg) = Jcfg(:,1) +!!! !alphasIcfg_list(:,2,idxI,NalphaIcfg) = IOR(Jdomo,ISHFT(1_8,n_core_orb)-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 +!!! NalphaIcfg_list(idxI) = NalphaIcfg +!!! !print *,i,j,"2| nalpha=",NalphaIcfg, " somo=",Jcfg(1,1),Jcfg(2,1) +!!! endif +!!! 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)) +!!! kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-2))) +!!! ndiffDOMO = 0 +!!! do k = kstart, idxI-1 +!!! do ii=1,N_int +!!! diffSOMO = IEOR(Icfg(ii,1),iand(act_bitmask(ii,1),psi_configuration(ii,1,k))) +!!! ndiffSOMO += POPCNT(diffSOMO) +!!! 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) +!!! nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO) +!!! end do +!!! if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4)) then +!!! ppExistsQ = .TRUE. +!!! EXIT +!!! endif +!!! end do +!!! ! Diagonal part (pp,qq) +!!! if(nholes > 0 .AND. (.NOT. ppExistsQ))then +!!! ! SOMO +!!! if(POPCNT(Jsomo) .ge. NSOMOMin) then +!!! NalphaIcfg += 1 +!!! alphasIcfg_list(:,1,idxI,NalphaIcfg) = Icfg(:,1) +!!! alphasIcfg_list(:,2,idxI,NalphaIcfg) = Icfg(:,2) +!!! NalphaIcfg_list(idxI) = NalphaIcfg +!!! endif +!!! endif +!!! +!!! NalphaIcfg = 0 +!!! enddo ! end loop idxI +!!! call wall_time(t1) +!!! print *, 'Preparation : ', t1 - t0 +!!! +!!!END_PROVIDER + + BEGIN_PROVIDER [ integer(bit_kind), alphasIcfg_list , (N_int,2,N_configuration,mo_num*mo_num)] &BEGIN_PROVIDER [ integer, NalphaIcfg_list, (N_configuration) ] implicit none !use bitmasks @@ -12,6 +536,7 @@ use bitmasks integer :: idxI ! The id of the Ith CFG integer(bit_kind) :: Icfg(N_int,2) + integer(bit_kind) :: Jcfg(N_int,2) integer :: NalphaIcfg logical,dimension(:,:),allocatable :: tableUniqueAlphas integer :: listholes(mo_num) @@ -20,31 +545,32 @@ use bitmasks 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(bit_kind) :: diffDOMO, xordiffSOMODOMO, diffSOMO integer :: ndiffSOMO integer :: ndiffDOMO integer :: nxordiffSOMODOMO integer :: ndiffAll - integer :: i,ii - integer :: j,jj + integer :: i,ii,iii, iint, jint, ipos, jpos + integer :: j,jj, 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 + integer :: countelec logical :: pqAlreadyGenQ logical :: pqExistsQ logical :: ppExistsQ integer*8 :: MS + integer :: listall(N_int*bit_kind_size), nelall double precision :: t0, t1 call wall_time(t0) @@ -57,45 +583,57 @@ use bitmasks do idxI = 1, N_configuration Icfg = psi_configuration(:,:,idxI) + Jcfg = psi_configuration(:,:,idxI) - Isomo = iand(act_bitmask(1,1),Icfg(1,1)) - Idomo = iand(act_bitmask(1,1),Icfg(1,2)) + !print *,"idxI=",idxI + do ii=1, N_int + Isomo(ii) = iand(act_bitmask(ii,1),psi_configuration(ii,1,idxI)) + Idomo(ii) = iand(act_bitmask(ii,2),psi_configuration(ii,2,idxI)) + !print *,Isomo(ii), Idomo(ii) + enddo ! 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 + listholes=-1 + 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) - if(IAND(Idomo,(IBSET(0_8,i-1))) .EQ. 0) then - if(IAND(Isomo,(IBSET(0_8,i-1))) .EQ. 0) then + + ! 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) = i + listvmos(nvmos) = iii vmotype(nvmos) = 1 - else if(POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))) .EQ. 1) then + else if(POPCNT(IAND(Isomo(iint),(IBSET(0_8,ipos)))) .EQ. 1) then nvmos += 1 - listvmos(nvmos) = i + listvmos(nvmos) = iii vmotype(nvmos) = 2 end if end if @@ -104,60 +642,63 @@ use bitmasks tableUniqueAlphas = .FALSE. ! 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) + do ii=1, N_int + !Isomo(ii) = iand(reunion_of_act_virt_bitmask(ii,1),psi_configuration(ii,1,idxI)) + !Idomo(ii) = iand(reunion_of_act_virt_bitmask(ii,2),psi_configuration(ii,2,idxI)) + Isomo(ii) = iand(act_bitmask(ii,1),psi_configuration(ii,1,idxI)) + Idomo(ii) = iand(act_bitmask(ii,2),psi_configuration(ii,2,idxI)) + !Isomo(ii) = psi_configuration(ii,1,idxI) + !Idomo(ii) = psi_configuration(ii,2,idxI) + Jsomo(ii) = Isomo(ii) + Jdomo(ii) = Idomo(ii) + enddo + if(Nsomo_I .EQ. 0) then kstart = 1 else kstart = cfg_seniority_index(max(NSOMOMin,Nsomo_I-2)) endif + kstart = 1 kend = idxI-1 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) + Jsomo(iint) = IBSET(Jsomo(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) - ! Check for Minimal alpha electrons (MS) - if(POPCNT(Jsomo).ge.MS)then - kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-4))) - kend = idxI-1 - else - cycle - endif - 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) + if(pp.eq.qq) cycle + jint = shiftr(qq-1,bit_kind_shift) + 1 + jpos = qq-shiftl((jint-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) + Jsomo(jint) = IBCLR(Jsomo(jint),jpos) endif + + Nsomo_J=0 + do ii=1, N_int + Jcfg(ii,1) = Jsomo(ii) + Jcfg(ii,2) = Jdomo(ii) + Nsomo_J += POPCNT(Jsomo(ii)) + enddo + ! Check for Minimal alpha electrons (MS) - if(POPCNT(Jsomo).lt.MS)then + 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) + Jsomo(jint) = IBSET(Jsomo(jint),jpos) + endif cycle endif @@ -169,15 +710,21 @@ use bitmasks pqAlreadyGenQ = .FALSE. ! First check if it can be generated before do k = kstart, kend - diffSOMO = IEOR(Jsomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k))) - ndiffSOMO = POPCNT(diffSOMO) - if((ndiffSOMO .NE. 0) .AND. (ndiffSOMO .NE. 2)) cycle - diffDOMO = IEOR(Jdomo,iand(reunion_of_act_virt_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 + ndiffSOMO = 0 + ndiffDOMO = 0 + nxordiffSOMODOMO = 0 + do ii = 1, N_int + diffSOMO = IEOR(Jcfg(ii,1),iand(act_bitmask(ii,1),psi_configuration(ii,1,k))) + ndiffSOMO += POPCNT(diffSOMO) + diffDOMO = IEOR(Jcfg(ii,2),iand(act_bitmask(ii,2),psi_configuration(ii,2,k))) + xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) + ndiffDOMO += POPCNT(diffDOMO) + nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO) + nxordiffSOMODOMO += POPCNT(diffSOMO) + POPCNT(diffDOMO) + end do + + if((ndiffSOMO .ne. 0) .and. (ndiffSOMO .ne. 2)) cycle + if((ndiffSOMO+ndiffDOMO) .EQ. 0) then pqAlreadyGenQ = .TRUE. ppExistsQ = .TRUE. @@ -189,86 +736,166 @@ 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) + Jsomo(jint) = IBSET(Jsomo(jint),jpos) + endif + cycle + endif pqExistsQ = .FALSE. + !print *, " ndiffSOMO=",ndiffSOMO, " ndiffDOMO=", ndiffDOMO, " nxordiffSOMODOMO=",nxordiffSOMODOMO, " p=",pp," q=",qq + if(.NOT. pqExistsQ) then - tableUniqueAlphas(p,q) = .TRUE. + tableUniqueAlphas(pp,qq) = .TRUE. + endif + + + if(vmotype(j) == 1)then + 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 !print *,tableUniqueAlphas(:,:) - ! prune list of alphas - Isomo = Icfg(1,1) - Idomo = Icfg(1,2) - Jsomo = Icfg(1,1) - Jdomo = Icfg(1,2) + do ii=1, N_int + Isomo(ii) = iand(act_bitmask(ii,1),psi_configuration(ii,1,idxI)) + Idomo(ii) = iand(act_bitmask(ii,2),psi_configuration(ii,2,idxI)) + !Isomo(ii) = psi_configuration(ii,1,idxI) + !Idomo(ii) = psi_configuration(ii,2,idxI) + Jsomo(ii) = Isomo(ii) + Jdomo(ii) = Idomo(ii) + enddo + !print *, " Isomo=",Isomo(1), " Idomo=", Idomo(1) + + !countelec=0 + !do ii=1, N_int + ! countelec += POPCNT(Icfg(ii,1))*1 + POPCNT(Icfg(ii,2))*2 + !enddo + !if(countelec .ne. 14)then + ! print *," idxI=",idxI, "00countelec=",countelec, " bit_kind_size=",bit_kind_size, " nvmo=",nvmos," mo_num=",mo_num + ! stop + !endif + 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) - if(POPCNT(Jsomo).ge.MS)then - kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-4))) - kend = idxI-1 - else - cycle - endif - 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) + Jsomo(iint) = IBSET(Jsomo(iint),ipos) + endif - ! SOMO - !print *,i,j,"|",NalphaIcfg, Jsomo, IOR(Jdomo,ISHFT(1_8,n_core_orb)-1) - if(POPCNT(Jsomo) .ge. NSOMOMin) then + do j = 1, nvmos + qq = listvmos(j) + if(pp.eq.qq) cycle + jint = shiftr(qq-1,bit_kind_shift) + 1 + jpos = qq-shiftl((jint-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) + Jsomo(jint) = IBCLR(Jsomo(jint),jpos) + endif + + if(tableUniqueAlphas(pp,qq)) then + + Nsomo_J = 0 + countelec = 0 + do ii=1, N_int + Jcfg(ii,1) = Jsomo(ii) + Jcfg(ii,2) = Jdomo(ii) + Nsomo_J += POPCNT(Jsomo(ii)) + countelec += POPCNT(Jsomo(ii))*1 + POPCNT(Jdomo(ii))*2 + enddo + + 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,1,idxI,NalphaIcfg) = Jsomo - alphasIcfg_list(1,2,idxI,NalphaIcfg) = IOR(Jdomo,ISHFT(1_8,n_core_orb)-1) + !if(idxI.eq.8)then + ! print *," 1 Idx = ",idxI, " Nalpha=",NalphaIcfg, " n_core_orb=",n_core_orb + !endif + alphasIcfg_list(:,1,idxI,NalphaIcfg) = Jcfg(:,1) + alphasIcfg_list(:,2,idxI,NalphaIcfg) = Jcfg(:,2) + if(n_core_orb .le. 64)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 NalphaIcfg_list(idxI) = NalphaIcfg endif + !print *," ", NalphaIcfg, Jsomo(1), Jsomo(2), "|", Jdomo(1), Jdomo(2) + endif + + if(vmotype(j) == 1)then + 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 ! 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)) - kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-2))) + do ii=1, N_int + Isomo(ii) = iand(act_bitmask(ii,1),psi_configuration(ii,1,idxI)) + Idomo(ii) = iand(act_bitmask(ii,2),psi_configuration(ii,2,idxI)) + enddo + !Icfg = psi_configuration(:,:,idxI) + + !kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-2))) + kstart = 1 + ndiffDOMO = 0 do k = kstart, idxI-1 - diffSOMO = IEOR(Isomo,iand(act_bitmask(1,1),psi_configuration(1,1,k))) - ndiffSOMO = POPCNT(diffSOMO) + ndiffSOMO = 0 + do ii=1,N_int + diffSOMO = IEOR(Icfg(ii,1),iand(act_bitmask(ii,1),psi_configuration(ii,1,k))) + ndiffSOMO += POPCNT(diffSOMO) + end do + ! ndiffSOMO cannot be 0 (I /= k) if idxI is a single ex + ! 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 - diffDOMO = IEOR(Idomo,iand(act_bitmask(1,1),psi_configuration(1,2,k))) - xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) - ndiffDOMO = POPCNT(diffDOMO) - nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO) + ndiffDOMO = 0 + nxordiffSOMODOMO = 0 + do ii=1,N_int + diffDOMO = IEOR(Icfg(ii,2),iand(act_bitmask(ii,2),psi_configuration(ii,2,k))) + xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) + ndiffDOMO += POPCNT(diffDOMO) + nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO) + end do if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4)) then ppExistsQ = .TRUE. EXIT @@ -277,15 +904,17 @@ use bitmasks ! Diagonal part (pp,qq) if(nholes > 0 .AND. (.NOT. ppExistsQ))then ! SOMO - if(POPCNT(Jsomo) .ge. NSOMOMin) then + if(Nsomo_I .ge. NSOMOMin) then NalphaIcfg += 1 - alphasIcfg_list(1,1,idxI,NalphaIcfg) = Icfg(1,1) - alphasIcfg_list(1,2,idxI,NalphaIcfg) = Icfg(1,2) + alphasIcfg_list(:,1,idxI,NalphaIcfg) = Icfg(:,1) + alphasIcfg_list(:,2,idxI,NalphaIcfg) = Icfg(:,2) NalphaIcfg_list(idxI) = NalphaIcfg endif + !print *," ---> ", NalphaIcfg, Icfg(1,1), Icfg(2,1), "|", Icfg(1,2), Icfg(2,2) endif NalphaIcfg = 0 + enddo ! end loop idxI call wall_time(t1) print *, 'Preparation : ', t1 - t0 @@ -303,6 +932,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 @@ -312,74 +942,86 @@ 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(bit_kind) :: diffDOMO, xordiffSOMODOMO, diffSOMO 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 ii=1, N_int + Isomo(ii) = iand(act_bitmask(ii,1),Icfg(ii,1)) + Idomo(ii) = iand(act_bitmask(ii,2),Icfg(ii,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 + nvmos=0 + 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 @@ -388,10 +1030,15 @@ END_PROVIDER allocate(tableUniqueAlphas(mo_num,mo_num)) tableUniqueAlphas = .FALSE. + ! Now find the allowed (p,q) excitations + do ii=1, N_int + Isomo(ii) = iand(act_bitmask(ii,1),Icfg(ii,1)) + Idomo(ii) = iand(act_bitmask(ii,2),Icfg(ii,2)) + Jsomo(ii) = Isomo(ii) + Jdomo(ii) = Idomo(ii) + 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 @@ -411,41 +1058,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((jint-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(ii,1) = Jsomo(ii) + Jcfg(ii,2) = Jdomo(ii) + 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 @@ -456,14 +1102,19 @@ 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 + diffSOMO = IEOR(Jcfg(ii,1),iand(reunion_of_act_virt_bitmask(ii,1),psi_configuration(ii,1,k))) + ndiffSOMO += POPCNT(diffSOMO) + diffDOMO = IEOR(Jcfg(ii,2),iand(reunion_of_act_virt_bitmask(ii,2),psi_configuration(ii,2,k))) + xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) + ndiffDOMO += POPCNT(diffDOMO) + nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO) + nxordiffSOMODOMO += POPCNT(diffSOMO) + POPCNT(diffDOMO) + 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. @@ -473,19 +1124,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 @@ -506,53 +1158,68 @@ 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 ii=1, N_int + Isomo(ii) = iand(act_bitmask(ii,1),Icfg(ii,1)) + Idomo(ii) = iand(act_bitmask(ii,2),Icfg(ii,2)) + Jsomo(ii) = Isomo(ii) + Jdomo(ii) = Idomo(ii) + 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((jint-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) + alphasIcfg_list(:,2,idxI,NalphaIcfg) = Jcfg(:,2) + if(n_core_orb .le. 64)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 @@ -560,15 +1227,25 @@ END_PROVIDER ! Check if this Icfg has been previously generated as a mono ppExistsQ = .False. - Isomo = iand(act_bitmask(1,1),Icfg(1,1)) - Idomo = iand(act_bitmask(1,1),Icfg(1,2)) + !Isomo = iand(act_bitmask(1,1),Icfg(1,1)) + !Idomo = iand(act_bitmask(1,2),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) + 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,2),psi_configuration(ii,2,k))) + xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) + ndiffDOMO += POPCNT(diffDOMO) + nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO) + end do if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then ppExistsQ = .TRUE. EXIT @@ -581,8 +1258,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 diff --git a/src/csf/conversion.irp.f b/src/csf/conversion.irp.f index 494c3bfa..92c8e669 100644 --- a/src/csf/conversion.irp.f +++ b/src/csf/conversion.irp.f @@ -114,6 +114,7 @@ subroutine convertWFfromCSFtoDET(N_st,psi_coef_cfg_in, psi_coef_det) integer :: idx integer MS MS = elec_alpha_num-elec_beta_num + !print *,"size=",size(tmp_psi_coef_det,1)," ",size(tmp_psi_coef_det,2) countcsf = 0 diff --git a/src/csf/obtain_I_foralpha.irp.f b/src/csf/obtain_I_foralpha.irp.f index 7d7ae09b..211d5af6 100644 --- a/src/csf/obtain_I_foralpha.irp.f +++ b/src/csf/obtain_I_foralpha.irp.f @@ -38,6 +38,7 @@ subroutine obtain_connected_J_givenI(idxI, givenI, connectedI, idxs_connectedI, integer :: holetype(mo_num) integer :: end_index integer :: Nsomo_I + integer :: listall(N_int*bit_kind_size), nelall ! ! 2 2 1 1 0 0 : 1 1 0 0 0 0 @@ -65,9 +66,12 @@ subroutine obtain_connected_J_givenI(idxI, givenI, connectedI, idxs_connectedI, ! Since CFGs are sorted wrt to seniority ! we don't have to search the full CFG list - Isomo = givenI(1,1) - Idomo = givenI(1,2) - Nsomo_I = POPCNT(Isomo) + Nsomo_I = 0 + do i=1,N_int + Isomo = givenI(i,1) + Idomo = givenI(i,2) + Nsomo_I += POPCNT(Isomo) + end do end_index = min(N_configuration,cfg_seniority_index(min(Nsomo_I+6,elec_num))-1) if(end_index .LT. 0) end_index= N_configuration !end_index = N_configuration @@ -83,17 +87,24 @@ subroutine obtain_connected_J_givenI(idxI, givenI, connectedI, idxs_connectedI, ! idxs_connectedI(nconnectedI)=i ! cycle !endif - Isomo = givenI(1,1) - Idomo = givenI(1,2) - Jsomo = psi_configuration(1,1,i) - Jdomo = psi_configuration(1,2,i) - diffSOMO = IEOR(Isomo,Jsomo) - ndiffSOMO = POPCNT(diffSOMO) - diffDOMO = IEOR(Idomo,Jdomo) - xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) - ndiffDOMO = POPCNT(diffDOMO) - nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO) - nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO + + ndiffSOMO = 0 + ndiffDOMO = 0 + nxordiffSOMODOMO = 0 + do ii=1,N_int + Isomo = givenI(ii,1) + Idomo = givenI(ii,2) + Jsomo = psi_configuration(ii,1,i) + Jdomo = psi_configuration(ii,2,i) + diffSOMO = IEOR(Isomo,Jsomo) + ndiffSOMO += POPCNT(diffSOMO) + diffDOMO = IEOR(Idomo,Jdomo) + xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) + ndiffDOMO += POPCNT(diffDOMO) + nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO) + nxordiffSOMODOMO += POPCNT(diffSOMO) + POPCNT(diffDOMO) + end do + if((nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then !------- ! MONO | @@ -144,25 +155,45 @@ subroutine obtain_connected_J_givenI(idxI, givenI, connectedI, idxs_connectedI, ! find out all pq holes possible nholes = 0 ! holes in SOMO - Isomo = psi_configuration(1,1,i) - Idomo = psi_configuration(1,2,i) - do iii = 1,n_act_orb - ii = list_act(iii) - if(POPCNT(IAND(Isomo,IBSET(0_8,ii-1))) .EQ. 1) then - nholes += 1 - listholes(nholes) = ii - holetype(nholes) = 1 - endif + !Isomo = psi_configuration(1,1,i) + !Idomo = psi_configuration(1,2,i) + !do iii = 1,n_act_orb + ! ii = list_act(iii) + ! if(POPCNT(IAND(Isomo,IBSET(0_8,ii-1))) .EQ. 1) then + ! nholes += 1 + ! listholes(nholes) = ii + ! holetype(nholes) = 1 + ! endif + !end do + + call bitstring_to_list(psi_configuration(1,1,i),listall,nelall,N_int) + + do iii=1,nelall + nholes += 1 + listholes(nholes) = listall(iii) + holetype(nholes) = 1 end do + ! holes in DOMO - do iii = 1,n_act_orb - ii = list_act(iii) - if(POPCNT(IAND(Idomo,IBSET(0_8,ii-1))) .EQ. 1) then - nholes += 1 - listholes(nholes) = ii - holetype(nholes) = 2 - endif + !do iii = 1,n_act_orb + ! ii = list_act(iii) + ! if(POPCNT(IAND(Idomo,IBSET(0_8,ii-1))) .EQ. 1) then + ! nholes += 1 + ! listholes(nholes) = ii + ! holetype(nholes) = 2 + ! endif + !end do + + call bitstring_to_list(psi_configuration(1,2,i),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 + ntotalconnectedI += max(1,(psi_config_data(i,2)-psi_config_data(i,1)+1)*nholes) endif end do @@ -199,6 +230,8 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI integer*8 :: Isomo integer*8 :: Jdomo integer*8 :: Jsomo + integer(bit_kind) :: Jcfg(N_int,2) + integer(bit_kind) :: Icfg(N_int,2) integer*8 :: IJsomo integer*8 :: diffSOMO integer*8 :: diffDOMO @@ -209,132 +242,261 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI integer :: iii,ii,i,j,k,l,p,q,nsomoJ,nsomoalpha,starti,endi,extyp,nholes integer :: listholes(mo_num) integer :: holetype(mo_num) - integer :: end_index - integer :: Nsomo_alpha + integer :: end_index, ishift + integer :: Nsomo_alpha, pp,qq, nperm, iint, ipos integer*8 :: MS + integer :: exc(0:2,2,2), tz, m, n, high, low + integer :: listall(N_int*bit_kind_size), nelall + integer :: nconnectedExtradiag, nconnectedDiag + integer(bit_kind) :: hole, particle, tmp MS = elec_alpha_num-elec_beta_num + nconnectedExtradiag=0 + nconnectedDiag=0 nconnectedI = 0 end_index = N_configuration ! Since CFGs are sorted wrt to seniority ! we don't have to search the full CFG list - Isomo = Ialpha(1,1) - Idomo = Ialpha(1,2) - Nsomo_alpha = POPCNT(Isomo) + !Isomo = Ialpha(1,1) + !Idomo = Ialpha(1,2) + !Nsomo_alpha = POPCNT(Isomo) + Icfg = Ialpha + Nsomo_alpha = 0 + !print *," Ialpha=" + do ii=1,N_int + Isomo = Ialpha(ii,1) + Idomo = Ialpha(ii,2) + Nsomo_alpha += POPCNT(Isomo) + !print *,Isomo, Idomo, "Nsomo=",Nsomo_alpha + end do end_index = min(N_configuration,cfg_seniority_index(min(Nsomo_alpha+4,elec_num))-1) - if(end_index .LT. 0) end_index= N_configuration + if(end_index .LT. 0 .OR. end_index .lt. idxI) end_index= N_configuration end_index = N_configuration p = 0 q = 0 - if (N_int > 1) stop 'obtain_connected_i_foralpha : N_int > 1' + !if (N_int > 1) stop 'obtain_connected_i_foralpha : N_int > 1' do i=idxI,end_index - Isomo = Ialpha(1,1) - Idomo = Ialpha(1,2) - Jsomo = psi_configuration(1,1,i) - Jdomo = psi_configuration(1,2,i) ! Check for Minimal alpha electrons (MS) - if(POPCNT(Isomo).lt.MS)then + if(Nsomo_alpha .lt. MS)then cycle endif - diffSOMO = IEOR(Isomo,Jsomo) - ndiffSOMO = POPCNT(diffSOMO) - !if(idxI.eq.1)then - ! print *," \t idxI=",i," diffS=",ndiffSOMO," popJs=", POPCNT(Jsomo)," popIs=",POPCNT(Isomo) + + ndiffSOMO = 0 + ndiffDOMO = 0 + nxordiffSOMODOMO = 0 + nsomoJ=0 + nsomoalpha=0 + do ii=1,N_int + Isomo = Ialpha(ii,1) + Idomo = Ialpha(ii,2) + Jsomo = psi_configuration(ii,1,i) + Jdomo = psi_configuration(ii,2,i) + nsomoJ += POPCNT(Jsomo) + nsomoalpha += POPCNT(Isomo) + diffSOMO = IEOR(Isomo,Jsomo) + ndiffSOMO += POPCNT(diffSOMO) + diffDOMO = IEOR(Idomo,Jdomo) + xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) + ndiffDOMO += POPCNT(diffDOMO) + nxordiffSOMODOMO += POPCNT(xordiffSOMODOMO) + nxordiffSOMODOMO += POPCNT(diffSOMO) + POPCNT(diffDOMO) + end do + !if(idxI.eq.218)then + ! print *,"I=",idxI,"Nsomo_alpha=",Nsomo_alpha,"nxordiffSOMODOMO(4)=",nxordiffSOMODOMO, " ndiffSOMO(2)=",ndiffSOMO, " ndiffDOMO=",ndiffDOMO !endif - diffDOMO = IEOR(Idomo,Jdomo) - xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) - ndiffDOMO = POPCNT(diffDOMO) - nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO) - nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO + !Jcfg = psi_configuration(:,:,i) + !print *,"nxordiffSOMODOMO(4)=",nxordiffSOMODOMO, " ndiffSOMO(2)=",ndiffSOMO + if((nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then select case(ndiffDOMO) case (0) ! SOMO -> VMO !print *,"obt SOMO -> VMO" extyp = 3 - IJsomo = IEOR(Isomo, Jsomo) -!IRP_IF WITHOUT_TRAILZ -! p = (popcnt(ieor( IAND(Isomo,IJsomo) , IAND(Isomo,IJsomo) -1))-1) + 1 -!IRP_ELSE - p = TRAILZ(IAND(Isomo,IJsomo)) + 1 -!IRP_ENDIF - IJsomo = IBCLR(IJsomo,p-1) -!IRP_IF WITHOUT_TRAILZ -! q = (popcnt(ieor(IJsomo,IJsomo-1))-1) + 1 -!IRP_ELSE - q = TRAILZ(IJsomo) + 1 -!IRP_ENDIF + !if(N_int .eq. 1) then + ! IJsomo = IEOR(Isomo, Jsomo) + ! p = TRAILZ(IAND(Isomo,IJsomo)) + 1 + ! IJsomo = IBCLR(IJsomo,p-1) + ! q = TRAILZ(IJsomo) + 1 + ! !print *," p=",p," q=",q + ! !call get_single_excitation_cfg(Jcfg, Icfg, p, q, N_int) + !else + ! Find p + do ii=1,N_int + Isomo = Ialpha(ii,1) + Jsomo = psi_configuration(ii,1,i) + IJsomo = IEOR(Isomo, Jsomo) + if(popcnt(IAND(Isomo,IJsomo)) > 0)then + p = TRAILZ(IAND(Isomo,IJsomo)) + 1 + (ii-1) * bit_kind_size + EXIT + endif + end do + ! Find q + do ii=1,N_int + Isomo = Ialpha(ii,1) + Jsomo = psi_configuration(ii,1,i) + IJsomo = IEOR(Isomo, Jsomo) + iint = shiftr(p-1,bit_kind_shift) + 1 + ipos = p-shiftl((iint-1),bit_kind_shift) + if(iint .eq. ii)then + IJsomo = IBCLR(IJsomo,ipos-1) + endif + if(popcnt(IJsomo) > 0)then + q = TRAILZ(IJsomo) + 1 + (ii-1) * bit_kind_size + EXIT + endif + enddo + !endif + !assert ( p == pp) + !assert ( q == qq) + !print *," 1--- p=",p," q=",q 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 -!IRP_IF WITHOUT_TRAILZ -! p = (popcnt(ieor( IEOR(Idomo,Jdomo),IEOR(Idomo,Jdomo) -1))-1) + 1 -!IRP_ELSE - p = TRAILZ(IEOR(Idomo,Jdomo)) + 1 -!IRP_ENDIF - Isomo = IEOR(Isomo, Jsomo) - Isomo = IBCLR(Isomo,p-1) -!IRP_IF WITHOUT_TRAILZ -! q = (popcnt(ieor(Isomo,Isomo-1))-1) + 1 -!IRP_ELSE - q = TRAILZ(Isomo) + 1 -!IRP_ENDIF + !if(N_int.eq.1)then + ! p = TRAILZ(IEOR(Idomo,Jdomo)) + 1 + ! Isomo = IEOR(Isomo, Jsomo) + ! Isomo = IBCLR(Isomo,p-1) + ! q = TRAILZ(Isomo) + 1 + !else + + ! Find p + do ii=1,N_int + Isomo = Ialpha(ii,1) + Jsomo = psi_configuration(ii,1,i) + Idomo = Ialpha(ii,2) + Jdomo = psi_configuration(ii,2,i) + if(popcnt(IEOR(Idomo,Jdomo)) > 0)then + p = TRAILZ(IEOR(Idomo,Jdomo)) + 1 + (ii-1) * bit_kind_size + EXIT + endif + end do + ! Find q + do ii=1,N_int + Isomo = Ialpha(ii,1) + Jsomo = psi_configuration(ii,1,i) + IJsomo = IEOR(Isomo, Jsomo) + iint = shiftr(p-1,bit_kind_shift) + 1 + ipos = p-shiftl((iint-1),bit_kind_shift) + if(iint .eq. ii)then + IJsomo = IBCLR(IJsomo,ipos-1) + endif + if(popcnt(IJsomo) > 0)then + q = TRAILZ(IJsomo) + 1 + (ii-1) * bit_kind_size + EXIT + endif + end do + !endif + !assert ( p == pp) + !assert ( q == qq) else ! SOMO -> SOMO !print *,"obt SOMO -> SOMO" extyp = 1 -!IRP_IF WITHOUT_TRAILZ -! q = (popcnt(ieor( IEOR(Idomo,Jdomo), IEOR(Idomo,Jdomo)-1))-1) + 1 -!IRP_ELSE - q = TRAILZ(IEOR(Idomo,Jdomo)) + 1 -!IRP_ENDIF - Isomo = IEOR(Isomo, Jsomo) - Isomo = IBCLR(Isomo,q-1) -!IRP_IF WITHOUT_TRAILZ -! p = (popcnt(ieor(Isomo,Isomo-1))-1) + 1 -!IRP_ELSE - p = TRAILZ(Isomo) + 1 -!IRP_ENDIF - ! Check for Minimal alpha electrons (MS) - !if(POPCNT(Isomo).lt.MS)then - ! cycle + !if(N_int.eq.1)then + ! q = TRAILZ(IEOR(Idomo,Jdomo)) + 1 + ! Isomo = IEOR(Isomo, Jsomo) + ! Isomo = IBCLR(Isomo,q-1) + ! p = TRAILZ(Isomo) + 1 + ! ! Check for Minimal alpha electrons (MS) + ! !if(POPCNT(Isomo).lt.MS)then + ! ! cycle + ! !endif + !else + ! Find p + !print *,"Ialpha somo=",Ialpha(1,1), Ialpha(2,1)," Ialpha domo=",Ialpha(1,2), Ialpha(2,2) + !print *,"J somo=",psi_configuration(1,1,i), psi_configuration(2,1,i)," J domo=",psi_configuration(1,2,i),& + !psi_configuration(2,2,i) + do ii=1,N_int + Isomo = Ialpha(ii,1) + Jsomo = psi_configuration(ii,1,i) + Idomo = Ialpha(ii,2) + Jdomo = psi_configuration(ii,2,i) + if(popcnt(IEOR(Idomo,Jdomo)) > 0)then + q = TRAILZ(IEOR(Idomo,Jdomo)) + 1 + (ii-1) * bit_kind_size + EXIT + endif + enddo + ! Find q + do ii=1,N_int + Isomo = Ialpha(ii,1) + Jsomo = psi_configuration(ii,1,i) + IJsomo = IEOR(Isomo, Jsomo) + iint = shiftr(q-1,bit_kind_shift) + 1 + ipos = q-shiftl((iint-1),bit_kind_shift) + if(iint .eq. ii)then + IJsomo = IBCLR(IJsomo,ipos-1) + endif + !print *,"ii=",ii," Isomo=",Isomo + if(popcnt(IJsomo) > 0)then + p = TRAILZ(IJsomo) + 1 + (ii-1) * bit_kind_size + EXIT + endif + enddo !endif - end if + !assert ( p == pp) + !assert ( q == qq) + endif + !print *," 2--- p=",p," q=",q case (2) ! DOMO -> SOMO !print *,"obt DOMO -> SOMO" extyp = 4 - IJsomo = IEOR(Isomo, Jsomo) -!IRP_IF WITHOUT_TRAILZ -! p = (popcnt(ieor( IAND(Jsomo,IJsomo), IAND(Jsomo,IJsomo)-1))-1) + 1 -!IRP_ELSE - p = TRAILZ(IAND(Jsomo,IJsomo)) + 1 -!IRP_ENDIF - IJsomo = IBCLR(IJsomo,p-1) -!IRP_IF WITHOUT_TRAILZ -! q = (popcnt(ieor( IJsomo , IJsomo -1))-1) + 1 -!IRP_ELSE - q = TRAILZ(IJsomo) + 1 -!IRP_ENDIF + !if(N_int.eq.1)then + ! IJsomo = IEOR(Isomo, Jsomo) + ! p = TRAILZ(IAND(Jsomo,IJsomo)) + 1 + ! IJsomo = IBCLR(IJsomo,p-1) + ! q = TRAILZ(IJsomo) + 1 + !else + ! Find p + do ii=1,N_int + Isomo = Ialpha(ii,1) + Jsomo = psi_configuration(ii,1,i) + Idomo = Ialpha(ii,2) + Jdomo = psi_configuration(ii,2,i) + IJsomo = IEOR(Isomo, Jsomo) + if(popcnt(IAND(Jsomo,IJsomo)) > 0)then + p = TRAILZ(IAND(Jsomo,IJsomo)) + 1 + (ii-1) * bit_kind_size + EXIT + endif + enddo + ! Find q + do ii=1,N_int + Isomo = Ialpha(ii,1) + Jsomo = psi_configuration(ii,1,i) + IJsomo = IEOR(Isomo, Jsomo) + iint = shiftr(p-1,bit_kind_shift) + 1 + ipos = p-shiftl((iint-1),bit_kind_shift) + if(iint .eq. ii)then + IJsomo = IBCLR(IJsomo,ipos-1) + endif + if(popcnt(IJsomo) > 0)then + q = TRAILZ(IJsomo) + 1 + (ii-1) * bit_kind_size + EXIT + endif + enddo + !endif + !assert ( p == pp) + !assert ( q == qq) + !print *," 3--- p=",p," q=",q case default print *,"something went wront in get connectedI" end select starti = psi_config_data(i,1) endi = psi_config_data(i,2) + nconnectedExtradiag+=1 nconnectedI += 1 - do k=1,N_int - connectedI(k,1,nconnectedI) = psi_configuration(k,1,i) - connectedI(k,2,nconnectedI) = psi_configuration(k,2,i) + do ii=1,N_int + connectedI(ii,1,nconnectedI) = psi_configuration(ii,1,i) + connectedI(ii,2,nconnectedI) = psi_configuration(ii,2,i) enddo idxs_connectedI(nconnectedI)=starti excitationIds(1,nconnectedI)=p @@ -343,28 +505,51 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI diagfactors(nconnectedI) = 1.0d0 else if((ndiffSOMO + ndiffDOMO) .EQ. 0) then ! find out all pq holes possible + !print *,"I = ",i + !print *,"I somo= ",psi_configuration(1,1,i), " domo=", psi_configuration(1,2,i) + !print *,"alp somo= ",Ialpha(1,1), " domo=", Ialpha(1,2) nholes = 0 ! holes in SOMO - Isomo = psi_configuration(1,1,i) - Idomo = psi_configuration(1,2,i) - do iii = 1,n_act_orb - ii = list_act(iii) - if(POPCNT(IAND(Isomo,IBSET(0_8,ii-1))) .EQ. 1) then - nholes += 1 - listholes(nholes) = ii - holetype(nholes) = 1 - endif + !Isomo = psi_configuration(1,1,i) + !Idomo = psi_configuration(1,2,i) + !do iii = 1,n_act_orb + ! ii = list_act(iii) + ! if(POPCNT(IAND(Isomo,IBSET(0_8,ii-1))) .EQ. 1) then + ! nholes += 1 + ! listholes(nholes) = ii + ! holetype(nholes) = 1 + ! endif + !end do + call bitstring_to_list(psi_configuration(1,1,i),listall,nelall,N_int) + + do iii=1,nelall + nholes += 1 + listholes(nholes) = listall(iii) + holetype(nholes) = 1 end do + ! holes in DOMO - do iii = 1,n_act_orb - ii = list_act(iii) - if(POPCNT(IAND(Idomo,IBSET(0_8,ii-1))) .EQ. 1) then - nholes += 1 - listholes(nholes) = ii - holetype(nholes) = 2 - endif + !do iii = 1,n_act_orb + ! ii = list_act(iii) + ! if(POPCNT(IAND(Idomo,IBSET(0_8,ii-1))) .EQ. 1) then + ! nholes += 1 + ! listholes(nholes) = ii + ! holetype(nholes) = 2 + ! endif + !end do + nelall=0 + listall=0 + call bitstring_to_list(psi_configuration(1,2,i),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 + do k=1,nholes p = listholes(k) q = p @@ -372,6 +557,7 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI if(holetype(k) .EQ. 1) then starti = psi_config_data(i,1) endi = psi_config_data(i,2) + nconnectedDiag+=1 nconnectedI += 1 connectedI(:,:,nconnectedI) = psi_configuration(:,:,i) idxs_connectedI(nconnectedI)=starti @@ -382,6 +568,7 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI else starti = psi_config_data(i,1) endi = psi_config_data(i,2) + nconnectedDiag+=1 nconnectedI += 1 connectedI(:,:,nconnectedI) = psi_configuration(:,:,i) idxs_connectedI(nconnectedI)=starti @@ -390,8 +577,10 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI excitationTypes(nconnectedI) = extyp diagfactors(nconnectedI) = 2.0d0 endif + !print *,excitationIds(1,nconnectedI), excitationIds(2,nconnectedI) enddo endif end do + !print *,"nconnectedExtradiag=",nconnectedExtradiag," nconnectedDiad=",nconnectedDiag end subroutine obtain_connected_I_foralpha diff --git a/src/csf/sigma_vector.irp.f b/src/csf/sigma_vector.irp.f index 2d6f2ee2..c6644c4c 100644 --- a/src/csf/sigma_vector.irp.f +++ b/src/csf/sigma_vector.irp.f @@ -146,7 +146,6 @@ ncfgprev = cfg_seniority_index(i+2) end do !print *," ^^^^^ N_CSF = ",n_CSF," N_CFG=",N_configuration - END_PROVIDER @@ -832,7 +831,7 @@ subroutine calculate_preconditioner_cfg(diag_energies) ! the configurations in psi_configuration ! returns : diag_energies : END_DOC - integer :: i,j,k,kk,l,p,q,noccp,noccq, ii, jj + integer :: i,j,k,kk,l,p,q,noccp,noccq, ii, jj, iii real*8,intent(out) :: diag_energies(n_CSF) integer :: nholes integer :: nvmos @@ -859,6 +858,7 @@ subroutine calculate_preconditioner_cfg(diag_energies) real*8 :: hpp real*8 :: meCC real*8 :: core_act_contrib + integer :: listall(N_int*bit_kind_size), nelall !PROVIDE h_core_ri PROVIDE core_fock_operator @@ -875,11 +875,11 @@ subroutine calculate_preconditioner_cfg(diag_energies) 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)) + !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) @@ -888,48 +888,63 @@ subroutine calculate_preconditioner_cfg(diag_energies) ! find out all pq holes possible nholes = 0 + listholes = -1 ! holes in SOMO - !do k = 1,mo_num - do kk = 1,n_act_orb - k = list_act(kk) - 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 - do kk = 1,n_act_orb - k = list_act(kk) - if(POPCNT(IAND(Idomo,IBSET(0_8,k-1))) .EQ. 1) then - nholes += 1 - listholes(nholes) = k - holetype(nholes) = 2 - endif - enddo + !do kk = 1,n_act_orb + ! k = list_act(kk) + ! if(POPCNT(IAND(Isomo,IBSET(0_8,k-1))) .EQ. 1) then + ! nholes += 1 + ! listholes(nholes) = k + ! holetype(nholes) = 1 + ! endif + !enddo + call bitstring_to_list(psi_configuration(1,1,i),listall,nelall,N_int) - ! 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 - do kk = 1,n_act_orb - k = list_act(kk) - !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 + do iii=1,nelall + nholes += 1 + listholes(nholes) = listall(iii) + holetype(nholes) = 1 + end do + + ! holes in DOMO + !do kk = 1,n_act_orb + ! k = list_act(kk) + ! if(POPCNT(IAND(Idomo,IBSET(0_8,k-1))) .EQ. 1) then + ! nholes += 1 + ! listholes(nholes) = k + ! holetype(nholes) = 2 + ! endif + !enddo + call bitstring_to_list(psi_configuration(1,2,i),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 k = n_core_orb+1,n_core_orb + n_act_orb + !!!do k = 1,mo_num + !!do kk = 1,n_act_orb + !! k = list_act(kk) + !! !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) @@ -1219,27 +1234,30 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod integer,intent(in) :: p,q integer,intent(in) :: extype integer,intent(out) :: pmodel,qmodel - !integer(bit_kind) :: Isomo(N_int) - !integer(bit_kind) :: Idomo(N_int) - !integer(bit_kind) :: Jsomo(N_int) - !integer(bit_kind) :: Jdomo(N_int) - integer*8 :: Isomo - integer*8 :: Idomo - integer*8 :: Jsomo - integer*8 :: Jdomo + integer(bit_kind) :: Isomo(N_int) + integer(bit_kind) :: Idomo(N_int) + integer(bit_kind) :: Jsomo(N_int) + integer(bit_kind) :: Jdomo(N_int) + !integer*8 :: Isomo + !integer*8 :: Idomo + !integer*8 :: Jsomo + !integer*8 :: Jdomo integer*8 :: mask - integer :: iint, ipos + integer :: iint, ipos, ii !integer(bit_kind) :: Isomotmp(N_int) !integer(bit_kind) :: Jsomotmp(N_int) integer*8 :: Isomotmp integer*8 :: Jsomotmp integer :: pos0,pos0prev + integer :: tmpp, tmpq ! 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) + do ii=1,N_int + Isomo(ii) = Ialpha(ii,1) + Idomo(ii) = Ialpha(ii,2) + Jsomo(ii) = Jcfg(ii,1) + Jdomo(ii) = Jcfg(ii,2) + end do pos0prev = 0 pmodel = p qmodel = q @@ -1253,40 +1271,155 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod ! 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)) + !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)) + + iint = shiftr(p-1,bit_kind_shift) + 1 + ipos = p-shiftl((iint-1),bit_kind_shift)-1 + tmpp = 0 + !print *,"iint=",iint, " p=",p + do ii=1,iint-1 + !mask = ISHFT(1_bit_kind,-1)-1_bit_kind + !Isomotmp = IAND(Isomo(ii),mask) + !tmpp += POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) + tmpp += POPCNT(Isomo(ii)) + end do + mask = ISHFT(1_bit_kind,ipos+1) - 1 + Isomotmp = IAND(Isomo(iint),mask) + !pmodel = tmpp + POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) + pmodel = tmpp + POPCNT(Isomotmp) + !print *,"iint=",iint, " ipos=",ipos,"pmodel=",pmodel, XOR(Isomotmp,mask),Isomo(iint) + + iint = shiftr(q-1,bit_kind_shift) + 1 + ipos = q-shiftl((iint-1),bit_kind_shift)-1 + tmpq = 0 + do ii=1,iint-1 + !mask = ISHFT(1_bit_kind,-1)-1_bit_kind + !Isomotmp = IAND(Isomo(ii),mask) + !tmpq += POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) + tmpq += POPCNT(Isomo(ii)) + end do + mask = ISHFT(1_bit_kind,ipos+1) - 1 + Isomotmp = IAND(Isomo(iint),mask) + !qmodel = tmpq + POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) + qmodel = tmpq + POPCNT(Isomotmp) + !print *,"iint=",iint, " ipos=",ipos,"qmodel=",qmodel 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)) + !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)) + + iint = shiftr(p-1,bit_kind_shift) + 1 + ipos = p-shiftl((iint-1),bit_kind_shift)-1 + tmpp = 0 + do ii=1,iint-1 + !mask = ISHFT(1_bit_kind,-1)-1_bit_kind + !Jsomotmp = IAND(Jsomo(ii),mask) + !tmpp += POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) + tmpp += POPCNT(Jsomo(ii)) + end do + mask = ISHFT(1_bit_kind,ipos+1) - 1 + Jsomotmp = IAND(Jsomo(iint),mask) + !pmodel = tmpp + POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) + pmodel = tmpp + POPCNT(Jsomotmp) + + iint = shiftr(q-1,bit_kind_shift) + 1 + ipos = q-shiftl((iint-1),bit_kind_shift)-1 + tmpq = 0 + do ii=1,iint-1 + !mask = ISHFT(1_bit_kind,-1)-1_bit_kind + !Jsomotmp = IAND(Jsomo(ii),mask) + !tmpq += POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) + tmpq += POPCNT(Jsomo(ii)) + end do + mask = ISHFT(1_bit_kind,ipos+1) - 1 + Jsomotmp = IAND(Jsomo(iint),mask) + !qmodel = tmpq + POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) + qmodel = tmpq + POPCNT(Jsomotmp) 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 + !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 + + iint = shiftr(p-1,bit_kind_shift) + 1 + ipos = p-shiftl((iint-1),bit_kind_shift)-1 + tmpp = 0 + do ii=1,iint-1 + !mask = ISHFT(1_bit_kind,-1)-1_bit_kind + !Isomotmp = IAND(Isomo(ii),mask) + !tmpp += POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) + tmpp += POPCNT(Isomo(ii)) + end do + mask = ISHFT(1_bit_kind,ipos+1) - 1 + Isomotmp = IAND(Isomo(iint),mask) + !pmodel = tmpp + POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) + pmodel = tmpp + POPCNT(Isomotmp) + + iint = shiftr(q-1,bit_kind_shift) + 1 + ipos = q-shiftl((iint-1),bit_kind_shift)-1 + tmpq = 0 + do ii=1,iint-1 + !mask = ISHFT(1_bit_kind,-1)-1_bit_kind + !Jsomotmp = IAND(Jsomo(ii),mask) + !tmpq += POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) + tmpq += POPCNT(Jsomo(ii)) + end do + mask = ISHFT(1_bit_kind,ipos+1) - 1 + Jsomotmp = IAND(Jsomo(iint),mask) + !qmodel = tmpq + POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) + 1 + qmodel = tmpq + POPCNT(Jsomotmp) + 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)) + !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)) + + iint = shiftr(p-1,bit_kind_shift) + 1 + ipos = p-shiftl((iint-1),bit_kind_shift)-1 + tmpp = 0 + do ii=1,iint-1 + !mask = ISHFT(1_bit_kind,-1)-1_bit_kind + !Isomotmp = IAND(Isomo(ii),mask) + !tmpp += POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) + tmpp += POPCNT(Isomo(ii)) + end do + mask = ISHFT(1_bit_kind,ipos+1) - 1 + Isomotmp = IAND(Isomo(iint),mask) + !pmodel = tmpp + POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) + 1 + pmodel = tmpp + POPCNT(Isomotmp) + 1 + + iint = shiftr(q-1,bit_kind_shift) + 1 + ipos = q-shiftl((iint-1),bit_kind_shift)-1 + tmpq = 0 + do ii=1,iint-1 + !mask = ISHFT(1_bit_kind,-1)-1_bit_kind + !Jsomotmp = IAND(Jsomo(ii),mask) + !tmpq += POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) + tmpq += POPCNT(Jsomo(ii)) + end do + mask = ISHFT(1_bit_kind,ipos+1) - 1 + Jsomotmp = IAND(Jsomo(iint),mask) + !qmodel = tmpq + POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) + qmodel = tmpq + POPCNT(Jsomotmp) endif case (4) ! DOMO -> SOMO @@ -1294,19 +1427,75 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod !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 + !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 + + iint = shiftr(p-1,bit_kind_shift) + 1 + ipos = p-shiftl((iint-1),bit_kind_shift)-1 + tmpp = 0 + do ii=1,iint-1 + !mask = ISHFT(1_bit_kind,-1)-1_bit_kind + !Jsomotmp = IAND(Jsomo(ii),mask) + !tmpp += POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) + tmpp += POPCNT(Jsomo(ii)) + end do + mask = ISHFT(1_bit_kind,ipos+1) - 1 + Jsomotmp = IAND(Jsomo(iint),mask) + !pmodel = tmpp + POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) + pmodel = tmpp + POPCNT(Jsomotmp) + + iint = shiftr(q-1,bit_kind_shift) + 1 + ipos = q-shiftl((iint-1),bit_kind_shift)-1 + tmpq = 0 + do ii=1,iint-1 + !mask = ISHFT(1_bit_kind,-1)-1_bit_kind + !Isomotmp = IAND(Isomo(ii),mask) + !tmpq += POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) + tmpq += POPCNT(Isomo(ii)) + end do + mask = ISHFT(1_bit_kind,ipos+1) - 1 + Isomotmp = IAND(Isomo(iint),mask) + !qmodel = tmpq + POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) + 1 + qmodel = tmpq + POPCNT(Isomotmp) + 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)) + !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)) + + iint = shiftr(p-1,bit_kind_shift) + 1 + ipos = p-shiftl((iint-1),bit_kind_shift)-1 + tmpp = 0 + do ii=1,iint-1 + !mask = ISHFT(1_bit_kind,-1)-1_bit_kind + !Jsomotmp = IAND(Jsomo(ii),mask) + !tmpp += POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) + tmpp += POPCNT(Jsomo(ii)) + end do + mask = ISHFT(1_bit_kind,ipos+1) - 1 + Jsomotmp = IAND(Jsomo(iint),mask) + !pmodel = tmpp + POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) + 1 + pmodel = tmpp + POPCNT(Jsomotmp) + 1 + + iint = shiftr(q-1,bit_kind_shift) + 1 + ipos = q-shiftl((iint-1),bit_kind_shift)-1 + tmpq = 0 + do ii=1,iint-1 + !mask = ISHFT(1_bit_kind,-1)-1_bit_kind + !Isomotmp = IAND(Isomo(ii),mask) + !tmpq += POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) + tmpq += POPCNT(Isomo(ii)) + end do + mask = ISHFT(1_bit_kind,ipos+1) - 1 + Isomotmp = IAND(Isomo(iint),mask) + !qmodel = tmpq + POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) + qmodel = tmpq + POPCNT(Isomotmp) endif case default print *,"something is wrong in convertOrbIdsToModelSpaceIds" @@ -1364,8 +1553,13 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze integer :: rowsTKI integer :: noccpp integer :: istart_cfg, iend_cfg, num_threads_max + integer :: iint, jint, ipos, jpos, Nsomo_I, iii integer :: nconnectedJ,nconnectedtotalmax,nconnectedmaxJ,maxnalphas,ntotJ - integer*8 :: MS, Isomo, Idomo, Jsomo, Jdomo, Ialpha, Ibeta + integer*8 :: MS,Ialpha, Ibeta + integer(bit_kind) :: Isomo(N_INT) + integer(bit_kind) :: Idomo(N_INT) + integer(bit_kind) :: Jsomo(N_INT) + integer(bit_kind) :: Jdomo(N_INT) integer :: moi, moj, mok, mol, starti, endi, startj, endj, cnti, cntj, cntk real*8 :: norm_coef_cfg, fac2eints real*8 :: norm_coef_det @@ -1380,6 +1574,8 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze real*8,dimension(:),allocatable:: diag_energies real*8 :: tmpvar, tmptot real*8 :: core_act_contrib + integer :: listall(N_int*bit_kind_size), nelall + integer :: countelec integer(omp_lock_kind), allocatable :: lock(:) call omp_set_max_active_levels(1) @@ -1408,8 +1604,8 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze !nconnectedtotalmax = 1000 !nconnectedmaxJ = 1000 maxnalphas = elec_num*mo_num - Icfg(1,1) = psi_configuration(1,1,1) - Icfg(1,2) = psi_configuration(1,2,1) + Icfg(:,1) = psi_configuration(:,1,1) + Icfg(:,2) = psi_configuration(:,2,1) allocate(listconnectedJ(N_INT,2,max(sze,10000))) allocate(idslistconnectedJ(max(sze,10000))) call obtain_connected_J_givenI(1, Icfg, listconnectedJ, idslistconnectedJ, nconnectedmaxJ, nconnectedtotalmax) @@ -1441,6 +1637,7 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze !$OMP shared(istart_cfg, iend_cfg, psi_configuration, mo_num, psi_config_data,& !$OMP N_int, N_st, psi_out, psi_in, h_core_ri, core_energy, h_act_ri, AIJpqContainer,& !$OMP pp, sze, NalphaIcfg_list,alphasIcfg_list, bit_tmp, & + !$OMP qq, iint, jint, ipos, jpos, nelall, listall, Nsomo_I, countelec,& !$OMP AIJpqMatrixDimsList, diag_energies, n_CSF, lock, NBFmax,nconnectedtotalmax, nconnectedmaxJ,maxnalphas,& !$OMP n_core_orb, n_act_orb, list_act, n, list_core, list_core_is_built,core_act_contrib, num_threads_max,& !$OMP n_core_orb_is_built, mo_integrals_map, mo_integrals_map_is_built) @@ -1463,11 +1660,13 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze ! else ! cycle - 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) + do ii=1,N_INT + Icfg(ii,1) = psi_configuration(ii,1,i) + Icfg(ii,2) = psi_configuration(ii,2,i) + Isomo(ii) = Icfg(ii,1) + Idomo(ii) = Icfg(ii,2) + enddo + NSOMOI = getNSOMO(Icfg) ! find out all pq holes possible nholes = 0 @@ -1477,42 +1676,86 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze ! list_core_inact ! bitmasks !do k = 1,mo_num - do kk = 1,n_act_orb - k = list_act(kk) - 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 = 1,mo_num - do kk = 1,n_act_orb - k = list_act(kk) - if(POPCNT(IAND(Idomo,IBSET(0_8,k-1))) .EQ. 1) then - nholes += 1 - listholes(nholes) = k - holetype(nholes) = 2 - endif - enddo + ! do kk = 1,n_act_orb + ! k = list_act(kk) + ! 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 = 1,mo_num + ! do kk = 1,n_act_orb + ! k = list_act(kk) + ! if(POPCNT(IAND(Idomo,IBSET(0_8,k-1))) .EQ. 1) then + ! nholes += 1 + ! listholes(nholes) = k + ! holetype(nholes) = 2 + ! endif + ! enddo + + ! ! find vmos + ! do kk = 1,n_act_orb + ! k = list_act(kk) + ! !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 + + ! find out all pq holes possible + 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 kk = 1,n_act_orb - k = list_act(kk) - !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 + ! find vmos + ! 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 + ! Icsf ids @@ -1531,16 +1774,31 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze extype = excitationTypes_single(j) ! Off diagonal terms call convertOrbIdsToModelSpaceIds(Icfg, singlesI(1,1,j), p, q, extype, pmodel, qmodel) - Jsomo = singlesI(1,1,j) - Jdomo = singlesI(1,2,j) + do ii=1,N_INT + Jsomo(ii) = singlesI(1,1,j) + Jdomo(ii) = singlesI(1,2,j) + enddo + + ! Get actual p pos + pp = p + iint = shiftr(pp-1,bit_kind_shift) + 1 + ipos = pp-shiftl((iint-1),bit_kind_shift)-1 + + ! Get actual q pos + qq = q + jint = shiftr(qq-1,bit_kind_shift) + 1 + jpos = qq-shiftl((jint-1),bit_kind_shift)-1 ! 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 + !if(POPCNT(IAND(Jsomo,IBSET(0_8,q-1))) .EQ. 1 .AND. POPCNT(IAND(Isomo,IBSET(0_8,q-1))) .EQ. 0) then + if(POPCNT(IAND(Jsomo(jint),IBSET(0_8,jpos))) .EQ. 1 .AND. POPCNT(IAND(Isomo(jint),IBSET(0_8,jpos))) .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 + !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 + if((POPCNT(IAND(Jdomo(jint),IBSET(0_8,jpos))) .EQ. 1 .AND. POPCNT(IAND(Idomo(jint),IBSET(0_8,jpos))) .EQ. 0) .AND.& + POPCNT(IAND(Isomo(jint),IBSET(0_8,jpos))) .EQ. 0) then nholes += 1 listholes(nholes) = q holetype(nholes) = 2 @@ -1576,10 +1834,12 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze 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 + !if(POPCNT(IAND(Jsomo,IBSET(0_8,q-1))) .EQ. 1 .AND. POPCNT(IAND(Isomo,IBSET(0_8,q-1))) .EQ. 0) then + if(POPCNT(IAND(Jsomo(jint),IBSET(0_8,jpos))) .EQ. 1 .AND. POPCNT(IAND(Isomo(jint),IBSET(0_8,jpos))) .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 + if((POPCNT(IAND(Jdomo(jint),IBSET(0_8,jpos))) .EQ. 1 .AND. POPCNT(IAND(Idomo(jint),IBSET(0_8,jpos))) .EQ. 0) .AND.& + POPCNT(IAND(Isomo(jint),IBSET(0_8,jpos))) .EQ. 0) then nholes -= 1 endif enddo @@ -1591,6 +1851,9 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze deallocate(excitationTypes_single) !print *," singles part psi(1,1)=",psi_out(1,1) + !do i=1,n_CSF + ! print *,"i=",i," psi(i)=",psi_out(1,i) + !enddo allocate(listconnectedJ(N_INT,2,max(sze,10000))) allocate(alphas_Icfg(N_INT,2,max(sze,10000))) @@ -1605,7 +1868,6 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze !!!====================!!! !!! Double Excitations !!! !!!====================!!! - ! Loop over all selected configurations !$OMP DO SCHEDULE(static) do i = istart_cfg,iend_cfg @@ -1615,8 +1877,10 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze ! else ! cycle - Icfg(1,1) = psi_configuration(1,1,i) - Icfg(1,2) = psi_configuration(1,2,i) + do ii=1,N_INT + Icfg(ii,1) = psi_configuration(ii,1,i) + Icfg(ii,2) = psi_configuration(ii,2,i) + enddo starti = psi_config_data(i,1) endi = psi_config_data(i,2) @@ -1627,14 +1891,15 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze Nalphas_Icfg = NalphaIcfg_list(i) alphas_Icfg(1:n_int,1:2,1:Nalphas_Icfg) = alphasIcfg_list(1:n_int,1:2,i,1:Nalphas_Icfg) - if(Nalphas_Icfg .GT. maxnalphas) then - print *,"Nalpha > maxnalpha" - endif + !if(Nalphas_Icfg .GT. maxnalphas) then + ! print *,"Nalpha > maxnalpha" + !endif - call obtain_connected_J_givenI(i, Icfg, listconnectedJ, idslistconnectedJ, nconnectedJ, ntotJ) + !call obtain_connected_J_givenI(i, Icfg, listconnectedJ, idslistconnectedJ, nconnectedJ, ntotJ) ! TODO : remove doubly excited for return - !print *,"I=",i," isomo=",psi_configuration(1,1,i)," idomo=",psi_configuration(1,2,i), " psiout=",psi_out(1,5) + !print *,"I=",i,"isomo=",psi_configuration(1,1,i),psi_configuration(2,1,i),POPCNT(psi_configuration(1,1,i)),POPCNT(psi_configuration(2,1,i)),& + !"idomo=",psi_configuration(1,2,i),psi_configuration(2,2,i),POPCNT(psi_configuration(1,2,i)),POPCNT(psi_configuration(2,2,i)), "Nalphas_Icfg=",Nalphas_Icfg do k = 1,Nalphas_Icfg ! Now generate all singly excited with respect to a given alpha CFG @@ -1645,15 +1910,18 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze call obtain_connected_I_foralpha(i, alphas_Icfg(1,1,k), connectedI_alpha, idxs_connectedI_alpha, & nconnectedI, excitationIds, excitationTypes, diagfactors) + !if(i .EQ. 218) then + ! print *,'k=',k,' kcfgSOMO=',alphas_Icfg(1,1,k),alphas_Icfg(2,1,k),' ',POPCNT(alphas_Icfg(1,1,k)),' & + ! kcfgDOMO=',alphas_Icfg(1,2,k),alphas_Icfg(2,2,k),' ',POPCNT(alphas_Icfg(1,2,k)), " NconnectedI=",nconnectedI + ! !print *,'k=',k,' kcfgSOMO=',alphas_Icfg(1,1,k),' ',POPCNT(alphas_Icfg(1,1,k)),' & + ! !kcfgDOMO=',alphas_Icfg(1,2,k),' ',POPCNT(alphas_Icfg(1,2,k)), " NconnectedI=",nconnectedI + !endif + if(nconnectedI .EQ. 0) then cycle endif - !if(i .EQ. 1) then - ! print *,'k=',k,' kcfgSOMO=',alphas_Icfg(1,1,k),' ',POPCNT(alphas_Icfg(1,1,k)),' kcfgDOMO=',alphas_Icfg(1,2,k),' ',POPCNT(alphas_Icfg(1,2,k)) - !endif - ! Here we do 2x the loop. One to count for the size of the matrix, then we compute. totcolsTKI = 0 rowsTKI = -1 @@ -1663,15 +1931,30 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze p = excitationIds(1,j) q = excitationIds(2,j) extype = excitationTypes(j) + !print *,"K=",k,"j=",j, "countelec=",countelec," p=",p," q=",q, " extype=",extype, "NSOMOalpha=",NSOMOalpha," NSOMOI=",NSOMOI, "alphas_Icfg(1,1,k)=",alphas_Icfg(1,1,k), & + !alphas_Icfg(2,1,k), " domo=",alphas_Icfg(1,2,k), alphas_Icfg(2,2,k), " connected somo=",connectedI_alpha(1,1,j), & + !connectedI_alpha(2,1,j), " domo=",connectedI_alpha(1,2,j), connectedI_alpha(2,2,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 rowsikpq = AIJpqMatrixDimsList(NSOMOalpha,extype,pmodel,qmodel,1) colsikpq = AIJpqMatrixDimsList(NSOMOalpha,extype,pmodel,qmodel,2) + !if(i.eq.218)then + ! print *,"j=",j," k=",k,"p=",p,"q=",q,"NSOMOalpha=",NSOMOalpha, "pmodel=",pmodel,"qmodel=",qmodel, "extype=",extype,& + ! "conn somo=",connectedI_alpha(1,1,j),connectedI_alpha(2,1,j),& + ! "conn domo=",connectedI_alpha(1,2,j),connectedI_alpha(2,2,j) + ! do m=1,colsikpq + ! print *,idxs_connectedI_alpha(j)+m-1 + ! enddo + !endif !print *,"j=",j," Nsomo=",NSOMOalpha," rowsikpq=",rowsikpq," colsikpq=",colsikpq, " p=",pmodel," q=",qmodel, " extyp=",extype totcolsTKI += colsikpq rowsTKI = rowsikpq enddo + !if(i.eq.1)then + ! print *,"n_st=",n_st,"rowsTKI=",rowsTKI, " nconnectedI=",nconnectedI, & + ! "totcolsTKI=",totcolsTKI + !endif allocate(TKI(n_st,rowsTKI,totcolsTKI)) ! coefficients of CSF ! Initialize the integral container ! dims : (totcolsTKI, nconnectedI) @@ -1701,10 +1984,10 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze TKI(kk,l,totcolsTKI+m) = AIJpqContainer(l,m,pmodel,qmodel,extype,NSOMOalpha) & * psi_in(kk,idxs_connectedI_alpha(j)+m-1) enddo - !if(i.eq.1) then - ! print *,AIJpqContainer(l,m,pmodel,qmodel,extype,NSOMOalpha) - !endif enddo + !if(i.eq.1) then + ! print *,"j=",j,"psi_in=",psi_in(1,idxs_connectedI_alpha(j)+m-1) + !endif enddo diagfactors_0 = diagfactors(j)*0.5d0 @@ -1743,16 +2026,24 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze rowsTKI = rowsikpq CCmattmp = 0.d0 + !if(i.eq.1)then + ! print *,"\t n_st=",n_st," colsikpq=",colsikpq," rowsTKI=",rowsTKI,& + ! " | ",size(TKIGIJ,1),size(AIJpqContainer,1),size(CCmattmp,1) + !endif call dgemm('N','N', n_st, colsikpq, rowsTKI, 1.d0, & TKIGIJ(1,1,j), size(TKIGIJ,1), & AIJpqContainer(1,1,pmodel,qmodel,extype,NSOMOalpha), & size(AIJpqContainer,1), 0.d0, & CCmattmp, size(CCmattmp,1) ) + !print *,"j=",j,"colsikpq=",colsikpq, "sizeTIG=",size(TKIGIJ,1),"sizeaijpq=",size(AIJpqContainer,1) do m = 1,colsikpq call omp_set_lock(lock(idxs_connectedI_alpha(j)+m-1)) do kk = 1,n_st psi_out(kk,idxs_connectedI_alpha(j)+m-1) += CCmattmp(kk,m) + !if(dabs(CCmattmp(kk,m)).gt.1e-10)then + ! print *, CCmattmp(kk,m), " | ",idxs_connectedI_alpha(j)+m-1 + !end if enddo call omp_unset_lock(lock(idxs_connectedI_alpha(j)+m-1)) enddo @@ -1787,6 +2078,10 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze !$OMP END DO !$OMP END PARALLEL + !print *," ----- " + !do i=1,sze + ! print *,"i=",i," psi_out(i)=",psi_out(1,i) + !end do call omp_set_max_active_levels(4) deallocate(diag_energies) diff --git a/src/davidson/diagonalization_hcfg.irp.f b/src/davidson/diagonalization_hcfg.irp.f index 659602a1..8e12b9c8 100644 --- a/src/davidson/diagonalization_hcfg.irp.f +++ b/src/davidson/diagonalization_hcfg.irp.f @@ -112,6 +112,8 @@ subroutine davidson_diag_cfg_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N double precision, allocatable :: U(:,:), U_csf(:,:), overlap(:,:) double precision, allocatable :: tmpU(:,:), tmpW(:,:) double precision, pointer :: W(:,:), W_csf(:,:) + !double precision, pointer :: W2(:,:), W_csf2(:,:) + !double precision, allocatable :: U2(:,:), U_csf2(:,:) logical :: disk_based double precision :: energy_shift(N_st_diag_in*davidson_sze_max) @@ -234,12 +236,15 @@ subroutine davidson_diag_cfg_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N call c_f_pointer(ptr_w, W_csf, (/sze_csf,N_st_diag*itermax/)) else allocate(W(sze,N_st_diag),W_csf(sze_csf,N_st_diag*itermax)) + !allocate(W2(sze,N_st_diag),W_csf2(sze_csf,N_st_diag*itermax)) endif allocate( & ! Large U(sze,N_st_diag), & + !U2(sze,N_st_diag), & U_csf(sze_csf,N_st_diag*itermax), & + !U_csf2(sze_csf,N_st_diag*itermax), & ! Small h(N_st_diag*itermax,N_st_diag*itermax), & @@ -325,7 +330,7 @@ subroutine davidson_diag_cfg_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N enddo enddo !tmpU =0.0d0 - !tmpU(1,2)=1.0d0 + !tmpU(1,1)=1.0d0 double precision :: irp_rdtsc double precision :: ticks_0, ticks_1 integer*8 :: irp_imax @@ -348,9 +353,9 @@ subroutine davidson_diag_cfg_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N !call convertWFfromDETtoCSF(N_st_diag,u_in(1,1),W_csf2(1,1)) !do i=1,sze_csf ! print *,"I=",i," qp=",W_csf2(i,1)," my=",W_csf(i,1)," diff=",dabs(W_csf2(i,1))-dabs(W_csf(i,1)) - ! if(dabs(dabs(W_csf2(i,1))-dabs(W_csf(i,1))) .gt. 1.0e-10)then - ! print *,"somo=",psi_configuration(1,1,i)," domo=",psi_configuration(1,2,i)," diff=",dabs(W_csf2(i,1))-dabs(W_csf(i,1)) - ! endif + ! !if(dabs(dabs(W_csf2(i,1))-dabs(W_csf(i,1))) .gt. 1.0e-10)then + ! ! print *,"somo=",psi_configuration(1,1,i)," domo=",psi_configuration(1,2,i)," diff=",dabs(W_csf2(i,1))-dabs(W_csf(i,1)) + ! !endif !end do !stop deallocate(tmpW) diff --git a/src/determinants/slater_rules.irp.f b/src/determinants/slater_rules.irp.f index a34608f9..b9710fd1 100644 --- a/src/determinants/slater_rules.irp.f +++ b/src/determinants/slater_rules.irp.f @@ -83,7 +83,7 @@ subroutine get_excitation(det1,det2,exc,degree,phase,Nint) ! exc(1,1,1) = q ! exc(1,2,1) = p - ! T^alpha_pq : exc(0,1,2) = 1 + ! T^beta_pq : exc(0,1,2) = 1 ! exc(0,2,2) = 1 ! exc(1,1,2) = q ! exc(1,2,2) = p @@ -434,6 +434,98 @@ subroutine get_single_excitation(det1,det2,exc,phase,Nint) end +subroutine get_single_excitation_cfg(cfg1,cfg2,p,q,Nint) + use bitmasks + implicit none + BEGIN_DOC + ! Returns the excitation operator between two singly excited configurations. + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: cfg1(Nint,2) + integer(bit_kind), intent(in) :: cfg2(Nint,2) + integer, intent(out) :: p, q + integer :: tz + integer :: l, ispin, idx_hole, idx_particle, ishift + integer :: nperm + integer :: i,j,k,m,n + integer :: high, low + integer :: a,b,c,d + integer(bit_kind) :: hole, particle, tmp + integer :: exc(0:2,2,2) + + ASSERT (Nint > 0) + nperm = 0 + p = 0 + q = 0 + exc(0,1,1) = 0 + exc(0,2,1) = 0 + exc(0,1,2) = 0 + exc(0,2,2) = 0 + do ispin = 1,2 + ishift = 1-bit_kind_size + do l=1,Nint + ishift = ishift + bit_kind_size + if (cfg1(l,ispin) == cfg2(l,ispin)) then + cycle + endif + tmp = xor( cfg1(l,ispin), cfg2(l,ispin) ) + particle = iand(tmp, cfg2(l,ispin)) + hole = iand(tmp, cfg1(l,ispin)) + if (particle /= 0_bit_kind) then + tz = trailz(particle) + exc(0,2,ispin) = 1 + exc(1,2,ispin) = tz+ishift + !print *,"part ",tz+ishift, " ispin=",ispin + endif + if (hole /= 0_bit_kind) then + tz = trailz(hole) + exc(0,1,ispin) = 1 + exc(1,1,ispin) = tz+ishift + !print *,"hole ",tz+ishift, " ispin=",ispin + endif + + if ( iand(exc(0,1,ispin),exc(0,2,ispin)) /= 1) then ! exc(0,1,ispin)/=1 and exc(0,2,ispin) /= 1 + cycle + endif + + high = max(exc(1,1,ispin), exc(1,2,ispin))-1 + low = min(exc(1,1,ispin), exc(1,2,ispin)) + + ASSERT (low >= 0) + ASSERT (high > 0) + + k = shiftr(high,bit_kind_shift)+1 + j = shiftr(low,bit_kind_shift)+1 + m = iand(high,bit_kind_size-1) + n = iand(low,bit_kind_size-1) + + if (j==k) then + nperm = nperm + popcnt(iand(cfg1(j,ispin), & + iand( shiftl(1_bit_kind,m)-1_bit_kind, & + not(shiftl(1_bit_kind,n))+1_bit_kind)) ) + else + nperm = nperm + popcnt( & + iand(cfg1(j,ispin), & + iand(not(0_bit_kind), & + (not(shiftl(1_bit_kind,n)) + 1_bit_kind) ))) & + + popcnt(iand(cfg1(k,ispin), & + (shiftl(1_bit_kind,m) - 1_bit_kind ) )) + + do i=j+1,k-1 + nperm = nperm + popcnt(cfg1(i,ispin)) + end do + + endif + + ! Set p and q + q = max(exc(1,1,1),exc(1,1,2)) + p = max(exc(1,2,1),exc(1,2,2)) + return + + enddo + enddo +end + subroutine bitstring_to_list_ab( string, list, n_elements, Nint) use bitmasks implicit none diff --git a/src/utils/linear_algebra.irp.f b/src/utils/linear_algebra.irp.f index 60c84a3b..1e33c7dc 100644 --- a/src/utils/linear_algebra.irp.f +++ b/src/utils/linear_algebra.irp.f @@ -1136,7 +1136,6 @@ subroutine ortho_svd(A,LDA,m,n) end -! QR to orthonormalize CSFs does not work :-( !subroutine ortho_qr_withB(A,LDA,B,m,n) ! implicit none ! BEGIN_DOC @@ -1223,7 +1222,7 @@ end ! ! !deallocate(WORK,TAU) !end - +! !subroutine ortho_qr_csf(A, LDA, B, m, n) bind(C, name="ortho_qr_csf") ! use iso_c_binding ! integer(c_int32_t), value :: LDA @@ -1234,6 +1233,7 @@ end ! call ortho_qr_withB(A,LDA,B,m,n) !end subroutine ortho_qr_csf + subroutine ortho_qr(A,LDA,m,n) implicit none BEGIN_DOC