From 13154e765b48517e94aa4e50c19dac5f417c8348 Mon Sep 17 00:00:00 2001 From: v1j4y Date: Wed, 10 Feb 2021 23:12:04 +0100 Subject: [PATCH] Working on fixing bugs. #143. --- .../configuration_CI_sigma_helpers.irp.f | 146 ++++++++++------- .../configuration_CI_sigma_helpers.org | 152 +++++++++++------- src/determinants/configurations.irp.f | 35 ++++ src/determinants/create_excitations.irp.f | 22 ++- src/mo_two_e_ints/core_quantities.irp.f | 2 +- 5 files changed, 233 insertions(+), 124 deletions(-) diff --git a/src/determinants/configuration_CI_sigma_helpers.irp.f b/src/determinants/configuration_CI_sigma_helpers.irp.f index 905470ca..bf6ec9bc 100644 --- a/src/determinants/configuration_CI_sigma_helpers.irp.f +++ b/src/determinants/configuration_CI_sigma_helpers.irp.f @@ -1,4 +1,4 @@ - subroutine obtain_associated_alphaI(idxI, Icfg, alphasIcfg, NalphaIcfg) + subroutine obtain_associated_alphaI(idxI, Icfg, alphasIcfg, NalphaIcfg, factor_alphaI) implicit none use bitmasks BEGIN_DOC @@ -10,6 +10,7 @@ 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) @@ -261,7 +262,30 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI integer*8 :: diffDOMO integer :: ndiffSOMO integer :: ndiffDOMO - integer :: i,j,k,l,p,q,nsomoJ,nsomoalpha,starti,endi,extyp + integer :: 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 @@ -281,17 +305,22 @@ 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 + if((ndiffSOMO + ndiffDOMO) .EQ. 0) cycle !print *,"-I--i=",i,diffSOMO,diffDOMO!Isomo,Jsomo,ndiffSOMO,ndiffDOMO !print *,POPCNT(IEOR(diffSOMO,diffDOMO)), ndiffDOMO if(POPCNT(IEOR(diffSOMO,diffDOMO)) .LE. 1 .AND. ndiffDOMO .LT. 3) 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(AND(Isomo,IJsomo)) + 1 + p = TRAILZ(IAND(Isomo,IJsomo)) + 1 IJsomo = IBCLR(IJsomo,p-1) q = TRAILZ(IJsomo) + 1 case (1) @@ -322,7 +351,7 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI !print *,"obt DOMO -> SOMO" extyp = 4 IJsomo = IEOR(Isomo, Jsomo) - p = TRAILZ(AND(Jsomo,IJsomo)) + 1 + p = TRAILZ(IAND(Jsomo,IJsomo)) + 1 IJsomo = IBCLR(IJsomo,p-1) q = TRAILZ(IJsomo) + 1 case default @@ -336,11 +365,10 @@ 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 + print *,"------ > output p,q in obt=",p,q endif end do - end subroutine obtain_connected_I_foralpha function getNSOMO(Icfg) result(NSOMO) @@ -390,55 +418,61 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod pos0prev = 0 pmodel = p qmodel = q - !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 = AND(Isomo,mask) - pmodel = POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) - mask = ISHFT(1_8,q) - 1 - Isomotmp = AND(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 = AND(Jsomo,mask) - pmodel = POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) - mask = ISHFT(1_8,q) - 1 - Jsomotmp = AND(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 = AND(Isomo,mask) - pmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) - mask = ISHFT(1_8,q) - 1 - Jsomo = AND(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 = AND(Jsomo,mask) - pmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) - mask = ISHFT(1_8,q) - 1 - Isomo = AND(Isomo,mask) - qmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) - case default - print *,"something is wrong in convertOrbIdsToModelSpaceIds" + + 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/determinants/configuration_CI_sigma_helpers.org b/src/determinants/configuration_CI_sigma_helpers.org index fc72077f..4fc34da7 100644 --- a/src/determinants/configuration_CI_sigma_helpers.org +++ b/src/determinants/configuration_CI_sigma_helpers.org @@ -16,7 +16,7 @@ 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) + subroutine obtain_associated_alphaI(idxI, Icfg, alphasIcfg, NalphaIcfg, factor_alphaI) implicit none use bitmasks BEGIN_DOC @@ -28,6 +28,7 @@ the input determinant \(|D_I\rangle\). 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) @@ -286,7 +287,30 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI integer*8 :: diffDOMO integer :: ndiffSOMO integer :: ndiffDOMO - integer :: i,j,k,l,p,q,nsomoJ,nsomoalpha,starti,endi,extyp + integer :: 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 @@ -306,17 +330,22 @@ 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 + if((ndiffSOMO + ndiffDOMO) .EQ. 0) cycle !print *,"-I--i=",i,diffSOMO,diffDOMO!Isomo,Jsomo,ndiffSOMO,ndiffDOMO !print *,POPCNT(IEOR(diffSOMO,diffDOMO)), ndiffDOMO if(POPCNT(IEOR(diffSOMO,diffDOMO)) .LE. 1 .AND. ndiffDOMO .LT. 3) 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(AND(Isomo,IJsomo)) + 1 + p = TRAILZ(IAND(Isomo,IJsomo)) + 1 IJsomo = IBCLR(IJsomo,p-1) q = TRAILZ(IJsomo) + 1 case (1) @@ -347,7 +376,7 @@ subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI !print *,"obt DOMO -> SOMO" extyp = 4 IJsomo = IEOR(Isomo, Jsomo) - p = TRAILZ(AND(Jsomo,IJsomo)) + 1 + p = TRAILZ(IAND(Jsomo,IJsomo)) + 1 IJsomo = IBCLR(IJsomo,p-1) q = TRAILZ(IJsomo) + 1 case default @@ -361,11 +390,10 @@ 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 + print *,"------ > output p,q in obt=",p,q endif end do - end subroutine obtain_connected_I_foralpha #+end_src @@ -438,56 +466,62 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod pos0prev = 0 pmodel = p qmodel = q - !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 = AND(Isomo,mask) - pmodel = POPCNT(mask) - POPCNT(XOR(Isomotmp,mask)) - mask = ISHFT(1_8,q) - 1 - Isomotmp = AND(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 = AND(Jsomo,mask) - pmodel = POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask)) - mask = ISHFT(1_8,q) - 1 - Jsomotmp = AND(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 = AND(Isomo,mask) - pmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) - mask = ISHFT(1_8,q) - 1 - Jsomo = AND(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 = AND(Jsomo,mask) - pmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask)) - mask = ISHFT(1_8,q) - 1 - Isomo = AND(Isomo,mask) - qmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask)) - case default - print *,"something is wrong in convertOrbIdsToModelSpaceIds" + + 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 #+end_src @@ -502,7 +536,7 @@ end subroutine convertOrbIdsToModelSpaceIds print *,mask print *,POPCNT(mask) isomo = 144 - isomo = AND(isomo,mask) + isomo = IAND(isomo,mask) print *,isomo print *,XOR(isomo,mask) print *,POPCNT(mask) - POPCNT(XOR(isomo,mask)) @@ -518,8 +552,8 @@ end subroutine convertOrbIdsToModelSpaceIds #+begin_src fortran print *,IBSET(0_8,4)-1 - print *,POPCNT(IBSET(0_8,4)-1) - POPCNT(AND(716,IBSET(0_8,4)-1)) - print *,POPCNT(IBSET(0_8,8)-1) - POPCNT(AND(716,IBSET(0_8,8)-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: diff --git a/src/determinants/configurations.irp.f b/src/determinants/configurations.irp.f index 0cb4b629..739b0829 100644 --- a/src/determinants/configurations.irp.f +++ b/src/determinants/configurations.irp.f @@ -729,6 +729,41 @@ BEGIN_PROVIDER [ integer(bit_kind), dominant_dets_of_cfgs, (N_int,2,N_dominant_d enddo END_PROVIDER +subroutine binary_search_cfg(cfgInp,addcfg) + use bitmasks + implicit none + BEGIN_DOC + ! Documentation for binary_search + ! + ! Does a binary search to find + ! the address of a configuration in a list of + ! configurations. + END_DOC + integer(bit_kind), intent(in) :: cfgInp(N_int,2) + integer , intent(out) :: addcfg + integer :: i,j,k,r,l + integer*8 :: key, key2 + logical :: found + !integer*8, allocatable :: bit_tmp(:) + !integer*8, external :: configuration_search_key + + !allocate(bit_tmp(0:N_configuration)) + !bit_tmp(0) = 0 + do i=1,N_configuration + !bit_tmp(i) = configuration_search_key(psi_configuration(1,1,i),N_int) + found = .True. + do k=1,N_int + found = found .and. (psi_configuration(k,1,i) == cfgInp(k,1)) & + .and. (psi_configuration(k,2,i) == cfgInp(k,2)) + enddo + if (found) then + addcfg = i + exit + endif + enddo + +end subroutine + BEGIN_PROVIDER [ integer, psi_configuration_to_psi_det, (2,N_configuration) ] &BEGIN_PROVIDER [ integer, psi_configuration_to_psi_det_data, (N_det) ] diff --git a/src/determinants/create_excitations.irp.f b/src/determinants/create_excitations.irp.f index c0b001c5..8611ca22 100644 --- a/src/determinants/create_excitations.irp.f +++ b/src/determinants/create_excitations.irp.f @@ -336,7 +336,7 @@ subroutine generate_all_singles_cfg(cfg,singles,n_singles,Nint) enddo end -subroutine generate_all_singles_cfg_with_type(cfg,singles,idxs_singles,pq_singles,ex_type_singles,n_singles,Nint) +subroutine generate_all_singles_cfg_with_type(cfgInp,singles,idxs_singles,pq_singles,ex_type_singles,n_singles,Nint) implicit none use bitmasks BEGIN_DOC @@ -352,11 +352,12 @@ subroutine generate_all_singles_cfg_with_type(cfg,singles,idxs_singles,pq_single integer, intent(inout) :: n_singles integer, intent(out) :: idxs_singles(*) integer, intent(out) :: ex_type_singles(*) - real*8 , intent(out) :: pq_singles(2,*) - integer(bit_kind), intent(in) :: cfg(Nint,2) + integer, intent(out) :: pq_singles(2,*) + integer(bit_kind), intent(in) :: cfgInp(Nint,2) integer(bit_kind), intent(out) :: singles(Nint,2,*) + integer(bit_kind) :: Jdet(Nint,2) - integer :: i,k, n_singles_ma, i_hole, i_particle, ex_type + integer :: i,k, n_singles_ma, i_hole, i_particle, ex_type, addcfg integer(bit_kind) :: single(Nint,2) logical :: i_ok @@ -364,17 +365,22 @@ subroutine generate_all_singles_cfg_with_type(cfg,singles,idxs_singles,pq_single !TODO !Make list of Somo and Domo for holes !Make list of Unocc and Somo for particles - do i_hole = 1, mo_num - do i_particle = 1, mo_num - call do_single_excitation_cfg_with_type(cfg,single,i_hole,i_particle,ex_type,i_ok) + do i_hole = 1+n_core_orb, n_core_orb + n_act_orb + do i_particle = 1+n_core_orb, n_core_orb + n_act_orb + if(i_hole .EQ. i_particle) cycle + addcfg = -1 + call do_single_excitation_cfg_with_type(cfgInp,single,i_hole,i_particle,ex_type,i_ok) if (i_ok) then + call binary_search_cfg(single,addcfg) + if(addcfg .EQ. -1) cycle n_singles = n_singles + 1 do k=1,Nint singles(k,1,n_singles) = single(k,1) singles(k,2,n_singles) = single(k,2) ex_type_singles(n_singles) = ex_type pq_singles(1,n_singles) = i_hole ! p - pq_singles(1,n_singles) = i_particle ! q + pq_singles(2,n_singles) = i_particle ! q + idxs_singles(n_singles) = addcfg enddo endif enddo diff --git a/src/mo_two_e_ints/core_quantities.irp.f b/src/mo_two_e_ints/core_quantities.irp.f index b0f48b6d..b764a1a6 100644 --- a/src/mo_two_e_ints/core_quantities.irp.f +++ b/src/mo_two_e_ints/core_quantities.irp.f @@ -53,7 +53,7 @@ BEGIN_PROVIDER [ double precision, h_core_ri, (mo_num, mo_num) ] enddo do k=1,mo_num do i=1,mo_num - h_core_ri(i,j) = h_core_ri(i,j) - big_array_exchange_integrals(i,k,j) + h_core_ri(i,j) = h_core_ri(i,j) - 0.5 * big_array_exchange_integrals(k,i,j) enddo enddo enddo