Moved print_e_components

ocaml/Input_determinants_by_hand.ml

@@ -81,9 +81,6 @@ end = struct
   ;;
 
   let write_n_det n =
-    let n_det_old =
-      Ezfio.get_determinants_n_det ()
-    in
     Det_number.to_int n
     |> Ezfio.set_determinants_n_det
   ;;

src/cipsi/cipsi.irp.f

@@ -65,7 +65,7 @@ subroutine run_cipsi
 
   do while (                                                         &
         (N_det < N_det_max) .and.                                    &
-        (maxval(abs(pt2(1:N_states))) > pt2_max) .and.               &
+        (maxval(abs(rpt2(1:N_states))) > pt2_max) .and.               &
         (correlation_energy_ratio <= correlation_energy_ratio_max)    &
         )
       write(*,'(A)')  '--------------------------------------------------------------------------------'
@@ -92,7 +92,7 @@ subroutine run_cipsi
     call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm,N_det,N_occ_pattern,N_states,psi_s2)
 
     do k=1,N_states
-      rpt2(:) = pt2(:)/(1.d0 + norm(k))
+      rpt2(k) = pt2(k)/(1.d0 + norm(k))
     enddo
     call save_energy(psi_energy_with_nucl_rep, rpt2)
 
@@ -103,9 +103,8 @@ subroutine run_cipsi
     if (qp_stop()) exit 
 
     n_det_before = N_det
-    to_select = N_det
+    to_select = N_det*int(sqrt(dble(N_states)))*selection_factor
     to_select = max(N_states_diag, to_select)
-!    to_select = min(to_select, N_det_max-n_det_before)
     call ZMQ_selection(to_select, pt2, variance, norm)
 
     PROVIDE  psi_coef
@@ -137,7 +136,7 @@ subroutine run_cipsi
         norm,0) ! Stochastic PT2
       SOFT_TOUCH threshold_generators
       do k=1,N_states
-        rpt2(:) = pt2(:)/(1.d0 + norm(k))
+        rpt2(k) = pt2(k)/(1.d0 + norm(k))
       enddo
       call save_energy(psi_energy_with_nucl_rep, pt2)
     endif
@@ -148,7 +147,7 @@ subroutine run_cipsi
 
 
     do k=1,N_states
-      rpt2(:) = pt2(:)/(1.d0 + norm(k))
+      rpt2(k) = pt2(k)/(1.d0 + norm(k))
     enddo
     call save_energy(psi_energy_with_nucl_rep, rpt2)
 

src/cipsi/pt2_stoch_routines.irp.f

@@ -129,13 +129,13 @@ subroutine ZMQ_pt2(E, pt2,relative_error, error, variance, norm, N_in)
   PROVIDE psi_bilinear_matrix_rows psi_det_sorted_order psi_bilinear_matrix_order
   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 N_generators_bitmask
+  PROVIDE psi_det_hii N_generators_bitmask selection_weight pseudo_sym
 
   if (h0_type == 'SOP') then
     PROVIDE psi_occ_pattern_hii det_to_occ_pattern
   endif
 
-  if (N_det < max(1000,N_states)) then
+  if (N_det < max(4,N_states)) then
     pt2=0.d0
     variance=0.d0
     norm=0.d0
@@ -333,13 +333,7 @@ subroutine ZMQ_pt2(E, pt2,relative_error, error, variance, norm, N_in)
     pt2(k) = 0.d0
   enddo
 
-  ! Adjust PT2 weights for next selection
-  double precision :: pt2_avg
-  pt2_avg = sum(pt2) / dble(N_states)
-  do k=1,N_states
-    pt2_match_weight(k) *= (pt2(k)/pt2_avg)**2
-  enddo
-  SOFT_TOUCH pt2_match_weight
+  call update_pt2_and_variance_weights(pt2, variance, norm, N_states)
 
 end subroutine
 

src/cipsi/run_selection_slave.irp.f

@@ -25,8 +25,8 @@ subroutine run_selection_slave(thread,iproc,energy)
   PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique
   PROVIDE psi_bilinear_matrix_rows psi_det_sorted_order psi_bilinear_matrix_order
   PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns
-  PROVIDE psi_bilinear_matrix_transp_order N_int pt2_F
-  PROVIDE psi_selectors_coef_transp psi_det_sorted
+  PROVIDE psi_bilinear_matrix_transp_order N_int pt2_F pseudo_sym
+  PROVIDE psi_selectors_coef_transp psi_det_sorted weight_selection
 
 
   zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
@@ -230,6 +230,8 @@ subroutine pull_selection_results(zmq_socket_pull, pt2, variance, norm, val, det
       endif
   else
       pt2(:) = 0.d0
+      variance(:) = 0.d0
+      norm(:) = 0.d0
   endif
 
   rc = f77_zmq_recv( zmq_socket_pull, ntask, 4, 0)

src/cipsi/selection.irp.f

@@ -6,15 +6,88 @@ BEGIN_PROVIDER [ double precision, pt2_match_weight, (N_states) ]
  ! Weights adjusted along the selection to make the PT2 contributions
  ! of each state coincide.
  END_DOC
- pt2_match_weight = 1.d0
+ pt2_match_weight(:) = 1.d0
 END_PROVIDER
 
+BEGIN_PROVIDER [ double precision, variance_match_weight, (N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! Weights adjusted along the selection to make the variances 
+ ! of each state coincide.
+ END_DOC
+ variance_match_weight(:) = 1.d0
+END_PROVIDER
+
+subroutine update_pt2_and_variance_weights(pt2, variance, norm, N_st)
+  implicit none
+  BEGIN_DOC
+! Updates the rPT2- and Variance- matching weights.
+  END_DOC
+  integer, intent(in)          :: N_st
+  double precision, intent(in) :: pt2(N_st)
+  double precision, intent(in) :: variance(N_st)
+  double precision, intent(in) :: norm(N_st)
+
+  double precision :: avg, rpt2(N_st), element, dt, x
+  integer          :: k
+
+  dt = n_iter * selection_factor
+
+  do k=1,N_st
+    rpt2(k) = pt2(k)/(1.d0 + norm(k))                                                     
+  enddo
+
+  avg = sum(rpt2(1:N_st)) / dble(N_st)
+  do k=1,N_st
+    element = exp(dt*(rpt2(k)/avg -1.d0))
+    element = min(1.2d0 , element)
+    pt2_match_weight(k) *= element
+  enddo
+
+  avg = sum(variance(1:N_st)) / dble(N_st)
+  do k=1,N_st
+    element = exp(dt*(variance(k)/avg -1.d0))
+    element = min(1.2d0 , element)
+    variance_match_weight(k) *= element
+  enddo
+
+  SOFT_TOUCH pt2_match_weight variance_match_weight
+end
+
+
 BEGIN_PROVIDER [ double precision, selection_weight, (N_states) ]
    implicit none
    BEGIN_DOC
    ! Weights used in the selection criterion
    END_DOC
-   selection_weight(1:N_states) = c0_weight(1:N_states) * pt2_match_weight(1:N_states)
+   select case (weight_selection)
+
+     case (0)
+      print *,  'Using input weights in selection'
+      selection_weight(1:N_states) = state_average_weight(1:N_states)
+
+     case (1)
+      print *,  'Using 1/c_max^2 weight in selection'
+      selection_weight(1:N_states) = c0_weight(1:N_states) 
+
+     case (2)
+      print *,  'Using pt2-matching weight in selection'
+      selection_weight(1:N_states) = c0_weight(1:N_states) * pt2_match_weight(1:N_states)
+
+     case (3)
+      print *,  'Using variance-matching weight in selection'
+      selection_weight(1:N_states) = c0_weight(1:N_states) * variance_match_weight(1:N_states)
+
+     case (4)
+      print *,  'Using variance- and pt2-matching weights in selection'
+      selection_weight(1:N_states) = c0_weight(1:N_states) * variance_match_weight(1:N_states) * pt2_match_weight(1:N_states)
+
+     case (5)
+      print *,  'Using variance-matching weight in selection'
+      selection_weight(1:N_states) = c0_weight(1:N_states) * variance_match_weight(1:N_states)
+
+    end select
+
 END_PROVIDER
 
 
@@ -621,11 +694,13 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
         variance(istate) = variance(istate) + alpha_h_psi * alpha_h_psi
         norm(istate) = norm(istate) + coef * coef
 
-!        if (h0_type == "Variance") then
-!          sum_e_pert = sum_e_pert - alpha_h_psi * alpha_h_psi * selection_weight(istate)
-!        else
+        if (weight_selection /= 5) then
+          ! Energy selection
           sum_e_pert = sum_e_pert + e_pert * selection_weight(istate)
-!        endif
+        else
+          ! Variance selection
+          sum_e_pert = sum_e_pert - alpha_h_psi * alpha_h_psi * selection_weight(istate)
+        endif
       end do
       if(pseudo_sym)then
        if(dabs(mat(1, p1, p2)).lt.thresh_sym)then 

src/cipsi/slave_cipsi.irp.f

@@ -17,7 +17,7 @@ subroutine provide_everything
   PROVIDE H_apply_buffer_allocated mo_two_e_integrals_in_map psi_det_generators psi_coef_generators psi_det_sorted_bit psi_selectors n_det_generators n_states generators_bitmask zmq_context N_states_diag
   PROVIDE pt2_e0_denominator mo_num N_int ci_energy mpi_master zmq_state zmq_context
   PROVIDE psi_det psi_coef threshold_generators state_average_weight
-  PROVIDE N_det_selectors pt2_stoch_istate N_det
+  PROVIDE N_det_selectors pt2_stoch_istate N_det selection_weight pseudo_sym
 end
 
 subroutine run_slave_main

src/cipsi/stochastic_cipsi.irp.f

@@ -70,7 +70,7 @@ subroutine run_stochastic_cipsi
       write(*,'(A)')  '--------------------------------------------------------------------------------'
 
 
-    to_select = N_det*int(sqrt(dble(N_states)))
+    to_select = N_det*int(sqrt(dble(N_states)))*selection_factor
     to_select = max(N_states_diag, to_select)
 
     pt2 = 0.d0
@@ -88,7 +88,7 @@ subroutine run_stochastic_cipsi
     call print_summary(psi_energy_with_nucl_rep,pt2,error,variance,norm,N_det,N_occ_pattern,N_states,psi_s2)
 
     do k=1,N_states
-      rpt2(:) = pt2(:)/(1.d0 + norm(k))
+      rpt2(k) = pt2(k)/(1.d0 + norm(k))
     enddo
     call save_energy(psi_energy_with_nucl_rep, rpt2)
 
@@ -128,7 +128,7 @@ subroutine run_stochastic_cipsi
       norm,0) ! Stochastic PT2
 
     do k=1,N_states
-      rpt2(:) = pt2(:)/(1.d0 + norm(k))
+      rpt2(k) = pt2(k)/(1.d0 + norm(k))
     enddo
     call save_energy(psi_energy_with_nucl_rep, rpt2)
 

src/cipsi/zmq_selection.irp.f

@@ -21,7 +21,8 @@ subroutine ZMQ_selection(N_in, pt2, variance, norm)
     PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique
     PROVIDE psi_bilinear_matrix_rows psi_det_sorted_order psi_bilinear_matrix_order
     PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns
-    PROVIDE psi_bilinear_matrix_transp_order
+    PROVIDE psi_bilinear_matrix_transp_order selection_weight pseudo_sym
+
 
     call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,'selection')
 
@@ -85,7 +86,11 @@ subroutine ZMQ_selection(N_in, pt2, variance, norm)
   endif
 
   integer :: nproc_target
-  nproc_target = nproc
+  if (N_det < 3*nproc) then
+    nproc_target = N_det/4
+  else
+    nproc_target = nproc
+  endif
   double precision :: mem
   mem = 8.d0 * N_det * (N_int * 2.d0 * 3.d0 +  3.d0 + 5.d0) / (1024.d0**3)
   call write_double(6,mem,'Estimated memory/thread (Gb)')
@@ -131,13 +136,7 @@ subroutine ZMQ_selection(N_in, pt2, variance, norm)
     norm(k) = norm(k) * f(k)
   enddo
 
-  ! Adjust PT2 weights for next selection
-  double precision :: pt2_avg
-  pt2_avg = sum(pt2) / dble(N_states)
-  do k=1,N_states
-    pt2_match_weight(k) *= (pt2(k)/pt2_avg)**2
-  enddo
-  SOFT_TOUCH pt2_match_weight
+  call update_pt2_and_variance_weights(pt2, variance, norm, N_states)
 
 end subroutine
 
@@ -159,9 +158,9 @@ subroutine selection_collector(zmq_socket_pull, b, N, pt2, variance, norm)
   integer(ZMQ_PTR), intent(in)   :: zmq_socket_pull
   type(selection_buffer), intent(inout) :: b
   integer, intent(in)            :: N
-  double precision, intent(inout)    :: pt2(N_states)
-  double precision, intent(inout)    :: variance(N_states)
-  double precision, intent(inout)    :: norm(N_states)
+  double precision, intent(out)      :: pt2(N_states)
+  double precision, intent(out)      :: variance(N_states)
+  double precision, intent(out)      :: norm(N_states)
   double precision                   :: pt2_mwen(N_states)
   double precision                   :: variance_mwen(N_states)
   double precision                   :: norm_mwen(N_states)

src/davidson/print_e_components.irp.f

@@ -0,0 +1,54 @@
+subroutine print_energy_components()
+  implicit none
+  BEGIN_DOC
+! Prints the different components of the energy.
+  END_DOC
+  integer, save :: ifirst = 0
+  double precision :: Vee, Ven, Vnn, Vecp, T, f
+  integer  :: i,j,k
+
+  Vnn = nuclear_repulsion
+
+  print *, 'Energy components'
+  print *, '================='
+  print *, ''
+  do k=1,N_states
+
+    Ven  = 0.d0
+    Vecp = 0.d0
+    T    = 0.d0
+
+    do j=1,mo_num
+      do i=1,mo_num
+        f = one_e_dm_mo_alpha(i,j,k) + one_e_dm_mo_beta(i,j,k)
+        Ven  = Ven  + f * mo_integrals_n_e(i,j)
+        Vecp = Vecp + f * mo_pseudo_integrals(i,j)
+        T    = T    + f * mo_kinetic_integrals(i,j)
+      enddo
+    enddo
+    Vee = psi_energy(k) - Ven - Vecp - T
+    
+    if (ifirst == 0) then
+      ifirst = 1
+      print *, 'Vnn  : Nucleus-Nucleus   potential energy'
+      print *, 'Ven  : Electron-Nucleus  potential energy'
+      print *, 'Vee  : Electron-Electron potential energy'
+      print *, 'Vecp : Potential energy of the pseudo-potentials'
+      print *, 'T    : Electronic kinetic energy'
+      print *, ''
+    endif
+
+    print *, 'State ', k
+    print *, '---------'
+    print *, ''
+    print *, 'Vnn  = ', Vnn
+    print *, 'Ven  = ', Ven
+    print *, 'Vee  = ', Vee
+    print *, 'Vecp = ', Vecp
+    print *, 'T    = ', T
+    print *, ''
+  enddo
+
+  print *, ''
+
+end

src/determinants/EZFIO.cfg

@@ -28,12 +28,18 @@ doc: Force the wave function to be an eigenfunction of |S^2|
 interface: ezfio,provider,ocaml
 default: True
 
-[used_weight]
+[weight_one_e_dm]
 type: integer
 doc: Weight used in the calculation of the one-electron density matrix. 0: 1./(c_0^2), 1: 1/N_states, 2: input state-average weight, 3: 1/(Norm_L3(Psi))
 interface: ezfio,provider,ocaml
 default: 1
 
+[weight_selection]
+type: integer
+doc: Weight used in the selection. 0: input state-average weight, 1: 1./(c_0^2), 2: rPT2 matching, 3: variance matching, 4: variance and rPT2 matching, 5: variance minimization and matching                                  
+interface: ezfio,provider,ocaml
+default: 2
+
 [threshold_generators]
 type: Threshold
 doc: Thresholds on generators (fraction of the square of the norm) 
@@ -89,6 +95,11 @@ doc: Weight of the states in state-average calculations.
 interface: ezfio
 size: (determinants.n_states)
 
+[selection_factor]
+type: double precision
+doc: f such that the number of determinants to add is f * N_det * sqrt(N_states)
+interface: ezfio,provider,ocaml
+default: 1.
 
 [thresh_sym]
 type: Threshold

src/determinants/density_matrix.irp.f

@@ -305,9 +305,9 @@ BEGIN_PROVIDER [ double precision, state_average_weight, (N_states) ]
    logical                        :: exists
 
    state_average_weight(:) = 1.d0
-   if (used_weight == 0) then
+   if (weight_one_e_dm == 0) then
      state_average_weight(:) = c0_weight(:)
-   else if (used_weight == 1) then
+   else if (weight_one_e_dm == 1) then
      state_average_weight(:) = 1./N_states
    else
      call ezfio_has_determinants_state_average_weight(exists)

src/determinants/utils.irp.f

@@ -43,4 +43,3 @@ BEGIN_PROVIDER [ double precision, S2_matrix_all_dets,(N_det,N_det) ]
  !$OMP END PARALLEL DO
 END_PROVIDER
 
-

src/fci/pt2.irp.f

@@ -46,7 +46,7 @@ subroutine run
   call ZMQ_pt2(psi_energy_with_nucl_rep,pt2,relative_error,error, variance, &
      norm,0) ! Stochastic PT2
   do k=1,N_states
-    rpt2(:) = pt2(:)/(1.d0 + norm(k))
+    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)

src/iterations/print_summary.irp.f

@@ -31,18 +31,19 @@ subroutine print_summary(e_,pt2_,error_,variance_,norm_,n_det_,n_occ_pattern_,n_
 
   write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
   write(*,fmt)
-  write(fmt,*) '(12X,', N_states_p, '(6X,A7,1X,I6,10X))'
+  write(fmt,*) '(13X,', N_states_p, '(6X,A7,1X,I6,10X))'
   write(*,fmt) ('State',k, k=1,N_states_p)
   write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
   write(*,fmt)
-  write(fmt,*) '(A12,', N_states_p, '(1X,F14.8,15X))'
+  write(fmt,*) '(A13,', N_states_p, '(1X,F14.8,15X))'
   write(*,fmt) '# E          ', e_(1:N_states_p)
   if (N_states_p > 1) then
     write(*,fmt) '# Excit. (au)', e_(1:N_states_p)-e_(1)
     write(*,fmt) '# Excit. (eV)', (e_(1:N_states_p)-e_(1))*27.211396641308d0
   endif
   write(fmt,*) '(A13,', 2*N_states_p, '(1X,F14.8))'
-  write(*,fmt) '# PT2'//pt2_string, (pt2_(k), error_(k), k=1,N_states_p)
+  write(*,fmt) '# PT2 '//pt2_string, (pt2_(k), error_(k), k=1,N_states_p)
+  write(*,fmt) '# rPT2'//pt2_string, (pt2_(k)*f(k), error_(k)*f(k), k=1,N_states_p)
   write(*,'(A)') '#'
   write(*,fmt) '# E+PT2      ', (e_(k)+pt2_(k),error_(k), k=1,N_states_p)
   write(*,fmt) '# E+rPT2     ', (e_(k)+pt2_(k)*f(k),error_(k)*f(k), k=1,N_states_p)
@@ -97,5 +98,7 @@ subroutine print_summary(e_,pt2_,error_,variance_,norm_,n_det_,n_occ_pattern_,n_
     enddo
   endif
 
+  call print_energy_components()
+
 end subroutine