Merge pull request #14 from scemama/master

Merger with scemama

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/CAS_SD/cas_sd.irp.f

@@ -1,6 +1,7 @@
 program full_ci
   implicit none
   integer                        :: i,k
+  integer                        :: N_det_old
 
   
   double precision, allocatable  :: pt2(:), norm_pert(:), H_pert_diag(:)
@@ -9,6 +10,7 @@ program full_ci
   allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st))
   character*(64)                 :: perturbation
   
+  N_det_old = 0
   pt2 = 1.d0
   diag_algorithm = "Lapack"
   if (N_det > n_det_max_cas_sd) then
@@ -29,6 +31,7 @@ program full_ci
   endif
 
   do while (N_det < n_det_max_cas_sd.and.maxval(abs(pt2(1:N_st))) > pt2_max)
+    N_det_old = N_det
     call H_apply_CAS_SD(pt2, norm_pert, H_pert_diag,  N_st)
 
     PROVIDE  psi_coef
@@ -53,10 +56,11 @@ program full_ci
     if (abort_all) then
       exit
     endif
+    if (N_det == N_det_old) then
+      exit
+    endif
   enddo
 
-  ! Check that it is a CAS-SD
-  logical :: in_cas
   integer :: exc_max, degree_min
   exc_max = 0
   print *,  'CAS determinants : ', N_det_generators
@@ -69,18 +73,4 @@ program full_ci
     print *,  ''
   enddo
   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)
-      degree_min = min(degree_min,degree)
-    enddo
-    if (degree_min > 2) then
-        print *,  'Error : This is not a CAS-SD : '
-        print *,  'Excited determinant:', degree_min
-        call debug_det(psi_det(1,1,k),N_int)
-        stop
-    endif
-  enddo
 end

src/CAS_SD/cas_sd_selected.irp.f

@@ -50,14 +50,30 @@ program full_ci
     print *,  'E        = ', CI_energy
     print *,  'E+PT2    = ', CI_energy+pt2
     print *,  '-----'
-    call ezfio_set_full_ci_energy(CI_energy)
+    call ezfio_set_cas_sd_energy(CI_energy(1))
     if (abort_all) then
       exit
     endif
   enddo
+  call diagonalize_CI
+
+   if(do_pt2_end)then
+    print*,'Last iteration only to compute the PT2'
+    threshold_selectors = 1.d0
+    threshold_generators = 0.999d0
+    call H_apply_CAS_SD_PT2(pt2, norm_pert, H_pert_diag,  N_st)
+
+    print *,  'Final step'
+    print *,  'N_det    = ', N_det
+    print *,  'N_states = ', N_states
+    print *,  'PT2      = ', pt2
+    print *,  'E        = ', CI_energy
+    print *,  'E+PT2    = ', CI_energy+pt2
+    print *,  '-----'
+    call ezfio_set_cas_sd_energy_pt2(CI_energy(1)+pt2(1))
+   endif
+
 
-  ! Check that it is a CAS-SD
-  logical :: in_cas
   integer :: exc_max, degree_min
   exc_max = 0
   print *,  'CAS determinants : ', N_det_cas
@@ -70,18 +86,4 @@ program full_ci
     print *,  ''
   enddo
   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_cas
-      call get_excitation_degree(psi_cas(1,1,k),psi_det(1,1,i),degree,N_int)
-      degree_min = min(degree_min,degree)
-    enddo
-    if (degree_min > 2) then
-        print *,  'Error : This is not a CAS-SD : '
-        print *,  'Excited determinant:', degree_min
-        call debug_det(psi_det(1,1,k),N_int)
-        stop
-    endif
-  enddo
 end

src/Determinants/H_apply.irp.f

@@ -53,13 +53,17 @@ subroutine resize_H_apply_buffer(new_size,iproc)
   double precision,  pointer     :: buffer_e2(:,:)
   integer                        :: i,j,k
   integer                        :: Ndet
+
+  BEGIN_DOC
+! Resizes the H_apply buffer of proc iproc. The buffer lock should
+! be set before calling this function.
+  END_DOC
   PROVIDE H_apply_buffer_allocated
   
   ASSERT (new_size > 0)
   ASSERT (iproc >= 0)
   ASSERT (iproc < nproc)
   
-  call omp_set_lock(H_apply_buffer_lock(1,iproc))
   allocate ( buffer_det(N_int,2,new_size),                           &
       buffer_coef(new_size,N_states),                                &
       buffer_e2(new_size,N_states) )
@@ -93,7 +97,6 @@ subroutine resize_H_apply_buffer(new_size,iproc)
   
   H_apply_buffer(iproc)%sze = new_size
   H_apply_buffer(iproc)%N_det = min(new_size,H_apply_buffer(iproc)%N_det)
-  call omp_unset_lock(H_apply_buffer_lock(1,iproc))
   
 end
 
@@ -101,8 +104,7 @@ subroutine copy_H_apply_buffer_to_wf
   use omp_lib
   implicit none
   BEGIN_DOC
-! Copies the H_apply buffer to psi_coef. You need to touch psi_det, psi_coef and N_det
+! Copies the H_apply buffer to psi_coef.
-! after calling this function.
 ! After calling this subroutine, N_det, psi_det and psi_coef need to be touched
   END_DOC
   integer(bit_kind), allocatable :: buffer_det(:,:,:)
@@ -181,42 +183,76 @@ subroutine copy_H_apply_buffer_to_wf
   call normalize(psi_coef,N_det)
   SOFT_TOUCH N_det psi_det psi_coef
   
-  call debug_unicity_of_determinants
+  logical :: found_duplicates
+  call remove_duplicates_in_psi_det(found_duplicates)
 end
 
-subroutine debug_unicity_of_determinants
+subroutine remove_duplicates_in_psi_det(found_duplicates)
   implicit none
+  logical, intent(out) :: found_duplicates
   BEGIN_DOC
-! This subroutine checks that there are no repetitions in the wave function
+! Removes duplicate determinants in the wave function.
   END_DOC
-  logical :: same, failed
+  integer                        :: i,j,k
-  integer :: i,k
+  integer(bit_kind), allocatable :: bit_tmp(:)
-  print *,  "======= DEBUG UNICITY ========="
+  logical,allocatable            :: duplicate(:)
-  failed = .False.
+
-  do i=2,N_det
+  allocate (duplicate(N_det), bit_tmp(N_det))
-    same = .True.
+
-    do k=1,N_int
+  do i=1,N_det
-      if ( psi_det_sorted_bit(k,1,i) /= psi_det_sorted_bit(k,1,i-1) ) then
+    integer, external            :: det_search_key
-        same = .False.
+    !$DIR FORCEINLINE
-        exit
+    bit_tmp(i) = det_search_key(psi_det_sorted_bit(1,1,i),N_int)
+    duplicate(i) = .False.
+  enddo
+
+  do i=1,N_det-1
+    if (duplicate(i)) then
+      cycle
     endif
-      if ( psi_det_sorted_bit(k,2,i) /= psi_det_sorted_bit(k,2,i-1) ) then
+    j = i+1
-        same = .False.
+    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_sorted_bit(k,1,i) /= psi_det_sorted_bit(k,1,j) ) &
+        .or. (psi_det_sorted_bit(k,2,i) /= psi_det_sorted_bit(k,2,j) ) ) then
+          duplicate(j) = .False.
           exit
         endif
       enddo
-    if (same) then
+      j += 1
-      failed = .True.
+      if (j > N_det) then
-      call debug_det(psi_det_sorted_bit(1,1,i))
+        exit
+      endif
+    enddo
+  enddo
+
+  found_duplicates = .False.
+  do i=1,N_det
+    if (duplicate(i)) then
+      found_duplicates = .True.
+      exit
     endif
   enddo
 
-  if (failed) then
+  if (found_duplicates) then
-    print *,  '======= Determinants not unique : Failed ! ========='
+    call write_bool(output_determinants,found_duplicates,'Found duplicate determinants')
-    stop
+    k=0
-  else
+    do i=1,N_det
-    print *,  '======= Determinants are unique : OK ! ========='
+      if (.not.duplicate(i)) then
+        k += 1
+        psi_det(:,:,k) = psi_det_sorted_bit (:,:,i)
+        psi_coef(k,:)  = psi_coef_sorted_bit(i,:)
       endif
+    enddo
+    N_det = k
+    TOUCH N_det psi_det psi_coef
+  endif
+  deallocate (duplicate,bit_tmp)
 end
 
 subroutine fill_H_apply_buffer_no_selection(n_selected,det_buffer,Nint,iproc)
@@ -231,11 +267,11 @@ subroutine fill_H_apply_buffer_no_selection(n_selected,det_buffer,Nint,iproc)
   integer                        :: i,j,k
   integer                        :: new_size
   PROVIDE H_apply_buffer_allocated
+  call omp_set_lock(H_apply_buffer_lock(1,iproc))
   new_size = H_apply_buffer(iproc)%N_det + n_selected
   if (new_size > H_apply_buffer(iproc)%sze) then
     call resize_h_apply_buffer(max(2*H_apply_buffer(iproc)%sze,new_size),iproc)
   endif
-  call omp_set_lock(H_apply_buffer_lock(1,iproc))
   do i=1,H_apply_buffer(iproc)%N_det
     ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,1,i)) )== elec_alpha_num)
     ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,2,i))) == elec_beta_num)
@@ -250,7 +286,7 @@ subroutine fill_H_apply_buffer_no_selection(n_selected,det_buffer,Nint,iproc)
   enddo
   do j=1,N_states
     do i=1,N_selected
-      H_apply_buffer(iproc)%coef(i,j) = 0.d0
+      H_apply_buffer(iproc)%coef(i+H_apply_buffer(iproc)%N_det,j) = 0.d0
     enddo
   enddo
   H_apply_buffer(iproc)%N_det = new_size

src/Determinants/README.rst

@@ -40,15 +40,11 @@ Documentation
 .. Do not edit this section. It was auto-generated from the
 .. NEEDED_MODULES file.
 
-`copy_h_apply_buffer_to_wf <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/H_apply.irp.f#L100>`_
+`copy_h_apply_buffer_to_wf <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/H_apply.irp.f#L103>`_
-  Copies the H_apply buffer to psi_coef. You need to touch psi_det, psi_coef and N_det
+  Copies the H_apply buffer to psi_coef.
-  after calling this function.
   After calling this subroutine, N_det, psi_det and psi_coef need to be touched
 
-`debug_unicity_of_determinants <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/H_apply.irp.f#L187>`_
+`fill_h_apply_buffer_no_selection <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/H_apply.irp.f#L258>`_
-  This subroutine checks that there are no repetitions in the wave function
-
-`fill_h_apply_buffer_no_selection <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/H_apply.irp.f#L222>`_
   Fill the H_apply buffer with determiants for CISD
 
 `h_apply_buffer_allocated <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/H_apply.irp.f#L15>`_
@@ -59,8 +55,12 @@ Documentation
   Buffer of determinants/coefficients/perturbative energy for H_apply.
   Uninitialized. Filled by H_apply subroutines.
 
+`remove_duplicates_in_psi_det <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/H_apply.irp.f#L190>`_
+  Removes duplicate determinants in the wave function.
+
 `resize_h_apply_buffer <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/H_apply.irp.f#L48>`_
-  Undocumented
+  Resizes the H_apply buffer of proc iproc. The buffer lock should
+  be set before calling this function.
 
 `cisd_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/SC2.irp.f#L1>`_
   CISD+SC2 method              :: take off all the disconnected terms of a CISD (selected or not)
@@ -187,10 +187,10 @@ Documentation
 `det_svd <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/det_svd.irp.f#L1>`_
   Computes the SVD of the Alpha x Beta determinant coefficient matrix
 
-`filter_3_highest_electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L426>`_
+`filter_3_highest_electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L423>`_
   Returns a determinant with only the 3 highest electrons
 
-`int_of_3_highest_electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L391>`_
+`int_of_3_highest_electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L388>`_
   Returns an integer*8 as :
   .br
   |_<--- 21 bits ---><--- 21 bits ---><--- 21 bits --->|
@@ -207,26 +207,26 @@ Documentation
 `n_det <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L3>`_
   Number of determinants in the wave function
 
-`psi_average_norm_contrib <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L276>`_
+`psi_average_norm_contrib <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L273>`_
   Contribution of determinants to the state-averaged density
 
-`psi_average_norm_contrib_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L306>`_
+`psi_average_norm_contrib_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L303>`_
   Wave function sorted by determinants contribution to the norm (state-averaged)
 
-`psi_coef <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L230>`_
+`psi_coef <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L227>`_
   The wave function coefficients. Initialized with Hartree-Fock if the EZFIO file
   is empty
 
-`psi_coef_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L305>`_
+`psi_coef_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L302>`_
   Wave function sorted by determinants contribution to the norm (state-averaged)
 
-`psi_coef_sorted_ab <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L453>`_
+`psi_coef_sorted_ab <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L450>`_
   Determinants on which we apply <i|H|j>.
   They are sorted by the 3 highest electrons in the alpha part,
   then by the 3 highest electrons in the beta part to accelerate
   the research of connected determinants.
 
-`psi_coef_sorted_bit <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L336>`_
+`psi_coef_sorted_bit <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L333>`_
   Determinants on which we apply <i|H|psi> for perturbation.
   They are sorted by determinants interpreted as integers. Useful
   to accelerate the search of a random determinant in the wave
@@ -239,46 +239,46 @@ Documentation
 `psi_det_size <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L47>`_
   Size of the psi_det/psi_coef arrays
 
-`psi_det_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L304>`_
+`psi_det_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L301>`_
   Wave function sorted by determinants contribution to the norm (state-averaged)
 
-`psi_det_sorted_ab <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L452>`_
+`psi_det_sorted_ab <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L449>`_
   Determinants on which we apply <i|H|j>.
   They are sorted by the 3 highest electrons in the alpha part,
   then by the 3 highest electrons in the beta part to accelerate
   the research of connected determinants.
 
-`psi_det_sorted_bit <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L335>`_
+`psi_det_sorted_bit <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L332>`_
   Determinants on which we apply <i|H|psi> for perturbation.
   They are sorted by determinants interpreted as integers. Useful
   to accelerate the search of a random determinant in the wave
   function.
 
-`psi_det_sorted_next_ab <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L454>`_
+`psi_det_sorted_next_ab <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L451>`_
   Determinants on which we apply <i|H|j>.
   They are sorted by the 3 highest electrons in the alpha part,
   then by the 3 highest electrons in the beta part to accelerate
   the research of connected determinants.
 
-`read_dets <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L583>`_
+`read_dets <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L580>`_
   Reads the determinants from the EZFIO file
 
-`save_wavefunction <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L630>`_
+`save_wavefunction <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L627>`_
   Save the wave function into the EZFIO file
 
-`save_wavefunction_general <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L649>`_
+`save_wavefunction_general <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L646>`_
   Save the wave function into the EZFIO file
 
-`save_wavefunction_unsorted <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L640>`_
+`save_wavefunction_unsorted <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L637>`_
   Save the wave function into the EZFIO file
 
-`sort_dets_by_3_highest_electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L474>`_
+`sort_dets_by_3_highest_electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L471>`_
   Determinants on which we apply <i|H|j>.
   They are sorted by the 3 highest electrons in the alpha part,
   then by the 3 highest electrons in the beta part to accelerate
   the research of connected determinants.
 
-`sort_dets_by_det_search_key <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L349>`_
+`sort_dets_by_det_search_key <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L346>`_
   Determinants are sorted are sorted according to their det_search_key.
   Useful to accelerate the search of a random determinant in the wave
   function.
@@ -316,7 +316,7 @@ Documentation
 `diag_algorithm <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/diagonalize_CI.irp.f#L1>`_
   Diagonalization algorithm (Davidson or Lapack)
 
-`diagonalize_ci <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/diagonalize_CI.irp.f#L96>`_
+`diagonalize_ci <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/diagonalize_CI.irp.f#L100>`_
   Replace the coefficients of the CI states by the coefficients of the
   eigenstates of the CI matrix
 
@@ -345,7 +345,7 @@ Documentation
 `ci_electronic_energy_mono <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/diagonalize_CI_mono.irp.f#L1>`_
   Eigenvectors/values of the CI matrix
 
-`diagonalize_ci_mono <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/diagonalize_CI_mono.irp.f#L59>`_
+`diagonalize_ci_mono <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/diagonalize_CI_mono.irp.f#L73>`_
   Replace the coefficients of the CI states by the coefficients of the
   eigenstates of the CI matrix
 
@@ -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_states_diag
-!     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_states_diag
+        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

@@ -88,19 +88,19 @@ subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2)
  double precision, intent(out) :: s2
  integer :: i,j,l
  double precision :: s2_tmp
- s2 = S_z2_Sz
+ s2 = 0.d0
  !$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
+ s2 +=  S_z2_Sz
 end
 

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/README.rst

@@ -51,7 +51,7 @@ Documentation
 `ci_electronic_energy_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_utils.irp.f#L78>`_
   Eigenvectors/values of the CI matrix
 
-`ci_energy_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_utils.irp.f#L132>`_
+`ci_energy_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_utils.irp.f#L144>`_
   N_states lowest eigenvalues of the dressed CI matrix
 
 `delta_ij <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_utils.irp.f#L43>`_
@@ -60,7 +60,7 @@ Documentation
 `delta_ij_non_cas <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_utils.irp.f#L34>`_
   Dressing matrix in SD basis
 
-`diagonalize_ci_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_utils.irp.f#L147>`_
+`diagonalize_ci_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC/mrcc_utils.irp.f#L159>`_
   Replace the coefficients of the CI states by the coefficients of the
   eigenstates of the CI matrix
 

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_states_diag
+        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
+

src/MonoInts/pot_ao_ints.irp.f

@@ -10,7 +10,7 @@
  integer           :: i,j,k,l,n_pt_in,m
  double precision  ::overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult
 
-  if (do_pseudo.eqv..TRUE.) then
+  if (do_pseudo) then
     ao_nucl_elec_integral = ao_nucl_elec_integral_pseudo
   else
     ao_nucl_elec_integral = 0.d0

src/Perturbation/selection.irp.f

@@ -19,6 +19,7 @@ subroutine fill_H_apply_buffer_selection(n_selected,det_buffer,e_2_pert_buffer,c
   ASSERT (Nint > 0)
   ASSERT (N_int == N_int)
   ASSERT (N_selected >= 0)
+  call omp_set_lock(H_apply_buffer_lock(1,iproc))
   smax = selection_criterion
   smin = selection_criterion_min
   new_size = H_apply_buffer(iproc)%N_det + n_selected
@@ -26,7 +27,6 @@ subroutine fill_H_apply_buffer_selection(n_selected,det_buffer,e_2_pert_buffer,c
   if (new_size > h_apply_buffer(iproc)%sze) then
     call resize_h_apply_buffer(max(h_apply_buffer(iproc)%sze*2,new_size),iproc)
   endif
-  call omp_set_lock(H_apply_buffer_lock(1,iproc))
   do i=1,H_apply_buffer(iproc)%N_det
     ASSERT (sum(popcnt(h_apply_buffer(iproc)%det(:,1,i)) )== elec_alpha_num)
     ASSERT (sum(popcnt(h_apply_buffer(iproc)%det(:,2,i))) == elec_beta_num)

src/Pseudo_integrals/int.f90

@@ -1,7 +1,15 @@
-!! Vps= <Phi_A|Vloc(C)+Cpp(C)| Phi_B>
+!! INFO : You can display equations using : http://www.codecogs.com/latex/eqneditor.php
+
+!!
+!! {\tt Vps}(C) = \langle \Phi_A|{\tt Vloc}(C)+{\tt Vpp}(C)| \Phi_B \rangle
+!!
+!! with :
+!!
+!! {\tt Vloc}(C)=\sum_{k=1}^{\tt klocmax} v_k r_C^{n_k} \exp(-dz_k r_C^2) \\
+!! 
+!! {\tt Vpp}(C)=\sum_{l=0}^{\tt lmax}\left( \sum_{k=1}^{\tt kmax} v_{kl}
+!!   r_C^{n_{kl}} \exp(-dz_{kl} r_C)^2 \right) |l\rangle \langle l|
 !!
-!! with: Vloc(C)=\sum_{k=1}^klocmax v_k rC**n_k exp(-dz_k rC**2)
-!!       Vpp(C)=\sum_{l=0}^lmax\sum_{k=1}^kmax v_kl rC**n_kl exp(-dz_kl rC**2)*|l><l|
 double precision function Vps &
 (a,n_a,g_a,b,n_b,g_b,c,klocmax,v_k,n_k,dz_k,lmax,kmax,v_kl,n_kl,dz_kl)
 implicit none