Corrected CAS-SD bug

scripts/generate_h_apply.py

@@ -220,7 +220,7 @@ class H_apply(object):
 
       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)
+          norm_psi(k) = norm_psi(k) + psi_coef_generators(i_generator,k)*psi_coef_generators(i_generator,k)
         delta_pt2(k) = pt2(k) - pt2_old(k)
       enddo
       """
@@ -275,7 +275,7 @@ class H_apply(object):
         if (select_max(i_generator) < selection_criterion_min*selection_criterion_factor) then
           !$ call omp_set_lock(lck)
           do k=1,N_st
-            norm_psi(k) = norm_psi(k) + psi_coef(i_generator,k)*psi_coef(i_generator,k)
+            norm_psi(k) = norm_psi(k) + psi_coef_generators(i_generator,k)*psi_coef_generators(i_generator,k)
 !            delta_pt2(k) = 0.d0
 !            pt2_old(k) = 0.d0
 !            pt2(k) = select_max(i_generator)

src/CAS_SD_selected/cas_sd.irp.f

@@ -57,9 +57,9 @@ program full_ci
 
   ! Check that it is a CAS-SD
   logical :: in_cas
-  integer :: exc_max
+  integer :: exc_max, degree_min
   exc_max = 0
-  print *,  'CAS determinants'
+  print *,  'CAS determinants : ', N_det_generators
   do i=1,N_det_generators
     do k=i,N_det_generators
       call get_excitation_degree(psi_det_generators(1,1,k),psi_det_generators(1,1,i),degree,N_int)
@@ -71,25 +71,16 @@ program full_ci
   print *,  'Max excitation degree in the CAS :', exc_max
   do i=1,N_det
     in_cas = .False.
+    degree_min = 1000
     do k=1,N_det_generators
       call get_excitation_degree(psi_det_generators(1,1,k),psi_det(1,1,i),degree,N_int)
-      if (degree == 0) then
+      degree_min = min(degree_min,degree)
-        in_cas = .True.
-        exit
-      endif
     enddo
-    if (.not.in_cas) then
+    if (degree_min > 2) then
-      do k=i,N_det
-        call get_excitation_degree(psi_det(1,1,k),psi_det(1,1,i),degree,N_int)
-        if (degree > exc_max+2) then
         print *,  'Error : This is not a CAS-SD : '
-          print *,  'CAS determinant:'
+        print *,  'Excited determinant:', degree_min
-          call debug_det(psi_det(1,1,i),N_int)
-          print *,  'Excited determinant:', degree
         call debug_det(psi_det(1,1,k),N_int)
         stop
     endif
   enddo
-    endif
-  enddo
 end

src/Generators_CAS/generators.irp.f

@@ -15,9 +15,9 @@ BEGIN_PROVIDER [ integer, N_det_generators ]
       do k=1,N_int
         good = good .and. (                                          &
             iand(not(cas_bitmask(k,1,l)), psi_det(k,1,i)) ==         &
-            iand(not(cas_bitmask(k,1,l)), psi_det(k,1,1)) ) .and. (  &
+            iand(not(cas_bitmask(k,1,l)), HF_bitmask(k,1)) ) .and. ( &
             iand(not(cas_bitmask(k,2,l)), psi_det(k,2,i)) ==         &
-            iand(not(cas_bitmask(k,2,l)), psi_det(k,2,1)) )
+            iand(not(cas_bitmask(k,2,l)), HF_bitmask(k,2)) )
       enddo
       if (good) then
         exit
@@ -31,7 +31,8 @@ BEGIN_PROVIDER [ integer, N_det_generators ]
   call write_int(output_dets,N_det_generators,'Number of generators')
 END_PROVIDER
 
-BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ]
+ BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ]
+&BEGIN_PROVIDER [ double precision, psi_coef_generators, (psi_det_size,N_states) ]
   implicit none
   BEGIN_DOC
   ! For Single reference wave functions, the generator is the
@@ -46,9 +47,9 @@ BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ]
       do k=1,N_int
         good = good .and. (                                         &
             iand(not(cas_bitmask(k,1,l)), psi_det(k,1,i)) ==         &
-            iand(not(cas_bitmask(k,1,l)), psi_det(k,1,1)) ) .and. (  &
+            iand(not(cas_bitmask(k,1,l)), HF_bitmask(k,1)) .and. (   &
             iand(not(cas_bitmask(k,2,l)), psi_det(k,2,i)) ==         &
-            iand(not(cas_bitmask(k,2,l)), psi_det(k,2,1)) )
+            iand(not(cas_bitmask(k,2,l)), HF_bitmask(k,2) )) )
       enddo
       if (good) then
         exit
@@ -60,6 +61,7 @@ BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ]
         psi_det_generators(k,1,m) = psi_det(k,1,i)
         psi_det_generators(k,2,m) = psi_det(k,2,i)
       enddo
+      psi_coef_generators(m,:) = psi_coef(m,:)
 !     call debug_det(psi_det_generators(1,1,m),N_int)
     endif
   enddo

src/Generators_full/generators.irp.f

@@ -34,7 +34,8 @@ BEGIN_PROVIDER [ integer, N_det_generators ]
  call write_int(output_dets,N_det_generators,'Number of generators')
 END_PROVIDER
 
-BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ]
+ BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ]
+&BEGIN_PROVIDER [ double precision, psi_coef_generators, (psi_det_size,N_states) ]
  implicit none
  BEGIN_DOC
  ! For Single reference wave functions, the generator is the
@@ -46,6 +47,7 @@ BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ]
      psi_det_generators(k,1,i) = psi_det_sorted(k,1,i)
      psi_det_generators(k,2,i) = psi_det_sorted(k,2,i)
    enddo
+   psi_coef_generators(i,:) = psi_coef_sorted(i,:)
  enddo
 
 END_PROVIDER

src/Generators_restart/generators.irp.f

@@ -17,8 +17,8 @@ BEGIN_PROVIDER [ integer, N_det_generators ]
 END_PROVIDER
 
 
- BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,N_det_generators) ]
+ BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ]
-&BEGIN_PROVIDER [ double precision, psi_coef_generators, (N_det_generators,N_states) ]
+&BEGIN_PROVIDER [ double precision, psi_coef_generators, (psi_det_size,N_states) ]
  implicit none
  BEGIN_DOC
  ! read wf

src/MRCC/mrcc_dress.irp.f

@@ -53,22 +53,6 @@ subroutine mrcc_dress(delta_ij_sd_,Ndet_sd,i_generator,n_selected,det_buffer,Nin
   double precision :: haj
   double precision :: f(N_states)
 
-
-!  call i_h_j(psi_det(1,1,1), psi_det(1,1,64),Nint,hka)
-!  call debug_det(psi_det(1,1,1), N_int)
-!  call debug_det(psi_det(1,1,64), N_int)
-  double precision :: phase
-  integer          :: exc(0:2,2,2)
-!  call get_excitation(psi_det(1,1,1),psi_det(1,1,64),exc,degree(1),phase,Nint)
-  integer :: h1, p1, h2, p2, s1, s2
-!  call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
-!  print *,  hka
-!  print *,  h1, p1, h2, p2
-!  print *,  s1, s2
-!  pause
-
-
-
   do i=1,N_tq
     call get_excitation_degree_vector(psi_sd,tq(1,1,i),degree,Nint,Ndet_sd,idx)
     call i_h_j(tq(1,1,i),tq(1,1,i),Nint,haa)
@@ -83,21 +67,6 @@ subroutine mrcc_dress(delta_ij_sd_,Ndet_sd,i_generator,n_selected,det_buffer,Nin
           delta_ij_sd_(idx(k), idx(j),m) += haj*hka* f(m)
           delta_ij_sd_(idx(j), idx(k),m) += haj*hka* f(m)
         enddo
-  call get_excitation(tq(1,1,i),psi_sd(1,1,idx(j)),exc,degree(1),phase,Nint)
-  call decode_exc(exc,degree(1),h1,p1,h2,p2,s1,s2)
-  if ( (h1 == 1).and. &
-       (p1 == 6).and. &
-       (h2 == 1).and. &
-       (p2 == 6).and. &
-       (s1 == 1).and. &
-       (s2 == 2) ) then
-    call debug_det(tq(1,1,i), N_int)
-    call debug_det(psi_sd(1,1,idx(j)), N_int)
-    print *,  haj
-    pause
-  endif
-
-
       enddo 
     enddo
   enddo
@@ -196,13 +165,6 @@ END_PROVIDER
         exit
       endif
     enddo
-!    if(i_state < min(N_states_diag,N_det))then
-!     print *, 'pb with the number of states'
-!     print *, 'i_state = ',i_state
-!     print *, 'N_states_diag ',N_states_diag
-!     print *,'stopping ...'
-!     stop
-!    endif
     deallocate(eigenvectors,eigenvalues)
   endif
   

src/Makefile.common

@@ -8,6 +8,9 @@ default:
 veryclean:
 	 clean_modules.sh
 
+clean:
+	IN_MAKE=1 make clean
+
 else  # Called by scripts/build_module.sh
 
 default: all .gitignore

src/SingleRefMethod/generators.irp.f

@@ -10,7 +10,7 @@ BEGIN_PROVIDER [ integer, N_det_generators ]
 END_PROVIDER
 
  BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ]
-&BEGIN_PROVIDER [ integer(bit_kind), psi_coef_generators, (1,N_states) ]
+&BEGIN_PROVIDER [ double precision, psi_coef_generators, (psi_det_size,N_states) ]
  implicit none
  BEGIN_DOC
  ! For Single reference wave functions, the generator is the