Working on fixing bugs. #143.

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,6 +418,11 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod
   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)
@@ -401,30 +434,30 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod
           ! remove all domos
           !print *,"type -> SOMO -> SOMO"
           mask = ISHFT(1_8,p) - 1
-        Isomotmp = AND(Isomo,mask)
+          Isomotmp = IAND(Isomo,mask)
           pmodel = POPCNT(mask) - POPCNT(XOR(Isomotmp,mask))
           mask = ISHFT(1_8,q) - 1
-        Isomotmp = AND(Isomo,mask)
+          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 = AND(Jsomo,mask)
+          Jsomotmp = IAND(Jsomo,mask)
           pmodel = POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask))
           mask = ISHFT(1_8,q) - 1
-        Jsomotmp = AND(Jsomo,mask)
+          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 = AND(Isomo,mask)
+          Isomo = IAND(Isomo,mask)
           pmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask))
           mask = ISHFT(1_8,q) - 1
-        Jsomo = AND(Jsomo,mask)
+          Jsomo = IAND(Jsomo,mask)
           qmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask))
        case (4)
           ! DOMO -> SOMO
@@ -432,13 +465,14 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod
           !print *,"type -> DOMO -> SOMO"
           !Isomo = IEOR(Isomo,Jsomo)
           mask = ISHFT(1_8,p) - 1
-        Jsomo = AND(Jsomo,mask)
+          Jsomo = IAND(Jsomo,mask)
           pmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask))
           mask = ISHFT(1_8,q) - 1
-        Isomo = AND(Isomo,mask)
+          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

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,6 +466,11 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod
   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)
@@ -449,30 +482,30 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod
           ! remove all domos
           !print *,"type -> SOMO -> SOMO"
           mask = ISHFT(1_8,p) - 1
-        Isomotmp = AND(Isomo,mask)
+          Isomotmp = IAND(Isomo,mask)
           pmodel = POPCNT(mask) - POPCNT(XOR(Isomotmp,mask))
           mask = ISHFT(1_8,q) - 1
-        Isomotmp = AND(Isomo,mask)
+          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 = AND(Jsomo,mask)
+          Jsomotmp = IAND(Jsomo,mask)
           pmodel = POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask))
           mask = ISHFT(1_8,q) - 1
-        Jsomotmp = AND(Jsomo,mask)
+          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 = AND(Isomo,mask)
+          Isomo = IAND(Isomo,mask)
           pmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask))
           mask = ISHFT(1_8,q) - 1
-        Jsomo = AND(Jsomo,mask)
+          Jsomo = IAND(Jsomo,mask)
           qmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask))
        case (4)
           ! DOMO -> SOMO
@@ -480,14 +513,15 @@ subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmod
           !print *,"type -> DOMO -> SOMO"
           !Isomo = IEOR(Isomo,Jsomo)
           mask = ISHFT(1_8,p) - 1
-        Jsomo = AND(Jsomo,mask)
+          Jsomo = IAND(Jsomo,mask)
           pmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask))
           mask = ISHFT(1_8,q) - 1
-        Isomo = AND(Isomo,mask)
+          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,4)-1) - POPCNT(IAND(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,8)-1) - POPCNT(IAND(716,IBSET(0_8,8)-1))
 #+end_src
 
 #+RESULTS:

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) ]
 

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, intent(out)           :: pq_singles(2,*)
-  integer(bit_kind), intent(in)  :: cfg(Nint,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_hole = 1+n_core_orb, n_core_orb + n_act_orb
-    do i_particle = 1, mo_num
+    do i_particle = 1+n_core_orb, n_core_orb + n_act_orb
-      call do_single_excitation_cfg_with_type(cfg,single,i_hole,i_particle,ex_type,i_ok)
+      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

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