Forgot file

src/Davidson/find_reference.irp.f

@@ -0,0 +1,41 @@
+subroutine find_reference(thresh,n_ref,result)
+  implicit none
+  double precision, intent(in) :: thresh
+  integer, intent(out) :: result(N_det),n_ref 
+  integer :: i,j,istate
+  double precision :: i_H_psi_array(1), E0, hii, norm
+  double precision :: de
+  integer(bit_kind), allocatable :: psi_ref_(:,:,:)
+  double precision, allocatable :: psi_ref_coef_(:,:)
+
+  allocate(psi_ref_coef_(N_det,1), psi_ref_(N_int,2,N_det))
+  n_ref = 1
+  result(1) = 1
+  istate = 1
+  psi_ref_coef_(1,1) = psi_coef(1,istate)
+  psi_ref_(:,:,1) = psi_det(:,:,1)
+  norm = psi_ref_coef_(1,1) * psi_ref_coef_(1,1)
+  call u_0_H_u_0(E0,psi_ref_coef_,n_ref,psi_ref_,N_int,1,size(psi_ref_coef_,1))
+  print *,  ''
+  print *,  'Reference determinants'
+  print *,  '======================'
+  print *,  ''
+  print *,  n_ref, ': E0 = ', E0 + nuclear_repulsion
+  call debug_det(psi_ref_(1,1,n_ref),N_int)
+  do i=2,N_det
+    call i_h_psi(psi_det(1,1,i),psi_ref_(1,1,1),psi_ref_coef_(1,istate),N_int, &
+      n_ref,size(psi_ref_coef_,1),1,i_H_psi_array)
+    call i_H_j(psi_det(1,1,i),psi_det(1,1,i),N_int,hii)
+    de = i_H_psi_array(istate)**2 / (E0 - hii)
+    if (dabs(de) > thresh) then
+      n_ref += 1
+      result(n_ref) = i
+      psi_ref_(:,:,n_ref) = psi_det(:,:,i)
+      psi_ref_coef_(n_ref,1) = psi_coef(i,istate)
+      call u_0_H_u_0(E0,psi_ref_coef_,n_ref,psi_ref_,N_int,1,size(psi_ref_coef_,1))
+      print *,  n_ref, ': E0 = ', E0 + nuclear_repulsion
+      call debug_det(psi_ref_(1,1,n_ref),N_int)
+    endif
+  enddo
+end
+