Davidson convergence criteria

scripts/generate_h_apply.py

@@ -163,6 +163,7 @@ class H_apply(object):
       call copy_h_apply_buffer_to_wf
       selection_criterion_min = selection_criterion_min*0.1d0
       selection_criterion = selection_criterion_min
+      call remove_small_contributions
       """
       self.data["keys_work"] = """
       e_2_pert_buffer = 0.d0

src/Dets/davidson.irp.f

@@ -112,8 +112,11 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun
   double precision               :: residual_norm(N_st)
   character*(16384)              :: write_buffer
   double precision               :: to_print(2,N_st)
+  double precision               :: cpu, wall
   
   call write_time(iunit)
+  call wall_time(wall)
+  call cpu_time(cpu)
   write(iunit,'(A)') ''
   write(iunit,'(A)') 'Davidson Diagonalization'
   write(iunit,'(A)') '------------------------'
@@ -272,12 +275,12 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun
       enddo
 
       write(iunit,'(X,I3,X,100(X,F16.10,X,E16.6))'), iter, to_print(:,1:N_st)
-      
-      call davidson_converged(lambda,residual_norm,N_st,converged)
+      call davidson_converged(lambda,residual_norm,wall,iter,cpu,N_st,converged)
       if (converged) then
         exit
       endif
       
+      
       ! Davidson step
       ! -------------
       
@@ -321,6 +324,7 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun
       !pause
       !END DEBUG
 
+
     enddo
 
     if (.not.converged) then
@@ -368,21 +372,22 @@ end
 &BEGIN_PROVIDER [ double precision, davidson_threshold ]
  implicit none
  BEGIN_DOC
- ! Can be : [  energy  | residual | both ]
+ ! Can be : [  energy  | residual | both | wall_time | cpu_time | iterations ]
  END_DOC
  davidson_criterion = 'both'
  davidson_threshold = 1.d-8
 END_PROVIDER
 
-subroutine davidson_converged(energy,residual,N_st,converged)
+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
+  integer, intent(in) :: N_st, iterations
   logical, intent(out) :: converged
   double precision, intent(in) :: energy(N_st), residual(N_st)
-  double precision :: E(N_st)
+  double precision, intent(in) :: wall, cpu
+  double precision :: E(N_st), time
   double precision, allocatable, save :: energy_old(:)
 
   if (.not.allocated(energy_old)) then
@@ -397,7 +402,16 @@ subroutine davidson_converged(energy,residual,N_st,converged)
   else if (davidson_criterion == 'residual') then
     converged = dabs(maxval(residual(1:N_st))) < davidson_threshold 
   else if (davidson_criterion == 'both') then
-    converged = dabs(maxval(residual(1:N_st))) + dabs(maxval(E(1:N_st)) ) < davidson_threshold  
+    converged = dabs(maxval(residual(1:N_st))) + dabs(maxval(E(1:N_st)) ) &
+       < davidson_threshold  
+  else if (davidson_criterion == 'wall_time') then
+    call wall_time(time)
+    converged = time - wall > davidson_threshold
+  else if (davidson_criterion == 'cpu_time') then
+    call cpu_time(time)
+    converged = time - cpu > davidson_threshold
+  else if (davidson_criterion == 'iterations') then
+    converged = iterations >= int(davidson_threshold)
   endif
   converged = converged.or.abort_here
 end

src/Dets/utils.irp.f

@@ -14,3 +14,36 @@ BEGIN_PROVIDER [ double precision, H_matrix_all_dets,(N_det,N_det) ]
  enddo
 END_PROVIDER
 
+subroutine remove_small_contributions
+  implicit none
+  BEGIN_DOC
+!  Remove determinants with small contributions
+  END_DOC
+  integer :: i,j,k, N_removed
+  logical keep
+  N_removed = 0
+  do i=N_det,1,-1
+    keep = .False.
+    do j=1,N_states
+      keep = keep .or. (dabs(psi_coef(i,j)) > selection_criterion_min)
+    enddo
+    if (.not.keep) then
+      do k=i+1,N_det
+        do j=1,N_int
+           psi_det(j,1,k-1) = psi_det(j,1,k)
+           psi_det(j,2,k-1) = psi_det(j,2,k)
+        enddo
+      enddo
+      do j=1,N_states
+        do k=i+1,N_det
+           psi_coef(k-1,j) = psi_coef(k,j)
+        enddo
+      enddo
+      N_removed += 1
+    endif
+  enddo
+  if (N_removed > 0) then
+    N_det -= N_removed
+    call write_int(output_dets,N_removed, 'Removed determinants')
+  endif
+end

src/Perturbation/cisd_selection.irp.f

@@ -10,6 +10,8 @@ program cisd
   allocate (pt2(N_st), norm_pert(N_st))
   
   pt2 = 1.d0
+! davidson_criterion = 'wall_time'
+! davidson_threshold = 1.d0
   do while (maxval(abs(pt2(1:N_st))) > 1.d-6)
     E_old = CI_energy(1)
     call H_apply_cisd_selection(pt2, norm_pert, H_pert_diag,  N_st)