From f2bbbfc2998f050d08c27a73b1b6cb909c87e573 Mon Sep 17 00:00:00 2001 From: v1j4y Date: Mon, 1 Mar 2021 12:28:45 +0100 Subject: [PATCH 1/7] Added my configuration_CI_sigma_helpers.org --- src/csf/configuration_CI_sigma_helpers.org | 103 ++- .../configuration_CI_sigma_helpers.org | 657 ++++++++++++++++++ 2 files changed, 742 insertions(+), 18 deletions(-) create mode 100644 src/determinants/configuration_CI_sigma_helpers.org diff --git a/src/csf/configuration_CI_sigma_helpers.org b/src/csf/configuration_CI_sigma_helpers.org index ede42b97..6677139a 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,9 +161,12 @@ 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. !print *,i,k,ndiffSOMO,ndiffDOMO !call debug_spindet(Jsomo,1) @@ -171,11 +177,13 @@ 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 + 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) @@ -236,8 +244,6 @@ the input determinant \(|D_I\rangle\). 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 @@ -287,6 +293,7 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI integer*8 :: diffDOMO integer :: ndiffSOMO integer :: ndiffDOMO + integer :: nxordiffSOMODOMO integer :: i,j,k,l,p,q,nsomoJ,nsomoalpha,starti,endi,extyp,nholes integer :: listholes(mo_num) integer :: holetype(mo_num) @@ -316,7 +323,7 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI 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) @@ -330,8 +337,7 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI diffDOMO = IEOR(Idomo,Jdomo) 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 *,POPCNT(IEOR(diffSOMO,diffDOMO)), ndiffDOMO if(POPCNT(IEOR(diffSOMO,diffDOMO)) .LE. 1 .AND. ndiffDOMO .LT. 3) then !call debug_spindet(Isomo,1) @@ -390,7 +396,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 +525,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,12 +559,12 @@ 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) case (4) ! DOMO -> SOMO ! remove all domos except one at p @@ -518,6 +576,15 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod mask = ISHFT(1_8,q) - 1 Isomo = IAND(Isomo,mask) qmodel = POPCNT(mask) - POPCNT(XOR(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 diff --git a/src/determinants/configuration_CI_sigma_helpers.org b/src/determinants/configuration_CI_sigma_helpers.org new file mode 100644 index 00000000..88b27a7a --- /dev/null +++ b/src/determinants/configuration_CI_sigma_helpers.org @@ -0,0 +1,657 @@ +#+title: Configuration Sigma Vector Helpers +#+author: Vijay Gopal Chilkuri +#+email: vijay.gopal.c@gmail.com + +* Generate the singly excited configurations on-the-fly + +The algorithm is based on the work by Garniron et. al. (see thesis Chap 5). +The basic idea is to generate \(|\alpha\rangle\)'s on-the-fly. + +The algorithm is based on the idea of splitting the list of \(|\alpha\rangle\)'s +into blocks associated with a selected determinant \(|D_I\rangle\). + +** Create a function to return a list of alphas + +Here we create a list of \(|\alpha\rangle\)'s associated with +the input determinant \(|D_I\rangle\). + +#+begin_src f90 :main no :tangle configuration_CI_sigma_helpers.irp.f + subroutine obtain_associated_alphaI(idxI, Icfg, alphasIcfg, NalphaIcfg, factor_alphaI) + implicit none + use bitmasks + BEGIN_DOC + ! Documentation for alphasI + ! Returns the associated alpha's for + ! the input configuration Icfg. + END_DOC + + integer,intent(in) :: idxI ! The id of the Ith CFG + integer(bit_kind),intent(in) :: Icfg(N_int,2) + integer,intent(out) :: NalphaIcfg + real*8 ,intent(out) :: factor_alphaI(*) + integer(bit_kind),intent(out) :: alphasIcfg(N_int,2,*) + logical,dimension(:,:),allocatable :: tableUniqueAlphas + integer :: listholes(mo_num) + integer :: holetype(mo_num) ! 1-> SOMO 2->DOMO + integer :: nholes + integer :: nvmos + integer :: listvmos(mo_num) + integer :: vmotype(mo_num) ! 1 -> VMO 2 -> SOMO + integer*8 :: Idomo + integer*8 :: Isomo + integer*8 :: Jdomo + integer*8 :: Jsomo + integer*8 :: diffSOMO + integer*8 :: diffDOMO + integer*8 :: xordiffSOMODOMO + integer :: ndiffSOMO + integer :: ndiffDOMO + integer :: nxordiffSOMODOMO + integer :: ndiffAll + integer :: i + integer :: j + integer :: k + integer :: hole + integer :: p + integer :: q + integer :: countalphas + logical :: pqAlreadyGenQ + logical :: pqExistsQ + logical :: ppExistsQ + Isomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,1)) + Idomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,2)) + !print*,"Input cfg" + !call debug_spindet(Isomo,1) + !call debug_spindet(Idomo,1) + + !print*,n_act_orb, "monum=",mo_num," n_core=",n_core_orb + + ! find out all pq holes possible + nholes = 0 + ! holes in SOMO + do i = 1,mo_num + if(POPCNT(IAND(Isomo,IBSET(0_8,i-1))) .EQ. 1) then + nholes += 1 + listholes(nholes) = i + holetype(nholes) = 1 + endif + end do + ! holes in DOMO + do i = 1,mo_num + if(POPCNT(IAND(Idomo,IBSET(0_8,i-1))) .EQ. 1) then + nholes += 1 + listholes(nholes) = i + holetype(nholes) = 2 + endif + end do + + ! find vmos + listvmos = -1 + vmotype = -1 + nvmos = 0 + do i = 1,mo_num + !print *,i,IBSET(0,i-1),POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))), POPCNT(IAND(Idomo,(IBSET(0_8,i-1)))) + if(POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))) .EQ. 0 .AND. POPCNT(IAND(Idomo,(IBSET(0_8,i-1)))) .EQ. 0) then + nvmos += 1 + listvmos(nvmos) = i + vmotype(nvmos) = 1 + else if(POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))) .EQ. 1 .AND. POPCNT(IAND(Idomo,(IBSET(0_8,i-1)))) .EQ. 0 ) then + nvmos += 1 + listvmos(nvmos) = i + vmotype(nvmos) = 2 + end if + end do + + !print *,"Nvmo=",nvmos + !print *,listvmos + !print *,vmotype + + allocate(tableUniqueAlphas(mo_num,mo_num)) + tableUniqueAlphas = .FALSE. + + ! Now find the allowed (p,q) excitations + Isomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,1)) + Idomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,2)) + !print *,"Isomo" + !call debug_spindet(Isomo,1) + !call debug_spindet(Idomo,1) + + !print *,"Nholes=",nholes," Nvmos=",nvmos, " idxi=",idxI + !do i = 1,nholes + ! print *,i,"->",listholes(i) + !enddo + !do i = 1,nvmos + ! print *,i,"->",listvmos(i) + !enddo + + ! TODO cfg_seniority_index + do i = 1,nholes + p = listholes(i) + do j = 1,nvmos + q = listvmos(j) + if(p .EQ. q) cycle + if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 1) then + ! SOMO -> VMO + Jsomo = IBCLR(Isomo,p-1) + Jsomo = IBSET(Jsomo,q-1) + Jdomo = Idomo + else if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 2) then + ! SOMO -> SOMO + Jsomo = IBCLR(Isomo,p-1) + Jsomo = IBCLR(Jsomo,q-1) + Jdomo = IBSET(Idomo,q-1) + else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then + ! DOMO -> VMO + Jsomo = IBSET(Isomo,p-1) + Jsomo = IBSET(Jsomo,q-1) + Jdomo = IBCLR(Idomo,p-1) + else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 2) then + ! DOMO -> SOMO + Jsomo = IBSET(Isomo,p-1) + Jsomo = IBCLR(Jsomo,q-1) + Jdomo = IBCLR(Idomo,p-1) + Jdomo = IBSET(Jdomo,q-1) + else + print*,"Something went wrong in obtain_associated_alphaI" + endif + + + pqAlreadyGenQ = .FALSE. + ! First check if it can be generated before + do k = 1, idxI-1 + diffSOMO = IEOR(Jsomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k))) + diffDOMO = IEOR(Jdomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k))) + xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) + ndiffSOMO = POPCNT(diffSOMO) + ndiffDOMO = POPCNT(diffDOMO) + nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO) + !if(POPCNT(IEOR(diffSOMO,diffDOMO)) .LE. 1 .AND. ndiffDOMO .LT. 3) then + if((ndiffSOMO+ndiffDOMO) .EQ. 0) then + pqAlreadyGenQ = .TRUE. + ppExistsQ = .TRUE. + EXIT + endif + if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then + pqAlreadyGenQ = .TRUE. + !ppExistsQ = .TRUE. + !print *,i,k,ndiffSOMO,ndiffDOMO + !call debug_spindet(Jsomo,1) + !call debug_spindet(Jdomo,1) + !call debug_spindet(iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k)),1) + !call debug_spindet(iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k)),1) + EXIT + endif + end do + + !print *,"(,",p,",",q,")",pqAlreadyGenQ + + if(pqAlreadyGenQ) cycle + + pqExistsQ = .FALSE. + ! now check if this exists in the selected list + !do k = idxI+1, N_configuration + ! diffSOMO = IEOR(OR(reunion_of_act_virt_bitmask(1,1),Jsomo),psi_configuration(1,1,k)) + ! diffDOMO = IEOR(OR(reunion_of_act_virt_bitmask(1,1),Jdomo),psi_configuration(1,2,k)) + ! ndiffSOMO = POPCNT(diffSOMO) + ! ndiffDOMO = POPCNT(diffDOMO) + ! if((ndiffSOMO + ndiffDOMO) .EQ. 0) then + ! pqExistsQ = .TRUE. + ! EXIT + ! endif + !end do + + if(.NOT. pqExistsQ) then + tableUniqueAlphas(p,q) = .TRUE. + !print *,p,q + !call debug_spindet(Jsomo,1) + !call debug_spindet(Jdomo,1) + endif + end do + end do + + !print *,tableUniqueAlphas(:,:) + + ! prune list of alphas + Isomo = Icfg(1,1) + Idomo = Icfg(1,2) + Jsomo = Icfg(1,1) + Jdomo = Icfg(1,2) + NalphaIcfg = 0 + do i = 1, nholes + p = listholes(i) + do j = 1, nvmos + q = listvmos(j) + if(p .EQ. q) cycle + if(tableUniqueAlphas(p,q)) then + if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 1) then + ! SOMO -> VMO + Jsomo = IBCLR(Isomo,p-1) + Jsomo = IBSET(Jsomo,q-1) + Jdomo = Idomo + else if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 2) then + ! SOMO -> SOMO + Jsomo = IBCLR(Isomo,p-1) + Jsomo = IBCLR(Jsomo,q-1) + Jdomo = IBSET(Idomo,q-1) + else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then + ! DOMO -> VMO + Jsomo = IBSET(Isomo,p-1) + Jsomo = IBSET(Jsomo,q-1) + Jdomo = IBCLR(Idomo,p-1) + else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 2) then + ! DOMO -> SOMO + Jsomo = IBSET(Isomo,p-1) + Jsomo = IBCLR(Jsomo,q-1) + Jdomo = IBCLR(Idomo,p-1) + Jdomo = IBSET(Jdomo,q-1) + else + print*,"Something went wrong in obtain_associated_alphaI" + endif + + ! SOMO + NalphaIcfg += 1 + !print *,p,q,"|",holetype(i),vmotype(j),NalphaIcfg + alphasIcfg(1,1,NalphaIcfg) = Jsomo + alphasIcfg(1,2,NalphaIcfg) = IOR(Jdomo,ISHFT(1_8,n_core_orb)-1) + endif + end do + end do + + ! Check if this Icfg has been previously generated as a mono + ppExistsQ = .False. + Isomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,1)) + Idomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,2)) + do k = 1, idxI-1 + diffSOMO = IEOR(Isomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k))) + diffDOMO = IEOR(Idomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k))) + xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) + ndiffSOMO = POPCNT(diffSOMO) + ndiffDOMO = POPCNT(diffDOMO) + nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO) + if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then + ppExistsQ = .TRUE. + EXIT + endif + end do + ! Diagonal part (pp,qq) + if(nholes > 0 .AND. (.NOT. ppExistsQ))then + ! SOMO + NalphaIcfg += 1 + !print *,p,q,"|",holetype(i),vmotype(j),NalphaIcfg + !call debug_spindet(Idomo,1) + !call debug_spindet(Jdomo,1) + alphasIcfg(1,1,NalphaIcfg) = Icfg(1,1) + alphasIcfg(1,2,NalphaIcfg) = Icfg(1,2) + endif + + end subroutine +#+end_src + +** Given an \(\alpha\) CFG, return all the \(|I\rangle\) CFGs + +Next step is to obtain the connected CFGs \(|I\rangle\) that belong to the selected space +given a RI configuration \(|\alpha\rangle\). + +#+begin_src f90 :main no :tangle ../cfgCI/obtain_I_foralpha.irp.f +subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI, nconnectedI, excitationIds, excitationTypes, diagfactors) + implicit none + use bitmasks + BEGIN_DOC + ! Documentation for obtain_connected_I_foralpha + ! This function returns all those selected configurations + ! which are connected to the input configuration + ! Ialpha by a single excitation. + ! + ! The type of excitations are ordered as follows: + ! Type 1 - SOMO -> SOMO + ! Type 2 - DOMO -> VMO + ! Type 3 - SOMO -> VMO + ! Type 4 - DOMO -> SOMO + ! + ! Order of operators + ! \alpha> = a^\dag_p a_q |I> = E_pq |I> + END_DOC + integer ,intent(in) :: idxI + integer(bit_kind),intent(in) :: Ialpha(N_int,2) + integer(bit_kind),intent(out) :: connectedI(N_int,2,*) + integer ,intent(out) :: idxs_connectedI(*) + integer,intent(out) :: nconnectedI + integer,intent(out) :: excitationIds(2,*) + integer,intent(out) :: excitationTypes(*) + real*8 ,intent(out) :: diagfactors(*) + integer*8 :: Idomo + integer*8 :: Isomo + integer*8 :: Jdomo + integer*8 :: Jsomo + integer*8 :: IJsomo + integer*8 :: diffSOMO + integer*8 :: diffDOMO + integer*8 :: xordiffSOMODOMO + integer :: ndiffSOMO + integer :: ndiffDOMO + integer :: nxordiffSOMODOMO + integer :: ii,i,j,k,l,p,q,nsomoJ,nsomoalpha,starti,endi,extyp,nholes + integer :: listholes(mo_num) + integer :: holetype(mo_num) + + nconnectedI = 0 + + p = 0 + q = 0 + do i=idxI,N_configuration + Isomo = Ialpha(1,1) + Idomo = Ialpha(1,2) + Jsomo = psi_configuration(1,1,i) + Jdomo = psi_configuration(1,2,i) + !call debug_spindet(Isomo,1) + !call debug_spindet(Idomo,1) + !print *,"-J--i=",i,Idomo,Jdomo,">",N_configuration + !call debug_spindet(Jsomo,1) + !call debug_spindet(Jdomo,1) + diffSOMO = IEOR(Isomo,Jsomo) + diffDOMO = IEOR(Idomo,Jdomo) + xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) + ndiffSOMO = POPCNT(diffSOMO) + ndiffDOMO = POPCNT(diffDOMO) + nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO) + !print *,"-I--i=",i,ndiffSOMO,ndiffDOMO,nxordiffSOMODOMO!Isomo,Jsomo,ndiffSOMO,ndiffDOMO + !if((ndiffSOMO + ndiffDOMO) .EQ. 0) cycle + !print *,POPCNT(IEOR(diffSOMO,diffDOMO)), ndiffDOMO + !if(POPCNT(IEOR(diffSOMO,diffDOMO)) .LE. 1 .AND. ndiffDOMO .LT. 3) then + if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then + !call debug_spindet(Isomo,1) + !call debug_spindet(Idomo,1) + !print *,"-J--i=",i,Idomo,Jdomo,">",N_configuration + !call debug_spindet(Jsomo,1) + !call debug_spindet(Jdomo,1) + select case(ndiffDOMO) + case (0) + ! SOMO -> VMO + !print *,"obt SOMO -> VMO" + extyp = 3 + IJsomo = IEOR(Isomo, Jsomo) + p = TRAILZ(IAND(Isomo,IJsomo)) + 1 + IJsomo = IBCLR(IJsomo,p-1) + q = TRAILZ(IJsomo) + 1 + case (1) + ! DOMO -> VMO + ! or + ! SOMO -> SOMO + nsomoJ = POPCNT(Jsomo) + nsomoalpha = POPCNT(Isomo) + if(nsomoJ .GT. nsomoalpha) then + ! DOMO -> VMO + !print *,"obt DOMO -> VMO" + extyp = 2 + p = TRAILZ(IEOR(Idomo,Jdomo)) + 1 + Isomo = IEOR(Isomo, Jsomo) + Isomo = IBCLR(Isomo,p-1) + q = TRAILZ(Isomo) + 1 + else + ! SOMO -> SOMO + !print *,"obt SOMO -> SOMO" + extyp = 1 + q = TRAILZ(IEOR(Idomo,Jdomo)) + 1 + Isomo = IEOR(Isomo, Jsomo) + Isomo = IBCLR(Isomo,q-1) + p = TRAILZ(Isomo) + 1 + end if + case (2) + ! DOMO -> SOMO + !print *,"obt DOMO -> SOMO" + extyp = 4 + IJsomo = IEOR(Isomo, Jsomo) + p = TRAILZ(IAND(Jsomo,IJsomo)) + 1 + IJsomo = IBCLR(IJsomo,p-1) + q = TRAILZ(IJsomo) + 1 + case default + print *,"something went wront in get connectedI" + end select + starti = psi_config_data(i,1) + endi = psi_config_data(i,2) + nconnectedI += 1 + connectedI(:,:,nconnectedI) = psi_configuration(:,:,i) + idxs_connectedI(nconnectedI)=starti + excitationIds(1,nconnectedI)=p + excitationIds(2,nconnectedI)=q + excitationTypes(nconnectedI) = extyp + diagfactors(nconnectedI) = 1.0d0 + !print *,"------ > output p,q in obt=",p,q + else if((ndiffSOMO + ndiffDOMO) .EQ. 0) then + ! find out all pq holes possible + nholes = 0 + ! holes in SOMO + Isomo = psi_configuration(1,1,i) + Idomo = psi_configuration(1,2,i) + do ii = 1,mo_num + if(POPCNT(IAND(Isomo,IBSET(0_8,ii-1))) .EQ. 1) then + nholes += 1 + listholes(nholes) = ii + holetype(nholes) = 1 + endif + end do + ! holes in DOMO + do ii = 1,mo_num + if(POPCNT(IAND(Idomo,IBSET(0_8,ii-1))) .EQ. 1) then + nholes += 1 + listholes(nholes) = ii + holetype(nholes) = 2 + endif + end do + + do k=1,nholes + p = listholes(k) + q = p + extyp = 1 + if(holetype(k) .EQ. 1) then + starti = psi_config_data(i,1) + endi = psi_config_data(i,2) + nconnectedI += 1 + connectedI(:,:,nconnectedI) = psi_configuration(:,:,i) + idxs_connectedI(nconnectedI)=starti + excitationIds(1,nconnectedI)=p + excitationIds(2,nconnectedI)=q + excitationTypes(nconnectedI) = extyp + diagfactors(nconnectedI) = 1.0d0 + !print *,"------ > output p,q in obt=",p,q + else + starti = psi_config_data(i,1) + endi = psi_config_data(i,2) + nconnectedI += 1 + connectedI(:,:,nconnectedI) = psi_configuration(:,:,i) + idxs_connectedI(nconnectedI)=starti + excitationIds(1,nconnectedI)=p + excitationIds(2,nconnectedI)=q + excitationTypes(nconnectedI) = extyp + diagfactors(nconnectedI) = 2.0d0 + !print *,"------ > output p,q in obt=",p,q + endif + enddo + endif + end do + +end subroutine obtain_connected_I_foralpha +#+end_src + +#+begin_src fortran + print *,TRAILZ(8) + print *,IBCLR(8,TRAILZ(9)) + print *,TRAILZ(IBCLR(8,TRAILZ(9))) + +#+end_src + +#+RESULTS: +| 3 | +| 8 | +| 3 | + +** Function to get the NSOMOs (seniority) + +#+begin_src f90 :main no :tangle configuration_CI_sigma_helpers.irp.f + function getNSOMO(Icfg) result(NSOMO) + implicit none + integer(bit_kind),intent(in) :: Icfg(N_int,2) + integer :: NSOMO + integer :: i + NSOMO = 0 + do i = 1,N_int + NSOMO += POPCNT(Icfg(i,1)) + enddo + end function getNSOMO +#+end_src + +** Function to convert p,q to model space ids + +This function converts the real orbital ids \(i,j\) to model +space ids \(p,q\) which depend only on the number of somos. + +#+begin_src f90 :main no :tangle configuration_CI_sigma_helpers.irp.f +subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmodel) + implicit none + BEGIN_DOC + ! This function converts the orbital ids + ! in real space to those used in model space + ! in order to identify the matrices required + ! for the calculation of MEs. + ! + ! The type of excitations are ordered as follows: + ! Type 1 - SOMO -> SOMO + ! Type 2 - DOMO -> VMO + ! Type 3 - SOMO -> VMO + ! Type 4 - DOMO -> SOMO + END_DOC + integer(bit_kind),intent(in) :: Ialpha(N_int,2) + integer(bit_kind),intent(in) :: Jcfg(N_int,2) + integer,intent(in) :: p,q + integer,intent(in) :: extype + integer,intent(out) :: pmodel,qmodel + integer*8 :: Isomo + integer*8 :: Idomo + integer*8 :: Jsomo + integer*8 :: Jdomo + integer*8 :: mask + integer*8 :: Isomotmp + integer*8 :: Jsomotmp + integer :: pos0,pos0prev + + ! TODO Flag (print) when model space indices is > 64 + Isomo = Ialpha(1,1) + Idomo = Ialpha(1,2) + Jsomo = Jcfg(1,1) + Jdomo = Jcfg(1,2) + pos0prev = 0 + pmodel = p + qmodel = q + + if(p .EQ. q) then + !print *,"input pq=",p,q,"extype=",extype + pmodel = 1 + qmodel = 1 + else + !print *,"input pq=",p,q,"extype=",extype + !call debug_spindet(Isomo,1) + !call debug_spindet(Idomo,1) + !call debug_spindet(Jsomo,1) + !call debug_spindet(Jdomo,1) + select case(extype) + case (1) + ! SOMO -> SOMO + ! remove all domos + !print *,"type -> SOMO -> SOMO" + mask = ISHFT(1_8,p) - 1 + Isomotmp = IAND(Isomo,mask) + pmodel = POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) + mask = ISHFT(1_8,q) - 1 + Isomotmp = IAND(Isomo,mask) + qmodel = POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) + case (2) + ! DOMO -> VMO + ! remove all domos except one at p + !print *,"type -> DOMO -> VMO" + mask = ISHFT(1_8,p) - 1 + Jsomotmp = IAND(Jsomo,mask) + pmodel = POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) + mask = ISHFT(1_8,q) - 1 + Jsomotmp = IAND(Jsomo,mask) + qmodel = POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) + case (3) + ! SOMO -> VMO + !print *,"type -> SOMO -> VMO" + !Isomo = IEOR(Isomo,Jsomo) + if(p.LT.q) then + mask = ISHFT(1_8,p) - 1 + Isomo = IAND(Isomo,mask) + pmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) + mask = ISHFT(1_8,q) - 1 + Jsomo = IAND(Jsomo,mask) + qmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) + 1 + else + mask = ISHFT(1_8,p) - 1 + Isomo = IAND(Isomo,mask) + pmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) + 1 + mask = ISHFT(1_8,q) - 1 + Jsomo = IAND(Jsomo,mask) + qmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) + endif + case (4) + ! DOMO -> SOMO + ! remove all domos except one at p + !print *,"type -> DOMO -> SOMO" + !Isomo = IEOR(Isomo,Jsomo) + if(p.LT.q) then + mask = ISHFT(1_8,p) - 1 + Jsomo = IAND(Jsomo,mask) + pmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) + mask = ISHFT(1_8,q) - 1 + Isomo = IAND(Isomo,mask) + qmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) + 1 + else + mask = ISHFT(1_8,p) - 1 + Jsomo = IAND(Jsomo,mask) + pmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) + 1 + mask = ISHFT(1_8,q) - 1 + Isomo = IAND(Isomo,mask) + qmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) + endif + case default + print *,"something is wrong in convertOrbIdsToModelSpaceIds" + end select + endif + !print *,p,q,"model ids=",pmodel,qmodel +end subroutine convertOrbIdsToModelSpaceIds +#+end_src + +#+begin_src fortran + integer :: i + integer :: count + integer :: mask + integer :: isomo + count = 0 + mask = ISHFT(1_8,5)-1 + print *,mask + print *,POPCNT(mask) + isomo = 144 + isomo = IAND(isomo,mask) + print *,isomo + print *,XOR(isomo,mask) + print *,POPCNT(mask) - POPCNT(XOR(isomo,mask)) + +#+end_src + +#+RESULTS: +| 31 | +| 5 | +| 16 | +| 15 | +| 1 | + +#+begin_src fortran + print *,IBSET(0_8,4)-1 + print *,POPCNT(IBSET(0_8,4)-1) - POPCNT(IAND(716,IBSET(0_8,4)-1)) + print *,POPCNT(IBSET(0_8,8)-1) - POPCNT(IAND(716,IBSET(0_8,8)-1)) +#+end_src + +#+RESULTS: +| 15 | +| 2 | +| 4 | From 4f5f612c5ddf2d2c2cbb157bd2fea68b36ce397f Mon Sep 17 00:00:00 2001 From: v1j4y Date: Mon, 1 Mar 2021 12:29:26 +0100 Subject: [PATCH 2/7] Added modifs to configuration_CI_sigma_helpres.org --- src/csf/configuration_CI_sigma_helpers.org | 110 +-- .../configuration_CI_sigma_helpers.org | 657 ------------------ 2 files changed, 69 insertions(+), 698 deletions(-) delete mode 100644 src/determinants/configuration_CI_sigma_helpers.org diff --git a/src/csf/configuration_CI_sigma_helpers.org b/src/csf/configuration_CI_sigma_helpers.org index 6677139a..88b27a7a 100644 --- a/src/csf/configuration_CI_sigma_helpers.org +++ b/src/csf/configuration_CI_sigma_helpers.org @@ -168,6 +168,12 @@ the input determinant \(|D_I\rangle\). !if(POPCNT(IEOR(diffSOMO,diffDOMO)) .LE. 1 .AND. ndiffDOMO .LT. 3) then if((ndiffSOMO+ndiffDOMO) .EQ. 0) then pqAlreadyGenQ = .TRUE. + ppExistsQ = .TRUE. + EXIT + endif + if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then + pqAlreadyGenQ = .TRUE. + !ppExistsQ = .TRUE. !print *,i,k,ndiffSOMO,ndiffDOMO !call debug_spindet(Jsomo,1) !call debug_spindet(Jdomo,1) @@ -183,16 +189,16 @@ the input determinant \(|D_I\rangle\). pqExistsQ = .FALSE. ! now check if this exists in the selected list - do k = idxI+1, N_configuration - diffSOMO = IEOR(OR(reunion_of_act_virt_bitmask(1,1),Jsomo),psi_configuration(1,1,k)) - diffDOMO = IEOR(OR(reunion_of_act_virt_bitmask(1,1),Jdomo),psi_configuration(1,2,k)) - ndiffSOMO = POPCNT(diffSOMO) - ndiffDOMO = POPCNT(diffDOMO) - if((ndiffSOMO + ndiffDOMO) .EQ. 0) then - pqExistsQ = .TRUE. - EXIT - endif - end do + !do k = idxI+1, N_configuration + ! diffSOMO = IEOR(OR(reunion_of_act_virt_bitmask(1,1),Jsomo),psi_configuration(1,1,k)) + ! diffDOMO = IEOR(OR(reunion_of_act_virt_bitmask(1,1),Jdomo),psi_configuration(1,2,k)) + ! ndiffSOMO = POPCNT(diffSOMO) + ! ndiffDOMO = POPCNT(diffDOMO) + ! if((ndiffSOMO + ndiffDOMO) .EQ. 0) then + ! pqExistsQ = .TRUE. + ! EXIT + ! endif + !end do if(.NOT. pqExistsQ) then tableUniqueAlphas(p,q) = .TRUE. @@ -242,6 +248,7 @@ the input determinant \(|D_I\rangle\). print*,"Something went wrong in obtain_associated_alphaI" endif + ! SOMO NalphaIcfg += 1 !print *,p,q,"|",holetype(i),vmotype(j),NalphaIcfg alphasIcfg(1,1,NalphaIcfg) = Jsomo @@ -250,6 +257,33 @@ the input determinant \(|D_I\rangle\). end do end do + ! Check if this Icfg has been previously generated as a mono + ppExistsQ = .False. + Isomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,1)) + Idomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,2)) + do k = 1, idxI-1 + diffSOMO = IEOR(Isomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k))) + diffDOMO = IEOR(Idomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k))) + xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) + ndiffSOMO = POPCNT(diffSOMO) + ndiffDOMO = POPCNT(diffDOMO) + nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO) + if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then + ppExistsQ = .TRUE. + EXIT + endif + end do + ! Diagonal part (pp,qq) + if(nholes > 0 .AND. (.NOT. ppExistsQ))then + ! SOMO + NalphaIcfg += 1 + !print *,p,q,"|",holetype(i),vmotype(j),NalphaIcfg + !call debug_spindet(Idomo,1) + !call debug_spindet(Jdomo,1) + alphasIcfg(1,1,NalphaIcfg) = Icfg(1,1) + alphasIcfg(1,2,NalphaIcfg) = Icfg(1,2) + endif + end subroutine #+end_src @@ -259,7 +293,7 @@ Next step is to obtain the connected CFGs \(|I\rangle\) that belong to the selec given a RI configuration \(|\alpha\rangle\). #+begin_src f90 :main no :tangle ../cfgCI/obtain_I_foralpha.irp.f -subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI, nconnectedI, excitationIds, excitationTypes) +subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI, nconnectedI, excitationIds, excitationTypes, diagfactors) implicit none use bitmasks BEGIN_DOC @@ -284,6 +318,7 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI integer,intent(out) :: nconnectedI integer,intent(out) :: excitationIds(2,*) integer,intent(out) :: excitationTypes(*) + real*8 ,intent(out) :: diagfactors(*) integer*8 :: Idomo integer*8 :: Isomo integer*8 :: Jdomo @@ -291,34 +326,14 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI integer*8 :: IJsomo integer*8 :: diffSOMO integer*8 :: diffDOMO + integer*8 :: xordiffSOMODOMO integer :: ndiffSOMO integer :: ndiffDOMO integer :: nxordiffSOMODOMO - integer :: i,j,k,l,p,q,nsomoJ,nsomoalpha,starti,endi,extyp,nholes + integer :: ii,i,j,k,l,p,q,nsomoJ,nsomoalpha,starti,endi,extyp,nholes integer :: listholes(mo_num) integer :: holetype(mo_num) - ! find out all pq holes possible - nholes = 0 - ! holes in SOMO - Isomo = psi_configuration(1,1,idxI) - Idomo = psi_configuration(1,2,idxI) - do i = n_core_orb+1,n_core_orb + n_act_orb - if(POPCNT(IAND(Isomo,IBSET(0_8,i-1))) .EQ. 1) then - nholes += 1 - listholes(nholes) = i - holetype(nholes) = 1 - endif - end do - ! holes in DOMO - do i = n_core_orb+1,n_core_orb + n_act_orb - if(POPCNT(IAND(Idomo,IBSET(0_8,i-1))) .EQ. 1) then - nholes += 1 - listholes(nholes) = i - holetype(nholes) = 2 - endif - end do - nconnectedI = 0 p = 0 @@ -335,11 +350,15 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI !call debug_spindet(Jdomo,1) diffSOMO = IEOR(Isomo,Jsomo) diffDOMO = IEOR(Idomo,Jdomo) + xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) ndiffSOMO = POPCNT(diffSOMO) ndiffDOMO = POPCNT(diffDOMO) nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO) + !print *,"-I--i=",i,ndiffSOMO,ndiffDOMO,nxordiffSOMODOMO!Isomo,Jsomo,ndiffSOMO,ndiffDOMO + !if((ndiffSOMO + ndiffDOMO) .EQ. 0) cycle !print *,POPCNT(IEOR(diffSOMO,diffDOMO)), ndiffDOMO - if(POPCNT(IEOR(diffSOMO,diffDOMO)) .LE. 1 .AND. ndiffDOMO .LT. 3) then + !if(POPCNT(IEOR(diffSOMO,diffDOMO)) .LE. 1 .AND. ndiffDOMO .LT. 3) then + if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then !call debug_spindet(Isomo,1) !call debug_spindet(Idomo,1) !print *,"-J--i=",i,Idomo,Jdomo,">",N_configuration @@ -565,17 +584,26 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod pmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) mask = ISHFT(1_8,q) - 1 Jsomo = IAND(Jsomo,mask) + qmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) + 1 + else + mask = ISHFT(1_8,p) - 1 + Isomo = IAND(Isomo,mask) + pmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) + 1 + mask = ISHFT(1_8,q) - 1 + Jsomo = IAND(Jsomo,mask) + qmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) + endif case (4) ! DOMO -> SOMO ! remove all domos except one at p !print *,"type -> DOMO -> SOMO" !Isomo = IEOR(Isomo,Jsomo) - mask = ISHFT(1_8,p) - 1 - Jsomo = IAND(Jsomo,mask) - pmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) - mask = ISHFT(1_8,q) - 1 - Isomo = IAND(Isomo,mask) - qmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) + if(p.LT.q) then + mask = ISHFT(1_8,p) - 1 + Jsomo = IAND(Jsomo,mask) + pmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) + mask = ISHFT(1_8,q) - 1 + Isomo = IAND(Isomo,mask) qmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) + 1 else mask = ISHFT(1_8,p) - 1 @@ -589,7 +617,7 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod print *,"something is wrong in convertOrbIdsToModelSpaceIds" end select endif - !print *,p,q,"model ids=",pmodel,qmodel + !print *,p,q,"model ids=",pmodel,qmodel end subroutine convertOrbIdsToModelSpaceIds #+end_src diff --git a/src/determinants/configuration_CI_sigma_helpers.org b/src/determinants/configuration_CI_sigma_helpers.org deleted file mode 100644 index 88b27a7a..00000000 --- a/src/determinants/configuration_CI_sigma_helpers.org +++ /dev/null @@ -1,657 +0,0 @@ -#+title: Configuration Sigma Vector Helpers -#+author: Vijay Gopal Chilkuri -#+email: vijay.gopal.c@gmail.com - -* Generate the singly excited configurations on-the-fly - -The algorithm is based on the work by Garniron et. al. (see thesis Chap 5). -The basic idea is to generate \(|\alpha\rangle\)'s on-the-fly. - -The algorithm is based on the idea of splitting the list of \(|\alpha\rangle\)'s -into blocks associated with a selected determinant \(|D_I\rangle\). - -** Create a function to return a list of alphas - -Here we create a list of \(|\alpha\rangle\)'s associated with -the input determinant \(|D_I\rangle\). - -#+begin_src f90 :main no :tangle configuration_CI_sigma_helpers.irp.f - subroutine obtain_associated_alphaI(idxI, Icfg, alphasIcfg, NalphaIcfg, factor_alphaI) - implicit none - use bitmasks - BEGIN_DOC - ! Documentation for alphasI - ! Returns the associated alpha's for - ! the input configuration Icfg. - END_DOC - - integer,intent(in) :: idxI ! The id of the Ith CFG - integer(bit_kind),intent(in) :: Icfg(N_int,2) - integer,intent(out) :: NalphaIcfg - real*8 ,intent(out) :: factor_alphaI(*) - integer(bit_kind),intent(out) :: alphasIcfg(N_int,2,*) - logical,dimension(:,:),allocatable :: tableUniqueAlphas - integer :: listholes(mo_num) - integer :: holetype(mo_num) ! 1-> SOMO 2->DOMO - integer :: nholes - integer :: nvmos - integer :: listvmos(mo_num) - integer :: vmotype(mo_num) ! 1 -> VMO 2 -> SOMO - integer*8 :: Idomo - integer*8 :: Isomo - integer*8 :: Jdomo - integer*8 :: Jsomo - integer*8 :: diffSOMO - integer*8 :: diffDOMO - integer*8 :: xordiffSOMODOMO - integer :: ndiffSOMO - integer :: ndiffDOMO - integer :: nxordiffSOMODOMO - integer :: ndiffAll - integer :: i - integer :: j - integer :: k - integer :: hole - integer :: p - integer :: q - integer :: countalphas - logical :: pqAlreadyGenQ - logical :: pqExistsQ - logical :: ppExistsQ - Isomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,1)) - Idomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,2)) - !print*,"Input cfg" - !call debug_spindet(Isomo,1) - !call debug_spindet(Idomo,1) - - !print*,n_act_orb, "monum=",mo_num," n_core=",n_core_orb - - ! find out all pq holes possible - nholes = 0 - ! holes in SOMO - do i = 1,mo_num - if(POPCNT(IAND(Isomo,IBSET(0_8,i-1))) .EQ. 1) then - nholes += 1 - listholes(nholes) = i - holetype(nholes) = 1 - endif - end do - ! holes in DOMO - do i = 1,mo_num - if(POPCNT(IAND(Idomo,IBSET(0_8,i-1))) .EQ. 1) then - nholes += 1 - listholes(nholes) = i - holetype(nholes) = 2 - endif - end do - - ! find vmos - listvmos = -1 - vmotype = -1 - nvmos = 0 - do i = 1,mo_num - !print *,i,IBSET(0,i-1),POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))), POPCNT(IAND(Idomo,(IBSET(0_8,i-1)))) - if(POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))) .EQ. 0 .AND. POPCNT(IAND(Idomo,(IBSET(0_8,i-1)))) .EQ. 0) then - nvmos += 1 - listvmos(nvmos) = i - vmotype(nvmos) = 1 - else if(POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))) .EQ. 1 .AND. POPCNT(IAND(Idomo,(IBSET(0_8,i-1)))) .EQ. 0 ) then - nvmos += 1 - listvmos(nvmos) = i - vmotype(nvmos) = 2 - end if - end do - - !print *,"Nvmo=",nvmos - !print *,listvmos - !print *,vmotype - - allocate(tableUniqueAlphas(mo_num,mo_num)) - tableUniqueAlphas = .FALSE. - - ! Now find the allowed (p,q) excitations - Isomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,1)) - Idomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,2)) - !print *,"Isomo" - !call debug_spindet(Isomo,1) - !call debug_spindet(Idomo,1) - - !print *,"Nholes=",nholes," Nvmos=",nvmos, " idxi=",idxI - !do i = 1,nholes - ! print *,i,"->",listholes(i) - !enddo - !do i = 1,nvmos - ! print *,i,"->",listvmos(i) - !enddo - - ! TODO cfg_seniority_index - do i = 1,nholes - p = listholes(i) - do j = 1,nvmos - q = listvmos(j) - if(p .EQ. q) cycle - if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 1) then - ! SOMO -> VMO - Jsomo = IBCLR(Isomo,p-1) - Jsomo = IBSET(Jsomo,q-1) - Jdomo = Idomo - else if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 2) then - ! SOMO -> SOMO - Jsomo = IBCLR(Isomo,p-1) - Jsomo = IBCLR(Jsomo,q-1) - Jdomo = IBSET(Idomo,q-1) - else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then - ! DOMO -> VMO - Jsomo = IBSET(Isomo,p-1) - Jsomo = IBSET(Jsomo,q-1) - Jdomo = IBCLR(Idomo,p-1) - else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 2) then - ! DOMO -> SOMO - Jsomo = IBSET(Isomo,p-1) - Jsomo = IBCLR(Jsomo,q-1) - Jdomo = IBCLR(Idomo,p-1) - Jdomo = IBSET(Jdomo,q-1) - else - print*,"Something went wrong in obtain_associated_alphaI" - endif - - - pqAlreadyGenQ = .FALSE. - ! First check if it can be generated before - do k = 1, idxI-1 - diffSOMO = IEOR(Jsomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k))) - diffDOMO = IEOR(Jdomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k))) - xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) - ndiffSOMO = POPCNT(diffSOMO) - ndiffDOMO = POPCNT(diffDOMO) - nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO) - !if(POPCNT(IEOR(diffSOMO,diffDOMO)) .LE. 1 .AND. ndiffDOMO .LT. 3) then - if((ndiffSOMO+ndiffDOMO) .EQ. 0) then - pqAlreadyGenQ = .TRUE. - ppExistsQ = .TRUE. - EXIT - endif - if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then - pqAlreadyGenQ = .TRUE. - !ppExistsQ = .TRUE. - !print *,i,k,ndiffSOMO,ndiffDOMO - !call debug_spindet(Jsomo,1) - !call debug_spindet(Jdomo,1) - !call debug_spindet(iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k)),1) - !call debug_spindet(iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k)),1) - EXIT - endif - end do - - !print *,"(,",p,",",q,")",pqAlreadyGenQ - - if(pqAlreadyGenQ) cycle - - pqExistsQ = .FALSE. - ! now check if this exists in the selected list - !do k = idxI+1, N_configuration - ! diffSOMO = IEOR(OR(reunion_of_act_virt_bitmask(1,1),Jsomo),psi_configuration(1,1,k)) - ! diffDOMO = IEOR(OR(reunion_of_act_virt_bitmask(1,1),Jdomo),psi_configuration(1,2,k)) - ! ndiffSOMO = POPCNT(diffSOMO) - ! ndiffDOMO = POPCNT(diffDOMO) - ! if((ndiffSOMO + ndiffDOMO) .EQ. 0) then - ! pqExistsQ = .TRUE. - ! EXIT - ! endif - !end do - - if(.NOT. pqExistsQ) then - tableUniqueAlphas(p,q) = .TRUE. - !print *,p,q - !call debug_spindet(Jsomo,1) - !call debug_spindet(Jdomo,1) - endif - end do - end do - - !print *,tableUniqueAlphas(:,:) - - ! prune list of alphas - Isomo = Icfg(1,1) - Idomo = Icfg(1,2) - Jsomo = Icfg(1,1) - Jdomo = Icfg(1,2) - NalphaIcfg = 0 - do i = 1, nholes - p = listholes(i) - do j = 1, nvmos - q = listvmos(j) - if(p .EQ. q) cycle - if(tableUniqueAlphas(p,q)) then - if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 1) then - ! SOMO -> VMO - Jsomo = IBCLR(Isomo,p-1) - Jsomo = IBSET(Jsomo,q-1) - Jdomo = Idomo - else if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 2) then - ! SOMO -> SOMO - Jsomo = IBCLR(Isomo,p-1) - Jsomo = IBCLR(Jsomo,q-1) - Jdomo = IBSET(Idomo,q-1) - else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then - ! DOMO -> VMO - Jsomo = IBSET(Isomo,p-1) - Jsomo = IBSET(Jsomo,q-1) - Jdomo = IBCLR(Idomo,p-1) - else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 2) then - ! DOMO -> SOMO - Jsomo = IBSET(Isomo,p-1) - Jsomo = IBCLR(Jsomo,q-1) - Jdomo = IBCLR(Idomo,p-1) - Jdomo = IBSET(Jdomo,q-1) - else - print*,"Something went wrong in obtain_associated_alphaI" - endif - - ! SOMO - NalphaIcfg += 1 - !print *,p,q,"|",holetype(i),vmotype(j),NalphaIcfg - alphasIcfg(1,1,NalphaIcfg) = Jsomo - alphasIcfg(1,2,NalphaIcfg) = IOR(Jdomo,ISHFT(1_8,n_core_orb)-1) - endif - end do - end do - - ! Check if this Icfg has been previously generated as a mono - ppExistsQ = .False. - Isomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,1)) - Idomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,2)) - do k = 1, idxI-1 - diffSOMO = IEOR(Isomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k))) - diffDOMO = IEOR(Idomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k))) - xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) - ndiffSOMO = POPCNT(diffSOMO) - ndiffDOMO = POPCNT(diffDOMO) - nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO) - if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then - ppExistsQ = .TRUE. - EXIT - endif - end do - ! Diagonal part (pp,qq) - if(nholes > 0 .AND. (.NOT. ppExistsQ))then - ! SOMO - NalphaIcfg += 1 - !print *,p,q,"|",holetype(i),vmotype(j),NalphaIcfg - !call debug_spindet(Idomo,1) - !call debug_spindet(Jdomo,1) - alphasIcfg(1,1,NalphaIcfg) = Icfg(1,1) - alphasIcfg(1,2,NalphaIcfg) = Icfg(1,2) - endif - - end subroutine -#+end_src - -** Given an \(\alpha\) CFG, return all the \(|I\rangle\) CFGs - -Next step is to obtain the connected CFGs \(|I\rangle\) that belong to the selected space -given a RI configuration \(|\alpha\rangle\). - -#+begin_src f90 :main no :tangle ../cfgCI/obtain_I_foralpha.irp.f -subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI, nconnectedI, excitationIds, excitationTypes, diagfactors) - implicit none - use bitmasks - BEGIN_DOC - ! Documentation for obtain_connected_I_foralpha - ! This function returns all those selected configurations - ! which are connected to the input configuration - ! Ialpha by a single excitation. - ! - ! The type of excitations are ordered as follows: - ! Type 1 - SOMO -> SOMO - ! Type 2 - DOMO -> VMO - ! Type 3 - SOMO -> VMO - ! Type 4 - DOMO -> SOMO - ! - ! Order of operators - ! \alpha> = a^\dag_p a_q |I> = E_pq |I> - END_DOC - integer ,intent(in) :: idxI - integer(bit_kind),intent(in) :: Ialpha(N_int,2) - integer(bit_kind),intent(out) :: connectedI(N_int,2,*) - integer ,intent(out) :: idxs_connectedI(*) - integer,intent(out) :: nconnectedI - integer,intent(out) :: excitationIds(2,*) - integer,intent(out) :: excitationTypes(*) - real*8 ,intent(out) :: diagfactors(*) - integer*8 :: Idomo - integer*8 :: Isomo - integer*8 :: Jdomo - integer*8 :: Jsomo - integer*8 :: IJsomo - integer*8 :: diffSOMO - integer*8 :: diffDOMO - integer*8 :: xordiffSOMODOMO - integer :: ndiffSOMO - integer :: ndiffDOMO - integer :: nxordiffSOMODOMO - integer :: ii,i,j,k,l,p,q,nsomoJ,nsomoalpha,starti,endi,extyp,nholes - integer :: listholes(mo_num) - integer :: holetype(mo_num) - - nconnectedI = 0 - - p = 0 - q = 0 - do i=idxI,N_configuration - Isomo = Ialpha(1,1) - Idomo = Ialpha(1,2) - Jsomo = psi_configuration(1,1,i) - Jdomo = psi_configuration(1,2,i) - !call debug_spindet(Isomo,1) - !call debug_spindet(Idomo,1) - !print *,"-J--i=",i,Idomo,Jdomo,">",N_configuration - !call debug_spindet(Jsomo,1) - !call debug_spindet(Jdomo,1) - diffSOMO = IEOR(Isomo,Jsomo) - diffDOMO = IEOR(Idomo,Jdomo) - xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) - ndiffSOMO = POPCNT(diffSOMO) - ndiffDOMO = POPCNT(diffDOMO) - nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO) - !print *,"-I--i=",i,ndiffSOMO,ndiffDOMO,nxordiffSOMODOMO!Isomo,Jsomo,ndiffSOMO,ndiffDOMO - !if((ndiffSOMO + ndiffDOMO) .EQ. 0) cycle - !print *,POPCNT(IEOR(diffSOMO,diffDOMO)), ndiffDOMO - !if(POPCNT(IEOR(diffSOMO,diffDOMO)) .LE. 1 .AND. ndiffDOMO .LT. 3) then - if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then - !call debug_spindet(Isomo,1) - !call debug_spindet(Idomo,1) - !print *,"-J--i=",i,Idomo,Jdomo,">",N_configuration - !call debug_spindet(Jsomo,1) - !call debug_spindet(Jdomo,1) - select case(ndiffDOMO) - case (0) - ! SOMO -> VMO - !print *,"obt SOMO -> VMO" - extyp = 3 - IJsomo = IEOR(Isomo, Jsomo) - p = TRAILZ(IAND(Isomo,IJsomo)) + 1 - IJsomo = IBCLR(IJsomo,p-1) - q = TRAILZ(IJsomo) + 1 - case (1) - ! DOMO -> VMO - ! or - ! SOMO -> SOMO - nsomoJ = POPCNT(Jsomo) - nsomoalpha = POPCNT(Isomo) - if(nsomoJ .GT. nsomoalpha) then - ! DOMO -> VMO - !print *,"obt DOMO -> VMO" - extyp = 2 - p = TRAILZ(IEOR(Idomo,Jdomo)) + 1 - Isomo = IEOR(Isomo, Jsomo) - Isomo = IBCLR(Isomo,p-1) - q = TRAILZ(Isomo) + 1 - else - ! SOMO -> SOMO - !print *,"obt SOMO -> SOMO" - extyp = 1 - q = TRAILZ(IEOR(Idomo,Jdomo)) + 1 - Isomo = IEOR(Isomo, Jsomo) - Isomo = IBCLR(Isomo,q-1) - p = TRAILZ(Isomo) + 1 - end if - case (2) - ! DOMO -> SOMO - !print *,"obt DOMO -> SOMO" - extyp = 4 - IJsomo = IEOR(Isomo, Jsomo) - p = TRAILZ(IAND(Jsomo,IJsomo)) + 1 - IJsomo = IBCLR(IJsomo,p-1) - q = TRAILZ(IJsomo) + 1 - case default - print *,"something went wront in get connectedI" - end select - starti = psi_config_data(i,1) - endi = psi_config_data(i,2) - nconnectedI += 1 - connectedI(:,:,nconnectedI) = psi_configuration(:,:,i) - idxs_connectedI(nconnectedI)=starti - excitationIds(1,nconnectedI)=p - excitationIds(2,nconnectedI)=q - excitationTypes(nconnectedI) = extyp - diagfactors(nconnectedI) = 1.0d0 - !print *,"------ > output p,q in obt=",p,q - else if((ndiffSOMO + ndiffDOMO) .EQ. 0) then - ! find out all pq holes possible - nholes = 0 - ! holes in SOMO - Isomo = psi_configuration(1,1,i) - Idomo = psi_configuration(1,2,i) - do ii = 1,mo_num - if(POPCNT(IAND(Isomo,IBSET(0_8,ii-1))) .EQ. 1) then - nholes += 1 - listholes(nholes) = ii - holetype(nholes) = 1 - endif - end do - ! holes in DOMO - do ii = 1,mo_num - if(POPCNT(IAND(Idomo,IBSET(0_8,ii-1))) .EQ. 1) then - nholes += 1 - listholes(nholes) = ii - holetype(nholes) = 2 - endif - end do - - do k=1,nholes - p = listholes(k) - q = p - extyp = 1 - if(holetype(k) .EQ. 1) then - starti = psi_config_data(i,1) - endi = psi_config_data(i,2) - nconnectedI += 1 - connectedI(:,:,nconnectedI) = psi_configuration(:,:,i) - idxs_connectedI(nconnectedI)=starti - excitationIds(1,nconnectedI)=p - excitationIds(2,nconnectedI)=q - excitationTypes(nconnectedI) = extyp - diagfactors(nconnectedI) = 1.0d0 - !print *,"------ > output p,q in obt=",p,q - else - starti = psi_config_data(i,1) - endi = psi_config_data(i,2) - nconnectedI += 1 - connectedI(:,:,nconnectedI) = psi_configuration(:,:,i) - idxs_connectedI(nconnectedI)=starti - excitationIds(1,nconnectedI)=p - excitationIds(2,nconnectedI)=q - excitationTypes(nconnectedI) = extyp - diagfactors(nconnectedI) = 2.0d0 - !print *,"------ > output p,q in obt=",p,q - endif - enddo - endif - end do - -end subroutine obtain_connected_I_foralpha -#+end_src - -#+begin_src fortran - print *,TRAILZ(8) - print *,IBCLR(8,TRAILZ(9)) - print *,TRAILZ(IBCLR(8,TRAILZ(9))) - -#+end_src - -#+RESULTS: -| 3 | -| 8 | -| 3 | - -** Function to get the NSOMOs (seniority) - -#+begin_src f90 :main no :tangle configuration_CI_sigma_helpers.irp.f - function getNSOMO(Icfg) result(NSOMO) - implicit none - integer(bit_kind),intent(in) :: Icfg(N_int,2) - integer :: NSOMO - integer :: i - NSOMO = 0 - do i = 1,N_int - NSOMO += POPCNT(Icfg(i,1)) - enddo - end function getNSOMO -#+end_src - -** Function to convert p,q to model space ids - -This function converts the real orbital ids \(i,j\) to model -space ids \(p,q\) which depend only on the number of somos. - -#+begin_src f90 :main no :tangle configuration_CI_sigma_helpers.irp.f -subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmodel) - implicit none - BEGIN_DOC - ! This function converts the orbital ids - ! in real space to those used in model space - ! in order to identify the matrices required - ! for the calculation of MEs. - ! - ! The type of excitations are ordered as follows: - ! Type 1 - SOMO -> SOMO - ! Type 2 - DOMO -> VMO - ! Type 3 - SOMO -> VMO - ! Type 4 - DOMO -> SOMO - END_DOC - integer(bit_kind),intent(in) :: Ialpha(N_int,2) - integer(bit_kind),intent(in) :: Jcfg(N_int,2) - integer,intent(in) :: p,q - integer,intent(in) :: extype - integer,intent(out) :: pmodel,qmodel - integer*8 :: Isomo - integer*8 :: Idomo - integer*8 :: Jsomo - integer*8 :: Jdomo - integer*8 :: mask - integer*8 :: Isomotmp - integer*8 :: Jsomotmp - integer :: pos0,pos0prev - - ! TODO Flag (print) when model space indices is > 64 - Isomo = Ialpha(1,1) - Idomo = Ialpha(1,2) - Jsomo = Jcfg(1,1) - Jdomo = Jcfg(1,2) - pos0prev = 0 - pmodel = p - qmodel = q - - if(p .EQ. q) then - !print *,"input pq=",p,q,"extype=",extype - pmodel = 1 - qmodel = 1 - else - !print *,"input pq=",p,q,"extype=",extype - !call debug_spindet(Isomo,1) - !call debug_spindet(Idomo,1) - !call debug_spindet(Jsomo,1) - !call debug_spindet(Jdomo,1) - select case(extype) - case (1) - ! SOMO -> SOMO - ! remove all domos - !print *,"type -> SOMO -> SOMO" - mask = ISHFT(1_8,p) - 1 - Isomotmp = IAND(Isomo,mask) - pmodel = POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) - mask = ISHFT(1_8,q) - 1 - Isomotmp = IAND(Isomo,mask) - qmodel = POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) - case (2) - ! DOMO -> VMO - ! remove all domos except one at p - !print *,"type -> DOMO -> VMO" - mask = ISHFT(1_8,p) - 1 - Jsomotmp = IAND(Jsomo,mask) - pmodel = POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) - mask = ISHFT(1_8,q) - 1 - Jsomotmp = IAND(Jsomo,mask) - qmodel = POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) - case (3) - ! SOMO -> VMO - !print *,"type -> SOMO -> VMO" - !Isomo = IEOR(Isomo,Jsomo) - if(p.LT.q) then - mask = ISHFT(1_8,p) - 1 - Isomo = IAND(Isomo,mask) - pmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) - mask = ISHFT(1_8,q) - 1 - Jsomo = IAND(Jsomo,mask) - qmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) + 1 - else - mask = ISHFT(1_8,p) - 1 - Isomo = IAND(Isomo,mask) - pmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) + 1 - mask = ISHFT(1_8,q) - 1 - Jsomo = IAND(Jsomo,mask) - qmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) - endif - case (4) - ! DOMO -> SOMO - ! remove all domos except one at p - !print *,"type -> DOMO -> SOMO" - !Isomo = IEOR(Isomo,Jsomo) - if(p.LT.q) then - mask = ISHFT(1_8,p) - 1 - Jsomo = IAND(Jsomo,mask) - pmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) - mask = ISHFT(1_8,q) - 1 - Isomo = IAND(Isomo,mask) - qmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) + 1 - else - mask = ISHFT(1_8,p) - 1 - Jsomo = IAND(Jsomo,mask) - pmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) + 1 - mask = ISHFT(1_8,q) - 1 - Isomo = IAND(Isomo,mask) - qmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) - endif - case default - print *,"something is wrong in convertOrbIdsToModelSpaceIds" - end select - endif - !print *,p,q,"model ids=",pmodel,qmodel -end subroutine convertOrbIdsToModelSpaceIds -#+end_src - -#+begin_src fortran - integer :: i - integer :: count - integer :: mask - integer :: isomo - count = 0 - mask = ISHFT(1_8,5)-1 - print *,mask - print *,POPCNT(mask) - isomo = 144 - isomo = IAND(isomo,mask) - print *,isomo - print *,XOR(isomo,mask) - print *,POPCNT(mask) - POPCNT(XOR(isomo,mask)) - -#+end_src - -#+RESULTS: -| 31 | -| 5 | -| 16 | -| 15 | -| 1 | - -#+begin_src fortran - print *,IBSET(0_8,4)-1 - print *,POPCNT(IBSET(0_8,4)-1) - POPCNT(IAND(716,IBSET(0_8,4)-1)) - print *,POPCNT(IBSET(0_8,8)-1) - POPCNT(IAND(716,IBSET(0_8,8)-1)) -#+end_src - -#+RESULTS: -| 15 | -| 2 | -| 4 | From 40a3a30ea17c01e7676f68eca99ac232258bd958 Mon Sep 17 00:00:00 2001 From: v1j4y Date: Mon, 1 Mar 2021 12:41:42 +0100 Subject: [PATCH 3/7] Removed deletion of domos in get_phase. --- src/csf/configuration_CI_sigma_helpers.irp.f | 315 ++++++++++++++++++- src/csf/sigma_vector.irp.f | 5 - 2 files changed, 302 insertions(+), 18 deletions(-) diff --git a/src/csf/configuration_CI_sigma_helpers.irp.f b/src/csf/configuration_CI_sigma_helpers.irp.f index f73362eb..eb532b03 100644 --- a/src/csf/configuration_CI_sigma_helpers.irp.f +++ b/src/csf/configuration_CI_sigma_helpers.irp.f @@ -1,3 +1,273 @@ + subroutine obtain_associated_alphaI(idxI, Icfg, alphasIcfg, NalphaIcfg, factor_alphaI) + implicit none + use bitmasks + BEGIN_DOC + ! Documentation for alphasI + ! Returns the associated alpha's for + ! the input configuration Icfg. + END_DOC + + integer,intent(in) :: idxI ! The id of the Ith CFG + integer(bit_kind),intent(in) :: Icfg(N_int,2) + integer,intent(out) :: NalphaIcfg + real*8 ,intent(out) :: factor_alphaI(*) + integer(bit_kind),intent(out) :: alphasIcfg(N_int,2,*) + logical,dimension(:,:),allocatable :: tableUniqueAlphas + integer :: listholes(mo_num) + integer :: holetype(mo_num) ! 1-> SOMO 2->DOMO + integer :: nholes + integer :: nvmos + integer :: listvmos(mo_num) + integer :: vmotype(mo_num) ! 1 -> VMO 2 -> SOMO + integer*8 :: Idomo + integer*8 :: Isomo + integer*8 :: Jdomo + integer*8 :: Jsomo + integer*8 :: diffSOMO + integer*8 :: diffDOMO + integer*8 :: xordiffSOMODOMO + integer :: ndiffSOMO + integer :: ndiffDOMO + integer :: nxordiffSOMODOMO + integer :: ndiffAll + integer :: i + integer :: j + integer :: k + integer :: hole + integer :: p + integer :: q + integer :: countalphas + logical :: pqAlreadyGenQ + logical :: pqExistsQ + logical :: ppExistsQ + Isomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,1)) + Idomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,2)) + !print*,"Input cfg" + !call debug_spindet(Isomo,1) + !call debug_spindet(Idomo,1) + + !print*,n_act_orb, "monum=",mo_num," n_core=",n_core_orb + + ! find out all pq holes possible + nholes = 0 + ! holes in SOMO + do i = 1,mo_num + if(POPCNT(IAND(Isomo,IBSET(0_8,i-1))) .EQ. 1) then + nholes += 1 + listholes(nholes) = i + holetype(nholes) = 1 + endif + end do + ! holes in DOMO + do i = 1,mo_num + if(POPCNT(IAND(Idomo,IBSET(0_8,i-1))) .EQ. 1) then + nholes += 1 + listholes(nholes) = i + holetype(nholes) = 2 + endif + end do + + ! find vmos + listvmos = -1 + vmotype = -1 + nvmos = 0 + do i = 1,mo_num + !print *,i,IBSET(0,i-1),POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))), POPCNT(IAND(Idomo,(IBSET(0_8,i-1)))) + if(POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))) .EQ. 0 .AND. POPCNT(IAND(Idomo,(IBSET(0_8,i-1)))) .EQ. 0) then + nvmos += 1 + listvmos(nvmos) = i + vmotype(nvmos) = 1 + else if(POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))) .EQ. 1 .AND. POPCNT(IAND(Idomo,(IBSET(0_8,i-1)))) .EQ. 0 ) then + nvmos += 1 + listvmos(nvmos) = i + vmotype(nvmos) = 2 + end if + end do + + !print *,"Nvmo=",nvmos + !print *,listvmos + !print *,vmotype + + allocate(tableUniqueAlphas(mo_num,mo_num)) + tableUniqueAlphas = .FALSE. + + ! Now find the allowed (p,q) excitations + Isomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,1)) + Idomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,2)) + !print *,"Isomo" + !call debug_spindet(Isomo,1) + !call debug_spindet(Idomo,1) + + !print *,"Nholes=",nholes," Nvmos=",nvmos, " idxi=",idxI + !do i = 1,nholes + ! print *,i,"->",listholes(i) + !enddo + !do i = 1,nvmos + ! print *,i,"->",listvmos(i) + !enddo + + ! TODO cfg_seniority_index + do i = 1,nholes + p = listholes(i) + do j = 1,nvmos + q = listvmos(j) + if(p .EQ. q) cycle + if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 1) then + ! SOMO -> VMO + Jsomo = IBCLR(Isomo,p-1) + Jsomo = IBSET(Jsomo,q-1) + Jdomo = Idomo + else if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 2) then + ! SOMO -> SOMO + Jsomo = IBCLR(Isomo,p-1) + Jsomo = IBCLR(Jsomo,q-1) + Jdomo = IBSET(Idomo,q-1) + else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then + ! DOMO -> VMO + Jsomo = IBSET(Isomo,p-1) + Jsomo = IBSET(Jsomo,q-1) + Jdomo = IBCLR(Idomo,p-1) + else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 2) then + ! DOMO -> SOMO + Jsomo = IBSET(Isomo,p-1) + Jsomo = IBCLR(Jsomo,q-1) + Jdomo = IBCLR(Idomo,p-1) + Jdomo = IBSET(Jdomo,q-1) + else + print*,"Something went wrong in obtain_associated_alphaI" + endif + + + pqAlreadyGenQ = .FALSE. + ! First check if it can be generated before + do k = 1, idxI-1 + diffSOMO = IEOR(Jsomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k))) + diffDOMO = IEOR(Jdomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k))) + xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) + ndiffSOMO = POPCNT(diffSOMO) + ndiffDOMO = POPCNT(diffDOMO) + nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO) + !if(POPCNT(IEOR(diffSOMO,diffDOMO)) .LE. 1 .AND. ndiffDOMO .LT. 3) then + if((ndiffSOMO+ndiffDOMO) .EQ. 0) then + pqAlreadyGenQ = .TRUE. + ppExistsQ = .TRUE. + EXIT + endif + if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then + pqAlreadyGenQ = .TRUE. + !ppExistsQ = .TRUE. + !print *,i,k,ndiffSOMO,ndiffDOMO + !call debug_spindet(Jsomo,1) + !call debug_spindet(Jdomo,1) + !call debug_spindet(iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k)),1) + !call debug_spindet(iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k)),1) + EXIT + endif + end do + + !print *,"(,",p,",",q,")",pqAlreadyGenQ + + if(pqAlreadyGenQ) cycle + + pqExistsQ = .FALSE. + ! now check if this exists in the selected list + !do k = idxI+1, N_configuration + ! diffSOMO = IEOR(OR(reunion_of_act_virt_bitmask(1,1),Jsomo),psi_configuration(1,1,k)) + ! diffDOMO = IEOR(OR(reunion_of_act_virt_bitmask(1,1),Jdomo),psi_configuration(1,2,k)) + ! ndiffSOMO = POPCNT(diffSOMO) + ! ndiffDOMO = POPCNT(diffDOMO) + ! if((ndiffSOMO + ndiffDOMO) .EQ. 0) then + ! pqExistsQ = .TRUE. + ! EXIT + ! endif + !end do + + if(.NOT. pqExistsQ) then + tableUniqueAlphas(p,q) = .TRUE. + !print *,p,q + !call debug_spindet(Jsomo,1) + !call debug_spindet(Jdomo,1) + endif + end do + end do + + !print *,tableUniqueAlphas(:,:) + + ! prune list of alphas + Isomo = Icfg(1,1) + Idomo = Icfg(1,2) + Jsomo = Icfg(1,1) + Jdomo = Icfg(1,2) + NalphaIcfg = 0 + do i = 1, nholes + p = listholes(i) + do j = 1, nvmos + q = listvmos(j) + if(p .EQ. q) cycle + if(tableUniqueAlphas(p,q)) then + if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 1) then + ! SOMO -> VMO + Jsomo = IBCLR(Isomo,p-1) + Jsomo = IBSET(Jsomo,q-1) + Jdomo = Idomo + else if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 2) then + ! SOMO -> SOMO + Jsomo = IBCLR(Isomo,p-1) + Jsomo = IBCLR(Jsomo,q-1) + Jdomo = IBSET(Idomo,q-1) + else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then + ! DOMO -> VMO + Jsomo = IBSET(Isomo,p-1) + Jsomo = IBSET(Jsomo,q-1) + Jdomo = IBCLR(Idomo,p-1) + else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 2) then + ! DOMO -> SOMO + Jsomo = IBSET(Isomo,p-1) + Jsomo = IBCLR(Jsomo,q-1) + Jdomo = IBCLR(Idomo,p-1) + Jdomo = IBSET(Jdomo,q-1) + else + print*,"Something went wrong in obtain_associated_alphaI" + endif + + ! SOMO + NalphaIcfg += 1 + !print *,p,q,"|",holetype(i),vmotype(j),NalphaIcfg + alphasIcfg(1,1,NalphaIcfg) = Jsomo + alphasIcfg(1,2,NalphaIcfg) = IOR(Jdomo,ISHFT(1_8,n_core_orb)-1) + endif + end do + end do + + ! Check if this Icfg has been previously generated as a mono + ppExistsQ = .False. + Isomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,1)) + Idomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,2)) + do k = 1, idxI-1 + diffSOMO = IEOR(Isomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k))) + diffDOMO = IEOR(Idomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k))) + xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) + ndiffSOMO = POPCNT(diffSOMO) + ndiffDOMO = POPCNT(diffDOMO) + nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO) + if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then + ppExistsQ = .TRUE. + EXIT + endif + end do + ! Diagonal part (pp,qq) + if(nholes > 0 .AND. (.NOT. ppExistsQ))then + ! SOMO + NalphaIcfg += 1 + !print *,p,q,"|",holetype(i),vmotype(j),NalphaIcfg + !call debug_spindet(Idomo,1) + !call debug_spindet(Jdomo,1) + alphasIcfg(1,1,NalphaIcfg) = Icfg(1,1) + alphasIcfg(1,2,NalphaIcfg) = Icfg(1,2) + endif + + end subroutine + function getNSOMO(Icfg) result(NSOMO) implicit none integer(bit_kind),intent(in) :: Icfg(N_int,2) @@ -47,6 +317,7 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod qmodel = q if(p .EQ. q) then + !print *,"input pq=",p,q,"extype=",extype pmodel = 1 qmodel = 1 else @@ -80,26 +351,44 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod ! SOMO -> VMO !print *,"type -> SOMO -> VMO" !Isomo = IEOR(Isomo,Jsomo) - mask = ISHFT(1_8,p) - 1 - Isomo = IAND(Isomo,mask) - pmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) - mask = ISHFT(1_8,q) - 1 - Jsomo = IAND(Jsomo,mask) - qmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) + if(p.LT.q) then + mask = ISHFT(1_8,p) - 1 + Isomo = IAND(Isomo,mask) + pmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) + mask = ISHFT(1_8,q) - 1 + Jsomo = IAND(Jsomo,mask) + qmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) + 1 + else + mask = ISHFT(1_8,p) - 1 + Isomo = IAND(Isomo,mask) + pmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) + 1 + mask = ISHFT(1_8,q) - 1 + Jsomo = IAND(Jsomo,mask) + qmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) + endif case (4) ! DOMO -> SOMO ! remove all domos except one at p !print *,"type -> DOMO -> SOMO" !Isomo = IEOR(Isomo,Jsomo) - mask = ISHFT(1_8,p) - 1 - Jsomo = IAND(Jsomo,mask) - pmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) - mask = ISHFT(1_8,q) - 1 - Isomo = IAND(Isomo,mask) - qmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) + if(p.LT.q) then + mask = ISHFT(1_8,p) - 1 + Jsomo = IAND(Jsomo,mask) + pmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) + mask = ISHFT(1_8,q) - 1 + Isomo = IAND(Isomo,mask) + qmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) + 1 + else + mask = ISHFT(1_8,p) - 1 + Jsomo = IAND(Jsomo,mask) + pmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) + 1 + mask = ISHFT(1_8,q) - 1 + Isomo = IAND(Isomo,mask) + qmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) + endif case default print *,"something is wrong in convertOrbIdsToModelSpaceIds" end select endif - !print *,p,q,"model ids=",pmodel,qmodel + !print *,p,q,"model ids=",pmodel,qmodel end subroutine convertOrbIdsToModelSpaceIds diff --git a/src/csf/sigma_vector.irp.f b/src/csf/sigma_vector.irp.f index 9c1ac56a..e14f11bb 100644 --- a/src/csf/sigma_vector.irp.f +++ b/src/csf/sigma_vector.irp.f @@ -65,11 +65,6 @@ subroutine get_phase_qp_to_cfg(Ialpha, Ibeta, phaseout) integer :: nbetas integer :: count, k - ! Remove the DOMOs - mask2 = IAND(Ialpha,Ibeta) - deta = IEOR(Ialpha,mask2) - detb = IEOR(Ibeta ,mask2) - ! Find how many alpha electrons there are in all the N_ints integer :: Na(N_int) do k=1,N_int From e166d0646aaf38e122bc77c08a68a3ea80d2dac8 Mon Sep 17 00:00:00 2001 From: v1j4y Date: Mon, 1 Mar 2021 12:50:50 +0100 Subject: [PATCH 4/7] UPdate configuration_sigma_vector.irp.f to sigma_vector.irp.f --- src/csf/sigma_vector.irp.f | 503 ++++++++++++++++++++++++++++--------- src/csf/tree_utils.h | 3 +- 2 files changed, 380 insertions(+), 126 deletions(-) diff --git a/src/csf/sigma_vector.irp.f b/src/csf/sigma_vector.irp.f index e14f11bb..d981da79 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, dimBasisCSF] + &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 + dimBasisCSF = 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)))) + dimBasisCSF += ncfg * dimcsfpercfg + !print *,i,">(",ncfg,ncfgprev,ncfgpersomo,")",",",detDimperBF,">",dimcsfpercfg, " | dimbas= ", dimBasisCSF + !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)))) + dimBasisCSF += ncfg * dimcsfpercfg + !print *,i,">(",ncfg,ncfgprev,ncfgpersomo,")",",",detDimperBF,">",dimcsfpercfg, " | dimbas= ", dimBasisCSF + endif + END_PROVIDER subroutine get_phase_qp_to_cfg(Ialpha, Ibeta, phaseout) use bitmasks @@ -99,7 +108,296 @@ 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, (dimBasisCSF,1)] + &BEGIN_PROVIDER [ integer, psi_config_data, (N_configuration,2)] + &BEGIN_PROVIDER [ integer, psi_csf_to_config_data, (dimBasisCSF)] + 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, mask, Ialpha,Ibeta + 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 *,"dimBasisCSF=",dimBasisCSF + 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,1,psi_configuration_to_psi_det_data(j)) + Ibeta = psi_det(1,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 + + subroutine convertWFfromDETtoCSF(psi_coef_det_in, psi_coef_cfg_out) + 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, mask, Ialpha,Ibeta + integer :: rows, cols, i, j, k + integer :: startdet, enddet + integer*8 MS + integer ndetI + integer :: getNSOMO + real*8,intent(in) :: psi_coef_det_in(n_det,1) + real*8,intent(out) :: psi_coef_cfg_out(dimBasisCSF,1) + 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 *,"dimBasisCSF=",dimBasisCSF + print *,"N_configurations=",N_configuration + print *,"n_core_orb=",n_core_orb + ! initialization + psi_coef_cfg_out(:,1) = 0.d0 + + 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,1,psi_configuration_to_psi_det_data(j)) + Ibeta = psi_det(1,2,psi_configuration_to_psi_det_data(j)) + 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_cfg_out(countcsf+1,1), size(psi_coef_cfg_out,1)) + + deallocate(tempCoeff) + deallocate(tempBuffer) + psi_config_data(i,1) = countcsf + 1 + countcsf += bfIcfg + psi_config_data(i,2) = countcsf + enddo + print *,"Norm det=",norm_det1, size(psi_coef_cfg_out,1), " Dim csf=", countcsf + + end subroutine convertWFfromDETtoCSF + + subroutine convertWFfromCSFtoDET(psi_coef_cfg_in, psi_coef_det) + implicit none + BEGIN_DOC + ! Documentation for convertCSFtoDET + ! This function converts the wavefunction + ! in CFG basis to DET basis using the + ! transformation matrix provided before. + END_DOC + real*8,intent(in) :: psi_coef_cfg_in(dimBasisCSF,1) + real*8,intent(out) :: psi_coef_det(N_det,1) + real*8 :: tmp_psi_coef_det(maxDetDimPerBF) + integer s, bfIcfg + integer countcsf + integer countdet + integer*8 :: Isomo, Idomo, Ialpha, Ibeta + 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 :: phasedet + ! number of states + integer istate + istate = 1 + countcsf = 1 + countdet = 1 + print *,"in function convertWFfromCSFtoDET()" + + + 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)) + + do j = 1,bfIcfg + tempCoeff(j) = psi_coef_cfg_in(countcsf,1) + countcsf += 1 + enddo + !print *,"dimcoef=",bfIcfg + !call printMatrix(tempCoeff,bfIcfg,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('T','N', ndetI, 1, bfIcfg, 1.d0, tempBuffer, size(tempBuffer,1), tempCoeff, size(tempCoeff,1), 0.d0, psi_coef_det(countdet,1), size(psi_coef_det,1)) + call dgemm('T','N', ndetI, 1, bfIcfg, 1.d0, tempBuffer, size(tempBuffer,1), tempCoeff, size(tempCoeff,1), 0.d0, tmp_psi_coef_det, size(tmp_psi_coef_det,1)) + + !call dgemv('N', NBFMax, maxDetDimPerBF, 1.d0, tempBuffer, size(tempBuffer,1), tempCoeff, 1, 0.d0, psi_coef_config(countcsf,1), 1) + + !print *,"result" + !call printMatrix(tmp_psi_coef_det,ndetI,1) + + countdet = 1 + do j=startdet,enddet + Ialpha = psi_det(1,1,psi_configuration_to_psi_det_data(j)) + Ibeta = psi_det(1,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 + psi_coef_det(psi_configuration_to_psi_det_data(j),1) = tmp_psi_coef_det(countdet)*phasedet + countdet += 1 + enddo + + deallocate(tempCoeff) + deallocate(tempBuffer) + !countdet += ndetI + enddo + + !countdet = 1 + !tmp_psi_coef_det = psi_coef_det(:,1) + !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 + ! print *,i,">>>",startdet,enddet + ! do k=1,ndetI + ! !psi_coef_det(startdet+k-1,1) = tmp_psi_coef_det(countdet) + ! psi_coef_det(countdet,1) = tmp_psi_coef_det(startdet+k-1) + ! countdet += 1 + ! enddo + !enddo + + print *,"End ncsfs=",countcsf + + end subroutine convertCSFtoDET + + BEGIN_PROVIDER [ integer, AIJpqMatrixDimsList, (0:NSOMOMax,0:NSOMOMax,4,NSOMOMax+1,NSOMOMax+1,2)] &BEGIN_PROVIDER [ integer, rowsmax] &BEGIN_PROVIDER [ integer, colsmax] use cfunctions @@ -122,9 +420,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 @@ -152,6 +454,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 @@ -167,6 +470,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 @@ -200,6 +506,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 @@ -216,6 +523,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 @@ -223,8 +532,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) @@ -239,6 +548,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 @@ -253,18 +563,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 @@ -276,6 +588,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 @@ -289,9 +602,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 @@ -327,13 +641,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 @@ -350,6 +669,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, & @@ -371,6 +694,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 @@ -384,6 +709,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 @@ -409,6 +736,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, & @@ -430,6 +761,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 @@ -443,13 +776,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) @@ -460,6 +794,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, & @@ -481,6 +819,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 @@ -494,18 +834,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 @@ -533,6 +875,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 @@ -545,94 +889,3 @@ end subroutine get_phase_qp_to_cfg end do 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 - 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(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 :: 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)))) - - allocate(tempBuffer(bfIcfg,ndetI)) - call getCSFtoDETTransformationMatrix(Isomo, MS, NBFMax, maxDetDimPerBF, tempBuffer) - DetToCSFTransformationMatrix(i,1:bfIcfg,1:ndetI) = tempBuffer(1:bfIcfg,1:ndetI) - deallocate(tempBuffer) - 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 - - 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 - - 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, 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) - - deallocate(tempCoeff) - deallocate(tempBuffer) - psi_config_data(i,1) = countcsf + 1 - countcsf += bfIcfg - psi_config_data(i,2) = countcsf - enddo - - END_PROVIDER diff --git a/src/csf/tree_utils.h b/src/csf/tree_utils.h index eae19c1b..049b09db 100644 --- a/src/csf/tree_utils.h +++ b/src/csf/tree_utils.h @@ -22,7 +22,8 @@ struct bin_tree { int NBF; }; -#include "/export/apps/pgi/linux86-64/18.10/include/cblas.h" +//#include "/export/apps/pgi/linux86-64/18.10/include/cblas.h" +#include "cblas.h" #define MAX_SOMO 32 From 67f59ba52a2bc0787ea3a5a6c55d5e134d0c7c78 Mon Sep 17 00:00:00 2001 From: v1j4y Date: Mon, 1 Mar 2021 13:41:10 +0100 Subject: [PATCH 5/7] Updated conversion. --- src/csf/conversion.irp.f | 51 ++++----- src/csf/sigma_vector.irp.f | 213 +++---------------------------------- 2 files changed, 36 insertions(+), 228 deletions(-) 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 d981da79..b0fdde6c 100644 --- a/src/csf/sigma_vector.irp.f +++ b/src/csf/sigma_vector.irp.f @@ -2,7 +2,7 @@ &BEGIN_PROVIDER [ integer, NCSFMax] &BEGIN_PROVIDER [ integer*8, NMO] &BEGIN_PROVIDER [ integer, NBFMax] - &BEGIN_PROVIDER [ integer, dimBasisCSF] + &BEGIN_PROVIDER [ integer, n_CSF] &BEGIN_PROVIDER [ integer, maxDetDimPerBF] implicit none BEGIN_DOC @@ -29,7 +29,7 @@ MS = elec_alpha_num-elec_beta_num !print *,"NSOMOMax=",NSOMOMax, cfg_seniority_index(0) ! number of cfgs = number of dets for 0 somos - dimBasisCSF = cfg_seniority_index(0)-1 + 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 @@ -40,8 +40,8 @@ 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)))) - dimBasisCSF += ncfg * dimcsfpercfg - !print *,i,">(",ncfg,ncfgprev,ncfgpersomo,")",",",detDimperBF,">",dimcsfpercfg, " | dimbas= ", dimBasisCSF + 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) @@ -50,8 +50,8 @@ 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)))) - dimBasisCSF += ncfg * dimcsfpercfg - !print *,i,">(",ncfg,ncfgprev,ncfgpersomo,")",",",detDimperBF,">",dimcsfpercfg, " | dimbas= ", dimBasisCSF + n_CSF += ncfg * dimcsfpercfg + !print *,i,">(",ncfg,ncfgprev,ncfgpersomo,")",",",detDimperBF,">",dimcsfpercfg, " | dimbas= ", n_CSF endif END_PROVIDER @@ -70,7 +70,7 @@ 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 @@ -109,9 +109,9 @@ end subroutine get_phase_qp_to_cfg BEGIN_PROVIDER [ real*8, DetToCSFTransformationMatrix, (0:NSOMOMax,NBFMax,maxDetDimPerBF)] - &BEGIN_PROVIDER [ real*8, psi_coef_config, (dimBasisCSF,1)] + &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, (dimBasisCSF)] + &BEGIN_PROVIDER [ integer, psi_csf_to_config_data, (n_CSF)] use cfunctions implicit none BEGIN_DOC @@ -119,7 +119,8 @@ end subroutine get_phase_qp_to_cfg ! Provides the matrix of transformatons for the ! conversion between determinant to CSF basis (in BFs) END_DOC - integer*8 :: Isomo, Idomo, mask, Ialpha,Ibeta + 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 @@ -132,7 +133,7 @@ end subroutine get_phase_qp_to_cfg MS = elec_alpha_num - elec_beta_num print *,"Maxbfdim=",NBFMax print *,"Maxdetdim=",maxDetDimPerBF - print *,"dimBasisCSF=",dimBasisCSF + print *,"n_CSF=",n_CSF print *,"N_configurations=",N_configuration print *,"n_core_orb=",n_core_orb ! initialization @@ -167,8 +168,8 @@ end subroutine get_phase_qp_to_cfg allocate(tempCoeff(ndetI)) countdet = 1 do j = startdet, enddet - Ialpha = psi_det(1,1,psi_configuration_to_psi_det_data(j)) - Ibeta = psi_det(1,2,psi_configuration_to_psi_det_data(j)) + 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) @@ -211,192 +212,6 @@ end subroutine get_phase_qp_to_cfg END_PROVIDER - subroutine convertWFfromDETtoCSF(psi_coef_det_in, psi_coef_cfg_out) - 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, mask, Ialpha,Ibeta - integer :: rows, cols, i, j, k - integer :: startdet, enddet - integer*8 MS - integer ndetI - integer :: getNSOMO - real*8,intent(in) :: psi_coef_det_in(n_det,1) - real*8,intent(out) :: psi_coef_cfg_out(dimBasisCSF,1) - 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 *,"dimBasisCSF=",dimBasisCSF - print *,"N_configurations=",N_configuration - print *,"n_core_orb=",n_core_orb - ! initialization - psi_coef_cfg_out(:,1) = 0.d0 - - 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,1,psi_configuration_to_psi_det_data(j)) - Ibeta = psi_det(1,2,psi_configuration_to_psi_det_data(j)) - 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_cfg_out(countcsf+1,1), size(psi_coef_cfg_out,1)) - - deallocate(tempCoeff) - deallocate(tempBuffer) - psi_config_data(i,1) = countcsf + 1 - countcsf += bfIcfg - psi_config_data(i,2) = countcsf - enddo - print *,"Norm det=",norm_det1, size(psi_coef_cfg_out,1), " Dim csf=", countcsf - - end subroutine convertWFfromDETtoCSF - - subroutine convertWFfromCSFtoDET(psi_coef_cfg_in, psi_coef_det) - implicit none - BEGIN_DOC - ! Documentation for convertCSFtoDET - ! This function converts the wavefunction - ! in CFG basis to DET basis using the - ! transformation matrix provided before. - END_DOC - real*8,intent(in) :: psi_coef_cfg_in(dimBasisCSF,1) - real*8,intent(out) :: psi_coef_det(N_det,1) - real*8 :: tmp_psi_coef_det(maxDetDimPerBF) - integer s, bfIcfg - integer countcsf - integer countdet - integer*8 :: Isomo, Idomo, Ialpha, Ibeta - 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 :: phasedet - ! number of states - integer istate - istate = 1 - countcsf = 1 - countdet = 1 - print *,"in function convertWFfromCSFtoDET()" - - - 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)) - - do j = 1,bfIcfg - tempCoeff(j) = psi_coef_cfg_in(countcsf,1) - countcsf += 1 - enddo - !print *,"dimcoef=",bfIcfg - !call printMatrix(tempCoeff,bfIcfg,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('T','N', ndetI, 1, bfIcfg, 1.d0, tempBuffer, size(tempBuffer,1), tempCoeff, size(tempCoeff,1), 0.d0, psi_coef_det(countdet,1), size(psi_coef_det,1)) - call dgemm('T','N', ndetI, 1, bfIcfg, 1.d0, tempBuffer, size(tempBuffer,1), tempCoeff, size(tempCoeff,1), 0.d0, tmp_psi_coef_det, size(tmp_psi_coef_det,1)) - - !call dgemv('N', NBFMax, maxDetDimPerBF, 1.d0, tempBuffer, size(tempBuffer,1), tempCoeff, 1, 0.d0, psi_coef_config(countcsf,1), 1) - - !print *,"result" - !call printMatrix(tmp_psi_coef_det,ndetI,1) - - countdet = 1 - do j=startdet,enddet - Ialpha = psi_det(1,1,psi_configuration_to_psi_det_data(j)) - Ibeta = psi_det(1,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 - psi_coef_det(psi_configuration_to_psi_det_data(j),1) = tmp_psi_coef_det(countdet)*phasedet - countdet += 1 - enddo - - deallocate(tempCoeff) - deallocate(tempBuffer) - !countdet += ndetI - enddo - - !countdet = 1 - !tmp_psi_coef_det = psi_coef_det(:,1) - !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 - ! print *,i,">>>",startdet,enddet - ! do k=1,ndetI - ! !psi_coef_det(startdet+k-1,1) = tmp_psi_coef_det(countdet) - ! psi_coef_det(countdet,1) = tmp_psi_coef_det(startdet+k-1) - ! countdet += 1 - ! enddo - !enddo - - print *,"End ncsfs=",countcsf - - end subroutine convertCSFtoDET - BEGIN_PROVIDER [ integer, AIJpqMatrixDimsList, (0:NSOMOMax,0:NSOMOMax,4,NSOMOMax+1,NSOMOMax+1,2)] &BEGIN_PROVIDER [ integer, rowsmax] &BEGIN_PROVIDER [ integer, colsmax] From c1ef95cd6053b761782690d8698a6ffb56d182aa Mon Sep 17 00:00:00 2001 From: v1j4y Date: Mon, 1 Mar 2021 14:18:13 +0100 Subject: [PATCH 6/7] Fixed initialization in get_phase. --- src/csf/sigma_vector.irp.f | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/csf/sigma_vector.irp.f b/src/csf/sigma_vector.irp.f index b0fdde6c..aa929a76 100644 --- a/src/csf/sigma_vector.irp.f +++ b/src/csf/sigma_vector.irp.f @@ -74,6 +74,10 @@ subroutine get_phase_qp_to_cfg(Ialpha, Ibeta, phaseout) integer :: nbetas integer :: count, k + ! Initliaze deta and detb + deta = Ialpha + detb = Ibeta + ! Find how many alpha electrons there are in all the N_ints integer :: Na(N_int) do k=1,N_int @@ -86,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)) ) From c38d25fdba61e6c773a7ba8b059edb6bb68cb4c7 Mon Sep 17 00:00:00 2001 From: v1j4y Date: Mon, 1 Mar 2021 18:28:45 +0100 Subject: [PATCH 7/7] Working on sigma_vector for CFG-CI. --- src/csf/configuration_CI_sigma_helpers.irp.f | 42 +- src/csf/sigma_vector.irp.f | 464 ++++++++++++++++++ .../diagonalization_hcsf_dressed.irp.f | 29 +- 3 files changed, 508 insertions(+), 27 deletions(-) diff --git a/src/csf/configuration_CI_sigma_helpers.irp.f b/src/csf/configuration_CI_sigma_helpers.irp.f index eb532b03..0290a5f7 100644 --- a/src/csf/configuration_CI_sigma_helpers.irp.f +++ b/src/csf/configuration_CI_sigma_helpers.irp.f @@ -1,4 +1,4 @@ - subroutine obtain_associated_alphaI(idxI, Icfg, alphasIcfg, NalphaIcfg, factor_alphaI) + subroutine obtain_associated_alphaI(idxI, Icfg, alphasIcfg, NalphaIcfg) implicit none use bitmasks BEGIN_DOC @@ -10,7 +10,6 @@ integer,intent(in) :: idxI ! The id of the Ith CFG integer(bit_kind),intent(in) :: Icfg(N_int,2) integer,intent(out) :: NalphaIcfg - real*8 ,intent(out) :: factor_alphaI(*) integer(bit_kind),intent(out) :: alphasIcfg(N_int,2,*) logical,dimension(:,:),allocatable :: tableUniqueAlphas integer :: listholes(mo_num) @@ -293,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) @@ -317,15 +323,9 @@ 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 - !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 diff --git a/src/csf/sigma_vector.irp.f b/src/csf/sigma_vector.irp.f index aa929a76..1a434aea 100644 --- a/src/csf/sigma_vector.irp.f +++ b/src/csf/sigma_vector.irp.f @@ -708,3 +708,467 @@ end subroutine get_phase_qp_to_cfg end do end do END_PROVIDER + +subroutine calculate_preconditioner_cfg(diag_energies) + implicit none + use bitmasks + BEGIN_DOC + ! Documentation for calculate_preconditioner + ! + ! Calculates the diagonal energies of + ! the configurations in psi_configuration + ! returns : diag_energies : + END_DOC + 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 + 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) + + 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 + + !!! Double Excitations !!! + + ! Loop over all selected configurations + do i = istart_cfg,iend_cfg + + 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) + + ! 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) + + 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 + + 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 subroutine calculate_sigma_vector_cfg_nst 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 = = ! -------------------------------------------