diff --git a/external/qp2-dependencies b/external/qp2-dependencies index bc856147..90ee61f5 160000 --- a/external/qp2-dependencies +++ b/external/qp2-dependencies @@ -1 +1 @@ -Subproject commit bc856147f6e626a6616b20344e5b8e3f30f44a92 +Subproject commit 90ee61f5041c7c94a0c605625a264860292813a0 diff --git a/src/csf/configuration_CI_sigma_helpers.irp.f b/src/csf/configuration_CI_sigma_helpers.irp.f index f73362eb..a4aa41ff 100644 --- a/src/csf/configuration_CI_sigma_helpers.irp.f +++ b/src/csf/configuration_CI_sigma_helpers.irp.f @@ -1,3 +1,566 @@ +use bitmasks + + BEGIN_PROVIDER [ integer(bit_kind), alphasIcfg_list , (N_int,2,N_configuration,mo_num*(mo_num))] +&BEGIN_PROVIDER [ integer, NalphaIcfg_list, (N_configuration) ] + implicit none + !use bitmasks + BEGIN_DOC + ! Documentation for alphasI + ! Returns the associated alpha's for + ! the input configuration Icfg. + END_DOC + + integer :: idxI ! The id of the Ith CFG + integer(bit_kind) :: Icfg(N_int,2) + integer :: NalphaIcfg + 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 :: kstart + integer :: kend + integer :: Nsomo_I + integer :: hole + integer :: p + integer :: q + integer :: countalphas + logical :: pqAlreadyGenQ + logical :: pqExistsQ + logical :: ppExistsQ + + double precision :: t0, t1 + call wall_time(t0) + + allocate(tableUniqueAlphas(mo_num,mo_num)) + NalphaIcfg_list = 0 + + do idxI = 1, N_configuration + + Icfg = psi_configuration(:,:,idxI) + + Isomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,1)) + Idomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,2)) + + ! 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 + if(IAND(Idomo,(IBSET(0_8,i-1))) .EQ. 0) then + if(IAND(Isomo,(IBSET(0_8,i-1))) .EQ. 0) then + nvmos += 1 + listvmos(nvmos) = i + vmotype(nvmos) = 1 + else if(POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))) .EQ. 1) then + nvmos += 1 + listvmos(nvmos) = i + vmotype(nvmos) = 2 + end if + end if + end do + + 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)) + Nsomo_I = POPCNT(Isomo) + if(Nsomo_I .EQ. 0) then + kstart = 1 + else + kstart = cfg_seniority_index(max(0,Nsomo_I-2)) + endif + kend = idxI-1 + + 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 + kstart = max(1,cfg_seniority_index(max(0,Nsomo_I-2))) + kend = idxI-1 + else if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 2) then + ! SOMO -> SOMO + Jsomo = IBCLR(Isomo,p-1) + Jsomo = IBCLR(Jsomo,q-1) + Jdomo = IBSET(Idomo,q-1) + kstart = max(1,cfg_seniority_index(max(0,Nsomo_I-4))) + kend = idxI-1 + else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then + ! DOMO -> VMO + Jsomo = IBSET(Isomo,p-1) + Jsomo = IBSET(Jsomo,q-1) + Jdomo = IBCLR(Idomo,p-1) + kstart = cfg_seniority_index(Nsomo_I) + kend = idxI-1 + else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 2) then + ! DOMO -> SOMO + Jsomo = IBSET(Isomo,p-1) + Jsomo = IBCLR(Jsomo,q-1) + Jdomo = IBCLR(Idomo,p-1) + Jdomo = IBSET(Jdomo,q-1) + kstart = max(1,cfg_seniority_index(max(0,Nsomo_I-2))) + kend = idxI-1 + else + print*,"Something went wrong in obtain_associated_alphaI" + endif + + ! Again, we don't have to search from 1 + ! we just use seniority to find the + ! first index with NSOMO - 2 to NSOMO + 2 + ! this is what is done in kstart, kend + + pqAlreadyGenQ = .FALSE. + ! First check if it can be generated before + do k = kstart, kend + diffSOMO = IEOR(Jsomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k))) + ndiffSOMO = POPCNT(diffSOMO) + if((ndiffSOMO .NE. 0) .AND. (ndiffSOMO .NE. 2)) cycle + diffDOMO = IEOR(Jdomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k))) + xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) + ndiffDOMO = POPCNT(diffDOMO) + nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO) + nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO + !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((nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then + pqAlreadyGenQ = .TRUE. + EXIT + endif + end do + + if(pqAlreadyGenQ) cycle + + pqExistsQ = .FALSE. + + if(.NOT. pqExistsQ) then + tableUniqueAlphas(p,q) = .TRUE. + endif + end do + end do + + !print *,tableUniqueAlphas(:,:) + + ! prune list of alphas + Isomo = Icfg(1,1) + Idomo = Icfg(1,2) + Jsomo = Icfg(1,1) + Jdomo = Icfg(1,2) + NalphaIcfg = 0 + do i = 1, nholes + p = listholes(i) + 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 *,i,j,"|",NalphaIcfg + alphasIcfg_list(1,1,idxI,NalphaIcfg) = Jsomo + alphasIcfg_list(1,2,idxI,NalphaIcfg) = IOR(Jdomo,ISHFT(1_8,n_core_orb)-1) + NalphaIcfg_list(idxI) = NalphaIcfg + endif + end do + end do + + ! Check if this Icfg has been previously generated as a mono + ppExistsQ = .False. + Isomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,1)) + Idomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,2)) + kstart = max(1,cfg_seniority_index(max(0,Nsomo_I-2))) + do k = kstart, idxI-1 + diffSOMO = IEOR(Isomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k))) + ndiffSOMO = POPCNT(diffSOMO) + if (ndiffSOMO /= 2) cycle + diffDOMO = IEOR(Idomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k))) + xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) + ndiffDOMO = POPCNT(diffDOMO) + nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO) + if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4)) then + ppExistsQ = .TRUE. + EXIT + endif + end do + ! Diagonal part (pp,qq) + if(nholes > 0 .AND. (.NOT. ppExistsQ))then + ! SOMO + NalphaIcfg += 1 + alphasIcfg_list(1,1,idxI,NalphaIcfg) = Icfg(1,1) + alphasIcfg_list(1,2,idxI,NalphaIcfg) = Icfg(1,2) + NalphaIcfg_list(idxI) = NalphaIcfg + endif + + NalphaIcfg = 0 + enddo ! end loop idxI + call wall_time(t1) + print *, 'Preparation : ', t1 - t0 + +END_PROVIDER + + subroutine obtain_associated_alphaI(idxI, Icfg, alphasIcfg, NalphaIcfg) + implicit none + use bitmasks + BEGIN_DOC + ! Documentation for alphasI + ! Returns the associated alpha's for + ! the input configuration Icfg. + END_DOC + + integer,intent(in) :: idxI ! The id of the Ith CFG + integer(bit_kind),intent(in) :: Icfg(N_int,2) + integer,intent(out) :: NalphaIcfg + integer(bit_kind),intent(out) :: alphasIcfg(N_int,2,*) + logical,dimension(:,:),allocatable :: tableUniqueAlphas + integer :: listholes(mo_num) + integer :: holetype(mo_num) ! 1-> SOMO 2->DOMO + integer :: nholes + integer :: nvmos + integer :: listvmos(mo_num) + integer :: vmotype(mo_num) ! 1 -> VMO 2 -> SOMO + integer*8 :: Idomo + integer*8 :: Isomo + integer*8 :: Jdomo + integer*8 :: Jsomo + integer*8 :: diffSOMO + integer*8 :: diffDOMO + integer*8 :: xordiffSOMODOMO + integer :: ndiffSOMO + integer :: ndiffDOMO + integer :: nxordiffSOMODOMO + integer :: ndiffAll + integer :: i + integer :: j + integer :: k + integer :: kstart + integer :: kend + integer :: Nsomo_I + 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)) + Nsomo_I = POPCNT(Isomo) + if(Nsomo_I .EQ. 0) then + kstart = 1 + else + kstart = cfg_seniority_index(max(0,Nsomo_I-2)) + endif + kend = idxI-1 + !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 + + 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 + kstart = max(1,cfg_seniority_index(max(0,Nsomo_I-2))) + kend = idxI-1 + else if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 2) then + ! SOMO -> SOMO + Jsomo = IBCLR(Isomo,p-1) + Jsomo = IBCLR(Jsomo,q-1) + Jdomo = IBSET(Idomo,q-1) + kstart = max(1,cfg_seniority_index(max(0,Nsomo_I-4))) + kend = idxI-1 + else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then + ! DOMO -> VMO + Jsomo = IBSET(Isomo,p-1) + Jsomo = IBSET(Jsomo,q-1) + Jdomo = IBCLR(Idomo,p-1) + kstart = cfg_seniority_index(Nsomo_I) + kend = idxI-1 + else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 2) then + ! DOMO -> SOMO + Jsomo = IBSET(Isomo,p-1) + Jsomo = IBCLR(Jsomo,q-1) + Jdomo = IBCLR(Idomo,p-1) + Jdomo = IBSET(Jdomo,q-1) + kstart = max(1,cfg_seniority_index(max(0,Nsomo_I-2))) + kend = idxI-1 + else + print*,"Something went wrong in obtain_associated_alphaI" + endif + + ! Again, we don't have to search from 1 + ! we just use seniortiy to find the + ! first index with NSOMO - 2 to NSOMO + 2 + ! this is what is done in kstart, kend + + pqAlreadyGenQ = .FALSE. + ! First check if it can be generated before + do k = kstart, kend + diffSOMO = IEOR(Jsomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k))) + ndiffSOMO = POPCNT(diffSOMO) + if((ndiffSOMO .NE. 0) .AND. (ndiffSOMO .NE. 2)) cycle + diffDOMO = IEOR(Jdomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k))) + xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) + ndiffDOMO = POPCNT(diffDOMO) + nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO) + nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO + !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((nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then + pqAlreadyGenQ = .TRUE. + !EXIT + !ppExistsQ = .TRUE. + !print *,i,k,ndiffSOMO,ndiffDOMO + !call debug_spindet(Jsomo,1) + !call debug_spindet(Jdomo,1) + !call debug_spindet(iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k)),1) + !call debug_spindet(iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k)),1) + EXIT + endif + end do + + !print *,"(,",p,",",q,")",pqAlreadyGenQ + + if(pqAlreadyGenQ) cycle + + 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 *,i,j,"|",NalphaIcfg + alphasIcfg(1,1,NalphaIcfg) = Jsomo + alphasIcfg(1,2,NalphaIcfg) = IOR(Jdomo,ISHFT(1_8,n_core_orb)-1) + !print *,"I = ",idxI, " Na=",NalphaIcfg," - ",Jsomo, 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) @@ -8,98 +571,3 @@ NSOMO += POPCNT(Icfg(i,1)) enddo end function getNSOMO - -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 - 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) - 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)) - 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)) - case default - print *,"something is wrong in convertOrbIdsToModelSpaceIds" - end select - endif - !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 4d409f50..680cdc39 100644 --- a/src/csf/sigma_vector.irp.f +++ b/src/csf/sigma_vector.irp.f @@ -8,6 +8,7 @@ end function logabsgamma BEGIN_PROVIDER [ integer, NSOMOMax] + &BEGIN_PROVIDER [ integer, NSOMOMin] &BEGIN_PROVIDER [ integer, NCSFMax] &BEGIN_PROVIDER [ integer*8, NMO] &BEGIN_PROVIDER [ integer, NBFMax] @@ -20,6 +21,7 @@ ! required for the calculation of prototype arrays. END_DOC NSOMOMax = min(elec_num, cfg_nsomo_max + 2) + NSOMOMin = 0 ! Note that here we need NSOMOMax + 2 sizes NCSFMax = max(1,nint((binom(NSOMOMax,(NSOMOMax+1)/2)-binom(NSOMOMax,((NSOMOMax+1)/2)+1)))) ! TODO: NCSFs for MS=0 NBFMax = NCSFMax @@ -101,9 +103,9 @@ subroutine get_phase_qp_to_cfg(Ialpha, Ibeta, phaseout) real*8,intent(out) :: phaseout integer(bit_kind) :: mask, deta(N_int), detb(N_int) integer :: nbetas - integer :: k + integer :: count, k - ! Initialize deta and detb + ! Initliaze deta and detb deta = Ialpha detb = Ibeta @@ -141,7 +143,7 @@ end subroutine get_phase_qp_to_cfg - BEGIN_PROVIDER [ integer, AIJpqMatrixDimsList, (0:NSOMOMax,0:NSOMOMax,4,NSOMOMax,NSOMOMax,2)] + BEGIN_PROVIDER [ integer, AIJpqMatrixDimsList, (NSOMOMin:NSOMOMax,4,NSOMOMax+1,NSOMOMax+1,2)] &BEGIN_PROVIDER [ integer, rowsmax] &BEGIN_PROVIDER [ integer, colsmax] use cfunctions @@ -160,14 +162,18 @@ end subroutine get_phase_qp_to_cfg cols = -1 integer*8 MS MS = elec_alpha_num-elec_beta_num - integer nsomomin nsomomin = elec_alpha_num-elec_beta_num rowsmax = 0 colsmax = 0 + print *,"NSOMOMax = ",NSOMOMax + print *,"NSOMOMin = ",NSOMOMin !allocate(AIJpqMatrixDimsList(NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,2)) ! Type ! 1. SOMO -> SOMO - do i = 2-iand(nsomomin,1), NSOMOMax, 2 + !print *,"Doing SOMO->SOMO" + AIJpqMatrixDimsList(NSOMOMin,1,1,1,1) = 1 + AIJpqMatrixDimsList(NSOMOMin,1,1,1,2) = 1 + do i = NSOMOMin, NSOMOMax, 2 Isomo = ISHFT(1_8,i)-1 do j = i-2,i-2, 2 Jsomo = ISHFT(1_8,j)-1 @@ -194,6 +200,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 @@ -201,15 +208,18 @@ end subroutine get_phase_qp_to_cfg colsmax = cols end if ! i -> j - AIJpqMatrixDimsList(nsomoi,nsomoj,1,k,l,1) = rows - AIJpqMatrixDimsList(nsomoi,nsomoj,1,k,l,2) = cols + AIJpqMatrixDimsList(nsomoi,1,k,l,1) = rows + AIJpqMatrixDimsList(nsomoi,1,k,l,2) = cols end do end do end do end do ! Type ! 2. DOMO -> VMO - do i = 0+iand(nsomomin,1), NSOMOMax, 2 + !print *,"Doing DOMO->VMO" + AIJpqMatrixDimsList(NSOMOMin,2,1,1,1) = 1 + AIJpqMatrixDimsList(NSOMOMin,2,1,1,2) = 1 + do i = NSOMOMin, NSOMOMax, 2 Isomo = ISHFT(1_8,i)-1 tmpsomo = ISHFT(1_8,i+2)-1 do j = i+2,i+2, 2 @@ -242,6 +252,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 @@ -249,8 +260,8 @@ end subroutine get_phase_qp_to_cfg colsmax = cols end if ! i -> j - AIJpqMatrixDimsList(nsomoi,nsomoj,2,k,l,1) = rows - AIJpqMatrixDimsList(nsomoi,nsomoj,2,k,l,2) = cols + AIJpqMatrixDimsList(nsomoi,2,k,l,1) = rows + AIJpqMatrixDimsList(nsomoi,2,k,l,2) = cols end do end do end do @@ -258,15 +269,17 @@ end subroutine get_phase_qp_to_cfg ! Type ! 3. SOMO -> VMO !print *,"Doing SOMO->VMO" - do i = 2-iand(nsomomin,1), NSOMOMax, 2 + AIJpqMatrixDimsList(NSOMOMin,3,1,1,1) = 1 + AIJpqMatrixDimsList(NSOMOMin,3,1,1,2) = 1 + do i = NSOMOMin, 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) 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) @@ -281,6 +294,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 @@ -288,25 +302,27 @@ end subroutine get_phase_qp_to_cfg colsmax = cols end if ! i -> j - AIJpqMatrixDimsList(i,j,3,k,l,1) = rows - AIJpqMatrixDimsList(i,j,3,k,l,2) = cols + AIJpqMatrixDimsList(i,3,k,l,1) = rows + AIJpqMatrixDimsList(i,3,k,l,2) = cols end do end do end do end do ! Type - ! 4. DOMO -> VMO + ! 4. DOMO -> SOMO !print *,"Doing DOMO->SOMO" - do i = 2-iand(nsomomin,1), NSOMOMax, 2 + AIJpqMatrixDimsList(NSOMOMin,4,1,1,1) = 1 + AIJpqMatrixDimsList(NSOMOMin,4,1,1,2) = 1 + do i = NSOMOMin, 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 @@ -318,6 +334,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 @@ -325,15 +342,16 @@ end subroutine get_phase_qp_to_cfg colsmax = cols end if ! i -> j - AIJpqMatrixDimsList(i,j,4,k,l,1) = rows - AIJpqMatrixDimsList(i,j,4,k,l,2) = cols + AIJpqMatrixDimsList(i,4,k,l,1) = rows + AIJpqMatrixDimsList(i,4,k,l,2) = cols end do 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, (NBFMax,NBFmax,NSOMOMax+1,NSOMOMax+1,4,NSOMOMin:NSOMOMax)] use cfunctions implicit none BEGIN_DOC @@ -364,69 +382,78 @@ end subroutine get_phase_qp_to_cfg cols = -1 integer*8 MS MS = 0 - touch AIJpqMatrixDimsList real*8,dimension(:,:),allocatable :: meMatrix integer maxdim - !maxdim = max(rowsmax,colsmax) - ! allocate matrix - !allocate(AIJpqMatrixDimsList(NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,2)) + ! Type ! 1. SOMO -> SOMO - do i = 2, NSOMOMax, 2 + AIJpqContainer = 0.d0 + AIJpqContainer(1,1,1,1,1,NSOMOMin) = 1.0d0 + integer :: rows_old, cols_old + rows_old = -1 + cols_old = -1 + allocate(meMatrix(1,1)) + do i = NSOMOMin+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 - do k = 1,i - do l = 1,i + j=i-2 + if(j .GT. NSOMOMax .OR. j .LT. 0) cycle + nsomoi = i + do k = 1,i + orbp = k + do l = 1,i - ! Define Jsomo - if(k .NE. l) then - Jsomo = IBCLR(Isomo, k-1) - Jsomo = IBCLR(Jsomo, l-1) - nsomoi = i - nsomoj = j - else - Isomo = ISHFT(1_8,i)-1 - Jsomo = ISHFT(1_8,i)-1 - nsomoi = i - nsomoj = i - endif + ! Define Jsomo + if(k .NE. l) then + Jsomo = IBCLR(Isomo, k-1) + Jsomo = IBCLR(Jsomo, l-1) + nsomoj = j + else + Isomo = ISHFT(1_8,i)-1 + Jsomo = ISHFT(1_8,i)-1 + nsomoj = i + endif - AIJpqContainer(nsomoi,nsomoj,1,k,l,:,:) = 0.0d0 - call getApqIJMatrixDims(Isomo, & - Jsomo, & - MS, & - rows, & - cols) + call getApqIJMatrixDims(Isomo, & + Jsomo, & + MS, & + rows, & + cols) - orbp = k - orbq = l - allocate(meMatrix(rows,cols)) - meMatrix = 0.0d0 - ! fill matrix - call getApqIJMatrixDriver(Isomo, & - Jsomo, & - orbp, & - orbq, & - MS, & - NMO, & - meMatrix, & - rows, & - cols) - ! i -> j - do ri = 1,rows - do ci = 1,cols - AIJpqContainer(nsomoi,nsomoj,1,k,l,ri,ci) = meMatrix(ri, ci) - end do + orbq = l + if ((rows /= rows_old).or.(cols /= cols_old)) then + deallocate(meMatrix) + allocate(meMatrix(rows,cols)) + rows_old = rows + cols_old = cols + endif + meMatrix = 0.0d0 + ! fill matrix + call getApqIJMatrixDriver(Isomo, & + Jsomo, & + orbp, & + orbq, & + MS, & + NMO, & + meMatrix, & + rows, & + cols) + ! i -> j + do ri = 1,rows + do ci = 1,cols + AIJpqContainer(ri,ci,k,l,1,nsomoi) = meMatrix(ri, ci) end do - deallocate(meMatrix) end do end do end do end do + deallocate(meMatrix) + ! Type ! 2. DOMO -> VMO - do i = 0, NSOMOMax, 2 + !print *,"Doing DOMO -> VMO" + !AIJpqContainer(NSOMOMin,2,1,1,1,1) = 1.0d0 + AIJpqContainer(1,1,1,1,2,NSOMOMin) = 1.0d0 + do i = NSOMOMin, NSOMOMax, 2 Isomo = ISHFT(1_8,i)-1 tmpsomo = ISHFT(1_8,i+2)-1 do j = i+2,i+2, 2 @@ -451,7 +478,12 @@ end subroutine get_phase_qp_to_cfg nsomoj = j endif - AIJpqContainer(nsomoi,nsomoj,2,k,l,:,:) = 0.0d0 + !print *,"k,l=",k,l + !call debug_spindet(Jsomo,1) + !call debug_spindet(Isomo,1) + + !AIJpqContainer(nsomoi,2,k,l,:,:) = 0.0d0 + AIJpqContainer(:,:,k,l,2,nsomoi) = 0.0d0 call getApqIJMatrixDims(Isomo, & Jsomo, & MS, & @@ -472,10 +504,13 @@ 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 - AIJpqContainer(nsomoi,nsomoj,2,k,l,ri,ci) = meMatrix(ri, ci) + !AIJpqContainer(nsomoi,2,k,l,ri,ci) = meMatrix(ri, ci) + AIJpqContainer(ri,ci,k,l,2,nsomoi) = meMatrix(ri, ci) end do end do deallocate(meMatrix) @@ -485,13 +520,16 @@ end subroutine get_phase_qp_to_cfg end do ! Type ! 3. SOMO -> VMO - do i = 2, NSOMOMax, 2 + !print *,"Doing SOMO -> VMO" + !AIJpqContainer(NSOMOMin,3,1,1,1,1) = 1.0d0 + AIJpqContainer(1,1,1,1,3,NSOMOMin) = 1.0d0 + do i = NSOMOMin, 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) @@ -502,7 +540,12 @@ end subroutine get_phase_qp_to_cfg Jsomo = ISHFT(1_8,j)-1 endif - AIJpqContainer(i,j,3,k,l,:,:) = 0.0d0 + !print *,"k,l=",k,l + !call debug_spindet(Jsomo,1) + !call debug_spindet(Isomo,1) + + !AIJpqContainer(i,3,k,l,:,:) = 0.0d0 + AIJpqContainer(:,:,k,l,3,i) = 0.0d0 call getApqIJMatrixDims(Isomo, & Jsomo, & MS, & @@ -523,10 +566,13 @@ 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 - AIJpqContainer(i,j,3,k,l,ri,ci) = meMatrix(ri, ci) + !AIJpqContainer(i,3,k,l,ri,ci) = meMatrix(ri, ci) + AIJpqContainer(ri,ci,k,l,3,i) = meMatrix(ri, ci) end do end do deallocate(meMatrix) @@ -536,24 +582,28 @@ end subroutine get_phase_qp_to_cfg end do ! Type ! 4. DOMO -> SOMO - do i = 2, NSOMOMax, 2 + !print *,"Doing DOMO -> SOMO" + !AIJpqContainer(NSOMOMin,4,1,1,1,1) = 1.0d0 + AIJpqContainer(1,1,1,1,4,NSOMOMin) = 1.0d0 + do i = NSOMOMin+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 endif - AIJpqContainer(i,j,4,k,l,:,:) = 0.0d0 + !AIJpqContainer(i,4,k,l,:,:) = 0.0d0 + AIJpqContainer(:,:,k,l,4,i) = 0.0d0 call getApqIJMatrixDims(Isomo, & Jsomo, & MS, & @@ -575,10 +625,13 @@ 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 - AIJpqContainer(i,j,4,k,l,ri,ci) = meMatrix(ri, ci) + !AIJpqContainer(i,4,k,l,ri,ci) = meMatrix(ri, ci) + AIJpqContainer(ri,ci,k,l,4,i) = meMatrix(ri, ci) end do end do deallocate(meMatrix) @@ -588,31 +641,178 @@ end subroutine get_phase_qp_to_cfg 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 -!!!!!! + PROVIDE h_core_ri + ! 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 BEGIN_PROVIDER [ real*8, DetToCSFTransformationMatrix, (0:NSOMOMax,NBFMax,maxDetDimPerBF)] - &BEGIN_PROVIDER [ real*8, psi_coef_config, (n_CSF)] + &BEGIN_PROVIDER [ real*8, psi_coef_config, (n_CSF,1)] &BEGIN_PROVIDER [ integer, psi_config_data, (N_configuration,2)] + &BEGIN_PROVIDER [ integer, psi_csf_to_config_data, (n_CSF)] use cfunctions - 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*8 :: Isomo, Idomo + integer(bit_kind) :: Ialpha(N_int),Ibeta(N_int) integer :: rows, cols, i, j, k - integer :: startdet, enddet - integer*8 MS, Isomo, Idomo + integer :: startdet, enddet, idx + integer*8 MS integer ndetI integer :: getNSOMO real*8,dimension(:,:),allocatable :: tempBuffer real*8,dimension(:),allocatable :: tempCoeff real*8 :: norm_det1, phasedet + + integer :: nt + + norm_det1 = 0.d0 MS = elec_alpha_num - elec_beta_num + print *,"Maxbfdim=",NBFMax + print *,"Maxdetdim=",maxDetDimPerBF + print *,"n_CSF=",n_CSF + print *,"N_configurations=",N_configuration + print *,"n_core_orb=",n_core_orb ! initialization psi_coef_config = 0.d0 DetToCSFTransformationMatrix(0,:,:) = 1.d0 @@ -634,10 +834,13 @@ end subroutine get_phase_qp_to_cfg countcsf = 0 integer countdet countdet = 0 - integer idx integer istate istate = 1 + psi_csf_to_config_data(1) = 1 phasedet = 1.0d0 + call omp_set_max_active_levels(1) + !$OMP PARALLEL + !$OMP MASTER do i = 1,N_configuration startdet = psi_configuration_to_psi_det(1,i) enddet = psi_configuration_to_psi_det(2,i) @@ -655,6 +858,9 @@ end subroutine get_phase_qp_to_cfg 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 @@ -667,14 +873,1023 @@ end subroutine get_phase_qp_to_cfg 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 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) deallocate(tempCoeff) deallocate(tempBuffer) psi_config_data(i,1) = countcsf + 1 + do k=1,bfIcfg + psi_csf_to_config_data(countcsf+k) = i + enddo countcsf += bfIcfg psi_config_data(i,2) = countcsf enddo + print *,"Norm det=",norm_det1, size(psi_coef_config,1), " Dim csf=", countcsf + !$OMP END MASTER + !$OMP END PARALLEL + call omp_set_max_active_levels(4) END_PROVIDER + +subroutine obtain_connected_I_foralpha_fromfilterdlist(idxI, nconnectedJ, idslistconnectedJ, listconnectedJ, 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 ,intent(in) :: nconnectedJ + integer(bit_kind),intent(in) :: listconnectedJ(N_int,2,*) + integer(bit_kind),intent(in) :: Ialpha(N_int,2) + integer(bit_kind),intent(out) :: connectedI(N_int,2,*) + integer ,intent(in) :: idslistconnectedJ(*) + 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,kk,l,p,q,nsomoJ,nsomoalpha,starti,endi,extyp,nholes, idxJ + integer :: listholes(mo_num) + integer :: holetype(mo_num) + integer :: end_index + integer :: Nsomo_alpha + logical :: isOKlistJ + + PROVIDE DetToCSFTransformationMatrix + + isOKlistJ = .False. + + nconnectedI = 0 + end_index = N_configuration + + ! Since CFGs are sorted wrt to seniority + ! we don't have to search the full CFG list + Isomo = Ialpha(1,1) + Idomo = Ialpha(1,2) + Nsomo_alpha = POPCNT(Isomo) + end_index = min(N_configuration,cfg_seniority_index(min(Nsomo_alpha+4,elec_num))-1) + if(end_index .LT. 0) end_index= N_configuration + !end_index = N_configuration + + + p = 0 + q = 0 + do i=1,nconnectedJ + idxJ = idslistconnectedJ(i) + Isomo = Ialpha(1,1) + Idomo = Ialpha(1,2) + Jsomo = listconnectedJ(1,1,i) + Jdomo = listconnectedJ(1,2,i) + diffSOMO = IEOR(Isomo,Jsomo) + ndiffSOMO = POPCNT(diffSOMO) + diffDOMO = IEOR(Idomo,Jdomo) + xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO) + ndiffDOMO = POPCNT(diffDOMO) + nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO) + nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO + if((nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then + select case(ndiffDOMO) + case (0) + ! SOMO -> VMO + !print *,"obt SOMO -> VMO" + extyp = 3 + IJsomo = IEOR(Isomo, Jsomo) +IRP_IF WITHOUT_TRAILZ + p = (popcnt(ieor( IAND(Isomo,IJsomo), IAND(Isomo,IJsomo)-1)) -1) + 1 +IRP_ELSE + p = TRAILZ(IAND(Isomo,IJsomo)) + 1 +IRP_ENDIF + IJsomo = IBCLR(IJsomo,p-1) +IRP_IF WITHOUT_TRAILZ + q = (popcnt(ieor(IJsomo,IJsomo-1))-1) + 1 +IRP_ELSE + q = TRAILZ(IJsomo) + 1 +IRP_ENDIF + case (1) + ! DOMO -> VMO + ! or + ! SOMO -> SOMO + nsomoJ = POPCNT(Jsomo) + nsomoalpha = POPCNT(Isomo) + if(nsomoJ .GT. nsomoalpha) then + ! DOMO -> VMO + !print *,"obt DOMO -> VMO" + extyp = 2 +IRP_IF WITHOUT_TRAILZ + p = (popcnt(ieor( IEOR(Idomo,Jdomo), IEOR(Idomo,Jdomo)-1))-1) + 1 +IRP_ELSE + p = TRAILZ(IEOR(Idomo,Jdomo)) + 1 +IRP_ENDIF + Isomo = IEOR(Isomo, Jsomo) + Isomo = IBCLR(Isomo,p-1) +IRP_IF WITHOUT_TRAILZ + q = (popcnt(ieor(Isomo,Isomo-1))-1) + 1 +IRP_ELSE + q = TRAILZ(Isomo) + 1 +IRP_ENDIF + else + ! SOMO -> SOMO + !print *,"obt SOMO -> SOMO" + extyp = 1 +IRP_IF WITHOUT_TRAILZ + q = (popcnt(ieor( IEOR(Idomo,Jdomo), IEOR(Idomo,Jdomo)-1))-1) + 1 +IRP_ELSE + q = TRAILZ(IEOR(Idomo,Jdomo)) + 1 +IRP_ENDIF + Isomo = IEOR(Isomo, Jsomo) + Isomo = IBCLR(Isomo,q-1) +IRP_IF WITHOUT_TRAILZ + p = (popcnt(ieor(Isomo,Isomo-1))-1) + 1 +IRP_ELSE + p = TRAILZ(Isomo) + 1 +IRP_ENDIF + end if + case (2) + ! DOMO -> SOMO + !print *,"obt DOMO -> SOMO" + extyp = 4 + IJsomo = IEOR(Isomo, Jsomo) +IRP_IF WITHOUT_TRAILZ + p = (popcnt(ieor(IAND(Jsomo,IJsomo) ,IAND(Jsomo,IJsomo) -1))-1) + 1 +IRP_ELSE + p = TRAILZ(IAND(Jsomo,IJsomo)) + 1 +IRP_ENDIF + IJsomo = IBCLR(IJsomo,p-1) +IRP_IF WITHOUT_TRAILZ + q = (popcnt(ieor(IJsomo,IJsomo-1))-1) + 1 +IRP_ELSE + q = TRAILZ(IJsomo) + 1 +IRP_ENDIF + case default + print *,"something went wront in get connectedI" + end select + starti = psi_config_data(idxJ,1) + endi = psi_config_data(idxJ,2) + nconnectedI += 1 + connectedI(:,:,nconnectedI) = listconnectedJ(:,:,i) + idxs_connectedI(nconnectedI)=starti + excitationIds(1,nconnectedI)=p + excitationIds(2,nconnectedI)=q + excitationTypes(nconnectedI) = extyp + diagfactors(nconnectedI) = 1.0d0 + else if((ndiffSOMO + ndiffDOMO) .EQ. 0) then + ! find out all pq holes possible + nholes = 0 + ! holes in SOMO + Isomo = listconnectedJ(1,1,i) + Idomo = listconnectedJ(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(idxJ,1) + endi = psi_config_data(idxJ,2) + nconnectedI += 1 + connectedI(:,:,nconnectedI) = listconnectedJ(:,:,i) + idxs_connectedI(nconnectedI)=starti + excitationIds(1,nconnectedI)=p + excitationIds(2,nconnectedI)=q + excitationTypes(nconnectedI) = extyp + diagfactors(nconnectedI) = 1.0d0 + else + starti = psi_config_data(idxJ,1) + endi = psi_config_data(idxJ,2) + nconnectedI += 1 + connectedI(:,:,nconnectedI) = listconnectedJ(:,:,i) + idxs_connectedI(nconnectedI)=starti + excitationIds(1,nconnectedI)=p + excitationIds(2,nconnectedI)=q + excitationTypes(nconnectedI) = extyp + diagfactors(nconnectedI) = 2.0d0 + endif + enddo + endif + end do + +end subroutine obtain_connected_I_foralpha_fromfilterdlist + + +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(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) + Idomo = Ialpha(1,2) + Jsomo = Jcfg(1,1) + Jdomo = Jcfg(1,2) + pos0prev = 0 + pmodel = p + qmodel = q + + if(p .EQ. q) then + pmodel = 1 + qmodel = 1 + else + 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 + +subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze, istart, iend, ishift, istep) + implicit none + use bitmasks + use omp_lib + 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(n_st,sze) + real*8,intent(out) :: psi_out(n_st,sze) + 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),dimension(:,:,:),allocatable :: listconnectedJ + integer(bit_kind),dimension(:,:,:),allocatable :: alphas_Icfg + integer(bit_kind),dimension(:,:,:),allocatable :: singlesI + integer(bit_kind),dimension(:,:,:),allocatable :: connectedI_alpha + integer,dimension(:),allocatable :: idxs_singlesI + integer,dimension(:),allocatable :: idxs_connectedI_alpha + integer,dimension(:,:),allocatable :: excitationIds_single + integer,dimension(:),allocatable :: excitationTypes_single + integer,dimension(:,:),allocatable :: excitationIds + integer,dimension(:),allocatable :: excitationTypes + integer,dimension(:),allocatable :: idslistconnectedJ + real*8,dimension(:),allocatable :: diagfactors + 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, num_threads_max + integer :: nconnectedJ,nconnectedtotalmax,nconnectedmaxJ,maxnalphas,ntotJ + 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,dimension(:),allocatable :: psi_out_tmp + real*8,dimension(:,:),allocatable :: CCmattmp + real*8, external :: mo_two_e_integral + real*8, external :: get_two_e_integral + real*8,dimension(:),allocatable:: diag_energies + real*8 :: tmpvar, tmptot + + integer(omp_lock_kind), allocatable :: lock(:) + call omp_set_max_active_levels(1) + + allocate(lock(sze)) + do i=1,sze + call omp_init_lock(lock(i)) + enddo + + !print *," sze = ",sze + allocate(diag_energies(n_CSF)) + call calculate_preconditioner_cfg(diag_energies) + + MS = 0 + norm_coef_cfg=0.d0 + + psi_out=0.d0 + + istart_cfg = psi_csf_to_config_data(istart) + iend_cfg = psi_csf_to_config_data(iend) + + !nconnectedtotalmax = 1000 + !nconnectedmaxJ = 1000 + maxnalphas = elec_num*mo_num + Icfg(1,1) = psi_configuration(1,1,1) + Icfg(1,2) = psi_configuration(1,2,1) + allocate(listconnectedJ(N_INT,2,max(sze,100))) + allocate(idslistconnectedJ(max(sze,100))) + call obtain_connected_J_givenI(1, Icfg, listconnectedJ, idslistconnectedJ, nconnectedmaxJ, nconnectedtotalmax) + deallocate(listconnectedJ) + deallocate(idslistconnectedJ) + + integer*8, allocatable :: bit_tmp(:) + integer*8, external :: configuration_search_key + double precision :: diagfactors_0 + allocate( bit_tmp(0:N_configuration+1)) + do j=1,N_configuration + bit_tmp(j) = configuration_search_key(psi_configuration(1,1,j),N_int) + enddo + + call omp_set_max_active_levels(1) + !$OMP PARALLEL & + !$OMP DEFAULT(NONE) & + !$OMP private(i,icfg, isomo, idomo, NSOMOI, NSOMOJ, nholes, k, listholes,& + !$OMP holetype, vmotype, nvmos, listvmos, starti, endi, & + !$OMP nsinglesI, singlesI,idxs_singlesI,excitationIds_single,& + !$OMP excitationTypes_single, idxI, p, q, extype, pmodel, qmodel,& + !$OMP Jsomo, Jdomo, startj, endj, kk, jj, ii, cnti, cntj, meCC1,& + !$OMP nconnectedJ,listconnectedJ,idslistconnectedJ,ntotJ, & + !$OMP Nalphas_Icfg,alphas_Icfg,connectedI_alpha, & + !$OMP idxs_connectedI_alpha,nconnectedI,excitationIds,excitationTypes,diagfactors,& + !$OMP totcolsTKI,rowsTKI,NSOMOalpha,rowsikpq, & + !$OMP colsikpq, GIJpqrs,TKIGIJ,j,l,m,TKI,CCmattmp, moi, moj, mok, mol,& + !$OMP diagfac, tmpvar, diagfactors_0) & + !$OMP shared(istart_cfg, iend_cfg, psi_configuration, mo_num, psi_config_data,& + !$OMP N_int, N_st, psi_out, psi_in, h_core_ri, AIJpqContainer,& + !$OMP sze, NalphaIcfg_list,alphasIcfg_list, bit_tmp, & + !$OMP AIJpqMatrixDimsList, diag_energies, n_CSF, lock, NBFmax,nconnectedtotalmax, nconnectedmaxJ,maxnalphas,& + !$OMP num_threads_max) + + allocate(singlesI(N_INT,2,max(sze,100))) + allocate(idxs_singlesI(max(sze,100))) + allocate(excitationIds_single(2,max(sze,100))) + allocate(excitationTypes_single(max(sze,100))) +! + + !!! Single Excitations !!! + + !$OMP DO SCHEDULE(dynamic,16) + do i=istart_cfg,iend_cfg + + ! if Seniority_range > 8 then + ! continue + ! else + ! cycle + + Icfg(1,1) = psi_configuration(1,1,i) + Icfg(1,2) = psi_configuration(1,2,i) + Isomo = Icfg(1,1) + Idomo = Icfg(1,2) + NSOMOI = getNSOMO(Icfg) + + ! 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(bit_tmp,Icfg,singlesI,idxs_singlesI,excitationIds_single,& + excitationTypes_single,nsinglesI,N_int) + + do j = 1,nsinglesI + idxI = idxs_singlesI(j) + NSOMOJ = getNSOMO(singlesI(1,1,j)) + p = excitationIds_single(1,j) + q = excitationIds_single(2,j) + extype = excitationTypes_single(j) + ! Off diagonal terms + call convertOrbIdsToModelSpaceIds(Icfg, singlesI(1,1,j), p, q, extype, pmodel, qmodel) + Jsomo = singlesI(1,1,j) + Jdomo = singlesI(1,2,j) + + ! 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 ii = starti, endi + cnti = ii-starti+1 + do jj = startj, endj + cntj = jj-startj+1 + meCC1 = AIJpqContainer(cnti,cntj,pmodel,qmodel,extype,NSOMOI)* h_core_ri(p,q) + call omp_set_lock(lock(jj)) + do kk = 1,n_st + psi_out(kk,jj) = psi_out(kk,jj) + meCC1 * psi_in(kk,ii) + enddo + call omp_unset_lock(lock(jj)) + 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 + !$OMP END DO + deallocate(singlesI) + deallocate(idxs_singlesI) + deallocate(excitationIds_single) + deallocate(excitationTypes_single) + + allocate(listconnectedJ(N_INT,2,max(sze,100))) + allocate(alphas_Icfg(N_INT,2,max(sze,100))) + allocate(connectedI_alpha(N_INT,2,max(sze,100))) + allocate(idxs_connectedI_alpha(max(sze,100))) + allocate(excitationIds(2,max(sze,100))) + allocate(excitationTypes(max(sze,100))) + allocate(diagfactors(max(sze,100))) + allocate(idslistconnectedJ(max(sze,100))) + allocate(CCmattmp(n_st,NBFmax)) + + ! Loop over all selected configurations + !$OMP DO SCHEDULE(dynamic,16) + do i = istart_cfg,iend_cfg +! TKI = 0.d0 +! GIJpqrs = 0.d0 +! TKIGIJ = 0.d0 + + ! if Seniority_range > 8 then + ! continue + ! else + ! cycle + + 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 + Nalphas_Icfg = 0 + ! TODO: + ! test if size(alphas_Icfg,1) < Nmo**2) then deallocate + allocate + + Nalphas_Icfg = NalphaIcfg_list(i) + alphas_Icfg(1:n_int,1:2,1:Nalphas_Icfg) = alphasIcfg_list(1:n_int,1:2,i,1:Nalphas_Icfg) + if(Nalphas_Icfg .GT. maxnalphas) then + print *,"Nalpha > maxnalpha" + endif + + call obtain_connected_J_givenI(i, Icfg, listconnectedJ, idslistconnectedJ, nconnectedJ, ntotJ) + + ! 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_fromfilterdlist(i,nconnectedJ, idslistconnectedJ, & + listconnectedJ, 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 + NSOMOalpha = getNSOMO(alphas_Icfg(:,:,k)) + do j = 1,nconnectedI + 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 + rowsikpq = AIJpqMatrixDimsList(NSOMOalpha,extype,pmodel,qmodel,1) + colsikpq = AIJpqMatrixDimsList(NSOMOalpha,extype,pmodel,qmodel,2) + totcolsTKI += colsikpq + rowsTKI = rowsikpq + enddo + + allocate(TKI(n_st,rowsTKI,totcolsTKI)) ! coefficients of CSF + ! Initialize the integral container + ! dims : (totcolsTKI, nconnectedI) + allocate(GIJpqrs(totcolsTKI,nconnectedI)) ! gpqrs + allocate(TKIGIJ(n_st,rowsTKI,nconnectedI)) ! TKI * gpqrs + !print *,"\t---rowsTKI=",rowsTKI," totCols=",totcolsTKI + + totcolsTKI = 0 + do j = 1,nconnectedI + 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,extype,pmodel,qmodel,1) + colsikpq = AIJpqMatrixDimsList(NSOMOalpha,extype,pmodel,qmodel,2) + rowsTKI = rowsikpq + do m = 1,colsikpq + do l = 1,rowsTKI + do kk = 1,n_st + TKI(kk,l,totcolsTKI+m) = AIJpqContainer(l,m,pmodel,qmodel,extype,NSOMOalpha) & + * psi_in(kk,idxs_connectedI_alpha(j)+m-1) + enddo + enddo + enddo + diagfactors_0 = diagfactors(j)*0.5d0 + moi = excitationIds(1,j) ! p + mok = excitationIds(2,j) ! q + do l=1,nconnectedI + moj = excitationIds(2,l) ! s + mol = excitationIds(1,l) ! r + diagfac = diagfactors_0 * diagfactors(l)* mo_two_e_integral(mok,mol,moi,moj)! g(pq,sr) = + do m = 1,colsikpq + ! = (ik|jl) + GIJpqrs(totcolsTKI+m,l) = diagfac + enddo + enddo + totcolsTKI += colsikpq + enddo + + + ! Do big BLAS + call dgemm('N','N', rowsTKI*n_st, nconnectedI, totcolsTKI, 1.d0, & + TKI, size(TKI,1)*size(TKI,2), GIJpqrs, size(GIJpqrs,1), 0.d0, & + TKIGIJ , size(TKIGIJ,1)*size(TKIGIJ,2) ) + + + ! Collect the result + totcolsTKI = 0 + do j = 1,nconnectedI + NSOMOI = getNSOMO(connectedI_alpha(1,1,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) + rowsikpq = AIJpqMatrixDimsList(NSOMOalpha,extype,pmodel,qmodel,1) + colsikpq = AIJpqMatrixDimsList(NSOMOalpha,extype,pmodel,qmodel,2) + rowsTKI = rowsikpq + CCmattmp = 0.d0 + + call dgemm('N','N', n_st, colsikpq, rowsTKI, 1.d0, & + TKIGIJ(1,1,j), size(TKIGIJ,1), & + AIJpqContainer(1,1,pmodel,qmodel,extype,NSOMOalpha), & + size(AIJpqContainer,1), 0.d0, & + CCmattmp, size(CCmattmp,1) ) + + do m = 1,colsikpq + call omp_set_lock(lock(idxs_connectedI_alpha(j)+m-1)) + do kk = 1,n_st + psi_out(kk,idxs_connectedI_alpha(j)+m-1) += CCmattmp(kk,m) + enddo + call omp_unset_lock(lock(idxs_connectedI_alpha(j)+m-1)) + enddo + totcolsTKI += colsikpq + enddo + + deallocate(TKI) ! coefficients of CSF + deallocate(GIJpqrs) ! gpqrs + deallocate(TKIGIJ) ! gpqrs + + enddo ! loop over alphas + enddo ! loop over I + !$OMP END DO + call omp_set_max_active_levels(4) + deallocate(CCmattmp) + deallocate(connectedI_alpha) + deallocate(idxs_connectedI_alpha) + deallocate(excitationIds) + deallocate(excitationTypes) + deallocate(diagfactors) + + + ! Add the diagonal contribution + !$OMP DO + do i = 1,n_CSF + do kk=1,n_st + psi_out(kk,i) += diag_energies(i)*psi_in(kk,i) + enddo + enddo + !$OMP END DO + + !$OMP END PARALLEL + call omp_set_max_active_levels(4) + + deallocate(diag_energies) + deallocate(bit_tmp) + +end subroutine calculate_sigma_vector_cfg_nst_naive_store + + + + +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),dimension(:,:,:),allocatable :: alphas_Icfg + integer(bit_kind),dimension(:,:,:),allocatable :: singlesI + integer(bit_kind),dimension(:,:,:),allocatable :: connectedI_alpha + integer,dimension(:),allocatable :: idxs_singlesI + integer,dimension(:),allocatable :: idxs_connectedI_alpha + integer,dimension(:,:),allocatable :: excitationIds_single + integer,dimension(:),allocatable :: excitationTypes_single + integer,dimension(:,:),allocatable :: excitationIds + integer,dimension(:),allocatable :: excitationTypes + real*8,dimension(:),allocatable :: diagfactors + 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) + + ! allocate + allocate(alphas_Icfg(N_INT,2,max(sze/2,100))) + allocate(singlesI(N_INT,2,max(sze/2,100))) + allocate(connectedI_alpha(N_INT,2,max(sze/2,100))) + allocate(idxs_singlesI(max(sze/2,100))) + allocate(idxs_connectedI_alpha(max(sze/2,100))) + allocate(excitationIds_single(2,max(sze/2,100))) + allocate(excitationTypes_single(max(sze/2,100))) + allocate(excitationIds(2,max(sze/2,100))) + allocate(excitationTypes(max(sze/2,100))) + allocate(diagfactors(max(sze/2,100))) + + + !print *," sze = ",sze + call calculate_preconditioner_cfg(diag_energies) + + MS = 0 + norm_coef_cfg=0.d0 + + psi_out=0.d0 + + istart_cfg = psi_csf_to_config_data(istart) + iend_cfg = psi_csf_to_config_data(iend) + + + !!! Single Excitations !!! + do i=istart_cfg,iend_cfg + print *,"I=",i + + ! if Seniority_range > 8 then + ! continue + ! else + ! cycle + + 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 + Nalphas_Icfg = 0 + ! TODO: + ! test if size(alphas_Icfg,1) < Nmo**2) then deallocate + allocate + !call obtain_associated_alphaI(i, Icfg, alphas_Icfg, Nalphas_Icfg) + Nalphas_Icfg = NalphaIcfg_list(i) + alphas_Icfg(1:N_int,1:2,1:Nalphas_Icfg) = alphasIcfg_list(1:n_int,1:2,i,1:Nalphas_Icfg) + + ! 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) + + totcolsTKI = 0 + rowsTKI = -1 + do j = 1,nconnectedI + NSOMOalpha = getNSOMO(alphas_Icfg(1,1,k)) + NSOMOI = getNSOMO(connectedI_alpha(1,1,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,extype,pmodel,qmodel,1) + colsikpq = AIJpqMatrixDimsList(NSOMOalpha,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(n_st,rowsTKI,totcolsTKI)) ! coefficients of CSF + ! Initialize the inegral container + ! dims : (totcolsTKI, nconnectedI) + allocate(GIJpqrs(totcolsTKI,nconnectedI)) ! gpqrs + allocate(TKIGIJ(n_st,rowsTKI,nconnectedI)) ! TKI * gpqrs + + totcolsTKI = 0 + do j = 1,nconnectedI + NSOMOalpha = getNSOMO(alphas_Icfg(1,1,k)) + NSOMOI = getNSOMO(connectedI_alpha(1,1,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) + rowsikpq = AIJpqMatrixDimsList(NSOMOalpha,extype,pmodel,qmodel,1) + colsikpq = AIJpqMatrixDimsList(NSOMOalpha,extype,pmodel,qmodel,2) + do m = 1,colsikpq + do l = 1,rowsTKI + do kk = 1,n_st + TKI(kk,l,totcolsTKI+m) = AIJpqContainer(l,m,pmodel,qmodel,extype,NSOMOalpha) * psi_in(kk,idxs_connectedI_alpha(j)+m-1) + 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)*size(TKI,2), GIJpqrs, size(GIJpqrs,1), 0.d0,& + TKIGIJ , size(TKIGIJ,1)*size(TKIGIJ,2) ) + + + ! Collect the result + totcolsTKI = 0 + do j = 1,nconnectedI + NSOMOalpha = getNSOMO(alphas_Icfg(1,1,k)) + NSOMOI = getNSOMO(connectedI_alpha(1,1,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,extype,pmodel,qmodel,1) + colsikpq = AIJpqMatrixDimsList(NSOMOalpha,extype,pmodel,qmodel,2) + do m = 1,colsikpq + do l = 1,rowsTKI + do kk = 1,n_st + psi_out(kk,idxs_connectedI_alpha(j)+m-1) = psi_out(kk,idxs_connectedI_alpha(j)+m-1) + & + AIJpqContainer(l,m,pmodel,qmodel,extype,NSOMOalpha) * TKIGIJ(kk,l,j) + 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 + deallocate(connectedI_alpha) + deallocate(idxs_connectedI_alpha) + deallocate(excitationIds) + deallocate(excitationTypes) + deallocate(diagfactors) + + + ! Add the diagonal contribution + do i = 1,n_CSF + do kk=1,n_st + psi_out(kk,i) += diag_energies(i)*psi_in(kk,i) + enddo + enddo + call omp_set_max_active_levels(4) + +end subroutine calculate_sigma_vector_cfg_nst_naive_store diff --git a/src/davidson/diagonalization_hcsf_dressed.irp.f b/src/davidson/diagonalization_hcsf_dressed.irp.f index 0c3c6f92..b6d5a3f8 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,ii logical, intent(inout) :: converged double precision, external :: u_dot_v, u_dot_u @@ -110,6 +110,7 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N integer :: order(N_st_diag_in) double precision :: cmax double precision, allocatable :: U(:,:), U_csf(:,:), overlap(:,:) + double precision, allocatable :: tmpU(:,:), tmpW(:,:) double precision, pointer :: W(:,:), W_csf(:,:) logical :: disk_based double precision :: energy_shift(N_st_diag_in*davidson_sze_max) @@ -303,11 +304,42 @@ 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) + PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals if ((sze > 100000).and.distributed_davidson) then - call H_u_0_nstates_zmq (W,U,N_st_diag,sze) + !call H_u_0_nstates_zmq (W,U,N_st_diag,sze) + allocate(tmpW(N_st_diag,sze_csf)) + allocate(tmpU(N_st_diag,sze_csf)) + do kk=1,N_st_diag + do ii=1,sze_csf + tmpU(kk,ii) = U_csf(ii,shift+kk) + enddo + enddo + call calculate_sigma_vector_cfg_nst_naive_store(tmpW,tmpU,N_st_diag,sze_csf,1,sze_csf,0,1) + do kk=1,N_st_diag + do ii=1,sze_csf + W_csf(ii,shift+kk)=tmpW(kk,ii) + enddo + enddo + deallocate(tmpW) + deallocate(tmpU) else - call H_u_0_nstates_openmp(W,U,N_st_diag,sze) + !call H_u_0_nstates_openmp(W,U,N_st_diag,sze) + allocate(tmpW(N_st_diag,sze_csf)) + allocate(tmpU(N_st_diag,sze_csf)) + do kk=1,N_st_diag + do ii=1,sze_csf + tmpU(kk,ii) = U_csf(ii,shift+kk) + enddo + enddo + call calculate_sigma_vector_cfg_nst_naive_store(tmpW,tmpU,N_st_diag,sze_csf,1,sze_csf,0,1) + do kk=1,N_st_diag + do ii=1,sze_csf + W_csf(ii,shift+kk)=tmpW(kk,ii) + enddo + enddo + deallocate(tmpW) + deallocate(tmpU) endif ! else ! ! Already computed in update below @@ -351,7 +383,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 = = ! -------------------------------------------