diff --git a/src/csf/configuration_CI_sigma_helpers.irp.f b/src/csf/configuration_CI_sigma_helpers.irp.f index f73362eb..0290a5f7 100644 --- a/src/csf/configuration_CI_sigma_helpers.irp.f +++ b/src/csf/configuration_CI_sigma_helpers.irp.f @@ -1,3 +1,272 @@ + subroutine obtain_associated_alphaI(idxI, Icfg, alphasIcfg, NalphaIcfg) + implicit none + use bitmasks + BEGIN_DOC + ! Documentation for alphasI + ! Returns the associated alpha's for + ! the input configuration Icfg. + END_DOC + + integer,intent(in) :: idxI ! The id of the Ith CFG + integer(bit_kind),intent(in) :: Icfg(N_int,2) + integer,intent(out) :: NalphaIcfg + integer(bit_kind),intent(out) :: alphasIcfg(N_int,2,*) + logical,dimension(:,:),allocatable :: tableUniqueAlphas + integer :: listholes(mo_num) + integer :: holetype(mo_num) ! 1-> SOMO 2->DOMO + integer :: nholes + integer :: nvmos + integer :: listvmos(mo_num) + integer :: vmotype(mo_num) ! 1 -> VMO 2 -> SOMO + integer*8 :: Idomo + integer*8 :: Isomo + integer*8 :: Jdomo + integer*8 :: Jsomo + integer*8 :: diffSOMO + integer*8 :: diffDOMO + integer*8 :: xordiffSOMODOMO + integer :: ndiffSOMO + integer :: ndiffDOMO + integer :: nxordiffSOMODOMO + integer :: ndiffAll + integer :: i + integer :: j + integer :: k + integer :: hole + integer :: p + integer :: q + integer :: countalphas + logical :: pqAlreadyGenQ + logical :: pqExistsQ + logical :: ppExistsQ + Isomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,1)) + Idomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,2)) + !print*,"Input cfg" + !call debug_spindet(Isomo,1) + !call debug_spindet(Idomo,1) + + !print*,n_act_orb, "monum=",mo_num," n_core=",n_core_orb + + ! find out all pq holes possible + nholes = 0 + ! holes in SOMO + do i = 1,mo_num + if(POPCNT(IAND(Isomo,IBSET(0_8,i-1))) .EQ. 1) then + nholes += 1 + listholes(nholes) = i + holetype(nholes) = 1 + endif + end do + ! holes in DOMO + do i = 1,mo_num + if(POPCNT(IAND(Idomo,IBSET(0_8,i-1))) .EQ. 1) then + nholes += 1 + listholes(nholes) = i + holetype(nholes) = 2 + endif + end do + + ! find vmos + listvmos = -1 + vmotype = -1 + nvmos = 0 + do i = 1,mo_num + !print *,i,IBSET(0,i-1),POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))), POPCNT(IAND(Idomo,(IBSET(0_8,i-1)))) + if(POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))) .EQ. 0 .AND. POPCNT(IAND(Idomo,(IBSET(0_8,i-1)))) .EQ. 0) then + nvmos += 1 + listvmos(nvmos) = i + vmotype(nvmos) = 1 + else if(POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))) .EQ. 1 .AND. POPCNT(IAND(Idomo,(IBSET(0_8,i-1)))) .EQ. 0 ) then + nvmos += 1 + listvmos(nvmos) = i + vmotype(nvmos) = 2 + end if + end do + + !print *,"Nvmo=",nvmos + !print *,listvmos + !print *,vmotype + + allocate(tableUniqueAlphas(mo_num,mo_num)) + tableUniqueAlphas = .FALSE. + + ! Now find the allowed (p,q) excitations + Isomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,1)) + Idomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,2)) + !print *,"Isomo" + !call debug_spindet(Isomo,1) + !call debug_spindet(Idomo,1) + + !print *,"Nholes=",nholes," Nvmos=",nvmos, " idxi=",idxI + !do i = 1,nholes + ! print *,i,"->",listholes(i) + !enddo + !do i = 1,nvmos + ! print *,i,"->",listvmos(i) + !enddo + + ! TODO cfg_seniority_index + do i = 1,nholes + p = listholes(i) + do j = 1,nvmos + q = listvmos(j) + if(p .EQ. q) cycle + if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 1) then + ! SOMO -> VMO + Jsomo = IBCLR(Isomo,p-1) + Jsomo = IBSET(Jsomo,q-1) + Jdomo = Idomo + else if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 2) then + ! SOMO -> SOMO + Jsomo = IBCLR(Isomo,p-1) + Jsomo = IBCLR(Jsomo,q-1) + Jdomo = IBSET(Idomo,q-1) + else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then + ! DOMO -> VMO + Jsomo = IBSET(Isomo,p-1) + Jsomo = IBSET(Jsomo,q-1) + Jdomo = IBCLR(Idomo,p-1) + else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 2) then + ! DOMO -> SOMO + Jsomo = IBSET(Isomo,p-1) + Jsomo = IBCLR(Jsomo,q-1) + Jdomo = IBCLR(Idomo,p-1) + Jdomo = IBSET(Jdomo,q-1) + else + print*,"Something went wrong in obtain_associated_alphaI" + endif + + + pqAlreadyGenQ = .FALSE. + ! First check if it can be generated before + do k = 1, idxI-1 + diffSOMO = IEOR(Jsomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k))) + diffDOMO = IEOR(Jdomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k))) + xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) + ndiffSOMO = POPCNT(diffSOMO) + ndiffDOMO = POPCNT(diffDOMO) + nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO) + !if(POPCNT(IEOR(diffSOMO,diffDOMO)) .LE. 1 .AND. ndiffDOMO .LT. 3) then + if((ndiffSOMO+ndiffDOMO) .EQ. 0) then + pqAlreadyGenQ = .TRUE. + ppExistsQ = .TRUE. + EXIT + endif + if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then + pqAlreadyGenQ = .TRUE. + !ppExistsQ = .TRUE. + !print *,i,k,ndiffSOMO,ndiffDOMO + !call debug_spindet(Jsomo,1) + !call debug_spindet(Jdomo,1) + !call debug_spindet(iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k)),1) + !call debug_spindet(iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k)),1) + EXIT + endif + end do + + !print *,"(,",p,",",q,")",pqAlreadyGenQ + + if(pqAlreadyGenQ) cycle + + pqExistsQ = .FALSE. + ! now check if this exists in the selected list + !do k = idxI+1, N_configuration + ! diffSOMO = IEOR(OR(reunion_of_act_virt_bitmask(1,1),Jsomo),psi_configuration(1,1,k)) + ! diffDOMO = IEOR(OR(reunion_of_act_virt_bitmask(1,1),Jdomo),psi_configuration(1,2,k)) + ! ndiffSOMO = POPCNT(diffSOMO) + ! ndiffDOMO = POPCNT(diffDOMO) + ! if((ndiffSOMO + ndiffDOMO) .EQ. 0) then + ! pqExistsQ = .TRUE. + ! EXIT + ! endif + !end do + + if(.NOT. pqExistsQ) then + tableUniqueAlphas(p,q) = .TRUE. + !print *,p,q + !call debug_spindet(Jsomo,1) + !call debug_spindet(Jdomo,1) + endif + end do + end do + + !print *,tableUniqueAlphas(:,:) + + ! prune list of alphas + Isomo = Icfg(1,1) + Idomo = Icfg(1,2) + Jsomo = Icfg(1,1) + Jdomo = Icfg(1,2) + NalphaIcfg = 0 + do i = 1, nholes + p = listholes(i) + do j = 1, nvmos + q = listvmos(j) + if(p .EQ. q) cycle + if(tableUniqueAlphas(p,q)) then + if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 1) then + ! SOMO -> VMO + Jsomo = IBCLR(Isomo,p-1) + Jsomo = IBSET(Jsomo,q-1) + Jdomo = Idomo + else if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 2) then + ! SOMO -> SOMO + Jsomo = IBCLR(Isomo,p-1) + Jsomo = IBCLR(Jsomo,q-1) + Jdomo = IBSET(Idomo,q-1) + else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then + ! DOMO -> VMO + Jsomo = IBSET(Isomo,p-1) + Jsomo = IBSET(Jsomo,q-1) + Jdomo = IBCLR(Idomo,p-1) + else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 2) then + ! DOMO -> SOMO + Jsomo = IBSET(Isomo,p-1) + Jsomo = IBCLR(Jsomo,q-1) + Jdomo = IBCLR(Idomo,p-1) + Jdomo = IBSET(Jdomo,q-1) + else + print*,"Something went wrong in obtain_associated_alphaI" + endif + + ! SOMO + NalphaIcfg += 1 + !print *,p,q,"|",holetype(i),vmotype(j),NalphaIcfg + alphasIcfg(1,1,NalphaIcfg) = Jsomo + alphasIcfg(1,2,NalphaIcfg) = IOR(Jdomo,ISHFT(1_8,n_core_orb)-1) + endif + end do + end do + + ! Check if this Icfg has been previously generated as a mono + ppExistsQ = .False. + Isomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,1)) + Idomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,2)) + do k = 1, idxI-1 + diffSOMO = IEOR(Isomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k))) + diffDOMO = IEOR(Idomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k))) + xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) + ndiffSOMO = POPCNT(diffSOMO) + ndiffDOMO = POPCNT(diffDOMO) + nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO) + if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then + ppExistsQ = .TRUE. + EXIT + endif + end do + ! Diagonal part (pp,qq) + if(nholes > 0 .AND. (.NOT. ppExistsQ))then + ! SOMO + NalphaIcfg += 1 + !print *,p,q,"|",holetype(i),vmotype(j),NalphaIcfg + !call debug_spindet(Idomo,1) + !call debug_spindet(Jdomo,1) + alphasIcfg(1,1,NalphaIcfg) = Icfg(1,1) + alphasIcfg(1,2,NalphaIcfg) = Icfg(1,2) + endif + + end subroutine + function getNSOMO(Icfg) result(NSOMO) implicit none integer(bit_kind),intent(in) :: Icfg(N_int,2) @@ -23,19 +292,26 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod ! Type 3 - SOMO -> VMO ! Type 4 - DOMO -> SOMO END_DOC - integer(bit_kind),intent(in) :: Ialpha(N_int,2) - integer(bit_kind),intent(in) :: Jcfg(N_int,2) - integer,intent(in) :: p,q - integer,intent(in) :: extype - integer,intent(out) :: pmodel,qmodel - integer*8 :: Isomo - integer*8 :: Idomo - integer*8 :: Jsomo - integer*8 :: Jdomo - integer*8 :: mask - integer*8 :: Isomotmp - integer*8 :: Jsomotmp - integer :: pos0,pos0prev + integer(bit_kind),intent(in) :: Ialpha(N_int,2) + integer(bit_kind),intent(in) :: Jcfg(N_int,2) + integer,intent(in) :: p,q + integer,intent(in) :: extype + integer,intent(out) :: pmodel,qmodel + !integer(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(bit_kind) :: Isomotmp(N_int) + !integer(bit_kind) :: Jsomotmp(N_int) + integer*8 :: Isomotmp + integer*8 :: Jsomotmp + integer :: pos0,pos0prev ! TODO Flag (print) when model space indices is > 64 Isomo = Ialpha(1,1) @@ -50,11 +326,6 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod pmodel = 1 qmodel = 1 else - !print *,"input pq=",p,q,"extype=",extype - !call debug_spindet(Isomo,1) - !call debug_spindet(Idomo,1) - !call debug_spindet(Jsomo,1) - !call debug_spindet(Jdomo,1) select case(extype) case (1) ! SOMO -> SOMO @@ -80,26 +351,44 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod ! SOMO -> VMO !print *,"type -> SOMO -> VMO" !Isomo = IEOR(Isomo,Jsomo) - mask = ISHFT(1_8,p) - 1 - Isomo = IAND(Isomo,mask) - pmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) - mask = ISHFT(1_8,q) - 1 - Jsomo = IAND(Jsomo,mask) - qmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) + if(p.LT.q) then + mask = ISHFT(1_8,p) - 1 + Isomo = IAND(Isomo,mask) + pmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) + mask = ISHFT(1_8,q) - 1 + Jsomo = IAND(Jsomo,mask) + qmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) + 1 + else + mask = ISHFT(1_8,p) - 1 + Isomo = IAND(Isomo,mask) + pmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) + 1 + mask = ISHFT(1_8,q) - 1 + Jsomo = IAND(Jsomo,mask) + qmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) + endif case (4) ! DOMO -> SOMO ! remove all domos except one at p !print *,"type -> DOMO -> SOMO" !Isomo = IEOR(Isomo,Jsomo) - mask = ISHFT(1_8,p) - 1 - Jsomo = IAND(Jsomo,mask) - pmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) - mask = ISHFT(1_8,q) - 1 - Isomo = IAND(Isomo,mask) - qmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) + if(p.LT.q) then + mask = ISHFT(1_8,p) - 1 + Jsomo = IAND(Jsomo,mask) + pmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) + mask = ISHFT(1_8,q) - 1 + Isomo = IAND(Isomo,mask) + qmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) + 1 + else + mask = ISHFT(1_8,p) - 1 + Jsomo = IAND(Jsomo,mask) + pmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) + 1 + mask = ISHFT(1_8,q) - 1 + Isomo = IAND(Isomo,mask) + qmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) + endif case default print *,"something is wrong in convertOrbIdsToModelSpaceIds" end select endif - !print *,p,q,"model ids=",pmodel,qmodel + !print *,p,q,"model ids=",pmodel,qmodel end subroutine convertOrbIdsToModelSpaceIds diff --git a/src/csf/configuration_CI_sigma_helpers.org b/src/csf/configuration_CI_sigma_helpers.org index ede42b97..88b27a7a 100644 --- a/src/csf/configuration_CI_sigma_helpers.org +++ b/src/csf/configuration_CI_sigma_helpers.org @@ -43,8 +43,10 @@ the input determinant \(|D_I\rangle\). integer*8 :: Jsomo integer*8 :: diffSOMO integer*8 :: diffDOMO + integer*8 :: xordiffSOMODOMO integer :: ndiffSOMO integer :: ndiffDOMO + integer :: nxordiffSOMODOMO integer :: ndiffAll integer :: i integer :: j @@ -55,6 +57,7 @@ the input determinant \(|D_I\rangle\). integer :: countalphas logical :: pqAlreadyGenQ logical :: pqExistsQ + logical :: ppExistsQ Isomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,1)) Idomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,2)) !print*,"Input cfg" @@ -66,7 +69,7 @@ the input determinant \(|D_I\rangle\). ! find out all pq holes possible nholes = 0 ! holes in SOMO - do i = n_core_orb+1,n_core_orb + n_act_orb + do i = 1,mo_num if(POPCNT(IAND(Isomo,IBSET(0_8,i-1))) .EQ. 1) then nholes += 1 listholes(nholes) = i @@ -74,7 +77,7 @@ the input determinant \(|D_I\rangle\). endif end do ! holes in DOMO - do i = n_core_orb+1,n_core_orb + n_act_orb + do i = 1,mo_num if(POPCNT(IAND(Idomo,IBSET(0_8,i-1))) .EQ. 1) then nholes += 1 listholes(nholes) = i @@ -86,7 +89,7 @@ the input determinant \(|D_I\rangle\). listvmos = -1 vmotype = -1 nvmos = 0 - do i = n_core_orb+1,n_core_orb + n_act_orb + do i = 1,mo_num !print *,i,IBSET(0,i-1),POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))), POPCNT(IAND(Idomo,(IBSET(0_8,i-1)))) if(POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))) .EQ. 0 .AND. POPCNT(IAND(Idomo,(IBSET(0_8,i-1)))) .EQ. 0) then nvmos += 1 @@ -126,7 +129,7 @@ the input determinant \(|D_I\rangle\). p = listholes(i) do j = 1,nvmos q = listvmos(j) - if(p == q) cycle + if(p .EQ. q) cycle if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 1) then ! SOMO -> VMO Jsomo = IBCLR(Isomo,p-1) @@ -158,10 +161,19 @@ the input determinant \(|D_I\rangle\). do k = 1, idxI-1 diffSOMO = IEOR(Jsomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k))) diffDOMO = IEOR(Jdomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k))) + xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) ndiffSOMO = POPCNT(diffSOMO) ndiffDOMO = POPCNT(diffDOMO) - if(POPCNT(IEOR(diffSOMO,diffDOMO)) .LE. 1 .AND. ndiffDOMO .LT. 3) then + nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO) + !if(POPCNT(IEOR(diffSOMO,diffDOMO)) .LE. 1 .AND. ndiffDOMO .LT. 3) then + if((ndiffSOMO+ndiffDOMO) .EQ. 0) then pqAlreadyGenQ = .TRUE. + ppExistsQ = .TRUE. + EXIT + endif + if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then + pqAlreadyGenQ = .TRUE. + !ppExistsQ = .TRUE. !print *,i,k,ndiffSOMO,ndiffDOMO !call debug_spindet(Jsomo,1) !call debug_spindet(Jdomo,1) @@ -171,20 +183,22 @@ the input determinant \(|D_I\rangle\). endif end do + !print *,"(,",p,",",q,")",pqAlreadyGenQ + if(pqAlreadyGenQ) cycle pqExistsQ = .FALSE. ! now check if this exists in the selected list - do k = idxI, N_configuration - diffSOMO = IEOR(OR(reunion_of_act_virt_bitmask(1,1),Jsomo),psi_configuration(1,1,k)) - diffDOMO = IEOR(OR(reunion_of_act_virt_bitmask(1,1),Jdomo),psi_configuration(1,2,k)) - ndiffSOMO = POPCNT(diffSOMO) - ndiffDOMO = POPCNT(diffDOMO) - if((ndiffSOMO + ndiffDOMO) .EQ. 0) then - pqExistsQ = .TRUE. - EXIT - endif - end do + !do k = idxI+1, N_configuration + ! diffSOMO = IEOR(OR(reunion_of_act_virt_bitmask(1,1),Jsomo),psi_configuration(1,1,k)) + ! diffDOMO = IEOR(OR(reunion_of_act_virt_bitmask(1,1),Jdomo),psi_configuration(1,2,k)) + ! ndiffSOMO = POPCNT(diffSOMO) + ! ndiffDOMO = POPCNT(diffDOMO) + ! if((ndiffSOMO + ndiffDOMO) .EQ. 0) then + ! pqExistsQ = .TRUE. + ! EXIT + ! endif + !end do if(.NOT. pqExistsQ) then tableUniqueAlphas(p,q) = .TRUE. @@ -234,16 +248,42 @@ the input determinant \(|D_I\rangle\). print*,"Something went wrong in obtain_associated_alphaI" endif + ! SOMO NalphaIcfg += 1 !print *,p,q,"|",holetype(i),vmotype(j),NalphaIcfg - !call debug_spindet(Idomo,1) - !call debug_spindet(Jdomo,1) alphasIcfg(1,1,NalphaIcfg) = Jsomo alphasIcfg(1,2,NalphaIcfg) = IOR(Jdomo,ISHFT(1_8,n_core_orb)-1) endif end do end do + ! Check if this Icfg has been previously generated as a mono + ppExistsQ = .False. + Isomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,1)) + Idomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,2)) + do k = 1, idxI-1 + diffSOMO = IEOR(Isomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k))) + diffDOMO = IEOR(Idomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k))) + xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) + ndiffSOMO = POPCNT(diffSOMO) + ndiffDOMO = POPCNT(diffDOMO) + nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO) + if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then + ppExistsQ = .TRUE. + EXIT + endif + end do + ! Diagonal part (pp,qq) + if(nholes > 0 .AND. (.NOT. ppExistsQ))then + ! SOMO + NalphaIcfg += 1 + !print *,p,q,"|",holetype(i),vmotype(j),NalphaIcfg + !call debug_spindet(Idomo,1) + !call debug_spindet(Jdomo,1) + alphasIcfg(1,1,NalphaIcfg) = Icfg(1,1) + alphasIcfg(1,2,NalphaIcfg) = Icfg(1,2) + endif + end subroutine #+end_src @@ -253,7 +293,7 @@ Next step is to obtain the connected CFGs \(|I\rangle\) that belong to the selec given a RI configuration \(|\alpha\rangle\). #+begin_src f90 :main no :tangle ../cfgCI/obtain_I_foralpha.irp.f -subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI, nconnectedI, excitationIds, excitationTypes) +subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI, nconnectedI, excitationIds, excitationTypes, diagfactors) implicit none use bitmasks BEGIN_DOC @@ -278,6 +318,7 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI integer,intent(out) :: nconnectedI integer,intent(out) :: excitationIds(2,*) integer,intent(out) :: excitationTypes(*) + real*8 ,intent(out) :: diagfactors(*) integer*8 :: Idomo integer*8 :: Isomo integer*8 :: Jdomo @@ -285,38 +326,19 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI integer*8 :: IJsomo integer*8 :: diffSOMO integer*8 :: diffDOMO + integer*8 :: xordiffSOMODOMO integer :: ndiffSOMO integer :: ndiffDOMO - integer :: i,j,k,l,p,q,nsomoJ,nsomoalpha,starti,endi,extyp,nholes + integer :: nxordiffSOMODOMO + integer :: ii,i,j,k,l,p,q,nsomoJ,nsomoalpha,starti,endi,extyp,nholes integer :: listholes(mo_num) integer :: holetype(mo_num) - ! find out all pq holes possible - nholes = 0 - ! holes in SOMO - Isomo = psi_configuration(1,1,idxI) - Idomo = psi_configuration(1,2,idxI) - do i = n_core_orb+1,n_core_orb + n_act_orb - if(POPCNT(IAND(Isomo,IBSET(0_8,i-1))) .EQ. 1) then - nholes += 1 - listholes(nholes) = i - holetype(nholes) = 1 - endif - end do - ! holes in DOMO - do i = n_core_orb+1,n_core_orb + n_act_orb - if(POPCNT(IAND(Idomo,IBSET(0_8,i-1))) .EQ. 1) then - nholes += 1 - listholes(nholes) = i - holetype(nholes) = 2 - endif - end do - nconnectedI = 0 p = 0 q = 0 - do i=idxI+1,N_configuration + do i=idxI,N_configuration Isomo = Ialpha(1,1) Idomo = Ialpha(1,2) Jsomo = psi_configuration(1,1,i) @@ -328,12 +350,15 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI !call debug_spindet(Jdomo,1) diffSOMO = IEOR(Isomo,Jsomo) diffDOMO = IEOR(Idomo,Jdomo) + xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) ndiffSOMO = POPCNT(diffSOMO) ndiffDOMO = POPCNT(diffDOMO) - if((ndiffSOMO + ndiffDOMO) .EQ. 0) cycle - !print *,"-I--i=",i,diffSOMO,diffDOMO!Isomo,Jsomo,ndiffSOMO,ndiffDOMO + nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO) + !print *,"-I--i=",i,ndiffSOMO,ndiffDOMO,nxordiffSOMODOMO!Isomo,Jsomo,ndiffSOMO,ndiffDOMO + !if((ndiffSOMO + ndiffDOMO) .EQ. 0) cycle !print *,POPCNT(IEOR(diffSOMO,diffDOMO)), ndiffDOMO - if(POPCNT(IEOR(diffSOMO,diffDOMO)) .LE. 1 .AND. ndiffDOMO .LT. 3) then + !if(POPCNT(IEOR(diffSOMO,diffDOMO)) .LE. 1 .AND. ndiffDOMO .LT. 3) then + if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then !call debug_spindet(Isomo,1) !call debug_spindet(Idomo,1) !print *,"-J--i=",i,Idomo,Jdomo,">",N_configuration @@ -390,7 +415,58 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI excitationIds(1,nconnectedI)=p excitationIds(2,nconnectedI)=q excitationTypes(nconnectedI) = extyp - print *,"------ > output p,q in obt=",p,q + diagfactors(nconnectedI) = 1.0d0 + !print *,"------ > output p,q in obt=",p,q + else if((ndiffSOMO + ndiffDOMO) .EQ. 0) then + ! find out all pq holes possible + nholes = 0 + ! holes in SOMO + Isomo = psi_configuration(1,1,i) + Idomo = psi_configuration(1,2,i) + do ii = 1,mo_num + if(POPCNT(IAND(Isomo,IBSET(0_8,ii-1))) .EQ. 1) then + nholes += 1 + listholes(nholes) = ii + holetype(nholes) = 1 + endif + end do + ! holes in DOMO + do ii = 1,mo_num + if(POPCNT(IAND(Idomo,IBSET(0_8,ii-1))) .EQ. 1) then + nholes += 1 + listholes(nholes) = ii + holetype(nholes) = 2 + endif + end do + + do k=1,nholes + p = listholes(k) + q = p + extyp = 1 + if(holetype(k) .EQ. 1) then + starti = psi_config_data(i,1) + endi = psi_config_data(i,2) + nconnectedI += 1 + connectedI(:,:,nconnectedI) = psi_configuration(:,:,i) + idxs_connectedI(nconnectedI)=starti + excitationIds(1,nconnectedI)=p + excitationIds(2,nconnectedI)=q + excitationTypes(nconnectedI) = extyp + diagfactors(nconnectedI) = 1.0d0 + !print *,"------ > output p,q in obt=",p,q + else + starti = psi_config_data(i,1) + endi = psi_config_data(i,2) + nconnectedI += 1 + connectedI(:,:,nconnectedI) = psi_configuration(:,:,i) + idxs_connectedI(nconnectedI)=starti + excitationIds(1,nconnectedI)=p + excitationIds(2,nconnectedI)=q + excitationTypes(nconnectedI) = extyp + diagfactors(nconnectedI) = 2.0d0 + !print *,"------ > output p,q in obt=",p,q + endif + enddo endif end do @@ -468,6 +544,7 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod qmodel = q if(p .EQ. q) then + !print *,"input pq=",p,q,"extype=",extype pmodel = 1 qmodel = 1 else @@ -501,28 +578,46 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod ! SOMO -> VMO !print *,"type -> SOMO -> VMO" !Isomo = IEOR(Isomo,Jsomo) - mask = ISHFT(1_8,p) - 1 - Isomo = IAND(Isomo,mask) - pmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) - mask = ISHFT(1_8,q) - 1 - Jsomo = IAND(Jsomo,mask) - qmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) + if(p.LT.q) then + mask = ISHFT(1_8,p) - 1 + Isomo = IAND(Isomo,mask) + pmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) + mask = ISHFT(1_8,q) - 1 + Jsomo = IAND(Jsomo,mask) + qmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) + 1 + else + mask = ISHFT(1_8,p) - 1 + Isomo = IAND(Isomo,mask) + pmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) + 1 + mask = ISHFT(1_8,q) - 1 + Jsomo = IAND(Jsomo,mask) + qmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) + endif case (4) ! DOMO -> SOMO ! remove all domos except one at p !print *,"type -> DOMO -> SOMO" !Isomo = IEOR(Isomo,Jsomo) - mask = ISHFT(1_8,p) - 1 - Jsomo = IAND(Jsomo,mask) - pmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) - mask = ISHFT(1_8,q) - 1 - Isomo = IAND(Isomo,mask) - qmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) + if(p.LT.q) then + mask = ISHFT(1_8,p) - 1 + Jsomo = IAND(Jsomo,mask) + pmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) + mask = ISHFT(1_8,q) - 1 + Isomo = IAND(Isomo,mask) + qmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) + 1 + else + mask = ISHFT(1_8,p) - 1 + Jsomo = IAND(Jsomo,mask) + pmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) + 1 + mask = ISHFT(1_8,q) - 1 + Isomo = IAND(Isomo,mask) + qmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) + endif case default print *,"something is wrong in convertOrbIdsToModelSpaceIds" end select endif - !print *,p,q,"model ids=",pmodel,qmodel + !print *,p,q,"model ids=",pmodel,qmodel end subroutine convertOrbIdsToModelSpaceIds #+end_src diff --git a/src/csf/conversion.irp.f b/src/csf/conversion.irp.f index fecc6123..7625fa8f 100644 --- a/src/csf/conversion.irp.f +++ b/src/csf/conversion.irp.f @@ -10,7 +10,7 @@ subroutine convertWFfromDETtoCSF(N_st,psi_coef_det_in, psi_coef_cfg_out) integer, intent(in) :: N_st double precision, intent(in) :: psi_coef_det_in(N_det,N_st) double precision, intent(out) :: psi_coef_cfg_out(n_CSF,N_st) - integer*8 :: Isomo, Idomo, mask + integer*8 :: Isomo, Idomo integer(bit_kind) :: Ialpha(N_int) ,Ibeta(N_int) integer :: rows, cols, i, j, k integer :: startdet, enddet @@ -20,10 +20,10 @@ subroutine convertWFfromDETtoCSF(N_st,psi_coef_det_in, psi_coef_cfg_out) double precision,allocatable :: tempCoeff(:,:) double precision :: phasedet integer :: idx - + ! initialization psi_coef_cfg_out(:,1) = 0.d0 - + integer s, bfIcfg integer countcsf countcsf = 0 @@ -32,7 +32,7 @@ subroutine convertWFfromDETtoCSF(N_st,psi_coef_det_in, psi_coef_cfg_out) startdet = psi_configuration_to_psi_det(1,i) enddet = psi_configuration_to_psi_det(2,i) ndetI = enddet-startdet+1 - + allocate(tempCoeff(ndetI,N_st)) do j = startdet, enddet idx = psi_configuration_to_psi_det_data(j) @@ -43,29 +43,29 @@ subroutine convertWFfromDETtoCSF(N_st,psi_coef_det_in, psi_coef_cfg_out) tempCoeff(j-startdet+1,k) = psi_coef_det_in(idx, k)*phasedet enddo enddo - + s = 0 do k=1,N_int if (psi_configuration(k,1,i) == 0_bit_kind) cycle s = s + popcnt(psi_configuration(k,1,i)) enddo bfIcfg = max(1,nint((binom(s,(s+1)/2)-binom(s,((s+1)/2)+1)))) - + ! perhaps blocking with CFGs of same seniority ! can be more efficient allocate(tempBuffer(bfIcfg,ndetI)) tempBuffer = DetToCSFTransformationMatrix(s,:bfIcfg,:ndetI) - + call dgemm('N','N', bfIcfg, N_st, ndetI, 1.d0, tempBuffer, size(tempBuffer,1),& tempCoeff, size(tempCoeff,1), 0.d0, psi_coef_cfg_out(countcsf+1,1),& size(psi_coef_cfg_out,1)) - + deallocate(tempCoeff) deallocate(tempBuffer) countcsf += bfIcfg enddo - -end + +end subroutine convertWFfromDETtoCSF subroutine convertWFfromCSFtoDET(N_st,psi_coef_cfg_in, psi_coef_det) @@ -88,42 +88,42 @@ subroutine convertWFfromCSFtoDET(N_st,psi_coef_cfg_in, psi_coef_det) integer :: ndetI integer :: getNSOMO double precision,allocatable :: tempBuffer(:,:) - double precision,allocatable :: tempCoeff (:,:) + double precision,allocatable :: tempCoeff(:,:) double precision :: phasedet integer :: idx - - countcsf = 0 - + + countcsf = 0 + do i = 1,N_configuration startdet = psi_configuration_to_psi_det(1,i) enddet = psi_configuration_to_psi_det(2,i) ndetI = enddet-startdet+1 - + s = 0 do k=1,N_int if (psi_configuration(k,1,i) == 0_bit_kind) cycle s = s + popcnt(psi_configuration(k,1,i)) enddo bfIcfg = max(1,nint((binom(s,(s+1)/2)-binom(s,((s+1)/2)+1)))) - + allocate(tempCoeff(bfIcfg,N_st)) - + do k=1,N_st do j = 1,bfIcfg tempCoeff(j,k) = psi_coef_cfg_in(countcsf+j,k) enddo enddo - + countcsf += bfIcfg ! perhaps blocking with CFGs of same seniority ! can be more efficient allocate(tempBuffer(bfIcfg,ndetI)) tempBuffer = DetToCSFTransformationMatrix(s,:bfIcfg,:ndetI) - + call dgemm('T','N', ndetI, N_st, bfIcfg, 1.d0, tempBuffer, size(tempBuffer,1),& tempCoeff, size(tempCoeff,1), 0.d0, tmp_psi_coef_det, & size(tmp_psi_coef_det,1)) - + do j=startdet,enddet idx = psi_configuration_to_psi_det_data(j) Ialpha(:) = psi_det(:,1,idx) @@ -133,16 +133,9 @@ subroutine convertWFfromCSFtoDET(N_st,psi_coef_cfg_in, psi_coef_det) psi_coef_det(idx,k) = tmp_psi_coef_det(j-startdet+1,k) * phasedet enddo enddo - + deallocate(tempCoeff) deallocate(tempBuffer) enddo -end - - - - - - - +end subroutine convertCSFtoDET diff --git a/src/csf/sigma_vector.irp.f b/src/csf/sigma_vector.irp.f index 9c1ac56a..1a434aea 100644 --- a/src/csf/sigma_vector.irp.f +++ b/src/csf/sigma_vector.irp.f @@ -1,9 +1,9 @@ - BEGIN_PROVIDER [ integer, NSOMOMax] -&BEGIN_PROVIDER [ integer, NCSFMax] -&BEGIN_PROVIDER [ integer*8, NMO] -&BEGIN_PROVIDER [ integer, NBFMax] -&BEGIN_PROVIDER [ integer, n_CSF] -&BEGIN_PROVIDER [ integer, maxDetDimPerBF] + BEGIN_PROVIDER [ integer, NSOMOMax] + &BEGIN_PROVIDER [ integer, NCSFMax] + &BEGIN_PROVIDER [ integer*8, NMO] + &BEGIN_PROVIDER [ integer, NBFMax] + &BEGIN_PROVIDER [ integer, n_CSF] + &BEGIN_PROVIDER [ integer, maxDetDimPerBF] implicit none BEGIN_DOC ! Documentation for NSOMOMax @@ -22,29 +22,38 @@ integer NSOMO integer dimcsfpercfg integer detDimperBF - real*8 :: coeff + real*8 :: coeff integer MS integer ncfgpersomo detDimperBF = 0 MS = elec_alpha_num-elec_beta_num + !print *,"NSOMOMax=",NSOMOMax, cfg_seniority_index(0) ! number of cfgs = number of dets for 0 somos n_CSF = cfg_seniority_index(0)-1 ncfgprev = cfg_seniority_index(0) do i = 0-iand(MS,1)+2, NSOMOMax,2 - if(cfg_seniority_index(i) .EQ. -1)then - ncfgpersomo = N_configuration + 1 - else - ncfgpersomo = cfg_seniority_index(i) - endif - ncfg = ncfgpersomo - ncfgprev - !detDimperBF = max(1,nint((binom(i,(i+1)/2)))) - dimcsfpercfg = max(1,nint((binom(i-2,(i-2+1)/2)-binom(i-2,((i-2+1)/2)+1)))) - n_CSF += ncfg * dimcsfpercfg - !if(cfg_seniority_index(i+2) == -1) EXIT - !if(detDimperBF > maxDetDimPerBF) maxDetDimPerBF = detDimperBF - ncfgprev = cfg_seniority_index(i) + if(cfg_seniority_index(i) .EQ. -1)then + ncfgpersomo = N_configuration + 1 + else + ncfgpersomo = cfg_seniority_index(i) + endif + ncfg = ncfgpersomo - ncfgprev + !detDimperBF = max(1,nint((binom(i,(i+1)/2)))) + dimcsfpercfg = max(1,nint((binom(i-2,(i-2+1)/2)-binom(i-2,((i-2+1)/2)+1)))) + n_CSF += ncfg * dimcsfpercfg + !print *,i,">(",ncfg,ncfgprev,ncfgpersomo,")",",",detDimperBF,">",dimcsfpercfg, " | dimbas= ", n_CSF + !if(cfg_seniority_index(i+2) == -1) EXIT + !if(detDimperBF > maxDetDimPerBF) maxDetDimPerBF = detDimperBF + ncfgprev = cfg_seniority_index(i) enddo -END_PROVIDER + if(NSOMOMax .EQ. elec_num)then + ncfgpersomo = N_configuration + 1 + ncfg = ncfgpersomo - ncfgprev + dimcsfpercfg = max(1,nint((binom(i-2,(i-2+1)/2)-binom(i-2,((i-2+1)/2)+1)))) + n_CSF += ncfg * dimcsfpercfg + !print *,i,">(",ncfg,ncfgprev,ncfgpersomo,")",",",detDimperBF,">",dimcsfpercfg, " | dimbas= ", n_CSF + endif + END_PROVIDER subroutine get_phase_qp_to_cfg(Ialpha, Ibeta, phaseout) use bitmasks @@ -61,14 +70,13 @@ subroutine get_phase_qp_to_cfg(Ialpha, Ibeta, phaseout) integer(bit_kind),intent(in) :: Ialpha(N_int) integer(bit_kind),intent(in) :: Ibeta(N_int) real*8,intent(out) :: phaseout - integer(bit_kind) :: mask, mask2(N_int), deta(N_int), detb(N_int) + integer(bit_kind) :: mask, deta(N_int), detb(N_int) integer :: nbetas integer :: count, k - ! Remove the DOMOs - mask2 = IAND(Ialpha,Ibeta) - deta = IEOR(Ialpha,mask2) - detb = IEOR(Ibeta ,mask2) + ! Initliaze deta and detb + deta = Ialpha + detb = Ibeta ! Find how many alpha electrons there are in all the N_ints integer :: Na(N_int) @@ -82,12 +90,12 @@ subroutine get_phase_qp_to_cfg(Ialpha, Ibeta, phaseout) do while(detb(k) /= 0_bit_kind) ! Find the lowest beta electron and clear it - ipos = trailz(detb(k)) - detb(k) = ibclr(detb(k),ipos) + ipos = trailz(detb(k)) + detb(k) = ibclr(detb(k),ipos) ! Create a mask will all MOs higher than the beta electron - mask = not(shiftl(1_bit_kind,ipos+1) - 1_bit_kind) - + mask = not(shiftl(1_bit_kind,ipos + 1) - 1_bit_kind) + ! Apply the mask to the alpha string to count how many electrons to cross nperm = popcnt( iand(mask, deta(k)) ) @@ -104,7 +112,111 @@ end subroutine get_phase_qp_to_cfg - BEGIN_PROVIDER [ integer, AIJpqMatrixDimsList, (0:NSOMOMax,0:NSOMOMax,4,NSOMOMax,NSOMOMax,2)] + BEGIN_PROVIDER [ real*8, DetToCSFTransformationMatrix, (0:NSOMOMax,NBFMax,maxDetDimPerBF)] + &BEGIN_PROVIDER [ real*8, psi_coef_config, (n_CSF,1)] + &BEGIN_PROVIDER [ integer, psi_config_data, (N_configuration,2)] + &BEGIN_PROVIDER [ integer, psi_csf_to_config_data, (n_CSF)] + use cfunctions + implicit none + BEGIN_DOC + ! Documentation for DetToCSFTransformationMatrix + ! Provides the matrix of transformatons for the + ! conversion between determinant to CSF basis (in BFs) + END_DOC + integer*8 :: Isomo, Idomo + integer(bit_kind) :: Ialpha(N_int),Ibeta(N_int) + integer :: rows, cols, i, j, k + integer :: startdet, enddet + integer*8 MS + integer ndetI + integer :: getNSOMO + real*8,dimension(:,:),allocatable :: tempBuffer + real*8,dimension(:),allocatable :: tempCoeff + real*8 :: norm_det1, phasedet + norm_det1 = 0.d0 + MS = elec_alpha_num - elec_beta_num + print *,"Maxbfdim=",NBFMax + print *,"Maxdetdim=",maxDetDimPerBF + print *,"n_CSF=",n_CSF + print *,"N_configurations=",N_configuration + print *,"n_core_orb=",n_core_orb + ! initialization + psi_coef_config = 0.d0 + DetToCSFTransformationMatrix(0,:,:) = 1.d0 + do i = 2-iand(elec_alpha_num-elec_beta_num,1), NSOMOMax,2 + Isomo = IBSET(0_8, i) - 1_8 + ! rows = Ncsfs + ! cols = Ndets + bfIcfg = max(1,nint((binom(i,(i+1)/2)-binom(i,((i+1)/2)+1)))) + ndetI = max(1,nint((binom(i,(i+1)/2)))) + + allocate(tempBuffer(bfIcfg,ndetI)) + call getCSFtoDETTransformationMatrix(Isomo, MS, NBFMax, maxDetDimPerBF, tempBuffer) + DetToCSFTransformationMatrix(i,:bfIcfg,:ndetI) = tempBuffer + deallocate(tempBuffer) + enddo + + integer s, bfIcfg + integer countcsf + countcsf = 0 + integer countdet + countdet = 0 + integer istate + istate = 1 + phasedet = 1.0d0 + do i = 1,N_configuration + startdet = psi_configuration_to_psi_det(1,i) + enddet = psi_configuration_to_psi_det(2,i) + ndetI = enddet-startdet+1 + + allocate(tempCoeff(ndetI)) + countdet = 1 + do j = startdet, enddet + Ialpha = psi_det(:,1,psi_configuration_to_psi_det_data(j)) + Ibeta = psi_det(:,2,psi_configuration_to_psi_det_data(j)) + !call debug_spindet(Ialpha,1,1) + !call debug_spindet(Ibeta ,1,1) + call get_phase_qp_to_cfg(Ialpha, Ibeta, phasedet) + !print *,">>",Ialpha,Ibeta,phasedet + tempCoeff(countdet) = psi_coef(psi_configuration_to_psi_det_data(j), istate)*phasedet + !tempCoeff(countdet) = psi_coef(psi_configuration_to_psi_det_data(j), istate) + norm_det1 += tempCoeff(countdet)*tempCoeff(countdet) + countdet += 1 + enddo + + !print *,"dimcoef=",bfIcfg,norm_det1 + !call printMatrix(tempCoeff,ndetI,1) + + s = 0 + do k=1,N_int + if (psi_configuration(k,1,i) == 0_bit_kind) cycle + s = s + popcnt(psi_configuration(k,1,i)) + enddo + bfIcfg = max(1,nint((binom(s,(s+1)/2)-binom(s,((s+1)/2)+1)))) + + ! perhaps blocking with CFGs of same seniority + ! can be more efficient + allocate(tempBuffer(bfIcfg,ndetI)) + tempBuffer = DetToCSFTransformationMatrix(s,:bfIcfg,:ndetI) + !print *,"csftodetdim=",bfIcfg,ndetI + !call printMatrix(tempBuffer,bfIcfg,ndetI) + + call dgemm('N','N', bfIcfg, 1, ndetI, 1.d0, tempBuffer, size(tempBuffer,1), tempCoeff, size(tempCoeff,1), 0.d0, psi_coef_config(countcsf+1,1), size(psi_coef_config,1)) + !call dgemv('N', NBFMax, maxDetDimPerBF, 1.d0, tempBuffer, size(tempBuffer,1), tempCoeff, 1, 0.d0, psi_coef_config(countcsf), 1) + + !call printMatrix(psi_coef_config(countcsf+1,1),bfIcfg,1) + deallocate(tempCoeff) + deallocate(tempBuffer) + psi_config_data(i,1) = countcsf + 1 + countcsf += bfIcfg + psi_config_data(i,2) = countcsf + psi_csf_to_config_data(countcsf) = i + enddo + print *,"Norm det=",norm_det1, size(psi_coef_config,1), " Dim csf=", countcsf + + END_PROVIDER + + BEGIN_PROVIDER [ integer, AIJpqMatrixDimsList, (0:NSOMOMax,0:NSOMOMax,4,NSOMOMax+1,NSOMOMax+1,2)] &BEGIN_PROVIDER [ integer, rowsmax] &BEGIN_PROVIDER [ integer, colsmax] use cfunctions @@ -127,9 +239,13 @@ end subroutine get_phase_qp_to_cfg nsomomin = elec_alpha_num-elec_beta_num rowsmax = 0 colsmax = 0 + print *,"NSOMOMax = ",NSOMOMax !allocate(AIJpqMatrixDimsList(NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,2)) ! Type ! 1. SOMO -> SOMO + !print *,"Doing SOMO->SOMO" + AIJpqMatrixDimsList(0,0,1,1,1,1) = 1 + AIJpqMatrixDimsList(0,0,1,1,1,2) = 1 do i = 2-iand(nsomomin,1), NSOMOMax, 2 Isomo = ISHFT(1_8,i)-1 do j = i-2,i-2, 2 @@ -157,6 +273,7 @@ end subroutine get_phase_qp_to_cfg MS, & rows, & cols) + !print *, "SOMO->SOMO \t",i,j,k,l,">",Isomo,Jsomo,">",rows, cols if(rowsmax .LT. rows) then rowsmax = rows end if @@ -172,6 +289,9 @@ end subroutine get_phase_qp_to_cfg end do ! Type ! 2. DOMO -> VMO + !print *,"Doing DOMO->VMO" + AIJpqMatrixDimsList(0,0,2,1,1,1) = 1 + AIJpqMatrixDimsList(0,0,2,1,1,2) = 1 do i = 0+iand(nsomomin,1), NSOMOMax, 2 Isomo = ISHFT(1_8,i)-1 tmpsomo = ISHFT(1_8,i+2)-1 @@ -205,6 +325,7 @@ end subroutine get_phase_qp_to_cfg MS, & rows, & cols) + !print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols if(rowsmax .LT. rows) then rowsmax = rows end if @@ -221,6 +342,8 @@ end subroutine get_phase_qp_to_cfg ! Type ! 3. SOMO -> VMO !print *,"Doing SOMO->VMO" + AIJpqMatrixDimsList(0,0,3,1,1,1) = 1 + AIJpqMatrixDimsList(0,0,3,1,1,2) = 1 do i = 2-iand(nsomomin,1), NSOMOMax, 2 Isomo = ISHFT(1_8,i)-1 do j = i,i, 2 @@ -228,8 +351,8 @@ end subroutine get_phase_qp_to_cfg if(j .GT. NSOMOMax .OR. j .LE. 0) then cycle end if - do k = 1,i - do l = 1,i + do k = 1,i+1 + do l = 1,i+1 if(k .NE. l) then Isomo = ISHFT(1_8,i+1)-1 Isomo = IBCLR(Isomo,l-1) @@ -244,6 +367,7 @@ end subroutine get_phase_qp_to_cfg MS, & rows, & cols) + !print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols if(rowsmax .LT. rows) then rowsmax = rows end if @@ -258,18 +382,20 @@ end subroutine get_phase_qp_to_cfg end do end do ! Type - ! 4. DOMO -> VMO + ! 4. DOMO -> SOMO !print *,"Doing DOMO->SOMO" + AIJpqMatrixDimsList(0,0,4,1,1,1) = 1 + AIJpqMatrixDimsList(0,0,4,1,1,2) = 1 do i = 2-iand(nsomomin,1), NSOMOMax, 2 do j = i,i, 2 if(j .GT. NSOMOMax .OR. j .LE. 0) then cycle end if - do k = 1,i - do l = 1,i + do k = 1,i+1 + do l = 1,i+1 if(k .NE. l) then Isomo = ISHFT(1_8,i+1)-1 - Isomo = IBCLR(Isomo,k+1-1) + Isomo = IBCLR(Isomo,k-1) Jsomo = ISHFT(1_8,j+1)-1 Jsomo = IBCLR(Jsomo,l-1) else @@ -281,6 +407,7 @@ end subroutine get_phase_qp_to_cfg MS, & rows, & cols) + !print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols if(rowsmax .LT. rows) then rowsmax = rows end if @@ -294,9 +421,10 @@ end subroutine get_phase_qp_to_cfg end do end do end do + print *,"Rowsmax=",rowsmax," Colsmax=",colsmax END_PROVIDER - BEGIN_PROVIDER [ real*8, AIJpqContainer, (0:NSOMOMax,0:NSOMOMax,4,NSOMOMax,NSOMOMax,NBFMax,NBFMax)] + BEGIN_PROVIDER [ real*8, AIJpqContainer, (0:NSOMOMax,0:NSOMOMax,4,NSOMOMax+1,NSOMOMax+1,NBFMax,NBFMax)] use cfunctions implicit none BEGIN_DOC @@ -332,13 +460,18 @@ end subroutine get_phase_qp_to_cfg integer maxdim !maxdim = max(rowsmax,colsmax) ! allocate matrix + !print *,"rowsmax =",rowsmax," colsmax=",colsmax + !print *,"NSOMOMax = ",NSOMOMax !allocate(AIJpqMatrixDimsList(NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,2)) ! Type ! 1. SOMO -> SOMO + !print *,"Doing SOMO -> SOMO" + AIJpqContainer(0,0,1,1,1,1,1) = 1.0d0 do i = 2, NSOMOMax, 2 Isomo = ISHFT(1_8,i)-1 do j = i-2,i-2, 2 if(j .GT. NSOMOMax .OR. j .LT. 0) cycle + !print *,"i,j=",i,j do k = 1,i do l = 1,i @@ -355,6 +488,10 @@ end subroutine get_phase_qp_to_cfg nsomoj = i endif + !print *,"k,l=",k,l + !call debug_spindet(Jsomo,1) + !call debug_spindet(Isomo,1) + AIJpqContainer(nsomoi,nsomoj,1,k,l,:,:) = 0.0d0 call getApqIJMatrixDims(Isomo, & Jsomo, & @@ -376,6 +513,8 @@ end subroutine get_phase_qp_to_cfg meMatrix, & rows, & cols) + !print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols,">",rowsmax,colsmax + !call printMatrix(meMatrix,rows,cols) ! i -> j do ri = 1,rows do ci = 1,cols @@ -389,6 +528,8 @@ end subroutine get_phase_qp_to_cfg end do ! Type ! 2. DOMO -> VMO + !print *,"Doing DOMO -> VMO" + AIJpqContainer(0,0,2,1,1,1,1) = 1.0d0 do i = 0, NSOMOMax, 2 Isomo = ISHFT(1_8,i)-1 tmpsomo = ISHFT(1_8,i+2)-1 @@ -414,6 +555,10 @@ end subroutine get_phase_qp_to_cfg nsomoj = j endif + !print *,"k,l=",k,l + !call debug_spindet(Jsomo,1) + !call debug_spindet(Isomo,1) + AIJpqContainer(nsomoi,nsomoj,2,k,l,:,:) = 0.0d0 call getApqIJMatrixDims(Isomo, & Jsomo, & @@ -435,6 +580,8 @@ end subroutine get_phase_qp_to_cfg meMatrix, & rows, & cols) + !print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols,">",rowsmax,colsmax + !call printMatrix(meMatrix,rows,cols) ! i -> j do ri = 1,rows do ci = 1,cols @@ -448,13 +595,14 @@ end subroutine get_phase_qp_to_cfg end do ! Type ! 3. SOMO -> VMO + !print *,"Doing SOMO -> VMO" do i = 2, NSOMOMax, 2 Isomo = ISHFT(1_8,i)-1 do j = i,i, 2 Jsomo = ISHFT(1_8,j)-1 if(j .GT. NSOMOMax .OR. j .LE. 0) cycle - do k = 1,i - do l = 1,i + do k = 1,i+1 + do l = 1,i+1 if(k .NE. l) then Isomo = ISHFT(1_8,i+1)-1 Isomo = IBCLR(Isomo,l-1) @@ -465,6 +613,10 @@ end subroutine get_phase_qp_to_cfg Jsomo = ISHFT(1_8,j)-1 endif + !print *,"k,l=",k,l + !call debug_spindet(Jsomo,1) + !call debug_spindet(Isomo,1) + AIJpqContainer(i,j,3,k,l,:,:) = 0.0d0 call getApqIJMatrixDims(Isomo, & Jsomo, & @@ -486,6 +638,8 @@ end subroutine get_phase_qp_to_cfg meMatrix, & rows, & cols) + !call printMatrix(meMatrix,rows,cols) + !print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols,">",rowsmax,colsmax ! i -> j do ri = 1,rows do ci = 1,cols @@ -499,18 +653,20 @@ end subroutine get_phase_qp_to_cfg end do ! Type ! 4. DOMO -> SOMO + !print *,"Doing DOMO -> SOMO" + AIJpqContainer(0,0,4,1,1,1,1) = 1.0d0 do i = 2, NSOMOMax, 2 Isomo = ISHFT(1_8,i)-1 do j = i,i, 2 Jsomo = ISHFT(1_8,i)-1 if(j .GT. NSOMOMax .OR. j .LE. 0) cycle - do k = 1,i - do l = 1,i + do k = 1,i+1 + do l = 1,i+1 if(k .NE. l) then Isomo = ISHFT(1_8,i+1)-1 Isomo = IBCLR(Isomo,k-1) Jsomo = ISHFT(1_8,j+1)-1 - Jsomo = IBCLR(Jsomo,l+1-1) + Jsomo = IBCLR(Jsomo,l-1) else Isomo = ISHFT(1_8,i)-1 Jsomo = ISHFT(1_8,j)-1 @@ -538,6 +694,8 @@ end subroutine get_phase_qp_to_cfg meMatrix, & rows, & cols) + !call printMatrix(meMatrix,rows,cols) + !print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols,">",rowsmax,colsmax ! i -> j do ri = 1,rows do ci = 1,cols @@ -551,93 +709,466 @@ end subroutine get_phase_qp_to_cfg end do END_PROVIDER - -!!!!!! - - BEGIN_PROVIDER [ real*8, DetToCSFTransformationMatrix, (0:NSOMOMax,NBFMax,maxDetDimPerBF)] - &BEGIN_PROVIDER [ real*8, psi_coef_config, (n_CSF)] - &BEGIN_PROVIDER [ integer, psi_config_data, (N_configuration,2)] - use cfunctions - use bitmasks +subroutine calculate_preconditioner_cfg(diag_energies) implicit none + use bitmasks BEGIN_DOC - ! Documentation for DetToCSFTransformationMatrix - ! Provides the matrix of transformatons for the - ! conversion between determinant to CSF basis (in BFs) + ! Documentation for calculate_preconditioner + ! + ! Calculates the diagonal energies of + ! the configurations in psi_configuration + ! returns : diag_energies : END_DOC - integer(bit_kind) :: mask(N_int), Ialpha(N_int),Ibeta(N_int) - integer :: rows, cols, i, j, k - integer :: startdet, enddet - integer*8 MS, Isomo, Idomo - integer ndetI + integer :: i,j,k,l,p,q,noccp,noccq, ii, jj + real*8,intent(out) :: diag_energies(n_CSF) + integer :: nholes + integer :: nvmos + integer :: listvmos(mo_num) + integer :: vmotype(mo_num) ! 1 -> VMO 2 -> SOMO + integer :: listholes(mo_num) + integer :: holetype(mo_num) ! 1-> SOMO 2->DOMO + integer*8 :: Idomo + integer*8 :: Isomo + integer*8 :: Jdomo + integer*8 :: Jsomo + integer*8 :: diffSOMO + integer*8 :: diffDOMO + integer :: NSOMOI + integer :: NSOMOJ + integer :: ndiffSOMO + integer :: ndiffDOMO + integer :: starti, endi, cnti, cntj, rows,cols + integer :: extype,pmodel,qmodel + integer(bit_kind) :: Icfg(N_INT,2) + integer(bit_kind) :: Jcfg(N_INT,2) + integer,external :: getNSOMO + real*8, external :: mo_two_e_integral + real*8 :: hpp + real*8 :: meCC + real*8 :: ecore + + ! initialize energies + diag_energies = 0.d0 + + ! calculate core energy + !call get_core_energy(ecore) + !diag_energies = ecore + + ! calculate the core energy + !print *,"Core energy=",ref_bitmask_energy + + 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)) + + starti = psi_config_data(i,1) + endi = psi_config_data(i,2) + + ! find out all pq holes possible + nholes = 0 + ! holes in SOMO + !do k = n_core_orb+1,n_core_orb + n_act_orb + do k = 1,mo_num + 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 + if(POPCNT(IAND(Idomo,IBSET(0_8,k-1))) .EQ. 1) then + nholes += 1 + listholes(nholes) = k + holetype(nholes) = 2 + endif + enddo + + ! 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 + !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) + + do k=1,nholes + p = listholes(k) + noccp = holetype(k) + + ! Calculate one-electron + ! and two-electron coulomb terms + do l=1,nholes + q = listholes(l) + noccq = holetype(l) + !print *,"--------------- K=",p," L=",q + + ! one-electron term + if(p.EQ.q) then + hpp = noccq * h_core_ri(p,q)!mo_one_e_integrals(q,q) + else + hpp = 0.d0 + endif + + + do j=starti,endi + ! coulomb term + ! (pp,qq) = + if(p.EQ.q) then + diag_energies(j) += hpp !+ 0.5d0 * (noccp * noccq * mo_two_e_integral(p,q,p,q)) + !print *,"hpp=",hpp,"diga= ",diag_energies(j) +! else +! diag_energies(j) += ! 0.5d0 * noccp * noccq * mo_two_e_integral(p,q,p,q) +! print *,"diga= ",diag_energies(j) + endif + enddo + enddo + + enddo + enddo + +end subroutine calculate_preconditioner_cfg + + +subroutine calculate_sigma_vector_cfg_nst(psi_out, psi_in, n_st, sze, istart, iend, ishift, istep) + implicit none + use bitmasks + BEGIN_DOC + ! Documentation for sigma-vector calculation + ! + ! Calculates the result of the + ! application of the hamiltonian to the + ! wavefunction in CFG basis once + ! TODO : Things prepare outside this routine + ! 1. Touch the providers for + ! a. ApqIJ containers + ! b. DET to CSF transformation matrices + ! 2. DET to CSF transcormation + ! 2. CSF to DET back transcormation + ! returns : psi_coef_out_det : + END_DOC + integer,intent(in) :: sze, istart,iend, istep, ishift, n_st + real*8,intent(in):: psi_in(sze,n_st) + real*8,intent(out):: psi_out(sze,n_st) + integer(bit_kind) :: Icfg(N_INT,2) + integer :: i,j,k,l,p,q,noccp,noccq, ii, jj, m, n, idxI, kk, nocck,orbk + integer(bit_kind) :: alphas_Icfg(N_INT,2,sze) + integer(bit_kind) :: singlesI(N_INT,2,sze) + integer(bit_kind) :: connectedI_alpha(N_INT,2,sze) + integer :: idxs_singlesI(sze) + integer :: idxs_connectedI_alpha(sze) + integer(bit_kind) :: psi_configuration_out(N_INT,2,sze) + real*8 :: psi_coef_out(n_CSF) + logical :: psi_coef_out_init(n_CSF) + integer :: excitationIds_single(2,sze) + integer :: excitationTypes_single(sze) + integer :: excitationIds(2,sze) + integer :: excitationTypes(sze) + real*8 :: diagfactors(sze) + integer :: nholes + integer :: nvmos + integer :: listvmos(mo_num) + integer :: vmotype(mo_num) ! 1 -> VMO 2 -> SOMO + integer :: listholes(mo_num) + integer :: holetype(mo_num) ! 1-> SOMO 2->DOMO + integer :: Nalphas_Icfg, nconnectedI, rowsikpq, colsikpq, nsinglesI + integer :: extype,NSOMOalpha,NSOMOI,NSOMOJ,pmodel,qmodel integer :: getNSOMO - real*8,dimension(:,:),allocatable :: tempBuffer - real*8,dimension(:),allocatable :: tempCoeff - real*8 :: norm_det1, phasedet - norm_det1 = 0.d0 - MS = elec_alpha_num - elec_beta_num - ! initialization - psi_coef_config = 0.d0 - DetToCSFTransformationMatrix(0,:,:) = 1.d0 - do i = 2-iand(elec_alpha_num-elec_beta_num,1), NSOMOMax,2 - Isomo = IBSET(0_8, i) - 1_8 - ! rows = Ncsfs - ! cols = Ndets - bfIcfg = max(1,nint((binom(i,(i+1)/2)-binom(i,((i+1)/2)+1)))) - ndetI = max(1,nint((binom(i,(i+1)/2)))) + integer :: totcolsTKI + integer :: rowsTKI + integer :: noccpp + integer :: istart_cfg, iend_cfg + integer*8 :: MS, Isomo, Idomo, Jsomo, Jdomo, Ialpha, Ibeta + integer :: moi, moj, mok, mol, starti, endi, startj, endj, cnti, cntj, cntk + real*8 :: norm_coef_cfg, fac2eints + real*8 :: norm_coef_det + real*8 :: meCC1, meCC2, diagfac + real*8,dimension(:,:,:),allocatable :: TKI + real*8,dimension(:,:),allocatable :: GIJpqrs + real*8,dimension(:,:,:),allocatable :: TKIGIJ + real*8, external :: mo_two_e_integral + real*8, external :: get_two_e_integral + real*8 :: diag_energies(n_CSF) + call calculate_preconditioner_cfg(diag_energies) - allocate(tempBuffer(bfIcfg,ndetI)) - call getCSFtoDETTransformationMatrix(Isomo, MS, NBFMax, maxDetDimPerBF, tempBuffer) - DetToCSFTransformationMatrix(i,1:bfIcfg,1:ndetI) = tempBuffer(1:bfIcfg,1:ndetI) - deallocate(tempBuffer) + MS = 0 + norm_coef_cfg=0.d0 + + psi_out=0.d0 + psi_coef_out_init = .False. + + istart_cfg = psi_csf_to_config_data(istart) + iend_cfg = psi_csf_to_config_data(iend) + + + !!! Single Excitations !!! + do i=istart_cfg,iend_cfg + + 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) + + ! find out all pq holes possible + nholes = 0 + ! holes in SOMO + ! list_act + ! list_core + ! list_core_inact + ! bitmasks + !do k = n_core_orb+1,n_core_orb + n_act_orb + do k = 1,mo_num + 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,mo_num + if(POPCNT(IAND(Idomo,IBSET(0_8,k-1))) .EQ. 1) then + nholes += 1 + listholes(nholes) = k + holetype(nholes) = 2 + endif + enddo + + ! 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 + !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 + + + ! Icsf ids + starti = psi_config_data(i,1) + endi = psi_config_data(i,2) + NSOMOI = getNSOMO(Icfg) + + call generate_all_singles_cfg_with_type(Icfg,singlesI,idxs_singlesI,excitationIds_single, & + excitationTypes_single,nsinglesI,N_int) + + do j = 1,nsinglesI + idxI = idxs_singlesI(j) + NSOMOJ = getNSOMO(singlesI(:,:,j)) + p = excitationIds_single(1,j) + q = excitationIds_single(2,j) + extype = excitationTypes_single(j) + ! Off diagonal terms + call convertOrbIdsToModelSpaceIds(Icfg, singlesI(:,:,j), p, q, extype, pmodel, qmodel) + Jsomo = singlesI(1,1,j) + Jdomo = singlesI(1,2,j) + + ! 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 + 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 + nholes += 1 + listholes(nholes) = q + holetype(nholes) = 2 + endif + + startj = psi_config_data(idxI,1) + endj = psi_config_data(idxI,2) + + !!! One-electron contribution !!! + do kk = 1,n_st + cnti = 0 + do ii = starti, endi + cnti += 1 + cntj = 0 + do jj = startj, endj + cntj += 1 + meCC1 = AIJpqContainer(NSOMOI,NSOMOJ,extype,pmodel,qmodel,cnti,cntj) + psi_out(jj,kk) += meCC1 * psi_in(ii,kk) * h_core_ri(p,q) + psi_coef_out_init(jj) = .True. + enddo + enddo + 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 + 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 + nholes -= 1 + endif + enddo enddo - integer s, bfIcfg - integer countcsf - countcsf = 0 - integer countdet - countdet = 0 - integer idx - integer istate - istate = 1 - phasedet = 1.0d0 - do i = 1,N_configuration - startdet = psi_configuration_to_psi_det(1,i) - enddet = psi_configuration_to_psi_det(2,i) - ndetI = enddet-startdet+1 + !!! Double Excitations !!! - allocate(tempCoeff(ndetI)) - countdet = 1 - do j = startdet, enddet - idx = psi_configuration_to_psi_det_data(j) - Ialpha(:) = psi_det(:,1,idx) - Ibeta(:) = psi_det(:,2,idx) - call get_phase_qp_to_cfg(Ialpha, Ibeta, phasedet) - tempCoeff(countdet) = psi_coef(idx, istate)*phasedet - norm_det1 += tempCoeff(countdet)*tempCoeff(countdet) - countdet += 1 - enddo + ! Loop over all selected configurations + do i = istart_cfg,iend_cfg - s = 0 - do k=1,N_int - if (psi_configuration(k,1,i) == 0_bit_kind) cycle - s = s + popcnt(psi_configuration(k,1,i)) - enddo - bfIcfg = max(1,nint((binom(s,(s+1)/2)-binom(s,((s+1)/2)+1)))) + Icfg(1,1) = psi_configuration(1,1,i) + Icfg(1,2) = psi_configuration(1,2,i) + starti = psi_config_data(i,1) + endi = psi_config_data(i,2) - ! perhaps blocking with CFGs of same seniority - ! can be more efficient - allocate(tempBuffer(bfIcfg,ndetI)) - tempBuffer = DetToCSFTransformationMatrix(s,:bfIcfg,:ndetI) + ! Returns all unique (checking the past) singly excited cfgs connected to I + call obtain_associated_alphaI(i, Icfg, alphas_Icfg, Nalphas_Icfg) + ! TODO : remove doubly excited for return + ! Here we do 2x the loop. One to count for the size of the matrix, then we compute. + do k = 1,Nalphas_Icfg + ! Now generate all singly excited with respect to a given alpha CFG + call obtain_connected_I_foralpha(i,alphas_Icfg(1,1,k),connectedI_alpha,idxs_connectedI_alpha,nconnectedI,excitationIds,excitationTypes,diagfactors) - call dgemm('N','N', bfIcfg, 1, ndetI, 1.d0, tempBuffer, size(tempBuffer,1), tempCoeff, size(tempCoeff,1), 0.d0, psi_coef_config(countcsf+1), size(psi_coef_config,1)) - !call dgemv('N', NBFMax, maxDetDimPerBF, 1.d0, tempBuffer, size(tempBuffer,1), tempCoeff, 1, 0.d0, psi_coef_config(countcsf), 1) + if(nconnectedI .EQ. 0) then + cycle + endif + totcolsTKI = 0 + rowsTKI = -1 + do j = 1,nconnectedI + NSOMOalpha = getNSOMO(alphas_Icfg(:,:,k)) + NSOMOI = getNSOMO(connectedI_alpha(:,:,j)) + p = excitationIds(1,j) + q = excitationIds(2,j) + extype = excitationTypes(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 + if(p.EQ.q) then + NSOMOalpha = NSOMOI + endif + rowsikpq = AIJpqMatrixDimsList(NSOMOalpha,NSOMOI,extype,pmodel,qmodel,1) + colsikpq = AIJpqMatrixDimsList(NSOMOalpha,NSOMOI,extype,pmodel,qmodel,2) + totcolsTKI += colsikpq + if(rowsTKI .LT. rowsikpq .AND. rowsTKI .NE. -1) then + print *,">",j,"Something is wrong in sigma-vector", rowsTKI, rowsikpq, "(p,q)=",pmodel,qmodel,"ex=",extype,"na=",NSOMOalpha," nI=",NSOMOI + !rowsTKI = rowsikpq + else + rowsTKI = rowsikpq + endif + enddo - deallocate(tempCoeff) - deallocate(tempBuffer) - psi_config_data(i,1) = countcsf + 1 - countcsf += bfIcfg - psi_config_data(i,2) = countcsf + allocate(TKI(rowsTKI,n_st,totcolsTKI)) ! coefficients of CSF + ! Initialize the inegral container + ! dims : (totcolsTKI, nconnectedI) + allocate(GIJpqrs(totcolsTKI,nconnectedI)) ! gpqrs + allocate(TKIGIJ(rowsTKI,n_st,nconnectedI)) ! gpqrs + + totcolsTKI = 0 + do j = 1,nconnectedI + NSOMOalpha = getNSOMO(alphas_Icfg(:,:,k)) + NSOMOI = getNSOMO(connectedI_alpha(:,:,j)) + p = excitationIds(1,j) + q = excitationIds(2,j) + extype = excitationTypes(j) + call convertOrbIdsToModelSpaceIds(alphas_Icfg(:,:,k), connectedI_alpha(:,:,j), p, q, extype, pmodel, qmodel) + rowsikpq = AIJpqMatrixDimsList(NSOMOalpha,NSOMOI,extype,pmodel,qmodel,1) + colsikpq = AIJpqMatrixDimsList(NSOMOalpha,NSOMOI,extype,pmodel,qmodel,2) + do kk = 1,n_st + do l = 1,rowsTKI + do m = 1,colsikpq + TKI(l,kk,totcolsTKI+m) = AIJpqContainer(NSOMOalpha,NSOMOI,extype,pmodel,qmodel,l,m) * psi_in(idxs_connectedI_alpha(j)+m-1,kk) + enddo + enddo + enddo + do m = 1,colsikpq + do l = 1,nconnectedI + ! = (ik|jl) + moi = excitationIds(1,j) ! p + mok = excitationIds(2,j) ! q + moj = excitationIds(2,l) ! s + mol = excitationIds(1,l) ! r + if(moi.EQ.mok .AND. moj.EQ.mol)then + diagfac = diagfactors(j) + diagfac *= diagfactors(l) + !print *,"integrals (",totcolsTKI+m,l,")",mok,moi,mol,moj, "|", diagfac + GIJpqrs(totcolsTKI+m,l) = diagfac*0.5d0*mo_two_e_integral(mok,mol,moi,moj) ! g(pq,sr) = + else + diagfac = diagfactors(j)*diagfactors(l) + !print *,"integrals (",totcolsTKI+m,l,")",mok,moi,mol,moj, "|", diagfac + GIJpqrs(totcolsTKI+m,l) = diagfac*0.5d0*mo_two_e_integral(mok,mol,moi,moj) ! g(pq,sr) = + !endif + endif + enddo + enddo + totcolsTKI += colsikpq + enddo + + + + ! Do big BLAS + ! TODO TKI, size(TKI,1)*size(TKI,2) + call dgemm('N','N', rowsTKI*n_st, nconnectedI, totcolsTKI, 1.d0, & + TKI, size(TKI,1)*n_st, GIJpqrs, size(GIJpqrs,1), 0.d0, & + TKIGIJ , size(TKIGIJ,1)*n_st ) + + + ! Collect the result + totcolsTKI = 0 + do j = 1,nconnectedI + NSOMOalpha = getNSOMO(alphas_Icfg(:,:,k)) + NSOMOI = getNSOMO(connectedI_alpha(:,:,j)) + p = excitationIds(1,j) + q = excitationIds(2,j) + extype = excitationTypes(j) + call convertOrbIdsToModelSpaceIds(alphas_Icfg(:,:,k), connectedI_alpha(:,:,j), p, q, extype, pmodel, qmodel) + rowsikpq = AIJpqMatrixDimsList(NSOMOalpha,NSOMOI,extype,pmodel,qmodel,1) + colsikpq = AIJpqMatrixDimsList(NSOMOalpha,NSOMOI,extype,pmodel,qmodel,2) + !print *,">j=",j,rowsikpq,colsikpq, ">>",totcolsTKI,",",idxs_connectedI_alpha(j) + do kk = 1,n_st + do m = 1,colsikpq + do l = 1,rowsTKI + psi_out(idxs_connectedI_alpha(j)+m-1,kk) += AIJpqContainer(NSOMOalpha,NSOMOI,extype,pmodel,qmodel,l,m) * TKIGIJ(l,kk,j) + psi_coef_out_init(idxs_connectedI_alpha(j)+m-1) = .True. + enddo + enddo + enddo + totcolsTKI += colsikpq + enddo + + deallocate(TKI) ! coefficients of CSF + ! Initialize the inegral container + ! dims : (totcolsTKI, nconnectedI) + deallocate(GIJpqrs) ! gpqrs + deallocate(TKIGIJ) ! gpqrs + + enddo ! loop over alphas + enddo ! loop over I + + + ! Add the diagonal contribution + do i = 1,n_CSF + psi_out(i,1) += 1.0d0*diag_energies(i)*psi_in(i,1) enddo - END_PROVIDER + +end subroutine calculate_sigma_vector_cfg_nst diff --git a/src/csf/tree_utils.h b/src/csf/tree_utils.h index a94bf005..332ddf69 100644 --- a/src/csf/tree_utils.h +++ b/src/csf/tree_utils.h @@ -23,6 +23,7 @@ struct bin_tree { }; #include "/opt/intel/oneapi/mkl/2021.1.1/include/mkl_cblas.h" +//#include "cblas.h" #define MAX_SOMO 32 diff --git a/src/davidson/diagonalization_hcsf_dressed.irp.f b/src/davidson/diagonalization_hcsf_dressed.irp.f index 08f6a9a0..00f901bb 100644 --- a/src/davidson/diagonalization_hcsf_dressed.irp.f +++ b/src/davidson/diagonalization_hcsf_dressed.irp.f @@ -88,7 +88,7 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N double precision, intent(out) :: energies(N_st_diag_in) integer :: iter, N_st_diag - integer :: i,j,k,l,m + integer :: i,j,k,l,m,kk logical, intent(inout) :: converged double precision, external :: u_dot_v, u_dot_u @@ -285,7 +285,7 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N ! Make random verctors eigenstates of S2 call convertWFfromDETtoCSF(N_st_diag,U,U_csf) - call convertWFfromCSFtoDET(N_st_diag,U_csf,U) + !call convertWFfromCSFtoDET(N_st_diag,U_csf,U) do while (.not.converged) itertot = itertot+1 @@ -302,11 +302,28 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N ! Compute |W_k> = \sum_i |i> ! ----------------------------------- - call convertWFfromCSFtoDET(N_st_diag,U_csf(1,shift+1),U) + !call convertWFfromCSFtoDET(N_st_diag,U_csf(1,shift+1),U) if ((sze > 100000).and.distributed_davidson) then - call H_u_0_nstates_zmq (W,U,N_st_diag,sze) + + !call convertWFfromCSFtoDET(N_st_diag,U_csf(1,shift+1),U) + !call convertWFfromCSFtoDET(N_st_diag,W_csf(1,shift+1),W) + !call H_u_0_nstates_zmq (W,U,N_st_diag,sze) + !call convertWFfromDETtoCSF(N_st_diag,U,U_csf(1,shift+1)) + !call convertWFfromDETtoCSF(N_st_diag,W,W_csf(1,shift+1)) + !call calculate_sigma_vector_cfg_nst(W_csf(1,shift+1),U_csf(1,shift+1),N_st_diag,sze_csf,1,sze_csf,0,1) + do kk=1,N_st_diag + call calculate_sigma_vector_cfg_nst(W_csf(1,shift+kk),U_csf(1,shift+kk),1,sze_csf,1,sze_csf,0,1) + enddo else - call H_u_0_nstates_openmp(W,U,N_st_diag,sze) + !call convertWFfromCSFtoDET(N_st_diag,U_csf(1,shift+1),U) + !call convertWFfromCSFtoDET(N_st_diag,W_csf(1,shift+1),W) + !call H_u_0_nstates_openmp(W,U,N_st_diag,sze) + !call convertWFfromDETtoCSF(N_st_diag,U,U_csf(1,shift+1)) + !call convertWFfromDETtoCSF(N_st_diag,W,W_csf(1,shift+1)) + !call calculate_sigma_vector_cfg_nst(W_csf(1,shift+1),U_csf(1,shift+1),N_st_diag,sze_csf,1,sze_csf,0,1) + do kk=1,N_st_diag + call calculate_sigma_vector_cfg_nst(W_csf(1,shift+kk),U_csf(1,shift+kk),1,sze_csf,1,sze_csf,0,1) + enddo endif else ! Already computed in update below @@ -350,7 +367,7 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N endif endif - call convertWFfromDETtoCSF(N_st_diag,W,W_csf(1,shift+1)) + !call convertWFfromDETtoCSF(N_st_diag,W,W_csf(1,shift+1)) ! Compute h_kl = = ! -------------------------------------------