Fixed energies of non-expected s2 (#9)

* Moved diag_algorithm in Davdison

configure

@@ -11,8 +11,6 @@ echo "QP_ROOT="$QP_ROOT
 
 # Check if the module to create new DFT functionals exists or not
 if [[ ! -d ${QP_ROOT}/src/new_functionals  ]] ; then
-   echo "${QP_ROOT}/src/new_functionals does NOT exists !!"  
-   echo "Creating a new one ..."
    cp -r ${QP_ROOT}/scripts/functionals/do_not_touch_func ${QP_ROOT}/src/new_functionals
 fi
 

ocaml/.gitignore

@@ -9,6 +9,7 @@ Input_ao_two_e_erf_ints.ml
 Input_ao_two_e_ints.ml
 Input_auto_generated.ml
 Input_becke_numerical_grid.ml
+Input_champ.ml
 Input_davidson.ml
 Input_density_for_dft.ml
 Input_determinants.ml
@@ -21,6 +22,7 @@ Input_nuclei.ml
 Input_perturbation.ml
 Input_pseudo.ml
 Input_scf_utils.ml
+Input_variance.ml
 qp_create_ezfio
 qp_create_ezfio.native
 qp_edit

src/davidson/EZFIO.cfg

@@ -40,3 +40,9 @@ doc: If |true|, use filter out all vectors with bad |S^2| values
 default: True
 interface: ezfio,provider,ocaml
 
+[n_det_max_full]
+type: Det_number_max
+doc: Maximum number of determinants where |H| is fully diagonalized 
+interface: ezfio,provider,ocaml
+default: 1000
+

src/davidson/diagonalization_hs2_dressed.irp.f

@@ -1,3 +1,20 @@
+BEGIN_PROVIDER [ character*(64), diag_algorithm ]
+  implicit none
+  BEGIN_DOC
+  ! Diagonalization algorithm (Davidson or Lapack)
+  END_DOC
+  if (N_det > N_det_max_full) then
+    diag_algorithm = "Davidson"
+  else
+    diag_algorithm = "Lapack"
+  endif
+
+  if (N_det < N_states) then
+    diag_algorithm = "Lapack"
+  endif
+END_PROVIDER
+
+
 BEGIN_PROVIDER [ integer, dressed_column_idx, (N_states) ]
  implicit none
  BEGIN_DOC
@@ -114,7 +131,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
   integer                        :: k_pairs, kl
 
   integer                        :: iter2, itertot
-  double precision, allocatable  :: y(:,:), h(:,:), lambda(:), s2(:)
+  double precision, allocatable  :: y(:,:), h(:,:), h_p(:,:), lambda(:), s2(:)
   real, allocatable              :: y_s(:,:)
   double precision, allocatable  :: s_(:,:), s_tmp(:,:)
   double precision               :: diag_h_mat_elem
@@ -264,6 +281,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
 
       ! Small
       h(N_st_diag*itermax,N_st_diag*itermax),                        &
+      h_p(N_st_diag*itermax,N_st_diag*itermax),                      &
       y(N_st_diag*itermax,N_st_diag*itermax),                        &
       s_(N_st_diag*itermax,N_st_diag*itermax),                       &
       s_tmp(N_st_diag*itermax,N_st_diag*itermax),                    &
@@ -408,27 +426,44 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
       ! Compute h_kl = <u_k | W_l> = <u_k| H |u_l>
       ! -------------------------------------------
 
+      call dgemm('T','N', shift2, shift2, sze,                       &
+          1.d0, U, size(U,1), W, size(W,1),                          &
+          0.d0, h, size(h_p,1))
+
       ! Penalty method
       ! --------------
 
       if (s2_eig) then
-        h = s_
+        h_p = s_
         do k=1,shift2
-          h(k,k) = h(k,k) + S_z2_Sz - expected_s2
+          h_p(k,k) = h_p(k,k) + S_z2_Sz - expected_s2
         enddo
         alpha = 0.1d0
+        h_p = h + alpha*h_p
       else
+        h_p = h
         alpha = 0.d0
       endif
 
-      call dgemm('T','N', shift2, shift2, sze,                       &
-          1.d0, U, size(U,1), W, size(W,1),                          &
-          alpha , h, size(h,1))
+      ! Diagonalize h_p
+      ! ---------------
 
-      ! Diagonalize h
-      ! -------------
+      call lapack_diag(lambda,y,h_p,size(h_p,1),shift2)
 
-      call lapack_diag(lambda,y,h,size(h,1),shift2)
+      ! Compute Energy for each eigenvector
+      ! -----------------------------------
+
+      call dgemm('N','N',shift2,shift2,shift2,                       &
+          1.d0, h, size(h,1), y, size(y,1),                          &
+          0.d0, s_tmp, size(s_tmp,1))
+
+      call dgemm('T','N',shift2,shift2,shift2,                       &
+          1.d0, y, size(y,1), s_tmp, size(s_tmp,1),                  &
+          0.d0, h, size(h,1))
+
+      do k=1,shift2
+        lambda(k) = h(k,k)
+      enddo
 
       ! Compute S2 for each eigenvector
       ! -------------------------------
@@ -441,8 +476,6 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
           1.d0, y, size(y,1), s_tmp, size(s_tmp,1),                  &
           0.d0, s_, size(s_,1))
 
-
-
       do k=1,shift2
         s2(k) = s_(k,k) + S_z2_Sz
       enddo

src/davidson/diagonalize_ci.irp.f

@@ -117,17 +117,24 @@ END_PROVIDER
        good_state_array = .False.
        call u_0_S2_u_0(s2_eigvalues,eigenvectors,N_det,psi_det,N_int,&
          N_det,size(eigenvectors,1))
-       do j=1,N_det
-         ! Select at least n_states states with S^2 values closed to "expected_s2"
-         if(dabs(s2_eigvalues(j)-expected_s2).le.0.5d0)then
-           i_state +=1
-           index_good_state_array(i_state) = j
-           good_state_array(j) = .True.
-         endif
-         if(i_state.eq.N_states) then
-           exit
-         endif
-       enddo
+       if (only_expected_s2) then
+          do j=1,N_det
+            ! Select at least n_states states with S^2 values closed to "expected_s2"
+            if(dabs(s2_eigvalues(j)-expected_s2).le.0.5d0)then
+              i_state +=1
+              index_good_state_array(i_state) = j
+              good_state_array(j) = .True.
+            endif
+            if(i_state.eq.N_states) then
+              exit
+            endif
+          enddo
+       else
+          do j=1,N_det
+             index_good_state_array(j) = j
+             good_state_array(j) = .True.
+          enddo
+       endif
        if(i_state .ne.0)then
          ! Fill the first "i_state" states that have a correct S^2 value
          do j = 1, i_state
@@ -185,6 +192,16 @@ END_PROVIDER
          CI_electronic_energy(j) = eigenvalues(j)
        enddo
      endif
+     do k=1,N_states_diag
+       CI_electronic_energy(k) = 0.d0
+       do j=1,N_det
+         do i=1,N_det
+           CI_electronic_energy(k) +=                                &
+               CI_eigenvectors(i,k) * CI_eigenvectors(j,k) *         &
+               H_matrix_all_dets(i,j)
+         enddo
+       enddo
+     enddo
      deallocate(eigenvectors,eigenvalues)
    endif
 

src/determinants/EZFIO.cfg

@@ -10,12 +10,6 @@ doc: Maximum number of determinants to be printed with the program print_wf
 interface: ezfio,provider,ocaml
 default: 10000
 
-[n_det_max_full]
-type: Det_number_max
-doc: Maximum number of determinants where |H| is fully diagonalized 
-interface: ezfio,provider,ocaml
-default: 1000
-
 [n_states]
 type: States_number
 doc: Number of states to consider

src/determinants/determinants.irp.f

@@ -1,22 +1,5 @@
 use bitmasks
 
-BEGIN_PROVIDER [ character*(64), diag_algorithm ]
-  implicit none
-  BEGIN_DOC
-  ! Diagonalization algorithm (Davidson or Lapack)
-  END_DOC
-  if (N_det > N_det_max_full) then
-    diag_algorithm = "Davidson"
-  else
-    diag_algorithm = "Lapack"
-  endif
-
-  if (N_det < N_states) then
-    diag_algorithm = "Lapack"
-  endif
-END_PROVIDER
-
-
 BEGIN_PROVIDER [ integer, N_det ]
   implicit none
   BEGIN_DOC

src/fci/40.fci.bats

@@ -59,7 +59,7 @@ function run_stoch() {
 @test "H2O2" { # 12.9214s
   qp set_file h2o2.ezfio
   qp set_mo_class --core="[1-2]" --act="[3-24]" --del="[25-38]"
-  run -151.004888189874 2.e-5
+  run -151.004888189874 4.e-5
 }
 
 @test "HBO" { # 13.3144s