Added EZFIO.cfg and s2_eig test in diag

src/Bitmask/README.rst

@@ -128,6 +128,10 @@ Documentation
   Subroutine to print the content of a determinant in '+-' notation and
   hexadecimal representation.
 
+`debug_spindet <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks_routines.irp.f#L155>`_
+  Subroutine to print the content of a determinant in '+-' notation and
+  hexadecimal representation.
+
 `list_to_bitstring <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks_routines.irp.f#L29>`_
   Returns the physical string "string(N_int,2)" from the array of
   occupations "list(N_int*bit_kind_size,2)
@@ -135,5 +139,8 @@ Documentation
 `print_det <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks_routines.irp.f#L138>`_
   Subroutine to print the content of a determinant using the '+-' notation
 
+`print_spindet <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks_routines.irp.f#L171>`_
+  Subroutine to print the content of a determinant using the '+-' notation
+
 
 

src/Bitmask/bitmasks_routines.irp.f

@@ -126,7 +126,7 @@ subroutine debug_det(string,Nint)
   END_DOC
   integer, intent(in)            :: Nint
   integer(bit_kind), intent(in)  :: string(Nint,2)
-  character*(512)                :: output(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))
@@ -143,7 +143,7 @@ subroutine print_det(string,Nint)
   END_DOC
   integer, intent(in)            :: Nint
   integer(bit_kind), intent(in)  :: string(Nint,2)
-  character*(512)                :: output(2)
+  character*(2048)                :: output(2)
 
   call bitstring_to_str( output(1), string(1,1), Nint )
   call bitstring_to_str( output(2), string(1,2), Nint )
@@ -151,3 +151,34 @@ subroutine print_det(string,Nint)
   print *,  trim(output(2))
 
 end
+
+subroutine debug_spindet(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(1)
+  call bitstring_to_hexa( output(1), string(1,1), Nint )
+  print *,  trim(output(1)) 
+  call print_spindet(string,Nint)
+
+end
+
+subroutine print_spindet(string,Nint)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+  ! Subroutine to print the content of a determinant using the '+-' notation
+  END_DOC
+  integer, intent(in)            :: Nint
+  integer(bit_kind), intent(in)  :: string(Nint,2)
+  character*(2048)                :: output(1)
+
+  call bitstring_to_str( output(1), string(1,1), Nint )
+  print *,  trim(output(1))
+
+end

src/Determinants/README.rst

@@ -532,7 +532,7 @@ Documentation
 `save_dets_qmcchem <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/save_for_qmcchem.irp.f#L1>`_
   Undocumented
 
-`save_for_qmc <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/save_for_qmcchem.irp.f#L48>`_
+`save_for_qmc <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/save_for_qmcchem.irp.f#L46>`_
   Undocumented
 
 `save_natorb <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/save_natorb.irp.f#L1>`_
@@ -623,61 +623,49 @@ Documentation
 `n_con_int <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L1131>`_
   Number of integers to represent the connections between determinants
 
-`create_wf_of_psi_svd_matrix <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L473>`_
+`create_wf_of_psi_svd_matrix <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L483>`_
   Matrix of wf coefficients. Outer product of alpha and beta determinants
 
-`generate_all_alpha_beta_det_products <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L528>`_
+`generate_all_alpha_beta_det_products <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L538>`_
   Create a wave function from all possible alpha x beta determinants
 
-`get_index_in_psi_det_alpha_unique <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L131>`_
+`get_index_in_psi_det_alpha_unique <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L139>`_
   Returns the index of the determinant in the ``psi_det_alpha_unique`` array
 
-`get_index_in_psi_det_beta_unique <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L212>`_
+`get_index_in_psi_det_beta_unique <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L221>`_
   Returns the index of the determinant in the ``psi_det_beta_unique`` array
 
-`n_det_alpha_unique <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L54>`_
-  Unique alpha determinants
-
-`n_det_beta_unique <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L91>`_
-  Unique beta determinants
-
 `psi_det_alpha <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L25>`_
   List of alpha determinants of psi_det
 
-`psi_det_alpha_unique <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L53>`_
-  Unique alpha determinants
-
 `psi_det_beta <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L39>`_
   List of beta determinants of psi_det
 
-`psi_det_beta_unique <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L90>`_
+`psi_svd_alpha <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L578>`_
-  Unique beta determinants
-
-`psi_svd_alpha <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L568>`_
   SVD wave function
 
-`psi_svd_beta <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L569>`_
+`psi_svd_beta <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L579>`_
   SVD wave function
 
-`psi_svd_coefs <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L570>`_
+`psi_svd_coefs <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L580>`_
   SVD wave function
 
-`psi_svd_matrix <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L457>`_
+`psi_svd_matrix <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L467>`_
   Matrix of wf coefficients. Outer product of alpha and beta determinants
 
-`psi_svd_matrix_columns <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L398>`_
+`psi_svd_matrix_columns <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L407>`_
   Matrix of wf coefficients. Outer product of alpha and beta determinants
 
-`psi_svd_matrix_rows <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L397>`_
+`psi_svd_matrix_rows <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L406>`_
   Matrix of wf coefficients. Outer product of alpha and beta determinants
 
-`psi_svd_matrix_values <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L396>`_
+`psi_svd_matrix_values <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L405>`_
   Matrix of wf coefficients. Outer product of alpha and beta determinants
 
 `spin_det_search_key <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L9>`_
   Return an integer*8 corresponding to a determinant index for searching
 
-`write_spindeterminants <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L294>`_
+`write_spindeterminants <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L303>`_
   Undocumented
 
 `cisd <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/truncate_wf.irp.f#L1>`_

src/Determinants/determinants.irp.f

@@ -147,9 +147,7 @@ END_PROVIDER
    !$DIR FORCEINLINE
    bit_tmp(i) = occ_pattern_search_key(psi_occ_pattern(1,1,i),N_int)
  enddo
- print*,'passed 1'
  call i8sort(bit_tmp,iorder,N_det)
- print*,'passed 2'
  !DIR$ IVDEP
  do i=1,N_det
   do k=1,N_int
@@ -189,7 +187,6 @@ END_PROVIDER
     endif
   enddo
  enddo
-  print*,'passed 3'
 
  N_occ_pattern=0
  do i=1,N_det

src/Determinants/diagonalize_CI.irp.f

@@ -66,6 +66,7 @@ END_PROVIDER
     enddo
     integer :: i_state
     double precision :: s2
+    if (s2_eig) then
       i_state = 0
       do j=1,N_det
         call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2)
@@ -81,13 +82,16 @@ END_PROVIDER
           exit
         endif
       enddo
-!    if(i_state < min(N_states_diag,N_det))then
+    else
-!     print *, 'pb with the number of states'
+      do j=1,N_det
-!     print *, 'i_state = ',i_state
+        call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2)
-!     print *, 'N_states_diag ',N_states_diag
+        do i=1,N_det
-!     print *,'stopping ...'
+          CI_eigenvectors(i,j) = eigenvectors(i,j)
-!     stop
+        enddo
-!    endif
+        CI_electronic_energy(j) = eigenvalues(j)
+        CI_eigenvectors_s2(j) = s2
+      enddo
+    endif
     deallocate(eigenvectors,eigenvalues)
   endif
   

src/Determinants/diagonalize_CI_mono.irp.f

@@ -32,6 +32,7 @@
     integer :: i_state
     double precision :: s2
     i_state = 0
+    if (s2_eig) then
       do j=1,N_det
         call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2)
         if(dabs(s2-expected_s2).le.0.3d0)then
@@ -51,6 +52,19 @@
           exit
         endif
       enddo
+    else
+      do j=1,N_det
+        call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2)
+        if(dabs(eigenvectors(1,j)).gt.0.9d0)then
+          i_state += 1
+          do i=1,N_det
+            CI_eigenvectors_mono(i,i_state) = eigenvectors(i,j)
+          enddo
+          CI_electronic_energy_mono(i_state) = eigenvalues(j)
+          CI_eigenvectors_s2_mono(i_state) = s2
+        endif
+      enddo
+    endif
     deallocate(eigenvectors,eigenvalues)
   endif
   

src/Determinants/s2.irp.f

@@ -92,13 +92,12 @@ subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2)
  !$OMP PARALLEL DO DEFAULT(NONE) &
  !$OMP PRIVATE(i,j,s2_tmp) SHARED(n,psi_coefs_tmp,psi_keys_tmp,N_int) &
  !$OMP REDUCTION(+:s2) SCHEDULE(dynamic)
- do i = 1, n
+ do i=1,n
    call get_s2(psi_keys_tmp(1,1,i),psi_keys_tmp(1,1,i),s2_tmp,N_int)
-!  print*,'s2_tmp = ',s2_tmp
+   s2 += psi_coefs_tmp(i)*psi_coefs_tmp(i)*s2_tmp
-  do j = 1, n
+   do j=i+1,n
      call get_s2(psi_keys_tmp(1,1,i),psi_keys_tmp(1,1,j),s2_tmp,N_int)
-    if (s2_tmp == 0.d0) cycle
+     s2 += (psi_coefs_tmp(i)+psi_coefs_tmp(i))*psi_coefs_tmp(j)*s2_tmp
-    s2 += psi_coefs_tmp(i)*psi_coefs_tmp(j)*s2_tmp
    enddo
  enddo
  !$OMP END PARALLEL DO

src/Determinants/save_for_qmcchem.irp.f

@@ -7,8 +7,6 @@ subroutine save_dets_qmcchem
  integer, allocatable :: occ(:,:,:), occ_tmp(:,:)
  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: occ, occ_tmp
 
- read_wf = .True.
- TOUCH read_wf
  call ezfio_set_determinants_det_num(N_det)
  call ezfio_set_determinants_det_coef(psi_coef_sorted(1,1))
 
@@ -46,6 +44,8 @@ subroutine save_dets_qmcchem
 end
 
 program save_for_qmc
-  call save_dets_qmcchem
+ read_wf = .True.
+ TOUCH read_wf
+!  call save_dets_qmcchem
  call write_spindeterminants
 end

src/Determinants/spindeterminants.irp.f

@@ -50,80 +50,88 @@ BEGIN_PROVIDER [ integer(bit_kind), psi_det_beta, (N_int,psi_det_size) ]
  enddo
 END_PROVIDER
 
- BEGIN_PROVIDER [ integer(bit_kind), psi_det_alpha_unique, (N_int,psi_det_size) ]
+
-&BEGIN_PROVIDER [ integer, N_det_alpha_unique ]
+BEGIN_TEMPLATE
+
+ BEGIN_PROVIDER [ integer(bit_kind), psi_det_$alpha_unique, (N_int,psi_det_size) ]
+&BEGIN_PROVIDER [ integer, N_det_$alpha_unique ]
  implicit none
  BEGIN_DOC
- ! Unique alpha determinants
+ ! Unique $alpha determinants
  END_DOC
 
- integer                        :: i,k
+ integer                        :: i,j,k
  integer, allocatable           :: iorder(:)
  integer*8, allocatable         :: bit_tmp(:)
  integer*8                      :: last_key
  integer*8, external            :: spin_det_search_key
+ logical,allocatable            :: duplicate(:)
 
- allocate ( iorder(N_det), bit_tmp(N_det))
+ allocate ( iorder(N_det), bit_tmp(N_det), duplicate(N_det) )
 
  do i=1,N_det
    iorder(i) = i
-   bit_tmp(i) = spin_det_search_key(psi_det_alpha(1,i),N_int)
+   bit_tmp(i) = spin_det_search_key(psi_det_$alpha(1,i),N_int)
  enddo
 
  call i8sort(bit_tmp,iorder,N_det)
 
- N_det_alpha_unique = 0
+ N_det_$alpha_unique = 0
  last_key = 0_8
  do i=1,N_det
-  if (bit_tmp(i) /= last_key) then
    last_key = bit_tmp(i)
-    N_det_alpha_unique += 1
+   N_det_$alpha_unique += 1
    do k=1,N_int
-      psi_det_alpha_unique(k,N_det_alpha_unique) = psi_det_alpha(k,iorder(i))
+     psi_det_$alpha_unique(k,N_det_$alpha_unique) = psi_det_$alpha(k,iorder(i))
    enddo
+   duplicate(i) = .False.
+ enddo
+
+ j=1
+ do i=1,N_det_$alpha_unique-1
+   if (duplicate(i)) then
+     cycle
+   endif
+   j = i+1
+   do while (bit_tmp(j)==bit_tmp(i))
+     if (duplicate(j)) then
+       j += 1
+       cycle
+     endif
+     duplicate(j) = .True.
+     do k=1,N_int
+       if (psi_det_$alpha_unique(k,i) /= psi_det_$alpha_unique(k,j)) then
+         duplicate(j) = .False.
+         exit
        endif
      enddo
-
+     j+=1
- deallocate (iorder, bit_tmp)
+     if (j > N_det_$alpha_unique) then
-END_PROVIDER
+       exit
-
- BEGIN_PROVIDER [ integer(bit_kind), psi_det_beta_unique, (N_int,psi_det_size) ]
-&BEGIN_PROVIDER [ integer, N_det_beta_unique ]
- implicit none
- BEGIN_DOC
- ! Unique beta determinants
- END_DOC
-
- integer                        :: i,k
- integer, allocatable           :: iorder(:)
- integer*8, allocatable         :: bit_tmp(:)
- integer*8                      :: last_key
- integer*8, external            :: spin_det_search_key
-
- allocate ( iorder(N_det), bit_tmp(N_det))
-
- do i=1,N_det
-   iorder(i) = i
-   bit_tmp(i) = spin_det_search_key(psi_det_beta(1,i),N_int)
- enddo
-
- call i8sort(bit_tmp,iorder,N_det)
-
- N_det_beta_unique = 0
- last_key = 0_8
- do i=1,N_det
-  if (bit_tmp(i) /= last_key) then
-    last_key = bit_tmp(i)
-    N_det_beta_unique += 1
-    do k=1,N_int
-      psi_det_beta_unique(k,N_det_beta_unique) = psi_det_beta(k,iorder(i))
-    enddo
      endif
    enddo
+ enddo
 
- deallocate (iorder, bit_tmp)
+ j=1
+ do i=2,N_det_$alpha_unique
+   if (duplicate(i)) then
+     cycle
+  else
+    j += 1
+    psi_det_$alpha_unique(:,j) = psi_det_$alpha_unique(:,i)
+  endif
+ enddo
+ N_det_$alpha_unique = j
+
+ deallocate (iorder, bit_tmp, duplicate)
 END_PROVIDER
 
+SUBST [ alpha ]
+
+alpha ;;
+beta ;;
+
+END_TEMPLATE
 
 
 
@@ -150,6 +158,7 @@ integer function get_index_in_psi_det_alpha_unique(key,Nint)
 
   !DIR$ FORCEINLINE
   det_ref = spin_det_search_key(key,Nint)
+
   !DIR$ FORCEINLINE
   det_search = spin_det_search_key(psi_det_alpha_unique(1,1),Nint)
 
@@ -439,6 +448,7 @@ BEGIN_PROVIDER  [ double precision, psi_svd_matrix_values, (N_det,N_states) ]
   do k=1,N_det
     i = get_index_in_psi_det_alpha_unique(psi_det(1,1,k),N_int)
     j = get_index_in_psi_det_beta_unique (psi_det(1,2,k),N_int)
+
     do l=1,N_states
       psi_svd_matrix_values(k,l) = psi_coef(k,l)
     enddo

src/MRCC/mrcc_utils.irp.f

@@ -110,6 +110,7 @@ END_PROVIDER
     integer :: i_state
     double precision :: s2
     i_state = 0
+    if (s2_eig) then
       do j=1,N_det
         call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2)
         if(dabs(s2-expected_s2).le.0.3d0)then
@@ -124,6 +125,17 @@ END_PROVIDER
           exit
         endif
       enddo
+    else
+      do j=1,N_det
+        call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2)
+        i_state += 1
+        do i=1,N_det
+          CI_eigenvectors_dressed(i,i_state) = eigenvectors(i,j)
+        enddo
+        CI_electronic_energy_dressed(i_state) = eigenvalues(j)
+        CI_eigenvectors_s2_dressed(i_state) = s2
+      enddo
+    endif
     deallocate(eigenvectors,eigenvalues)
   endif
   

src/MonoInts/EZFIO.cfg

@@ -0,0 +1,6 @@
+[do_pseudo]
+type: logical
+doc:  Using pseudo potential integral of not
+interface: input
+default: False
+