Removed diagonalize_s2

plugins/MRCC_Utils/mrcc_utils.irp.f

@@ -167,7 +167,7 @@ END_PROVIDER
      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
+         if(dabs(s2-expected_s2).le.0.5d0)then
            i_state += 1
            do i=1,N_det
              CI_eigenvectors_dressed(i,i_state) = eigenvectors(i,j)
@@ -193,7 +193,7 @@ END_PROVIDER
      deallocate(eigenvectors,eigenvalues)
    endif
    
-   if(diagonalize_s2.and.n_states_diag > 1.and. n_det >= n_states_diag)then
+   if(s2_eig.and.n_states_diag > 1.and. n_det >= n_states_diag)then
       ! Diagonalizing S^2 within the "n_states_diag" states found
       allocate(s2_eigvalues(N_states_diag))
       call diagonalize_s2_betweenstates(psi_det,CI_eigenvectors_dressed,n_det,size(psi_det,3),size(CI_eigenvectors_dressed,1),min(n_states_diag,n_det),s2_eigvalues)
@@ -204,8 +204,6 @@ END_PROVIDER
         enddo
       enddo
 
-     if(s2_eig)then
-
       ! Browsing the "n_states_diag" states and getting the lowest in energy "n_states" ones that have the S^2 value
       ! closer to the "expected_s2" set as input
 
@@ -213,7 +211,7 @@ END_PROVIDER
       good_state_array = .False.
       i_state = 0
       do j = 1, N_states_diag
-         if(dabs(s2_eigvalues(j)-expected_s2).le.0.3d0)then
+        if(dabs(s2_eigvalues(j)-expected_s2).le.0.5d0)then
           good_state_array(j) = .True.
           i_state +=1
           index_good_state_array(i_state) = j
@@ -261,26 +259,6 @@ END_PROVIDER
       enddo
       deallocate(index_good_state_array,good_state_array)
 
-     else
-
-       ! Sorting the N_states_diag by energy, whatever the S^2 value is
-
-       allocate(e_array(n_states_diag),iorder(n_states_diag))
-       do j = 1, N_states_diag
-         call u0_H_u_0_mrcc(e_0,CI_eigenvectors_dressed(1,j),n_det,psi_det,N_int,mrcc_state)
-         e_array(j) = e_0
-         iorder(j) = j
-       enddo
-       call dsort(e_array,iorder,n_states_diag)
-       do j = 1, N_states_diag
-         CI_electronic_energy_dressed(j) = e_array(j)
-         do i = 1, N_det
-           CI_eigenvectors_dressed(i,j) = psi_coef(i,iorder(j))
-         enddo
-         CI_eigenvectors_s2_dressed(j) = s2_eigvalues(iorder(j))
-       enddo
-       deallocate(e_array,iorder)
-     endif
       deallocate(s2_eigvalues)
 
    endif
@@ -297,6 +275,7 @@ BEGIN_PROVIDER [ double precision, CI_energy_dressed, (N_states_diag) ]
   character*(8)                  :: st
   call write_time(output_determinants)
   do j=1,N_states_diag
+    write(st,'(I4)') j
     CI_energy_dressed(j) = CI_electronic_energy_dressed(j) + nuclear_repulsion
     call write_double(output_determinants,CI_energy(j),'Energy of state '//trim(st))
     call write_double(output_determinants,CI_eigenvectors_s2(j),'S^2 of state '//trim(st))

src/Determinants/EZFIO.cfg

@@ -40,12 +40,6 @@ doc: Force the wave function to be an eigenfunction of S^2
 interface: ezfio,provider,ocaml
 default: False
 
-[diagonalize_s2]
-type: logical
-doc: Diagonalize the S^2 operator within the n_states_diag states required. Notice : the vectors are sorted by increasing S^2 values. 
-interface: ezfio,provider,ocaml
-default: True 
-
 [threshold_davidson]
 type: Threshold
 doc: Thresholds of Davidson's algorithm

src/Determinants/diagonalize_CI.irp.f

@@ -92,7 +92,7 @@ END_PROVIDER
          call get_s2_u0(psi_det,eigenvectors(1,j),N_det,size(eigenvectors,1),s2)
          s2_eigvalues(j) = s2
          ! Select at least n_states states with S^2 values closed to "expected_s2"
-         if(dabs(s2-expected_s2).le.0.3d0)then
+         if(dabs(s2-expected_s2).le.0.5d0)then
            i_state +=1
            index_good_state_array(i_state) = j
            good_state_array(j) = .True.
@@ -133,7 +133,7 @@ END_PROVIDER
          print*,'  We did not find any state with S^2 values close to ',expected_s2
          print*,'  We will then set the first N_states eigenvectors of the H matrix'
          print*,'  as the CI_eigenvectors'
-         print*,'  You should consider more states and maybe ask for diagonalize_s2 to be .True. or just enlarge the CI space'
+         print*,'  You should consider more states and maybe ask for s2_eig to be .True. or just enlarge the CI space'
          print*,''
          do j=1,min(N_states_diag,N_det)
            do i=1,N_det
@@ -159,7 +159,8 @@ END_PROVIDER
      deallocate(eigenvectors,eigenvalues)
    endif
    
-   if(diagonalize_s2.and.n_states_diag > 1.and. n_det >= n_states_diag)then
+
+   if( s2_eig.and.(n_states_diag > 1).and.(n_det >= n_states_diag) )then
       ! Diagonalizing S^2 within the "n_states_diag" states found
       allocate(s2_eigvalues(N_states_diag))
       call diagonalize_s2_betweenstates(psi_det,CI_eigenvectors,n_det,size(psi_det,3),size(CI_eigenvectors,1),min(n_states_diag,n_det),s2_eigvalues)
@@ -170,8 +171,6 @@ END_PROVIDER
         enddo
       enddo
      
-     if(s2_eig)then
-       
       ! Browsing the "n_states_diag" states and getting the lowest in energy "n_states" ones that have the S^2 value
       ! closer to the "expected_s2" set as input
       
@@ -179,7 +178,7 @@ END_PROVIDER
       good_state_array = .False.
       i_state = 0
       do j = 1, N_states_diag
-         if(dabs(s2_eigvalues(j)-expected_s2).le.0.3d0)then
+        if(dabs(s2_eigvalues(j)-expected_s2).le.0.5d0)then
           good_state_array(j) = .True.
           i_state +=1
           index_good_state_array(i_state) = j
@@ -227,26 +226,6 @@ END_PROVIDER
       enddo
       deallocate(index_good_state_array,good_state_array)
        
-     else
-       
-       ! Sorting the N_states_diag by energy, whatever the S^2 value is
-       
-       allocate(e_array(n_states_diag),iorder(n_states_diag))
-       do j = 1, N_states_diag
-         call u0_H_u_0(e_0,CI_eigenvectors(1,j),n_det,psi_det,N_int)
-         e_array(j) = e_0
-         iorder(j) = j
-       enddo
-       call dsort(e_array,iorder,n_states_diag)
-       do j = 1, N_states_diag
-         CI_electronic_energy(j) = e_array(j)
-         do i = 1, N_det
-           CI_eigenvectors(i,j) = psi_coef(i,iorder(j))
-         enddo
-         CI_eigenvectors_s2(j) = s2_eigvalues(iorder(j))
-       enddo
-       deallocate(e_array,iorder)
-     endif
      deallocate(s2_eigvalues)
      
    endif