Fixed determinstic PT2

src/cipsi/cipsi.irp.f

@@ -6,22 +6,27 @@ subroutine run_cipsi
   END_DOC
   integer                        :: i,j,k
   double precision, allocatable  :: pt2(:), variance(:), norm(:), rpt2(:), zeros(:)
-  integer                        :: n_det_before, to_select
+  integer                        :: to_select
+  logical, external :: qp_stop
 
+  double precision :: threshold_generators_save
   double precision :: rss
   double precision, external :: memory_of_double
+  PROVIDE H_apply_buffer_allocated 
+
+  N_iter = 1
+  threshold_generators = 1.d0
+  SOFT_TOUCH threshold_generators
+
   rss = memory_of_double(N_states)*4.d0
   call check_mem(rss,irp_here)
 
-  N_iter = 1
   allocate (pt2(N_states), zeros(N_states), rpt2(N_states), norm(N_states), variance(N_states))
 
   double precision               :: hf_energy_ref
   logical                        :: has
   double precision               :: relative_error
 
-  PROVIDE H_apply_buffer_allocated
-
   relative_error=PT2_relative_error
 
   zeros = 0.d0
@@ -55,35 +60,35 @@ subroutine run_cipsi
     call save_wavefunction
   endif
 
-  n_det_before = 0
-
   double precision :: correlation_energy_ratio
-  double precision :: threshold_generators_save
-  threshold_generators_save = threshold_generators
   double precision :: error(N_states)
-  logical, external :: qp_stop
 
   correlation_energy_ratio = 0.d0
 
   do while (                                                         &
         (N_det < N_det_max) .and.                                    &
         (maxval(abs(rpt2(1:N_states))) > pt2_max) .and.               &
-        (maxval(variance(1:N_states)) > variance_max) .and.               &
+        (maxval(abs(variance(1:N_states))) > variance_max) .and.     &
-        (correlation_energy_ratio <= correlation_energy_ratio_max)    &
+        (correlation_energy_ratio <= correlation_energy_ratio_max)   &
         )
       write(*,'(A)')  '--------------------------------------------------------------------------------'
 
 
+    to_select = int(sqrt(dble(N_states))*dble(N_det)*selection_factor)
+    to_select = max(N_states_diag, to_select)
     if (do_pt2) then
       pt2 = 0.d0
       variance = 0.d0
       norm = 0.d0
+      threshold_generators_save = threshold_generators
       threshold_generators = 1.d0
       SOFT_TOUCH threshold_generators
       call ZMQ_pt2(psi_energy_with_nucl_rep,pt2,relative_error,error, variance, &
         norm, 0) ! Stochastic PT2
       threshold_generators = threshold_generators_save
       SOFT_TOUCH threshold_generators
+    else
+    call ZMQ_selection(to_select, pt2, variance, norm)
     endif
 
     do k=1,N_states
@@ -95,7 +100,7 @@ subroutine run_cipsi
     correlation_energy_ratio = min(1.d0,correlation_energy_ratio)
 
     call write_double(6,correlation_energy_ratio, 'Correlation ratio')
-    call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm,N_det,N_occ_pattern,N_states,psi_s2)
+    call print_summary(psi_energy_with_nucl_rep,pt2,error,variance,norm,N_det,N_occ_pattern,N_states,psi_s2)
 
     call save_energy(psi_energy_with_nucl_rep, rpt2)
 
@@ -105,11 +110,9 @@ subroutine run_cipsi
 
     if (qp_stop()) exit
 
-    n_det_before = N_det
+    ! Add selected determinants
-    to_select = int(sqrt(dble(N_states))*dble(N_det)*selection_factor)
+    call copy_H_apply_buffer_to_wf()
-    to_select = max(N_states_diag, to_select)
+!    call save_wavefunction
-    to_select = max(to_select,1)
-    call ZMQ_selection(to_select, pt2, variance, norm)
 
     PROVIDE  psi_coef
     PROVIDE  psi_det
@@ -119,10 +122,6 @@ subroutine run_cipsi
     call save_wavefunction
     call save_energy(psi_energy_with_nucl_rep, zeros)
     if (qp_stop()) exit 
-print *,  (N_det < N_det_max)
-print *,  (maxval(abs(rpt2(1:N_states))) > pt2_max)
-print *,  (maxval(variance(1:N_states)) > variance_max)
-print *,  (correlation_energy_ratio <= correlation_energy_ratio_max)
   enddo
 
   if (.not.qp_stop()) then
@@ -152,8 +151,8 @@ print *,  (correlation_energy_ratio <= correlation_energy_ratio_max)
       rpt2(k) = pt2(k)/(1.d0 + norm(k))
     enddo
 
-    call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm,N_det,N_occ_pattern,N_states,psi_s2)
     call save_energy(psi_energy_with_nucl_rep, rpt2)
+    call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm,N_det,N_occ_pattern,N_states,psi_s2)
     call save_iterations(psi_energy_with_nucl_rep(1:N_states),rpt2,N_det)
     call print_extrapolated_energy()
   endif

src/cipsi/run_selection_slave.irp.f

@@ -63,7 +63,12 @@ subroutine run_selection_slave(thread,iproc,energy)
         call create_selection_buffer(bsize, bsize*2, buf)
         buffer_ready = .True.
       else
-        ASSERT (N == buf%N)
+        if (N /= buf%N) then
+          print *, 'N=', N
+          print *, 'buf%N=', buf%N
+          print *, 'bug in ', irp_here
+          stop '-1'
+        end if
       end if
       call select_connected(i_generator,energy,pt2,variance,norm,buf,subset,pt2_F(i_generator))
     endif

src/cipsi/zmq_selection.irp.f

@@ -129,8 +129,6 @@ subroutine ZMQ_selection(N_in, pt2, variance, norm)
       call make_selection_buffer_s2(b)
     endif
     call fill_H_apply_buffer_no_selection(b%cur,b%det,N_int,0)
-    call copy_H_apply_buffer_to_wf()
-    call save_wavefunction
   endif
   call delete_selection_buffer(b)
   do k=1,N_states