Merge branch 'csf' of github.com:QuantumPackage/qp2 into csf

RELEASE_NOTES.org

@@ -51,6 +51,7 @@
     - Added ~print_energy~
     - Added ~print_hamiltonian~
     - Added input for two body RDM
+    - Added keyword ~save_wf_after_selection~ 
 
 *** Code
 
@@ -81,3 +82,4 @@
       
 
 
+    

config/ifort_avx.cfg

@@ -32,7 +32,7 @@ OPENMP  : 1          ; Append OpenMP flags
 #
 [OPT]
 FC       : -traceback
-FCFLAGS  : -mavx -O2 -ip -ftz -g
+FCFLAGS  : -xAVX -O2 -ip -ftz -g
 
 # Profiling flags
 #################

docs/source/users_guide/quickstart.rst

@@ -128,6 +128,12 @@ and the atomic basis set:
   ao_two_e_erf_ints  density_for_dft       electrons     mo_two_e_ints      scf_utils
   ao_two_e_ints      determinants          ezfio         nuclei             work
 
+If you need to run using an already existing EZFIO database, use
+
+.. code::  bash
+
+  qp set_file hcn
+
 
 Run a Hartree-Fock calculation
 ------------------------------

src/bitmask/bitmasks_routines.irp.f

@@ -125,6 +125,41 @@ subroutine bitstring_to_str( output, string, Nint )
   output(ibuf:ibuf) = '|'
 end
 
+subroutine configuration_to_str( output, string, Nint )
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+! Transform the bit string of a configuration to a string for printing
+  END_DOC
+  character*(*), intent(out)     :: output
+  integer, intent(in)            :: Nint
+  integer(bit_kind), intent(in)  :: string(Nint,2)
+
+  integer                        :: i, j, ibuf
+  integer(bit_kind)              :: itemp
+
+  ibuf = 1
+  output = ''
+  output(ibuf:ibuf) = '|'
+  ibuf = ibuf+1
+  do i=1,Nint
+    itemp = 1_bit_kind
+    do j=1,bit_kind_size
+      if (iand(itemp,string(i,2)) == itemp) then
+        output(ibuf:ibuf) = '2'
+      else if (iand(itemp,string(i,1)) == itemp) then
+        output(ibuf:ibuf) = '1'
+      else
+        output(ibuf:ibuf) = '0'
+      endif
+      ibuf = ibuf+1
+      itemp = shiftl(itemp,1)
+    enddo
+  enddo
+  output(ibuf:ibuf) = '|'
+end
+
+
 
 subroutine bitstring_to_hexa( output, string, Nint )
   use bitmasks
@@ -166,6 +201,25 @@ subroutine debug_det(string,Nint)
 
 end
 
+subroutine debug_cfg(string,Nint)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+  ! Subroutine to print the content of a determinant in '+-' notation and
+  ! hexadecimal representation.
+  END_DOC
+  integer, intent(in)            :: Nint
+  integer(bit_kind), intent(in)  :: string(Nint,2)
+  character*(2048)                :: output(2)
+  call bitstring_to_hexa( output(1), string(1,1), Nint )
+  call bitstring_to_hexa( output(2), string(1,2), Nint )
+  print *,  trim(output(1)) , '|', trim(output(2))
+
+  call configuration_to_str( output(1), string, Nint )
+  print *,  trim(output(1))
+
+end
+
 subroutine print_det(string,Nint)
   use bitmasks
   implicit none

src/cipsi/EZFIO.cfg

@@ -4,6 +4,12 @@ doc: If true, computes the one- and two-body rdms with perturbation theory
 interface: ezfio,provider,ocaml
 default: False
 
+[save_wf_after_selection]
+type: logical
+doc: If true, saves the wave function after the selection, before the diagonalization
+interface: ezfio,provider,ocaml
+default: False
+
 [seniority_max]
 type: integer
 doc: Maximum number of allowed open shells. Using -1 selects all determinants

src/cipsi/cipsi.irp.f

@@ -114,7 +114,10 @@ subroutine run_cipsi
 
     ! Add selected determinants
     call copy_H_apply_buffer_to_wf()
-!    call save_wavefunction
+
+    if (save_wf_after_selection) then
+      call save_wavefunction
+    endif
 
     PROVIDE  psi_coef
     PROVIDE  psi_det

src/cipsi/selection.irp.f

@@ -166,6 +166,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
 
   integer                        :: l_a, nmax, idx
   integer, allocatable           :: indices(:), exc_degree(:), iorder(:)
+  double precision, parameter :: norm_thr = 1.d-16
   allocate (indices(N_det),                                          &
       exc_degree(max(N_det_alpha_unique,N_det_beta_unique)))
 
@@ -185,7 +186,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
       i = psi_bilinear_matrix_rows(l_a)
       if (nt + exc_degree(i) <= 4) then
         idx = psi_det_sorted_order(psi_bilinear_matrix_order(l_a))
-        if (psi_average_norm_contrib_sorted(idx) > 1.d-20) then
+        if (psi_average_norm_contrib_sorted(idx) > norm_thr) then
           indices(k) = idx
           k=k+1
         endif
@@ -212,7 +213,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
         idx = psi_det_sorted_order(                                  &
             psi_bilinear_matrix_order(                               &
             psi_bilinear_matrix_transp_order(l_a)))
-        if (psi_average_norm_contrib_sorted(idx) > 1.d-20) then
+        if (psi_average_norm_contrib_sorted(idx) > norm_thr) then
           indices(k) = idx
           k=k+1
         endif
@@ -742,12 +743,9 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
 
         alpha_h_psi = mat(istate, p1, p2)
 
-        do jstate=1,N_states
+        pt2_data % overlap(:,istate) = pt2_data % overlap(:,istate) + coef(:) * coef(istate)
-          pt2_data % overlap(jstate,istate) += coef(jstate) * coef(istate)
+        pt2_data % variance(istate)  = pt2_data % variance(istate) + alpha_h_psi * alpha_h_psi
-        enddo
+        pt2_data % pt2(istate)       = pt2_data % pt2(istate)      + e_pert(istate)
-
-        pt2_data % variance(istate) += alpha_h_psi * alpha_h_psi
-        pt2_data % pt2(istate)      += e_pert(istate)
 
 !!!DEBUG
 !        delta_E = E0(istate) - Hii + E_shift
@@ -1578,7 +1576,7 @@ subroutine get_d0_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
           phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int)
           hij = mo_two_e_integral(p1, p2, h1, h2) * phase
         end if
-        mat(:, p1, p2) += coefs(:) * hij
+        mat(:, p1, p2) = mat(:, p1, p2) + coefs(:) * hij
       end do
     end do
   else ! AA BB
@@ -1595,7 +1593,7 @@ subroutine get_d0_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
         else
           hij = (mo_two_e_integral(p1, p2, puti, putj) -  mo_two_e_integral(p2, p1, puti, putj))* get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2, N_int)
         end if
-        mat(:, puti, putj) += coefs(:) * hij
+        mat(:, puti, putj) = mat(:, puti, putj) + coefs(:) * hij
       end do
     end do
   end if
@@ -1654,18 +1652,18 @@ subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
       do putj=1, hfix-1
         if(lbanned(putj, ma) .or. banned(putj, puti,bant)) cycle
         hij = (mo_two_e_integral(p1, p2, putj, hfix)-mo_two_e_integral(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int)
-        tmp_row(1:N_states,putj) += hij * coefs(1:N_states)
+        tmp_row(1:N_states,putj) = tmp_row(1:N_states,putj) + hij * coefs(1:N_states)
       end do
       do putj=hfix+1, mo_num
         if(lbanned(putj, ma) .or. banned(putj, puti,bant)) cycle
         hij = (mo_two_e_integral(p1, p2, hfix, putj)-mo_two_e_integral(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int)
-        tmp_row(1:N_states,putj) += hij * coefs(1:N_states)
+        tmp_row(1:N_states,putj) = tmp_row(1:N_states,putj) + hij * coefs(1:N_states)
       end do
 
       if(ma == 1) then
-        mat(1:N_states,1:mo_num,puti) += tmp_row(1:N_states,1:mo_num)
+        mat(1:N_states,1:mo_num,puti) = mat(1:N_states,1:mo_num,puti) + tmp_row(1:N_states,1:mo_num)
       else
-        mat(1:N_states,puti,1:mo_num) += tmp_row(1:N_states,1:mo_num)
+        mat(1:N_states,puti,1:mo_num) = mat(1:N_states,puti,1:mo_num) + tmp_row(1:N_states,1:mo_num)
       end if
     end if
 
@@ -1679,22 +1677,22 @@ subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
       putj = p1
       if(.not. banned(putj,puti,bant)) then
         hij = mo_two_e_integral(p2,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix, N_int)
-        tmp_row(:,puti) += hij * coefs(:)
+        tmp_row(:,puti) = tmp_row(:,puti) + hij * coefs(:)
       end if
 
       putj = p2
       if(.not. banned(putj,puti,bant)) then
         hij = mo_two_e_integral(p1,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix, N_int)
-        tmp_row2(:,puti) += hij * coefs(:)
+        tmp_row2(:,puti) = tmp_row2(:,puti) + hij * coefs(:)
       end if
     end do
 
     if(mi == 1) then
-      mat(:,:,p1) += tmp_row(:,:)
+      mat(:,:,p1) = mat(:,:,p1) + tmp_row(:,:)
-      mat(:,:,p2) += tmp_row2(:,:)
+      mat(:,:,p2) = mat(:,:,p2) + tmp_row2(:,:)
     else
-      mat(:,p1,:) += tmp_row(:,:)
+      mat(:,p1,:) = mat(:,p1,:) + tmp_row(:,:)
-      mat(:,p2,:) += tmp_row2(:,:)
+      mat(:,p2,:) = mat(:,p2,:) + tmp_row2(:,:)
     end if
   else
     if(p(0,ma) == 3) then
@@ -1707,16 +1705,16 @@ subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
         do putj=1,hfix-1
           if(lbanned(putj,ma) .or. banned(puti,putj,1)) cycle
           hij = (mo_two_e_integral(p1, p2, putj, hfix)-mo_two_e_integral(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int)
-          tmp_row(:,putj) += hij * coefs(:)
+          tmp_row(:,putj) = tmp_row(:,putj) + hij * coefs(:)
         end do
         do putj=hfix+1,mo_num
           if(lbanned(putj,ma) .or. banned(puti,putj,1)) cycle
           hij = (mo_two_e_integral(p1, p2, hfix, putj)-mo_two_e_integral(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int)
-          tmp_row(:,putj) += hij * coefs(:)
+          tmp_row(:,putj) = tmp_row(:,putj) + hij * coefs(:)
         end do
 
-        mat(:, :puti-1, puti) += tmp_row(:,:puti-1)
+        mat(:, :puti-1, puti) = mat(:, :puti-1, puti) + tmp_row(:,:puti-1)
-        mat(:, puti, puti:) += tmp_row(:,puti:)
+        mat(:, puti, puti:) = mat(:, puti, puti:) + tmp_row(:,puti:)
       end do
     else
       hfix = h(1,mi)
@@ -1730,19 +1728,19 @@ subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
         putj = p2
         if(.not. banned(puti,putj,1)) then
           hij = mo_two_e_integral(pfix, p1, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1, N_int)
-          tmp_row(:,puti) += hij * coefs(:)
+          tmp_row(:,puti) = tmp_row(:,puti) + hij * coefs(:)
         end if
 
         putj = p1
         if(.not. banned(puti,putj,1)) then
           hij = mo_two_e_integral(pfix, p2, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2, N_int)
-          tmp_row2(:,puti) += hij * coefs(:)
+          tmp_row2(:,puti) = tmp_row2(:,puti) + hij * coefs(:)
         end if
       end do
-      mat(:,:p2-1,p2) += tmp_row(:,:p2-1)
+      mat(:,:p2-1,p2) = mat(:,:p2-1,p2) + tmp_row(:,:p2-1)
-      mat(:,p2,p2:) += tmp_row(:,p2:)
+      mat(:,p2,p2:) = mat(:,p2,p2:) + tmp_row(:,p2:)
-      mat(:,:p1-1,p1) += tmp_row2(:,:p1-1)
+      mat(:,:p1-1,p1) = mat(:,:p1-1,p1) + tmp_row2(:,:p1-1)
-      mat(:,p1,p1:) += tmp_row2(:,p1:)
+      mat(:,p1,p1:) = mat(:,p1,p1:) + tmp_row2(:,p1:)
     end if
   end if
   deallocate(lbanned)
@@ -1765,7 +1763,7 @@ subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
         if(bannedOrb(p1, s1) .or. bannedOrb(p2, s2) .or. banned(p1, p2, 1)) cycle
         call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int)
         call i_h_j(gen, det, N_int, hij)
-        mat(:, p1, p2) += coefs(:) * hij
+        mat(:, p1, p2) = mat(:, p1, p2) + coefs(:) * hij
       end do
     end do
 end
@@ -1818,9 +1816,9 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
 
         hij = (mo_two_e_integral(p1, p2, h1, h2) - mo_two_e_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int)
         if(ma == 1) then
-          mat(:, putj, puti) += coefs(:) * hij
+          mat(:, putj, puti) = mat(:, putj, puti) + coefs(:) * hij
         else
-          mat(:, puti, putj) += coefs(:) * hij
+          mat(:, puti, putj) = mat(:, puti, putj) + coefs(:) * hij
         end if
       end do
     else
@@ -1836,7 +1834,7 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
           p1 = p(turn2(i), 1)
 
           hij = mo_two_e_integral(p1, p2, h1, h2) * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2,N_int)
-          mat(:, puti, putj) += coefs(:) * hij
+          mat(:, puti, putj) = mat(:, puti, putj) + coefs(:) * hij
         end do
       end do
     end if
@@ -1856,7 +1854,7 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
         p1 = p(i1, ma)
         p2 = p(i2, ma)
         hij = (mo_two_e_integral(p1, p2, h1, h2) - mo_two_e_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2,N_int)
-        mat(:, puti, putj) += coefs(:) * hij
+        mat(:, puti, putj) = mat(:, puti, putj) + coefs(:) * hij
       end do
       end do
     else if(tip == 3) then
@@ -1870,7 +1868,7 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
         p2 = p(i, ma)
 
         hij = mo_two_e_integral(p1, p2, h1, h2) * get_phase_bi(phasemask, mi, ma, h1, p1, h2, p2,N_int)
-        mat(:, min(puti, putj), max(puti, putj)) += coefs(:) * hij
+        mat(:, min(puti, putj), max(puti, putj)) = mat(:, min(puti, putj), max(puti, putj)) + coefs(:) * hij
       end do
     else ! tip == 4
       puti = p(1, sp)
@@ -1881,7 +1879,7 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
         h1 = h(1, mi)
         h2 = h(2, mi)
         hij = (mo_two_e_integral(p1, p2, h1, h2) - mo_two_e_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2,N_int)
-        mat(:, puti, putj) += coefs(:) * hij
+        mat(:, puti, putj) = mat(:, puti, putj) + coefs(:) * hij
       end if
     end if
   end if

src/cipsi/stochastic_cipsi.irp.f

@@ -104,7 +104,9 @@ subroutine run_stochastic_cipsi
 
     ! Add selected determinants
     call copy_H_apply_buffer_to_wf()
-!    call save_wavefunction
+    if (save_wf_after_selection) then
+      call save_wavefunction
+    endif
 
     PROVIDE  psi_coef
     PROVIDE  psi_det

src/davidson/davidson_parallel.irp.f

@@ -428,7 +428,7 @@ subroutine H_S2_u_0_nstates_zmq(v_0,s_0,u_0,N_st,sze)
 
   integer :: istep, imin, imax, ishift, ipos
   integer, external :: add_task_to_taskserver
-  integer, parameter :: tasksize=10000
+  integer, parameter :: tasksize=20000
   character*(100000) :: task
   istep=1
   ishift=0

src/davidson/u0_h_u0.irp.f

@@ -211,6 +211,7 @@ subroutine H_S2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend,
   double precision               :: rss, mem, ratio
   double precision, allocatable  :: utl(:,:)
   integer, parameter             :: block_size=128
+  logical                        :: u_is_sparse
 
 !  call resident_memory(rss)
 !  mem = dble(singles_beta_csc_size) / 1024.d0**3
@@ -222,6 +223,7 @@ subroutine H_S2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend,
 !  endif
 compute_singles=.True.
 
+
   maxab = max(N_det_alpha_unique, N_det_beta_unique)+1
   allocate(idx0(maxab))
 
@@ -249,8 +251,8 @@ compute_singles=.True.
       !$OMP          singles_alpha_csc,singles_alpha_csc_idx,        &
       !$OMP          singles_beta_csc,singles_beta_csc_idx)          &
       !$OMP   PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i,     &
-      !$OMP          lcol, lrow, l_a, l_b, utl, kk,                  &
+      !$OMP          lcol, lrow, l_a, l_b, utl, kk, u_is_sparse,     &
-      !$OMP          buffer, doubles, n_doubles,                     &
+      !$OMP          buffer, doubles, n_doubles, umax,               &
       !$OMP          tmp_det2, hij, sij, idx, l, kcol_prev,          &
       !$OMP          singles_a, n_singles_a, singles_b, ratio,       &
       !$OMP          n_singles_b, k8, last_found,left,right,right_max)
@@ -266,6 +268,22 @@ compute_singles=.True.
 
   kcol_prev=-1
 
+  ! Check if u has multiple zeros
+  kk=1 ! Avoid division by zero
+  !$OMP DO
+  do k=1,N_det
+    umax = 0.d0
+    do l=1,N_st
+      umax = max(umax, dabs(u_t(l,k)))
+    enddo
+    if (umax < 1.d-20) then
+      !$OMP ATOMIC
+      kk = kk+1
+    endif
+  enddo
+  !$OMP END DO
+  u_is_sparse = N_det / kk < 20  ! 5%
+
   ASSERT (iend <= N_det)
   ASSERT (istart > 0)
   ASSERT (istep  > 0)
@@ -399,9 +417,13 @@ compute_singles=.True.
       ! Loop over alpha singles
       ! -----------------------
 
+      double precision :: umax
+
       !DIR$ LOOP COUNT avg(1000)
       do k = 1,n_singles_a,block_size
+        umax = 0.d0
         ! Prefetch u_t(:,l_a)
+        if (u_is_sparse) then
           do kk=0,block_size-1
             if (k+kk > n_singles_a) exit
             l_a = singles_a(k+kk)
@@ -409,8 +431,19 @@ compute_singles=.True.
 
             do l=1,N_st
               utl(l,kk+1) = u_t(l,l_a)
+              umax = max(umax, dabs(utl(l,kk+1)))
             enddo
           enddo
+        else
+          do kk=0,block_size-1
+            if (k+kk > n_singles_a) exit
+            l_a = singles_a(k+kk)
+            ASSERT (l_a <= N_det)
+            utl(:,kk+1) = u_t(:,l_a)
+          enddo
+          umax = 1.d0
+        endif
+        if (umax < 1.d-20) cycle
 
         do kk=0,block_size-1
           if (k+kk > n_singles_a) exit
@@ -490,7 +523,9 @@ compute_singles=.True.
     tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
     !DIR$ LOOP COUNT avg(1000)
     do i=1,n_singles_a,block_size
+      umax = 0.d0
       ! Prefetch u_t(:,l_a)
+      if (u_is_sparse) then
         do kk=0,block_size-1
           if (i+kk > n_singles_a) exit
           l_a = singles_a(i+kk)
@@ -498,8 +533,19 @@ compute_singles=.True.
 
           do l=1,N_st
             utl(l,kk+1) = u_t(l,l_a)
+            umax = max(umax, dabs(utl(l,kk+1)))
           enddo
         enddo
+      else
+        do kk=0,block_size-1
+          if (i+kk > n_singles_a) exit
+          l_a = singles_a(i+kk)
+          ASSERT (l_a <= N_det)
+          utl(:,kk+1) = u_t(:,l_a)
+        enddo
+        umax = 1.d0
+      endif
+      if (umax < 1.d-20) cycle
 
       do kk=0,block_size-1
         if (i+kk > n_singles_a) exit
@@ -524,7 +570,9 @@ compute_singles=.True.
 
     !DIR$ LOOP COUNT avg(50000)
     do i=1,n_doubles,block_size
+      umax = 0.d0
       ! Prefetch u_t(:,l_a)
+      if (u_is_sparse) then
         do kk=0,block_size-1
           if (i+kk > n_doubles) exit
           l_a = doubles(i+kk)
@@ -532,8 +580,19 @@ compute_singles=.True.
 
           do l=1,N_st
             utl(l,kk+1) = u_t(l,l_a)
+            umax = max(umax, dabs(utl(l,kk+1)))
           enddo
         enddo
+      else
+        do kk=0,block_size-1
+          if (i+kk > n_doubles) exit
+          l_a = doubles(i+kk)
+          ASSERT (l_a <= N_det)
+          utl(:,kk+1) = u_t(:,l_a)
+        enddo
+        umax = 1.d0
+      endif
+      if (umax < 1.d-20) cycle
 
       do kk=0,block_size-1
         if (i+kk > n_doubles) exit
@@ -599,18 +658,32 @@ compute_singles=.True.
     tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
     !DIR$ LOOP COUNT avg(1000)
     do i=1,n_singles_b,block_size
+      umax = 0.d0
+      if (u_is_sparse) then
         do kk=0,block_size-1
           if (i+kk > n_singles_b) exit
           l_b = singles_b(i+kk)
-        ASSERT (l_b <= N_det)
-
           l_a = psi_bilinear_matrix_transp_order(l_b)
+          ASSERT (l_b <= N_det)
           ASSERT (l_a <= N_det)
 
           do l=1,N_st
             utl(l,kk+1) = u_t(l,l_a)
+            umax = max(umax, dabs(utl(l,kk+1)))
           enddo
         enddo
+      else
+        do kk=0,block_size-1
+          if (i+kk > n_singles_b) exit
+          l_b = singles_b(i+kk)
+          l_a = psi_bilinear_matrix_transp_order(l_b)
+          ASSERT (l_b <= N_det)
+          ASSERT (l_a <= N_det)
+          utl(:,kk+1) = u_t(:,l_a)
+        enddo
+        umax = 1.d0
+      endif
+      if (umax < 1.d-20) cycle
 
       do kk=0,block_size-1
         if (i+kk > n_singles_b) exit
@@ -634,17 +707,31 @@ compute_singles=.True.
 
     !DIR$ LOOP COUNT avg(50000)
     do i=1,n_doubles,block_size
+      umax = 0.d0
+      if (u_is_sparse) then
         do kk=0,block_size-1
           if (i+kk > n_doubles) exit
           l_b = doubles(i+kk)
-        ASSERT (l_b <= N_det)
           l_a = psi_bilinear_matrix_transp_order(l_b)
+          ASSERT (l_b <= N_det)
           ASSERT (l_a <= N_det)
-
           do l=1,N_st
             utl(l,kk+1) = u_t(l,l_a)
+            umax = max(umax, dabs(utl(l,kk+1)))
           enddo
         enddo
+      else
+        do kk=0,block_size-1
+          if (i+kk > n_doubles) exit
+          l_b = doubles(i+kk)
+          l_a = psi_bilinear_matrix_transp_order(l_b)
+          ASSERT (l_b <= N_det)
+          ASSERT (l_a <= N_det)
+          utl(:,kk+1) = u_t(:,l_a)
+        enddo
+        umax = 1.d0
+      endif
+      if (umax < 1.d-20) cycle
 
       do kk=0,block_size-1
         if (i+kk > n_doubles) exit
@@ -671,6 +758,16 @@ compute_singles=.True.
     ! Initial determinant is at k_a in alpha-major representation
     ! -----------------------------------------------------------------------
 
+    if (u_is_sparse) then
+      umax = 0.d0
+      do l=1,N_st
+        umax = max(umax, dabs(u_t(l,k_a)))
+      enddo
+    else
+      umax = 1.d0
+    endif
+    if (umax < 1.d-20) cycle
+
     krow = psi_bilinear_matrix_rows(k_a)
     ASSERT (krow <= N_det_alpha_unique)
 

src/determinants/configurations.irp.f

@@ -42,6 +42,7 @@ subroutine configuration_to_dets_size(o,sze,n_alpha,Nint)
     amax -= popcnt( o(k,2) )
   enddo
   if (binom_int(bmax, amax) > huge(1)) then
+    print *, bmax, amax
     print *,  irp_here, ': Too many determinants to generate'
     stop 1
   endif
@@ -254,16 +255,13 @@ end
     endif
     dup = .True.
     do k=1,N_int
-      if ( (tmp_array(k,1,i) /= tmp_array(k,1,j)) &
+      dup = dup .and. (tmp_array(k,1,i) == tmp_array(k,1,j)) &
-      .or. (tmp_array(k,2,i) /= tmp_array(k,2,j)) ) then
+                .and. (tmp_array(k,2,i) == tmp_array(k,2,j))
-         dup = .False.
-         exit
-      endif
     enddo
     if (dup) then
       duplicate(j) = .True.
     endif
-    j+=1
+    j = j+1
     if (j>N_det) then
       exit
     endif
@@ -335,18 +333,20 @@ BEGIN_PROVIDER [ integer, cfg_seniority_index, (0:elec_num) ]
  enddo
 END_PROVIDER
 
-
 BEGIN_PROVIDER [ integer, det_to_configuration, (N_det) ]
+ implicit none
+ BEGIN_DOC
  ! Returns the index of the configuration for each determinant
  END_DOC
  integer :: i,j,k,r,l
- integer*8 :: key
+ integer*8 :: key, key2
  integer(bit_kind) :: occ(N_int,2)
  logical :: found
  integer*8, allocatable :: bit_tmp(:)
  integer*8, external            :: configuration_search_key
 
- allocate(bit_tmp(N_configuration))
+ 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)
  enddo
@@ -361,36 +361,37 @@ BEGIN_PROVIDER [ integer, det_to_configuration, (N_det) ]
 
     key = configuration_search_key(occ,N_int)
 
-    ! TODO: Binary search
+    ! Binary search
-    l = 1
+    l = 0
-    r = N_configuration
+    r = N_configuration+1
-!    do while(r-l > 32)
+    j = shiftr(r-l,1)
-!      j = shiftr(r+l,1)
+    do while (j>=1)
-!      if (bit_tmp(j) < key) then
+      j = j+l
-!        l = j
+      if (bit_tmp(j) == key) then
-!      else
+        do while (bit_tmp(j) == bit_tmp(j-1))
-!        r = j
+          j = j-1
-!      endif
+        enddo
-!    enddo
+        do while (bit_tmp(j) == key)
-    do j=l,r
           found = .True.
           do k=1,N_int
-        if ( (occ(k,1) /= psi_configuration(k,1,j)) &
+            found = found .and. (psi_configuration(k,1,j) == occ(k,1)) &
-        .or. (occ(k,2) /= psi_configuration(k,2,j)) ) then
+                          .and. (psi_configuration(k,2,j) == occ(k,2))
-          found = .False.
-          exit
-        endif
           enddo
           if (found) then
             det_to_configuration(i) = j
             exit
           endif
+          j = j+1
+        enddo
+        if (found) exit
+      else if (bit_tmp(j) > key) then
+        r = j
+      else
+        l = j
+      endif
+      j = shiftr(r-l,1)
     enddo
 
-    if (.not.found) then
-      print *,  '3 bug in ',  irp_here
-      stop -1
-    endif
  enddo
  !$OMP END PARALLEL DO
  deallocate(bit_tmp)
@@ -584,3 +585,39 @@ BEGIN_PROVIDER [ integer(bit_kind), dominant_dets_of_cfgs, (N_int,2,N_dominant_d
    i += sze
  enddo
 END_PROVIDER
+
+ BEGIN_PROVIDER [ integer, psi_configuration_to_psi_det, (2,N_configuration) ]
+&BEGIN_PROVIDER [ integer, psi_configuration_to_psi_det_data, (N_det) ]
+
+ implicit none
+ BEGIN_DOC
+ ! psi_configuration_to_psi_det_data(k) -> i : i is the index of the
+ ! determinant in psi_det
+ !
+ ! psi_configuration_to_psi_det(1:2,k) gives the first and last index of the
+ ! determinants of configuration k in array psi_configuration_to_psi_det_data.
+ END_DOC
+
+ integer :: i, k, iorder
+ integer, allocatable :: confs(:)
+ allocate (confs(N_det))
+
+ do i=1,N_det
+   psi_configuration_to_psi_det_data(i) = i
+   confs(i) = det_to_configuration(i)
+ enddo
+
+ call isort(confs, psi_configuration_to_psi_det_data, N_det)
+ k=1
+ psi_configuration_to_psi_det(1,1) = 1
+ do i=2,N_det
+   if (confs(i) /= confs(i-1)) then
+     psi_configuration_to_psi_det(2,k) = i-1
+     k = k+1
+     psi_configuration_to_psi_det(1,k) = i
+   endif
+ enddo
+ psi_configuration_to_psi_det(2,k) = N_det
+
+END_PROVIDER
+

src/determinants/configurations.org

@@ -0,0 +1,120 @@
+# -*- mode:org -*-
+#+TITLE: CFG-CI
+#+AUTHOR: Vijay Gopal Chilkuri
+#+FILE: configurations.org
+#+EMAIL: vijay.gopal.c@gmail.com
+#+OPTIONS: toc:t
+#+LATEX_CLASS: article
+#+LATEX_HEADER: \usepackage{tabularx}
+#+LATEX_HEADER: \usepackage{braket}
+#+LATEX_HEADER: \usepackage{minted}
+
+* Configuration based CI
+
+Here we write the main functions that perform the functions necessary for
+the Configuration based CI.
+
+There are three main functions required for doing the CI
+
+- Convert wavefunction from determinant basis to configuration state function (CSF) basis
+
+- Apply the Hamiltonian to the wavefunction in CSF basis
+
+- Convert the converged wavefunction back to determinant basis
+
+** TODO[0/3] Convert basis from DET to CSF
+
+The conversion of basis is done by going via bonded functions (BFs).
+Importantly, all the CSFs of a chosen configuration (CFG) are kept.
+
+The advantage is that the sigma-vector can be performed efficiently
+via BLAS level 3 operations.
+
+
+- [ ] Write a function to calculate the maximum dimensions required
+
+  Prototype array contains the \( <I|\hat{E}_{pq}|J> \) for all possible
+  CFGs \( I, J\) and all \(4\) types of excitations for all possible model
+  orbitals \(p,q\). Note that the orbital ids \(p,q\) here do not refer to
+  the actual MO ids, they simply refer to the orbitals involved in that spefcific
+  SOMO, for e.g. an excitation of the type [2 2 2 1 1 1 1 0] -> [ 2 2 1 1 1 1 1]
+  implies an excitation from orbital \(3\) to orbital \(8\) which are the real MO ids.
+  However, the prototype only concerns the SOMOs like so [2 1 1 1 1 0] -> [ 1 1 1 1 1 1]
+  therefore \(p,q\) are model space ids \(1,6\).
+
+  #+begin_src f90 :main no :tangle configurations_sigma_vector.irp.f
+
+  BEGIN_PROVIDER [ integer, NSOMOMax]
+ &BEGIN_PROVIDER [ integer, NCSFMax]
+  implicit none
+  BEGIN_DOC
+  ! Documentation for NSOMOMax
+  ! The maximum number of SOMOs for the current calculation.
+  ! required for the calculation of prototype arrays.
+  END_DOC
+  NSOMOMax = 8
+  NCSFMax = 14 ! TODO: NCSFs for MS=0
+  END_PROVIDER
+  #+end_src
+
+  The prototype matrix AIJpqMatrixList has the following dimensions
+  \(\left(NSOMOMax, NSOMOMax, 4, NSOMOMax, NSOMOMax,NCSFMAx,NCSFMax\right)\) where the first two
+  indices represent the somos in \(I,J\) followed by the type of excitation and
+  finally the two model space orbitals \(p,q\).
+
+- [ ] Read the transformation matrix based on the number of SOMOs
+
+  #+begin_src f90 :main no :tangle configurations_sigma_vector.irp.f
+  BEGIN_PROVIDER [ double precision, AIJpqMatrixDimsList, (NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,2)]
+  use cfunctions
+  implicit none
+  BEGIN_DOC
+  ! Documentation for AIJpqMatrixList
+  ! The prototype matrix containing the <I|E_{pq}|J>
+  ! matrices for each I,J somo pair and orb ids.
+  END_DOC
+  integer i,j,k,l
+  integer*8 Isomo, Jsomo
+  Isomo = 0
+  Jsomo = 0
+  integer*8 rows, cols
+  rows = -1
+  cols = -1
+  integer*8 MS
+  MS = 0
+  print *,"NSOMOMax = ",NSOMOMax
+  !allocate(AIJpqMatrixDimsList(NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,2))
+  do i = 2, NSOMOMax, 2
+     Isomo = ISHFT(1,i)-1
+     do j = i-2,i+2, 2
+        Jsomo = ISHFT(1,j)-1
+        if(j .GT. NSOMOMax .OR. j .LE. 0) then
+           cycle
+        end if
+        do k = 1,NSOMOMax
+           do l = k,NSOMOMax
+              call getApqIJMatrixDims(Isomo,           &
+                   Jsomo, &
+                   MS,                       &
+                   rows,                     &
+                   cols)
+              print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols
+              ! i -> j
+              AIJpqMatrixDimsList(i,j,1,k,l,1) = rows
+              AIJpqMatrixDimsList(i,j,1,k,l,2) = cols
+              AIJpqMatrixDimsList(i,j,1,l,k,1) = rows
+              AIJpqMatrixDimsList(i,j,1,l,k,2) = cols
+              ! j -> i
+              AIJpqMatrixDimsList(j,i,1,k,l,1) = rows
+              AIJpqMatrixDimsList(j,i,1,k,l,2) = cols
+              AIJpqMatrixDimsList(j,i,1,l,k,1) = rows
+              AIJpqMatrixDimsList(j,i,1,l,k,2) = cols
+           end do
+        end do
+     end do
+  end do
+  END_PROVIDER
+
+  #+end_src
+
+- [ ] Perform the conversion by matrix-vector BLAS level 2 call

src/determinants/configurations_sigma_vector.irp.f

@@ -0,0 +1,61 @@
+  BEGIN_PROVIDER [ integer, NSOMOMax]
+ &BEGIN_PROVIDER [ integer, NCSFMax]
+  implicit none
+  BEGIN_DOC
+  ! Documentation for NSOMOMax
+  ! The maximum number of SOMOs for the current calculation.
+  ! required for the calculation of prototype arrays.
+  END_DOC
+  NSOMOMax = 8
+  NCSFMax = 14 ! TODO: NCSFs for MS=0
+  END_PROVIDER
+
+  BEGIN_PROVIDER [ double precision, AIJpqMatrixDimsList, (NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,2)]
+  use cfunctions
+  implicit none
+  BEGIN_DOC
+  ! Documentation for AIJpqMatrixList
+  ! The prototype matrix containing the <I|E_{pq}|J>
+  ! matrices for each I,J somo pair and orb ids.
+  END_DOC
+  integer i,j,k,l
+  integer*8 Isomo, Jsomo
+  Isomo = 0
+  Jsomo = 0
+  integer*8 rows, cols
+  rows = -1
+  cols = -1
+  integer*8 MS
+  MS = 0
+  print *,"NSOMOMax = ",NSOMOMax
+  !allocate(AIJpqMatrixDimsList(NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,2))
+  do i = 2, NSOMOMax, 2
+     Isomo = ISHFT(1,i)-1
+     do j = i-2,i+2, 2
+        Jsomo = ISHFT(1,j)-1
+        if(j .GT. NSOMOMax .OR. j .LE. 0) then
+           cycle
+        end if
+        do k = 1,NSOMOMax
+           do l = k,NSOMOMax
+              call getApqIJMatrixDims(Isomo,           &
+                   Jsomo, &
+                   MS,                       &
+                   rows,                     &
+                   cols)
+              print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols
+              ! i -> j
+              AIJpqMatrixDimsList(i,j,1,k,l,1) = rows
+              AIJpqMatrixDimsList(i,j,1,k,l,2) = cols
+              AIJpqMatrixDimsList(i,j,1,l,k,1) = rows
+              AIJpqMatrixDimsList(i,j,1,l,k,2) = cols
+              ! j -> i
+              AIJpqMatrixDimsList(j,i,1,k,l,1) = rows
+              AIJpqMatrixDimsList(j,i,1,k,l,2) = cols
+              AIJpqMatrixDimsList(j,i,1,l,k,1) = rows
+              AIJpqMatrixDimsList(j,i,1,l,k,2) = cols
+           end do
+        end do
+     end do
+  end do
+  END_PROVIDER

src/determinants/create_excitations.irp.f

@@ -99,7 +99,8 @@ logical function is_spin_flip_possible(key_in,i_flip,ispin)
   other_spin(1) = 2
   other_spin(2) = 1
   if(popcnt(iand(key_tmp(k,1),key_in(k,ispin))) == 1 .and. popcnt(iand(key_tmp(k,1),key_in(k,other_spin(ispin)))) == 0  )then
-    ! There is a spin "ispin" in the orbital i_flip   AND  There is no electron of opposit spin in the same orbital "i_flip"
+    ! There is a spin "ispin" in the orbital i_flip   AND
+    ! There is no electron of opposit spin in the same orbital "i_flip"
     is_spin_flip_possible = .True.
     return
   else
@@ -107,3 +108,115 @@ logical function is_spin_flip_possible(key_in,i_flip,ispin)
   endif
 end
 
+subroutine do_single_excitation_cfg(key_in,key_out,i_hole,i_particle,ok)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+  ! Applies the single excitation operator to a configuration
+  ! If the excitation is possible, ok is True
+  END_DOC
+  integer, intent(in)            :: i_hole,i_particle
+  integer(bit_kind), intent(in)  :: key_in(N_int,2)
+  logical , intent(out)          :: ok
+  integer                        :: k,j,i
+  integer(bit_kind)              :: mask
+  integer(bit_kind)              :: key_out(N_int,2)
+
+  ASSERT (i_hole > 0)
+  ASSERT (i_particle <= mo_num)
+
+  ok = .True.
+  key_out(:,:) = key_in(:,:)
+
+  ! hole
+  k = shiftr(i_hole-1,bit_kind_shift)+1
+  j = i_hole-shiftl(k-1,bit_kind_shift)-1
+  mask = ibset(0_bit_kind,j)
+
+  ! Check if the position j is singly occupied
+  ! 1  ->  0  (SOMO)
+  ! 0      0  (DOMO)
+  if (iand(key_out(k,1),mask) /= 0_bit_kind) then
+     key_out(k,1) = ibclr(key_out(k,1),j)
+
+  ! Check if the position j is doubly occupied
+  ! 0  ->  1  (SOMO)
+  ! 1      0  (DOMO)
+  else if (iand(key_out(k,2),mask) /= 0_bit_kind) then
+     key_out(k,1) = ibset(key_out(k,1),j)
+     key_out(k,2) = ibclr(key_out(k,2),j)
+
+  ! The position j is unoccupied: Not OK
+  ! 0  ->  0  (SOMO)
+  ! 0      0  (DOMO)
+  else
+     ok =.False.
+     return
+  endif
+
+
+  ! particle
+  k = shiftr(i_particle-1,bit_kind_shift)+1
+  j = i_particle-shiftl(k-1,bit_kind_shift)-1
+  mask = ibset(0_bit_kind,j)
+
+  ! Check if the position j is singly occupied
+  ! 1  ->  0  (SOMO)
+  ! 0      1  (DOMO)
+  if (iand(key_out(k,1),mask) /= 0_bit_kind) then
+     key_out(k,1) = ibclr(key_out(k,1),j)
+     key_out(k,2) = ibset(key_out(k,2),j)
+
+  ! Check if the position j is doubly occupied : Not OK
+  ! 0  ->  1  (SOMO)
+  ! 1      0  (DOMO)
+  else if (iand(key_out(k,2),mask) /= 0_bit_kind) then
+     ok = .False.
+     return
+
+  ! Position at j is unoccupied
+  ! 0  ->  0  (SOMO)
+  ! 0      0  (DOMO)
+  else
+     key_out(k,1) = ibset(key_out(k,1),j)
+  endif
+
+end
+
+subroutine generate_all_singles_cfg(cfg,singles,n_singles,Nint)
+  implicit none
+  use bitmasks
+  BEGIN_DOC
+  ! Generate all single excitation wrt a configuration
+  !
+  ! n_singles : on input, max number of singles :
+  !   elec_alpha_num * (mo_num - elec_beta_num)
+  !   on output, number of generated singles
+  END_DOC
+  integer, intent(in)            :: Nint
+  integer, intent(inout)         :: n_singles
+  integer(bit_kind), intent(in)  :: cfg(Nint,2)
+  integer(bit_kind), intent(out) :: singles(Nint,2,*)
+
+  integer           :: i,k, n_singles_ma, i_hole, i_particle
+  integer(bit_kind) :: single(Nint,2)
+  logical           :: i_ok
+
+  n_singles = 0
+  !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(cfg,single,i_hole,i_particle,i_ok)
+      if (i_ok) then
+        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)
+        enddo
+      endif
+    enddo
+  enddo
+end
+

src/mo_two_e_ints/map_integrals.irp.f

@@ -99,10 +99,16 @@ double precision function get_two_e_integral(i,j,k,l,map)
   type(map_type), intent(inout)  :: map
   real(integral_kind)            :: tmp
   PROVIDE mo_two_e_integrals_in_map mo_integrals_cache
-  if (banned_excitation(i,k) .or. banned_excitation(j,l)) then
+  if (use_banned_excitation) then
+    if (banned_excitation(i,k)) then
       get_two_e_integral = 0.d0
       return
     endif
+    if (banned_excitation(j,l)) then
+      get_two_e_integral = 0.d0
+      return
+    endif
+  endif
   ii = l-mo_integrals_cache_min
   ii = ior(ii, k-mo_integrals_cache_min)
   ii = ior(ii, j-mo_integrals_cache_min)
@@ -283,16 +289,18 @@ subroutine get_mo_two_e_integrals_exch_ii(k,l,sze,out_val,map)
 end
 
  BEGIN_PROVIDER [ logical, banned_excitation, (mo_num,mo_num) ]
+&BEGIN_PROVIDER [ logical, use_banned_excitation  ]
  implicit none
  use map_module
  BEGIN_DOC
  ! If true, the excitation is banned in the selection. Useful with local MOs.
  END_DOC
  banned_excitation = .False.
- integer :: i,j
+ integer :: i,j, icount
  integer(key_kind)              :: idx
  double precision :: tmp
-! double precision :: buffer(mo_num)
+
+ icount = 1 ! Avoid division by zero
  do j=1,mo_num
    do i=1,j-1
     call two_e_integrals_index(i,j,j,i,idx)
@@ -300,8 +308,14 @@ BEGIN_PROVIDER [ logical, banned_excitation, (mo_num,mo_num) ]
     call map_get(mo_integrals_map,idx,tmp)
     banned_excitation(i,j) = dabs(tmp) < 1.d-14
     banned_excitation(j,i) = banned_excitation(i,j)
+    if (banned_excitation(i,j)) icount = icount+2
   enddo
  enddo
+ use_banned_excitation =  (mo_num*mo_num) / icount <= 100  !1%
+ if (use_banned_excitation) then
+   print *, 'Using sparsity of exchange integrals'
+ endif
+
 END_PROVIDER
 
 

src/two_body_rdm/example.irp.f

@@ -284,3 +284,4 @@ subroutine routine_full_mos
    print*,'wee_tot_st_av_3         = ',wee_tot_st_av_3
 
 end
+

travis/compilation.sh

@@ -8,7 +8,7 @@ tar -zxf $HOME/cache/config.tgz
 # Configure QP2
 cd qp2
 source ./quantum_package.rc
-ninja -j 1 -v
+ninja -j 1 -v || exit -1
 
 # Create cache
 cd ..

travis/configuration.sh

@@ -2,7 +2,7 @@
 # Stage 1
 
 # Configure QP2
-./configure --install all --config ./config/travis.cfg
+./configure --install all --config ./config/travis.cfg || exit -1
 
 # Create cache
 cd ../

travis/testing.sh

@@ -8,7 +8,7 @@ tar -zxf $HOME/cache/compil.tgz
 # Configure QP2
 cd qp2
 source ./quantum_package.rc
-qp_test -a && rm $HOME/cache/compil.tgz
+exec qp_test -a && rm $HOME/cache/compil.tgz