Merge branch 'master' of github.com:LCPQ/quantum_package

scripts/generate_h_apply.py

@@ -17,6 +17,8 @@ copy_buffer
 finalization
 generate_psi_guess
 init_thread
+printout_now
+printout_always
 deinit_thread
 """.split()
 
@@ -84,6 +86,8 @@ class H_apply(object):
         s[k] = ""
     s["size_max"] = str(1024*128) 
     s["copy_buffer"] = "call copy_h_apply_buffer_to_wf"
+    s["printout_now"]   = """write(output_Dets,*)  &
+       100.*float(i_generator)/float(N_det_generators), '% in ', wall_2-wall_1, 's'"""
     self.data = s
 
   def __setitem__(self,key,value):
@@ -154,14 +158,38 @@ class H_apply(object):
   double precision, intent(inout):: pt2(N_st) 
   double precision, intent(inout):: norm_pert(N_st) 
   double precision, intent(inout):: H_pert_diag(N_st)
+  double precision               :: delta_pt2(N_st), norm_psi(N_st), pt2_old(N_st)
   PROVIDE CI_electronic_energy N_det_generators key_pattern_not_in_ref                                                                                                      
   do k=1,N_st
     pt2(k) = 0.d0
     norm_pert(k) = 0.d0
     H_pert_diag(k) = 0.d0
+    norm_psi(k) = 0.d0
+    delta_pt2(k) = 0.d0
+    pt2_old(k) = 0.d0
   enddo
+        write(output_Dets,'(A12, X, A8, 3(2X, A9), 2X, A8, 2X, A8, 2X, A8)') &
+                 'N_generators', 'Norm', 'Delta PT2', 'PT2', 'Est. PT2', 'secs'
+        write(output_Dets,'(A12, X, A8, 3(2X, A9), 2X, A8, 2X, A8, 2X, A8)') &
+                 '============', '========', '=========', '=========', '=========', &
+                 '========='
       """ 
 
+      self.data["printout_always"] = """
+      do k=1,N_st
+        norm_psi(k) = norm_psi(k) + psi_coef(i_generator,k)*psi_coef(i_generator,k)
+        delta_pt2(k) = pt2(k) - pt2_old(k)
+      enddo
+      """
+      self.data["printout_now"] = """
+      do k=1,N_st
+        write(output_Dets,'(I10, 4(2X, F9.6), 2X, F8.1)') &
+                 i_generator, norm_psi(k), delta_pt2(k), pt2(k), &
+                 pt2(k)/norm_psi(k), &
+                 wall_2-wall_1
+         pt2_old(k) = pt2(k)
+      enddo
+      """
       if self.openmp:
         self.data["omp_parallel"]    += """&
  !$OMP SHARED(N_st) PRIVATE(e_2_pert_buffer,coef_pert_buffer) &

src/CISD/cisd.irp.f

@@ -5,6 +5,9 @@ program cisd
   print *,   'HF      = ', HF_energy
   print *,  'N_states = ', N_states
   call H_apply_cisd
+! do i=1,N_det
+!   print '(100(X,O32))',  det_connections(:,i)
+! enddo
   print *,  'N_det = ', N_det
   do i = 1,N_states
    print *,  'energy  = ',CI_energy(i) 

src/Dets/H_apply_template.f

@@ -214,7 +214,6 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene
       occ_hole_tmp)
   $omp_end_parallel
   $finalization
-  abort_here = abort_all
 end
 
 subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator $parameters )
@@ -344,7 +343,9 @@ subroutine $subroutine($params_main)
   
   PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map N_det_reference psi_generators
   integer                        :: i_generator, k
+  double precision               :: wall_0, wall_1, wall_2
   
+  call wall_time(wall_1)
   do i_generator=1,N_det_generators
     call $subroutine_diexc(psi_generators(1,1,i_generator),          &
         generators_bitmask(1,1,d_hole1,i_bitmask_gen),               &
@@ -359,6 +360,12 @@ subroutine $subroutine($params_main)
     if (abort_here) then
       exit
     endif
+    call wall_time(wall_2)
+    $printout_always
+    if (wall_2 - wall_0 > 2.d0) then
+        wall_0 = wall_2
+        $printout_now
+    endif
   enddo
   
   $copy_buffer

src/Dets/README.rst

@@ -77,10 +77,10 @@ Documentation
 `key_pattern_not_in_ref <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/connected_to_ref.irp.f#L222>`_
   Min and max values of the integers of the keys of the reference
 
-`davidson_converged <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L381>`_
+`davidson_converged <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L383>`_
   True if the Davidson algorithm is converged
 
-`davidson_criterion <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L371>`_
+`davidson_criterion <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L373>`_
   Can be : [  energy  | residual | both | wall_time | cpu_time | iterations ]
 
 `davidson_diag <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L18>`_
@@ -127,7 +127,7 @@ Documentation
 `davidson_sze_max <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L9>`_
   Max number of Davidson sizes
 
-`davidson_threshold <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L372>`_
+`davidson_threshold <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L374>`_
   Can be : [  energy  | residual | both | wall_time | cpu_time | iterations ]
 
 `n_det <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L20>`_
@@ -207,7 +207,9 @@ Documentation
   .br
   idx(0) is the number of determinants that interact with key1
 
-`filter_connected_i_h_psi0 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/filter_connected.irp.f#L94>`_
+`filter_connected_davidson <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/filter_connected.irp.f#L101>`_
+  Filters out the determinants that are not connected by H
+  .br
   returns the array idx which contains the index of the
   .br
   determinants in the array key1 that interact
@@ -216,7 +218,16 @@ Documentation
   .br
   idx(0) is the number of determinants that interact with key1
 
-`filter_connected_i_h_psi0_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/filter_connected.irp.f#L193>`_
+`filter_connected_i_h_psi0 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/filter_connected.irp.f#L233>`_
+  returns the array idx which contains the index of the
+  .br
+  determinants in the array key1 that interact
+  .br
+  via the H operator with key2.
+  .br
+  idx(0) is the number of determinants that interact with key1
+
+`filter_connected_i_h_psi0_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/filter_connected.irp.f#L332>`_
   standard filter_connected_i_H_psi but returns in addition
   .br
   the array of the index of the non connected determinants to key1
@@ -252,6 +263,9 @@ Documentation
   s1,s2 : Spins (1:alpha, 2:beta)
   degree : Degree of excitation
 
+`det_connections <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L898>`_
+  .br
+
 `diag_h_mat_elem <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L659>`_
   Computes <i|H|i>
 
@@ -299,6 +313,9 @@ Documentation
   .br
   to repeat the excitations
 
+`n_con_int <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L890>`_
+  Number of integers to represent the connections between determinants
+
 `h_matrix_all_dets <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/utils.irp.f#L1>`_
   H matrix on the basis of the slater deter;inants defined by psi_det
 

src/Dets/davidson.irp.f

@@ -114,6 +114,8 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun
   double precision               :: to_print(2,N_st)
   double precision               :: cpu, wall
   
+  PROVIDE N_con_int det_connections
+
   call write_time(iunit)
   call wall_time(wall)
   call cpu_time(cpu)

src/Dets/filter_connected.irp.f

@@ -32,7 +32,9 @@ subroutine filter_connected(key1,key2,Nint,sze,idx)
     do i=1,sze
       degree_x2 = popcnt(    xor( key1(1,1,i), key2(1,1))) &
                 + popcnt(    xor( key1(1,2,i), key2(1,2)))
-      if (degree_x2 < 5) then
+      if (degree_x2 > 4) then
+        cycle
+      else
         idx(l) = i
         l = l+1
       endif
@@ -46,7 +48,9 @@ subroutine filter_connected(key1,key2,Nint,sze,idx)
           popcnt(xor( key1(2,1,i), key2(2,1))) +                     &
           popcnt(xor( key1(1,2,i), key2(1,2))) +                     &
           popcnt(xor( key1(2,2,i), key2(2,2)))
-      if (degree_x2 < 5) then
+      if (degree_x2 > 4) then
+        cycle
+      else
         idx(l) = i
         l = l+1
       endif
@@ -62,7 +66,9 @@ subroutine filter_connected(key1,key2,Nint,sze,idx)
           popcnt(xor( key1(2,2,i), key2(2,2))) +                     &
           popcnt(xor( key1(3,1,i), key2(3,1))) +                     &
           popcnt(xor( key1(3,2,i), key2(3,2)))
-      if (degree_x2 < 5) then
+      if (degree_x2 > 4) then
+        cycle
+      else
         idx(l) = i
         l = l+1
       endif
@@ -91,6 +97,139 @@ subroutine filter_connected(key1,key2,Nint,sze,idx)
   idx(0) = l-1
 end
 
+
+subroutine filter_connected_davidson(key1,key2,Nint,sze,idx)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+  ! Filters out the determinants that are not connected by H
+  !
+  ! returns the array idx which contains the index of the 
+  !
+  ! determinants in the array key1 that interact 
+  !
+  ! via the H operator with key2.
+  !
+  ! idx(0) is the number of determinants that interact with key1
+  END_DOC
+  integer, intent(in)            :: Nint, sze
+  integer(bit_kind), intent(in)  :: key1(Nint,2,sze)
+  integer(bit_kind), intent(in)  :: key2(Nint,2)
+  integer, intent(out)           :: idx(0:sze)
+  
+  integer                        :: i,j,k,l
+  integer                        :: degree_x2
+  integer :: j_int, j_start
+  integer*8 :: itmp
+  
+  PROVIDE N_con_int det_connections
+  ASSERT (Nint > 0)
+  ASSERT (sze >= 0)
+
+  l=1
+  
+  if (Nint==1) then
+
+    i = idx(0)
+    do j_int=1,N_con_int
+      itmp = det_connections(j_int,i)
+      do while (itmp /= 0_8)
+        j_start = ishft(j_int-1,11) + ishft(trailz(itmp),5)
+        do j = j_start+1, min(j_start+32,i-1)
+          degree_x2 = popcnt(xor( key1(1,1,j), key2(1,1))) +             &
+              popcnt(xor( key1(1,2,j), key2(1,2)))
+          if (degree_x2 > 4) then
+            cycle
+          else
+            idx(l) = j
+            l = l+1
+          endif
+        enddo
+        itmp = iand(itmp-1_8,itmp)
+      enddo
+    enddo
+    
+  else if (Nint==2) then
+    
+    
+    i = idx(0)
+    do j_int=1,N_con_int
+      itmp = det_connections(j_int,i)
+      do while (itmp /= 0_8)
+        j_start = ishft(j_int-1,11) + ishft(trailz(itmp),5)
+        do j = j_start+1, min(j_start+32,i-1)
+          degree_x2 = popcnt(xor( key1(1,1,j), key2(1,1))) +         &
+              popcnt(xor( key1(2,1,j), key2(2,1))) +                 &
+              popcnt(xor( key1(1,2,j), key2(1,2))) +                 &
+              popcnt(xor( key1(2,2,j), key2(2,2)))
+          if (degree_x2 > 4) then
+            cycle
+          else
+            idx(l) = j
+            l = l+1
+          endif
+        enddo
+        itmp = iand(itmp-1_8,itmp)
+      enddo
+    enddo
+    
+  else if (Nint==3) then
+    
+    !DIR$ LOOP COUNT (1000)
+    i = idx(0)
+    do j_int=1,N_con_int
+      itmp = det_connections(j_int,i)
+      do while (itmp /= 0_8)
+        j_start = ishft(j_int-1,11) + ishft(trailz(itmp),5)
+        do j = j_start+1, min(j_start+32,i-1)
+          degree_x2 = popcnt(xor( key1(1,1,j), key2(1,1))) +         &
+              popcnt(xor( key1(1,2,j), key2(1,2))) +                 &
+              popcnt(xor( key1(2,1,j), key2(2,1))) +                 &
+              popcnt(xor( key1(2,2,j), key2(2,2))) +                 &
+              popcnt(xor( key1(3,1,j), key2(3,1))) +                 &
+              popcnt(xor( key1(3,2,j), key2(3,2)))
+          if (degree_x2 > 4) then
+            cycle
+          else
+            idx(l) = j
+            l = l+1
+          endif
+        enddo
+        itmp = iand(itmp-1_8,itmp)
+      enddo
+    enddo
+    
+  else
+    
+    !DIR$ LOOP COUNT (1000)
+    i = idx(0)
+    do j_int=1,N_con_int
+      itmp = det_connections(j_int,i)
+      do while (itmp /= 0_8)
+        j_start = ishft(j_int-1,11) + ishft(trailz(itmp),5)
+        do j = j_start+1, min(j_start+32,i-1)
+          degree_x2 = 0
+          !DEC$ LOOP COUNT MIN(4)
+          do k=1,Nint
+            degree_x2 = degree_x2+ popcnt(xor( key1(k,1,j), key2(k,1))) +&
+                popcnt(xor( key1(k,2,j), key2(k,2)))
+            if (degree_x2 > 4) then
+              exit
+            endif
+          enddo
+          if (degree_x2 <= 5) then
+            idx(l) = j
+            l = l+1
+          endif
+        enddo
+        itmp = iand(itmp-1_8,itmp)
+      enddo
+    enddo
+
+  endif
+  idx(0) = l-1
+end
+
 subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx)
   use bitmasks
   BEGIN_DOC
@@ -123,12 +262,12 @@ subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx)
     do i=1,sze
       degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) +             &
           popcnt(xor( key1(1,2,i), key2(1,2)))
-      if (degree_x2 < 5) then
+      if (degree_x2 > 4) then
-        if(degree_x2 .ne. 0)then
+        cycle
+      else if(degree_x2 .ne. 0)then
         idx(l) = i
         l = l+1
       endif
-      endif
     enddo
     
   else if (Nint==2) then
@@ -139,12 +278,12 @@ subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx)
           popcnt(xor( key1(2,1,i), key2(2,1))) +                     &
           popcnt(xor( key1(1,2,i), key2(1,2))) +                     &
           popcnt(xor( key1(2,2,i), key2(2,2)))
-      if (degree_x2 < 5) then
+      if (degree_x2 > 4) then
-        if(degree_x2 .ne. 0)then
+        cycle
+      else if(degree_x2 .ne. 0)then
         idx(l) = i
         l = l+1
       endif
-      endif
     enddo
     
   else if (Nint==3) then
@@ -157,12 +296,12 @@ subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx)
           popcnt(xor( key1(2,2,i), key2(2,2))) +                     &
           popcnt(xor( key1(3,1,i), key2(3,1))) +                     &
           popcnt(xor( key1(3,2,i), key2(3,2)))
-      if (degree_x2 < 5) then
+      if (degree_x2 > 4) then
-        if(degree_x2 .ne. 0)then
+        cycle
+      else if(degree_x2 .ne. 0)then
         idx(l) = i
         l = l+1
       endif
-      endif
     enddo
     
   else
@@ -178,12 +317,12 @@ subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx)
           exit
         endif
       enddo
-      if (degree_x2 <= 5) then
+      if (degree_x2 > 4) then
-        if(degree_x2 .ne. 0)then
+        cycle
+      else if(degree_x2 .ne. 0)then
           idx(l) = i
           l = l+1
       endif
-      endif
     enddo
     
   endif

src/Dets/slater_rules.irp.f

@@ -503,7 +503,7 @@ subroutine i_H_psi(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array)
   double precision               :: hij
   integer                        :: idx(0:Ndet)
   BEGIN_DOC
-  ! <key|H|psi> for the various Nstate
+  ! <key|H|psi> for the various Nstates
   END_DOC
   
   ASSERT (Nint > 0)
@@ -864,12 +864,15 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)
   Vt = 0.d0
   !$OMP DO SCHEDULE(guided)
   do i=1,n
-        call filter_connected(keys_tmp,keys_tmp(1,1,i),Nint,i-1,idx)
+        idx(0) = i
+        call filter_connected_davidson(keys_tmp,keys_tmp(1,1,i),Nint,i-1,idx)
         do jj=1,idx(0)
           j = idx(jj)
+          if ( (dabs(u_0(j)) > 1.d-7).or.((dabs(u_0(i)) > 1.d-7)) ) then
             call i_H_j(keys_tmp(1,1,j),keys_tmp(1,1,i),Nint,hij)
             vt (i) = vt (i) + hij*u_0(j)
             vt (j) = vt (j) + hij*u_0(i)
+          endif
         enddo
   enddo
   !$OMP END DO
@@ -884,3 +887,47 @@ end
 
 
 
+BEGIN_PROVIDER [ integer, N_con_int ]
+  implicit none
+  BEGIN_DOC
+  ! Number of integers to represent the connections between determinants
+  END_DOC
+  N_con_int = 1 + ishft(N_det-1,-11)
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer*8, det_connections, (N_con_int,N_det) ]
+  implicit none
+  BEGIN_DOC
+  !
+  END_DOC
+  integer                        :: i,j
+  integer                        :: degree
+  integer                        :: j_int, j_k, j_l
+  integer, allocatable           :: idx(:)
+  !$OMP PARALLEL DEFAULT (NONE) &
+  !$OMP SHARED(N_det, N_con_int, psi_det,N_int, det_connections) &
+  !$OMP PRIVATE(i,j_int,j_k,j_l,j,degree,idx)
+  allocate (idx(0:N_det))
+  !$OMP DO SCHEDULE(guided)
+  do i=1,N_det
+    do j_int=1,N_con_int
+      det_connections(j_int,i) = 0_8
+      j_k = ishft(j_int-1,11)
+      do j_l = j_k,min(j_k+2047,N_det), 32 
+        do j = j_l+1,min(j_l+32,i)
+          !DIR$ FORCEINLINE
+          call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int)
+          if (degree < 3) then
+            det_connections(j_int,i) = ibset( det_connections(j_int,i), iand(63,ishft(j_l,-5)) )
+            exit
+          endif
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP ENDDO
+  deallocate(idx)
+  !$OMP ENDPARALLEL
+
+END_PROVIDER
+