OpenMP in save_MRCC_wf

plugins/CAS_SD_ZMQ/cassd_zmq.irp.f

@@ -5,16 +5,28 @@ program cassd_zmq
   integer                        :: degree
   integer                        :: n_det_before, to_select
   double precision               :: threshold_davidson_in
+  double precision               :: error(N_states)
   
   allocate (pt2(N_states))
 
   double precision               :: hf_energy_ref
   logical                        :: has
+  integer                        :: N_states_p
+  character*(512)                :: fmt
+  character*(8) :: pt2_string
+
   pt2 = -huge(1.d0)
+  error = 0.d0
   threshold_davidson_in = threshold_davidson
   threshold_davidson = threshold_davidson_in * 100.d0
   SOFT_TOUCH threshold_davidson
 
+  if (do_pt2) then
+    pt2_string = '        '
+  else
+    pt2_string = '(approx)'
+  endif
+
   call diagonalize_CI
   call save_wavefunction
   
@@ -45,7 +57,6 @@ program cassd_zmq
   double precision               :: E_CI_before(N_states)
   
   
-  print*,'Beginning the selection ...'
   if (.True.) then ! Avoid pre-calculation of CI_energy
     E_CI_before(1:N_states) = CI_energy(1:N_states)
   endif
@@ -60,6 +71,8 @@ program cassd_zmq
           (maxval(abs(pt2(1:N_states))) > pt2_max) .and.               &
           (correlation_energy_ratio <= correlation_energy_ratio_max)    &
           )
+       write(*,'(A)')  '--------------------------------------------------------------------------------'
+
 
       correlation_energy_ratio = (CI_energy(1) - hf_energy_ref)  /     &
                       (E_CI_before(1) + pt2(1) - hf_energy_ref)
@@ -98,6 +111,67 @@ program cassd_zmq
       to_select = min(to_select, N_det_max-n_det_before)
       call ZMQ_selection(to_select, pt2)
       
+      N_states_p = min(N_det,N_states)
+
+      print *, ''
+      print '(A,I12)',  'Summary at N_det = ', N_det
+      print '(A)',      '-----------------------------------'
+      print *, ''
+      call write_double(6,correlation_energy_ratio, 'Correlation ratio')
+      print *, ''
+
+      write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
+      write(*,fmt)
+      write(fmt,*) '(12X,', 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) '# E          ', E_CI_before(1:N_states_p)
+      if (N_states_p > 1) then
+        write(*,fmt) '# Excit. (au)', E_CI_before(1:N_states_p)-E_CI_before(1)
+        write(*,fmt) '# Excit. (eV)', (E_CI_before(1:N_states_p)-E_CI_before(1))*27.211396641308d0
+      endif
+      write(fmt,*) '(A12,', 2*N_states_p, '(1X,F14.8))'
+      write(*,fmt) '# PT2'//pt2_string, (pt2(k), error(k), k=1,N_states_p)
+      write(*,'(A)') '#'
+      write(*,fmt) '# E+PT2      ', (E_CI_before(k)+pt2(k),error(k), k=1,N_states_p)
+      if (N_states_p > 1) then
+        write(*,fmt) '# Excit. (au)', ( (E_CI_before(k)+pt2(k)-E_CI_before(1)-pt2(1)), &
+          dsqrt(error(k)*error(k)+error(1)*error(1)), k=1,N_states_p)
+        write(*,fmt) '# Excit. (eV)', ( (E_CI_before(k)+pt2(k)-E_CI_before(1)-pt2(1))*27.211396641308d0, &
+          dsqrt(error(k)*error(k)+error(1)*error(1))*27.211396641308d0, k=1,N_states_p)
+      endif
+      write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
+      write(*,fmt)
+      print *,  ''
+
+      print *,  'N_det             = ', N_det
+      print *,  'N_states          = ', N_states
+      print*,   'correlation_ratio = ', correlation_energy_ratio
+
+      do k=1, N_states_p
+        print*,'State ',k
+        print *,  'PT2             = ', pt2(k)
+        print *,  'E               = ', E_CI_before(k)
+        print *,  'E+PT2'//pt2_string//'   = ', E_CI_before(k)+pt2(k), ' +/- ', error(k)
+      enddo
+
+      print *,  '-----'
+      if(N_states.gt.1)then
+        print *, 'Variational Energy difference (au | eV)'
+        do i=2, N_states_p
+          print*,'Delta E = ', (E_CI_before(i) - E_CI_before(1)), &
+            (E_CI_before(i) - E_CI_before(1)) * 27.211396641308d0
+        enddo
+        print *,  '-----'
+        print*, 'Variational + perturbative Energy difference (au | eV)'
+        do i=2, N_states_p
+          print*,'Delta E = ', (E_CI_before(i)+ pt2(i) - (E_CI_before(1) + pt2(1))), &
+            (E_CI_before(i)+ pt2(i) - (E_CI_before(1) + pt2(1))) * 27.211396641308d0
+        enddo
+      endif
+
       PROVIDE  psi_coef
       PROVIDE  psi_det
       PROVIDE  psi_det_sorted

plugins/Full_CI_ZMQ/fci_zmq.irp.f

@@ -22,8 +22,8 @@ program fci_zmq
   threshold_davidson = threshold_davidson_in * 100.d0
   SOFT_TOUCH threshold_davidson
 
-  call diagonalize_CI
+!  call diagonalize_CI
-  call save_wavefunction
+!  call save_wavefunction
   
   call ezfio_has_hartree_fock_energy(has)
   if (has) then

plugins/Generators_CAS/generators.irp.f

@@ -35,10 +35,10 @@ END_PROVIDER
   integer, external :: number_of_holes,number_of_particles
   integer, allocatable :: nongen(:)
   integer :: inongen
-  inongen = 0
 
   allocate(nongen(N_det))
 
+  inongen = 0
   m=0
   do i=1,N_det
     good = ( number_of_holes(psi_det_sorted(1,1,i)) ==0).and.(number_of_particles(psi_det_sorted(1,1,i))==0 )

plugins/dress_zmq/alpha_factory.irp.f

@@ -66,7 +66,7 @@ subroutine generate_singles_and_doubles(delta_ij_loc, i_generator, bitmask_index
   integer(bit_kind), allocatable:: preinteresting_det(:,:,:)
   integer ,allocatable :: abuf(:), labuf(:)
   
-  allocate(abuf(N_det*6), labuf(N_det))
+  allocate(abuf(0:N_det*6), labuf(0:N_det))
   allocate(preinteresting_det(N_int,2,N_det))
   
   PROVIDE fragment_count
@@ -387,7 +387,7 @@ subroutine alpha_callback_mask(delta_ij_loc, i_gen, sp, mask, bannedOrb, banned,
   integer(bit_kind), allocatable :: det_minilist(:,:,:)
 
 
-  allocate(abuf(siz), labuf(N_det), putten(N_det), det_minilist(N_int, 2, N_det))
+  allocate(abuf(0:siz), labuf(0:N_det), putten(N_det), det_minilist(N_int, 2, N_det))
   
   do i=1,siz
     abuf(i) = psi_from_sorted_gen(rabuf(i))
@@ -638,7 +638,7 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, indexes, ab
   integer(bit_kind),intent(in)   :: mask(N_int, 2), det(N_int, 2, N_sel)
   logical, intent(inout)         :: bannedOrb(mo_tot_num, 2), banned(mo_tot_num, mo_tot_num, 2)
   integer, intent(inout)         :: indexes(0:mo_tot_num, 0:mo_tot_num)
-  integer, intent(inout) :: abuf(*)
+  integer, intent(inout) :: abuf(0:*)
   integer                        :: i, ii, j, k, l, h(0:2,2), p(0:4,2), nt, s
   integer(bit_kind)              :: perMask(N_int, 2), mobMask(N_int, 2), negMask(N_int, 2)
   integer :: phasemask(2,N_int*bit_kind_size)
@@ -704,7 +704,7 @@ subroutine get_d2(i_gen, gen, banned, bannedOrb, indexes, abuf, mask, h, p, sp)
   implicit none
 
   integer(bit_kind), intent(in) :: mask(N_int, 2), gen(N_int, 2)
-  integer, intent(inout) :: abuf(*)
+  integer, intent(inout) :: abuf(0:*)
   integer, intent(in) :: i_gen
   logical, intent(in) :: bannedOrb(mo_tot_num, 2), banned(mo_tot_num, mo_tot_num,2)
   integer, intent(inout) :: indexes(0:mo_tot_num, 0:mo_tot_num)
@@ -832,7 +832,7 @@ subroutine get_d1(i_gen, gen, banned, bannedOrb, indexes, abuf, mask, h, p, sp)
   implicit none
 
   integer(bit_kind), intent(in)  :: mask(N_int, 2), gen(N_int, 2)
-  integer, intent(inout)         :: abuf(*)
+  integer, intent(inout)         :: abuf(0:*)
   integer,intent(in)             :: i_gen
   logical, intent(in)            :: bannedOrb(mo_tot_num, 2), banned(mo_tot_num, mo_tot_num,2)
   integer(bit_kind)              :: det(N_int, 2)

plugins/mrcepa0/dressing.irp.f

@@ -236,6 +236,7 @@ subroutine mrcc_part_dress(delta_ij_, delta_ij_s2_, i_generator,n_selected,det_b
   double precision, intent(inout) :: contrib(N_states)
   double precision :: sdress, hdress
 
+  PROVIDE n_act_orb elec_num
 
   if (perturbative_triples) then
     PROVIDE one_anhil fock_virt_total fock_core_inactive_total one_creat
@@ -1025,7 +1026,7 @@ subroutine filter_tq(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_m
     if (good) then
       if (.not. is_in_wavefunction(det_buffer(1,1,i),Nint)) then
         N_tq += 1
-        do k=1,N_int
+        do k=1,Nint
           tq(k,1,N_tq) = det_buffer(k,1,i)
           tq(k,2,N_tq) = det_buffer(k,2,i)
         enddo
@@ -1131,6 +1132,7 @@ subroutine get_cc_coef(tq,c_alpha)
   integer                        :: i_state, k_sd, l_sd, i_I
   logical                        :: ok
   
+  PROVIDE n_act_orb elec_num
   if (perturbative_triples) then
     PROVIDE one_anhil fock_virt_total fock_core_inactive_total one_creat
   endif

plugins/mrcepa0/save_mrcc_wavefunction.irp.f

@@ -3,6 +3,7 @@ program save_mrcc_wf
   
   threshold_generators = 1.d0
   threshold_selectors = 1.d0
+  PROVIDE N_int psi_det
   TOUCH threshold_generators threshold_selectors
 
   mrmode=5
@@ -21,12 +22,19 @@ subroutine run1
   double precision :: c_alpha(N_states)
   call set_generators_bitmasks_as_holes_and_particles
 
+  call get_cc_coef(psi_det(1,1,1), c_alpha)
+  !$OMP PARALLEL DO DEFAULT(SHARED) &
+  !$OMP PRIVATE(k,c_alpha) SCHEDULE(static,64)
   do k=1,N_det
     if (maxval(abs(psi_coef(k,1:N_states))) == 0.d0) then
+      if (iand(k,1023) == 0) then
+        print *,  k, '/', N_det
+      endif
       call get_cc_coef(psi_det(1,1,k), c_alpha)
       psi_coef(k,1:N_states) = c_alpha(1:N_states)
     endif
   enddo
+  !$OMP END PARALLEL DO
   SOFT_TOUCH psi_coef
 end
 

src/Determinants/slater_rules.irp.f

@@ -568,7 +568,6 @@ subroutine i_H_j_s2(key_i,key_j,Nint,hij,s2)
       call get_mono_excitation_from_fock(key_i,key_j,p,m,spin,phase,hij)
       
     case (0)
-      print *,irp_here,": ZERO"
       double precision, external :: diag_S_mat_elem
       s2 = diag_S_mat_elem(key_i,Nint)
       hij = diag_H_mat_elem(key_i,Nint)