Merge pull request #8 from QuantumPackage/dev-stable-tc-scf

Dev stable tc scf

src/davidson/diagonalization_hs2_dressed.irp.f

@@ -461,6 +461,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
        integer :: lwork, info
        double precision, allocatable :: work(:)
 
+
 !       y = h
        y = h_p
        lwork = -1

src/tc_bi_ortho/dav_h_tc_s2.irp.f

@@ -1,7 +1,7 @@
 
 ! ---
 
-subroutine davidson_hs2_nonsym_b1space(u_in, H_jj, s2_out,energies, sze, N_st, N_st_diag_in, converged, hcalc)
+subroutine davidson_hs2_nonsym_b1space(u_in, H_jj, s2_out,energies, sze, N_st, N_st_diag_in, n_it_max_dav, converged, hcalc)
 
   use mmap_module
 
@@ -21,11 +21,17 @@ subroutine davidson_hs2_nonsym_b1space(u_in, H_jj, s2_out,energies, sze, N_st, N
   ! Initial guess vectors are not necessarily orthonormal
   !
   ! hcalc subroutine to compute W = H U (see routine hcalc_template for template of input/output)
+  ! 
+  ! !!! WARNING !!! IT SEEMS THAT IF THE NUMBER OF MACRO ITERATIONS EXCEEDS n_it_max_dav, 
+  ! 
+  !                 THE RECONTRACTION IS WRONG. YOU SHOULD CONSIDER CALLING MULTIPLE TIME THE ROUTINE
+  !
+  !                 SEE FOR INSTANCE IN tc_bi_ortho/tc_h_eigvectors.irp.f
   END_DOC
 
   implicit none
 
-  integer,           intent(in)   :: sze, N_st, N_st_diag_in
+  integer,           intent(in)   :: sze, N_st, N_st_diag_in, n_it_max_dav
   double precision,  intent(in)   :: H_jj(sze)
   logical,          intent(inout) :: converged
   double precision, intent(inout) :: u_in(sze,N_st_diag_in)
@@ -246,7 +252,9 @@ subroutine davidson_hs2_nonsym_b1space(u_in, H_jj, s2_out,energies, sze, N_st, N
 
   itertot = 0
 
-  do while (.not.converged)
+!  do while (.not.converged.or.itertot.le.n_it_max_dav)
+  integer :: iiii
+  do iiii = 1, n_it_max_dav
 
     itertot = itertot + 1
     if(itertot == 8) then
@@ -522,7 +530,7 @@ subroutine davidson_hs2_nonsym_b1space(u_in, H_jj, s2_out,energies, sze, N_st, N
         enddo
       endif
     enddo
-
+   if(converged)exit
   enddo ! loop over while
 
   ! ---

src/tc_bi_ortho/tc_h_eigvectors.irp.f

@@ -223,7 +223,11 @@ end
       call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,i), psi_det(1,1,i), N_int, H_jj(i))
     enddo
  !!!! Preparing the left-eigenvector
+    print*,'---------------------------------'
+    print*,'---------------------------------'
     print*,'Computing the left-eigenvector '
+    print*,'---------------------------------'
+    print*,'---------------------------------'
     vec_tmp = 0.d0
     do istate = 1, N_states
      vec_tmp(1:N_det,istate) = psi_l_coef_bi_ortho(1:N_det,istate)
@@ -233,12 +237,21 @@ end
     enddo
 !    call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_left_tc_bi_orth, N_det, n_states, n_states_diag, converged, htcdag_bi_ortho_calc_tdav)
 !    call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_left_tc_bi_orth, N_det, n_states, n_states_diag, converged, H_tc_dagger_u_0_opt)
-    call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2_eigvec_tc_bi_orth, eigval_left_tc_bi_orth, N_det, n_states, n_states_diag, converged, H_tc_s2_dagger_u_0_opt)
+    integer :: n_it_max
+    n_it_max = 1
+    converged = .False.
+    do while (.not.converged)
+     call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2_eigvec_tc_bi_orth, eigval_left_tc_bi_orth, N_det, n_states, n_states_diag, n_it_max, converged, H_tc_s2_dagger_u_0_opt)
+    enddo
     do istate = 1, N_states
      leigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate)
     enddo
 
+    print*,'---------------------------------'
+    print*,'---------------------------------'
     print*,'Computing the right-eigenvector '
+    print*,'---------------------------------'
+    print*,'---------------------------------'
  !!!! Preparing the right-eigenvector
     vec_tmp = 0.d0
     do istate = 1, N_states
@@ -249,7 +262,10 @@ end
     enddo
 !    call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_right_tc_bi_orth, N_det, n_states, n_states_diag, converged, htc_bi_ortho_calc_tdav)
 !    call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_right_tc_bi_orth, N_det, n_states, n_states_diag, converged, H_tc_u_0_opt)
-    call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2_eigvec_tc_bi_orth, eigval_right_tc_bi_orth, N_det, n_states, n_states_diag, converged, H_tc_s2_dagger_u_0_opt)
+    converged = .False.
+    do while (.not.converged)
+     call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2_eigvec_tc_bi_orth, eigval_right_tc_bi_orth, N_det, n_states, n_states_diag, n_it_max, converged, H_tc_s2_dagger_u_0_opt)
+    enddo
     do istate = 1, N_states
      reigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate)
     enddo

src/tc_bi_ortho/test_s2_tc.irp.f

@@ -105,7 +105,9 @@ subroutine routine_test_s2_davidson
  do istate = 1, N_states
   leigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate)
  enddo
- call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2, energies, N_det, n_states, n_states_diag, converged, H_tc_s2_dagger_u_0_opt)
+ integer :: n_it_max
+ n_it_max = 1
+ call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2, energies, N_det, n_states, n_states_diag, n_it_max, converged, H_tc_s2_dagger_u_0_opt)
  double precision, allocatable :: v_0_new(:,:),s_0_new(:,:)
  integer :: sze,N_st
  logical           :: do_right 
@@ -138,7 +140,8 @@ subroutine routine_test_s2_davidson
  do istate = 1, N_states
   leigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate)
  enddo
- call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2, energies, N_det, n_states, n_states_diag, converged, H_tc_s2_u_0_opt)
+ n_it_max = 1
+ call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2, energies, N_det, n_states, n_states_diag, n_it_max, converged, H_tc_s2_u_0_opt)
  sze = N_det
  N_st = 1
  do_right = .True.