Update print in PT2

src/cipsi/pt2_stoch_routines.irp.f

@@ -131,7 +131,7 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in)
   PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns
   PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp psi_det_sorted
   PROVIDE psi_det_hii selection_weight pseudo_sym
-  PROVIDE list_act list_inact list_core list_virt list_del seniority_max
+  PROVIDE n_act_orb n_inact_orb n_core_orb n_virt_orb n_del_orb seniority_max
   PROVIDE excitation_beta_max  excitation_alpha_max excitation_max
 
   if (h0_type == 'CFG') then
@@ -290,9 +290,9 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in)
       call set_multiple_levels_omp(.False.)
 
 
-      print '(A)', '========== ======================= ===================== ===================== ==========='
+      print '(A)', '========== ==================== ================ ================ ================ ============= ==========='
-      print '(A)', ' Samples          Energy                Variance               Norm^2            Seconds'
+      print '(A)', ' Samples           Energy             PT2            Variance          Norm^2       Convergence    Seconds'
-      print '(A)', '========== ======================= ===================== ===================== ==========='
+      print '(A)', '========== ==================== ================ ================ ================ ============= ==========='
 
       PROVIDE global_selection_buffer
 
@@ -316,7 +316,8 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in)
       call end_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'pt2')
       call set_multiple_levels_omp(.True.)
 
-      print '(A)', '========== ======================= ===================== ===================== ==========='
+      print '(A)', '========== ==================== ================ ================ ================ ============= ==========='
+
 
       do k=1,N_states
         pt2_overlap(pt2_stoch_istate,k) = pt2_data % overlap(k,pt2_stoch_istate)
@@ -414,6 +415,17 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_
   double precision :: rss
   double precision, external :: memory_of_double, memory_of_int
 
+  character(len=20) :: format_str1, str_error1, format_str2, str_error2
+  character(len=20) :: format_str3, str_error3, format_str4, str_error4
+  character(len=20) :: format_value1, format_value2, format_value3, format_value4
+  character(len=20) :: str_value1, str_value2, str_value3, str_value4
+  character(len=20) :: str_conv
+  double precision  :: value1, value2, value3, value4
+  double precision  :: error1, error2, error3, error4
+  integer           :: size1,size2,size3,size4
+
+  double precision :: conv_crit
+
   sending =.False.
 
   rss  = memory_of_int(pt2_n_tasks_max*2+N_det_generators*2)
@@ -537,14 +549,60 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_
 
           if ((time - time1 > 1.d0) .or. (n==N_det_generators)) then
             time1 = time
-            print '(I10, X, F12.6, X, G10.3, X, F10.6, X, G10.3, X, F10.6, X, G10.3, X, F10.1)', c, &
+
-              pt2_data     % pt2(pt2_stoch_istate) +E, &
+            value1 =  pt2_data     % pt2(pt2_stoch_istate) + E
-              pt2_data_err % pt2(pt2_stoch_istate), &
+            error1 =  pt2_data_err % pt2(pt2_stoch_istate)           
-              pt2_data     % variance(pt2_stoch_istate), &
+            value2 =  pt2_data     % pt2(pt2_stoch_istate)
-              pt2_data_err % variance(pt2_stoch_istate), &
+            error2 =  pt2_data_err % pt2(pt2_stoch_istate)           
-              pt2_data     % overlap(pt2_stoch_istate,pt2_stoch_istate), &
+            value3 =  pt2_data     % variance(pt2_stoch_istate)           
-              pt2_data_err % overlap(pt2_stoch_istate,pt2_stoch_istate), &
+            error3 =  pt2_data_err % variance(pt2_stoch_istate)           
-              time-time0
+            value4 =  pt2_data     % overlap(pt2_stoch_istate,pt2_stoch_istate)           
+            error4 =  pt2_data_err % overlap(pt2_stoch_istate,pt2_stoch_istate)           
+            
+            ! Max size of the values (FX.Y) with X=size
+            size1 = 15
+            size2 = 12
+            size3 = 12
+            size4 = 12
+
+            ! To generate the format: number(error)
+            call format_w_error(value1,error1,size1,8,format_value1,str_error1)
+            call format_w_error(value2,error2,size2,8,format_value2,str_error2)
+            call format_w_error(value3,error3,size3,8,format_value3,str_error3)
+            call format_w_error(value4,error4,size4,8,format_value4,str_error4)
+
+            ! value > string with the right format
+            write(str_value1,'('//format_value1//')') value1
+            write(str_value2,'('//format_value2//')') value2
+            write(str_value3,'('//format_value3//')') value3
+            write(str_value4,'('//format_value4//')') value4
+
+            ! Convergence criterion
+            conv_crit = dabs(pt2_data_err % pt2(pt2_stoch_istate)) /    &
+                  (1.d-20 + dabs(pt2_data % pt2(pt2_stoch_istate)) )
+            write(str_conv,'(G10.3)') conv_crit
+
+            write(*,'(I10,X,X,A20,X,A16,X,A16,X,A16,X,A12,X,F10.1)') c,&
+            adjustl(adjustr(str_value1)//'('//str_error1(1:1)//')'),&
+            adjustl(adjustr(str_value2)//'('//str_error2(1:1)//')'),&
+            adjustl(adjustr(str_value3)//'('//str_error3(1:1)//')'),&
+            adjustl(adjustr(str_value4)//'('//str_error4(1:1)//')'),&
+            adjustl(str_conv),& 
+            time-time0
+
+            ! Old print
+            !print '(I10, X, F12.6, X, G10.3, X, F10.6, X, G10.3, X, F10.6, X, G10.3, X, F10.1,1pE16.6,1pE16.6)', c, &
+            !  pt2_data     % pt2(pt2_stoch_istate) +E, &
+            !  pt2_data_err % pt2(pt2_stoch_istate), &
+            !  pt2_data     % variance(pt2_stoch_istate), &
+            !  pt2_data_err % variance(pt2_stoch_istate), &
+            !  pt2_data     % overlap(pt2_stoch_istate,pt2_stoch_istate), &
+            !  pt2_data_err % overlap(pt2_stoch_istate,pt2_stoch_istate), &
+            !  time-time0, &
+            !  pt2_data     % pt2(pt2_stoch_istate), &
+            !  dabs(pt2_data_err % pt2(pt2_stoch_istate)) /    &
+            !      (1.d-20 + dabs(pt2_data % pt2(pt2_stoch_istate)) )
+
             if (stop_now .or. (                                      &
                   (do_exit .and. (dabs(pt2_data_err % pt2(pt2_stoch_istate)) /    &
                   (1.d-20 + dabs(pt2_data % pt2(pt2_stoch_istate)) ) <= relative_error))) ) then
@@ -844,6 +902,7 @@ END_PROVIDER
      if (tooth_width == 0.d0) then
        tooth_width = max(1.d-15,sum(tilde_w(pt2_n_0(t):pt2_n_0(t+1))))
      endif
+     ASSERT(tooth_width > 0.d0)
      do i=pt2_n_0(t)+1, pt2_n_0(t+1)
        pt2_w(i) = tilde_w(i) * pt2_W_T / tooth_width
      end do

src/cipsi/run_pt2_slave.irp.f

@@ -83,6 +83,7 @@ subroutine run_pt2_slave_small(thread,iproc,energy)
   buffer_ready = .False.
   n_tasks = 1
 
+!  sending = .False.
   done = .False.
   do while (.not.done)
 
@@ -116,13 +117,14 @@ subroutine run_pt2_slave_small(thread,iproc,energy)
     do k=1,n_tasks
       call pt2_alloc(pt2_data(k),N_states)
       b%cur = 0
-!      double precision :: time2
+!double precision :: time2
-!      call wall_time(time2)
+!call wall_time(time2)
       call select_connected(i_generator(k),energy,pt2_data(k),b,subset(k),pt2_F(i_generator(k)))
-!      call wall_time(time1)
+!call wall_time(time1)
-!      print *,  i_generator(1), time1-time2, n_tasks, pt2_F(i_generator(1))
+!print *,  i_generator(1), time1-time2, n_tasks, pt2_F(i_generator(1))
     enddo
     call wall_time(time1)
+!print *,  '-->', i_generator(1), time1-time0, n_tasks
 
     integer, external :: tasks_done_to_taskserver
     if (tasks_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id,n_tasks) == -1) then

src/cipsi/selection.irp.f

@@ -571,7 +571,6 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
   double precision, external :: diag_H_mat_elem_fock
   double precision :: E_shift
   double precision :: s_weight(N_states,N_states)
-  logical, external :: is_in_wavefunction
   PROVIDE dominant_dets_of_cfgs N_dominant_dets_of_cfgs
   do jstate=1,N_states
     do istate=1,N_states
@@ -751,7 +750,10 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
         if (delta_E < 0.d0) then
             tmp = -tmp
         endif
+
+        !e_pert(istate) = alpha_h_psi * alpha_h_psi / (E0(istate) - Hii)
         e_pert(istate) = 0.5d0 * (tmp - delta_E)
+
         if (dabs(alpha_h_psi) > 1.d-4) then
           coef(istate) = e_pert(istate) / alpha_h_psi
         else
@@ -864,6 +866,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
 !!!BEGIN_DEBUG
 !      ! To check if the pt2 is taking determinants already in the wf
 !      if (is_in_wavefunction(det(N_int,1),N_int)) then
+!        logical, external :: is_in_wavefunction
 !        print*, 'A determinant contributing to the pt2 is already in'
 !        print*, 'the wave function:'
 !        call  print_det(det(N_int,1),N_int)

src/cipsi/slave_cipsi.irp.f

@@ -311,7 +311,7 @@ subroutine run_slave_main
           if (mpi_master) then
             print *,  'Running PT2'
           endif
-          !$OMP PARALLEL PRIVATE(i) NUM_THREADS(nproc_target)
+          !$OMP PARALLEL PRIVATE(i) NUM_THREADS(nproc_target+1)
           i = omp_get_thread_num()
           call run_pt2_slave(0,i,pt2_e0_denominator)
           !$OMP END PARALLEL

src/utils/format_w_error.irp.f

@@ -0,0 +1,71 @@
+subroutine format_w_error(value,error,size_nb,max_nb_digits,format_value,str_error)
+
+  implicit none
+ 
+  BEGIN_DOC
+  ! Format for double precision, value(error)
+  END_DOC
+
+  ! in
+  ! | value         | double precision | value...      |
+  ! | error         | double precision | error...      |
+  ! | size_nb       | integer          | X in FX.Y     |
+  ! | max_nb_digits | integer          | Max Y in FX.Y |
+
+  ! out
+  ! | format_value  | character | string FX.Y for the format |
+  ! | str_error     | character | string of the error |
+
+  ! internal
+  ! | str_size      | character | size in string format                             |
+  ! | nb_digits     | integer   | number of digits Y in FX.Y depending of the error |
+  ! | str_nb_digits | character | nb_digits in string format                        |
+  ! | str_exp       | character | string of the value in exponential format         |
+
+  ! in
+  double precision, intent(in)   :: error, value
+  integer, intent(in)            :: size_nb, max_nb_digits
+
+  ! out
+  character(len=20), intent(out) :: str_error, format_value
+
+  ! internal
+  character(len=20)              :: str_size, str_nb_digits, str_exp
+  integer                        :: nb_digits
+
+  ! max_nb_digit: Y max
+  ! size_nb = Size of the double: X (FX.Y)
+  write(str_size,'(I3)') size_nb
+
+  ! Error
+  write(str_exp,'(1pE20.0)') error
+  str_error = trim(adjustl(str_exp))
+  
+  ! Number of digit: Y (FX.Y) from the exponent
+  str_nb_digits = str_exp(19:20)
+  read(str_nb_digits,*) nb_digits
+ 
+  ! If the error is 0d0
+  if (error <= 1d-16) then 
+    write(str_nb_digits,*) max_nb_digits
+  endif
+
+  ! If the error is too small 
+  if (nb_digits > max_nb_digits) then
+      write(str_nb_digits,*) max_nb_digits
+      str_error(1:1) = '0'
+  endif
+
+  ! If the error is too big (>= 0.5)
+  if (error >= 0.5d0) then
+    str_nb_digits = '1'
+    str_error(1:1) = '*'
+  endif
+
+  ! FX.Y,A1,A1,A1 for value(str_error)
+  !string = 'F'//trim(adjustl(str_size))//'.'//trim(adjustl(str_nb_digits))//',A1,A1,A1'
+
+  ! FX.Y just for the value 
+  format_value = 'F'//trim(adjustl(str_size))//'.'//trim(adjustl(str_nb_digits))
+
+end