Dynamic Davidson threshold

src/cipsi/selection.irp.f

@@ -70,7 +70,10 @@ subroutine update_pt2_and_variance_weights(pt2, variance, norm, N_st)
     variance_match_weight(k) = product(memo_variance(k,:))
   enddo
 
-  SOFT_TOUCH pt2_match_weight variance_match_weight
+  threshold_davidson_pt2 = min(1.d-6, &
+     max(threshold_davidson, 1.e-1 * PT2_relative_error * minval(abs(rpt2(1:N_states)))) )
+
+  SOFT_TOUCH pt2_match_weight variance_match_weight threshold_davidson_pt2
 end
 
 

src/davidson/EZFIO.cfg

@@ -1,9 +1,15 @@
 [threshold_davidson]
 type: Threshold
-doc: Thresholds of Davidson's algorithm
+doc: Thresholds of Davidson's algorithm if threshold_davidson_from_pt2 is false.
 interface: ezfio,provider,ocaml
 default: 1.e-10
 
+[threshold_davidson_from_pt2]
+type: logical
+doc: Thresholds of Davidson's algorithm is set to E(rPT2)*threshold_davidson_from_pt2
+interface: ezfio,provider,ocaml
+default: true
+
 [n_states_diag]
 type: States_number
 doc: Controls the number of states to consider during the Davdison diagonalization. The number of states is n_states * n_states_diag

src/davidson/diagonalization_hs2_dressed.irp.f

@@ -14,6 +14,16 @@ BEGIN_PROVIDER [ character*(64), diag_algorithm ]
   endif
 END_PROVIDER
 
+BEGIN_PROVIDER [ double precision, threshold_davidson_pt2 ]
+ implicit none
+ BEGIN_DOC
+ ! Threshold of Davidson's algorithm, using PT2 as a guide
+ END_DOC
+ threshold_davidson_pt2 = threshold_davidson
+
+END_PROVIDER
+
+
 
 BEGIN_PROVIDER [ integer, dressed_column_idx, (N_states) ]
  implicit none
@@ -171,11 +181,9 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
   endif
   overlap = 0.d0
 
-  PROVIDE nuclear_repulsion expected_s2 psi_bilinear_matrix_order psi_bilinear_matrix_order_reverse
+  PROVIDE nuclear_repulsion expected_s2 psi_bilinear_matrix_order psi_bilinear_matrix_order_reverse threshold_davidson_pt2
 
   call write_time(6)
-  call wall_time(wall)
-  call cpu_time(cpu)
   write(6,'(A)') ''
   write(6,'(A)') 'Davidson Diagonalization'
   write(6,'(A)') '------------------------'
@@ -224,6 +232,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
     else if (m==1.and.disk_based_davidson) then
       m=0
       disk_based = .True.
+      itermax = 6
     else
       nproc_target = nproc_target - 1
     endif
@@ -596,11 +605,17 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
       else
         write(*,'(1X,I3,1X,100(1X,F16.10,1X,F11.6,1X,E11.3))') iter-1, to_print(1:3,1:N_st)
       endif
-      call davidson_converged(lambda,residual_norm,wall,iter,cpu,N_st,converged)
+
+      ! Check convergence
+      if (iter > 1) then
+          converged = dabs(maxval(residual_norm(1:N_st))) < threshold_davidson_pt2
+      endif   
+      
+
       do k=1,N_st
         if (residual_norm(k) > 1.e8) then
-        print *,  ''
+          print *, 'Davidson failed'
-          stop 'Davidson failed'
+          stop -1
         endif
       enddo
       if (converged) then

src/davidson/parameters.irp.f

@@ -1,49 +0,0 @@
-BEGIN_PROVIDER [ character(64), davidson_criterion ]
- implicit none
- BEGIN_DOC
- ! Can be : [  energy  | residual | both | wall_time | cpu_time | iterations ]
- END_DOC
- davidson_criterion = 'residual'
-END_PROVIDER
-
-subroutine davidson_converged(energy,residual,wall,iterations,cpu,N_st,converged)
-  implicit none
-  BEGIN_DOC
-! True if the Davidson algorithm is converged
-  END_DOC
-  integer, intent(in) :: N_st, iterations
-  logical, intent(out) :: converged
-  double precision, intent(in) :: energy(N_st), residual(N_st)
-  double precision, intent(in) :: wall, cpu
-  double precision :: E(N_st), time
-  double precision, allocatable, save :: energy_old(:)
-
-  if (iterations < 2) then
-    converged = .False.
-    return
-  endif
-
-  if (.not.allocated(energy_old)) then
-    allocate(energy_old(N_st))
-    energy_old = 0.d0
-  endif
-
-  E = energy - energy_old
-  energy_old = energy
-  if (davidson_criterion == 'energy') then
-    converged = dabs(maxval(E(1:N_st))) < threshold_davidson
-  else if (davidson_criterion == 'residual') then
-    converged = dabs(maxval(residual(1:N_st))) < threshold_davidson
-  else if (davidson_criterion == 'both') then
-    converged = dabs(maxval(residual(1:N_st))) + dabs(maxval(E(1:N_st)) ) &
-       < threshold_davidson
-  else if (davidson_criterion == 'wall_time') then
-    call wall_time(time)
-    converged = time - wall > threshold_davidson
-  else if (davidson_criterion == 'cpu_time') then
-    call cpu_time(time)
-    converged = time - cpu > threshold_davidson
-  else if (davidson_criterion == 'iterations') then
-    converged = iterations >= int(threshold_davidson)
-  endif
-end

src/ezfio_files/qp_stop.irp.f

@@ -23,6 +23,7 @@ logical function qp_stop()
     INQUIRE(FILE=trim(qp_kill_filename), EXIST=qp_stop)
     if (qp_stop) then
       qp_stop_variable = 1
+      ! qp_stop is true
       return
     endif