Merge pull request #202 from Ydrnan/qp2_dev

update binom_func to avoid .999...

src/cipsi/pt2_stoch_routines.irp.f

@@ -290,9 +290,13 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in)
       call set_multiple_levels_omp(.False.)
 
 
-      print '(A)', '========== ======================= ===================== ===================== ==========='
+      ! old
-      print '(A)', ' Samples          Energy                Variance               Norm^2            Seconds'
+      !print '(A)', '========== ======================= ===================== ===================== ==========='
-      print '(A)', '========== ======================= ===================== ===================== ==========='
+      !print '(A)', ' Samples          Energy                Variance               Norm^2            Seconds'
+      !print '(A)', '========== ======================= ===================== ===================== ==========='
+      print '(A)', '========== ==================== ================ ================ ================ ============= ==========='
+      print '(A)', ' Samples           Energy             PT2            Variance          Norm^2       Convergence    Seconds'
+      print '(A)', '========== ==================== ================ ================ ================ ============= ==========='
 
       PROVIDE global_selection_buffer
 
@@ -316,7 +320,10 @@ 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)', '========== ======================= ===================== ===================== ==========='
+      ! old
+      !print '(A)', '========== ======================= ===================== ===================== ==========='
+      print '(A)', '========== ==================== ================ ================ ================ ============= ==========='
+
 
       do k=1,N_states
         pt2_overlap(pt2_stoch_istate,k) = pt2_data % overlap(k,pt2_stoch_istate)
@@ -414,6 +421,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 +555,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

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

src/utils/util.irp.f

@@ -37,6 +37,10 @@ double precision function binom_func(i,j)
   else
     binom_func = dexp( logfact(i)-logfact(j)-logfact(i-j) )
   endif
+
+  ! To avoid .999999 numbers
+  binom_func = binom_func + 0.5d0
+
 end