Working on mrcc

plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f

@@ -275,6 +275,8 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2, error)
   
   integer, allocatable :: f(:)
   logical, allocatable :: d(:) 
+  logical :: do_exit
+
 
   zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
   allocate(task_id(pt2_n_tasks_max), index(pt2_n_tasks_max), f(N_det_generators))
@@ -296,6 +298,7 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2, error)
   more = 1
   time0 = omp_get_wtime()
 
+  do_exit = .false.
   do while (n <= N_det_generators)
     if(f(pt2_J(n)) == 0) then
       d(pt2_J(n)) = .true.
@@ -328,6 +331,9 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2, error)
           S2(p) += x**2
         end do
         avg = E0 + S(t) / dble(c)
+        if ((avg /= 0.d0) .or. (n == N_det_generators) ) then 
+          do_exit = .true.
+        endif
         pt2(pt2_stoch_istate) = avg
         ! 1/(N-1.5) : see  Brugger, The American Statistician (23) 4 p. 32 (1969)
         if(c > 2) then
@@ -336,7 +342,7 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2, error)
           error(pt2_stoch_istate) = eqt
           if(mod(c,10)==0 .or. n==N_det_generators) then
             print '(G10.3, 2X, F16.10, 2X, G16.3, 2X, F16.4, A20)', c, avg+E, eqt, time-time0, ''
-            if( dabs(error(pt2_stoch_istate) / pt2(pt2_stoch_istate)) < relative_error) then
+            if(do_exit .and. (dabs(error(pt2_stoch_istate)) / (1.d-20 + dabs(pt2(pt2_stoch_istate)) ) <= relative_error)) then 
               if (zmq_abort(zmq_to_qp_run_socket) == -1) then
                 call sleep(10)
                 if (zmq_abort(zmq_to_qp_run_socket) == -1) then

plugins/dress_zmq/alpha_factory.irp.f

@@ -95,6 +95,7 @@ subroutine generate_singles_and_doubles(delta_ij_loc, i_generator, bitmask_index
   allocate (indices(N_det),  &
             exc_degree(max(N_det_alpha_unique,N_det_beta_unique)))
 
+  PROVIDE psi_det_sorted_gen_order
   !PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique
   !PROVIDE psi_bilinear_matrix_rows psi_det_sorted_gen_order psi_bilinear_matrix_order
   !PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns

plugins/dress_zmq/dress_stoch_routines.irp.f

@@ -25,7 +25,7 @@ END_PROVIDER
   enddo
 
   
-  if(N_det_generators < 1024) then
+  if(N_det_generators < 128) then
     pt2_minDetInFirstTeeth = 1
     pt2_N_teeth = 1
   else
@@ -485,13 +485,13 @@ subroutine dress_collector(zmq_socket_pull, E, relative_error, delta, delta_s2,
   double precision, external :: omp_get_wtime
   integer, allocatable :: dot_f(:)
   integer, external :: zmq_delete_tasks, dress_find_sample
-  logical :: found
+  logical :: do_exit
   integer :: worker_id
   worker_id=1
   zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
 
   
-  found = .false.
+  do_exit = .false.
   delta = 0d0
   delta_s2 = 0d0
   allocate(task_id(pt2_n_tasks_max))
@@ -513,7 +513,7 @@ subroutine dress_collector(zmq_socket_pull, E, relative_error, delta, delta_s2,
   time0 = -1d0 ! omp_get_wtime()
   more = 1
 
-  do while (.not. found)
+  do 
     if(dot_f(m) == 0) then
       E0 = 0
       do i=dress_dot_n_0(m),1,-1
@@ -532,6 +532,9 @@ subroutine dress_collector(zmq_socket_pull, E, relative_error, delta, delta_s2,
       end do
       t = dress_dot_t(m)
       avg = E0 + S(t) / dble(c)
+      if ((avg /= 0.d0) .or. (m == dress_N_cp) ) then
+        do_exit = .true.
+      endif
       if (c > 2) then
         eqt = dabs((S2(t) / c) - (S(t)/c)**2)
         error = sqrt(eqt / (dble(c)-1.5d0))
@@ -543,12 +546,11 @@ subroutine dress_collector(zmq_socket_pull, E, relative_error, delta, delta_s2,
         error =1.d0
       endif
       m += 1
-! /!\ avg == 0.d0 needs to be handled here
+      if(do_exit .and. (dabs(error) / (1.d-20 + dabs(avg) ) <= relative_error)) then 
-      if(dabs(error / avg) <= relative_error) then
         integer, external :: zmq_put_dvector
         integer, external :: zmq_put_int
         i= zmq_put_int(zmq_to_qp_run_socket, worker_id, 'ending', (m-1))
-        found = .true.
+        exit
       end if
     else
       do

plugins/mrcc/mrcc.irp.f

@@ -10,7 +10,80 @@ program shifted_bk
   PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns
   PROVIDE psi_bilinear_matrix_transp_order
 
+  read_wf = .True.
+  SOFT_TOUCH read_wf
+  call set_generators_bitmasks_as_holes_and_particles
+  if (.True.) then
+    integer :: i,j
+    do j=1,N_states
+      do i=1,N_det
+        psi_coef(i,j) = CI_eigenvectors(i,j)
+      enddo
+    enddo
+    TOUCH psi_coef
+  endif
 
   call dress_zmq()
+  if (.true.) then
+    call ezfio_set_mrcc_energy(ci_energy_dressed(1))
+  endif
+  if (do_pt2) then
+    call run_pt2(N_states, ci_energy_dressed)
+  endif
+end
+
+
+subroutine run_pt2(N_st,energy) 
+  implicit none 
+  integer :: i,j,k 
+  integer, intent(in)          :: N_st 
+  double precision, intent(in) :: energy(N_st) 
+  double precision :: pt2(N_st)
+  double precision :: norm_pert(N_st),H_pert_diag(N_st)
+  
+  pt2 = 0d0
+  !if(lambda_mrcc_pt2(0)  == 0) return
+  
+  print*,'Last iteration only to compute the PT2' 
+  
+  N_det_generators = N_det_cas
+  N_det_selectors = N_det_non_ref
+
+  do i=1,N_det_generators
+    do k=1,N_int
+      psi_det_generators(k,1,i) = psi_ref(k,1,i)
+      psi_det_generators(k,2,i) = psi_ref(k,2,i)
+    enddo
+    do k=1,N_st
+      psi_coef_generators(i,k) = psi_ref_coef(i,k)
+    enddo
+  enddo
+  do i=1,N_det
+    do k=1,N_int
+      psi_selectors(k,1,i) = psi_det_sorted(k,1,i)
+      psi_selectors(k,2,i) = psi_det_sorted(k,2,i)
+    enddo
+    do k=1,N_st
+      psi_selectors_coef(i,k) = psi_coef_sorted(i,k)
+    enddo
+  enddo
+
+  SOFT_TOUCH N_det_selectors psi_selectors_coef psi_selectors N_det_generators psi_det_generators psi_coef_generators ci_eigenvectors_dressed ci_eigenvectors_s2_dressed ci_electronic_energy_dressed
+  SOFT_TOUCH psi_ref_coef_diagonalized psi_ref_energy_diagonalized
+
+  call H_apply_mrcepa_PT2(pt2, norm_pert, H_pert_diag,  N_st)
+
+!  call ezfio_set_full_ci_energy_pt2(energy+pt2)
+
+  print *,  'Final step'
+  print *,  'N_det    = ', N_det
+  print *,  'N_states = ', N_states
+  print *,  'PT2      = ', pt2
+  print *,  'E        = ', energy
+  print *,  'E+PT2    = ', energy+pt2
+  print *,  '-----'
+
+  call ezfio_set_mrcc_energy_pt2(energy(1)+pt2(1))
+
 end 
 

plugins/mrcc/mrcc_routines.irp.f

@@ -8,12 +8,12 @@ subroutine generator_start(i_gen, iproc, interesting)
   logical, external :: deteq
   PROVIDE dij
   interesting = .true.
-  do i=1,N_det_ref
+!  do i=1,N_det_ref
-    if(deteq(psi_det_generators(1,1,i_gen), psi_ref(1,1,i), N_int)) then
+!    if(deteq(psi_det_generators(1,1,i_gen), psi_ref(1,1,i), N_int)) then
-      interesting = .false.
+!      interesting = .false.
-      exit
+!      exit
-    end if
+!    end if
-  end do
+!  end do
 end subroutine
 
  BEGIN_PROVIDER [ double precision, hij_cache_, (N_det,Nproc) ]
@@ -65,6 +65,11 @@ subroutine dress_with_alpha_buffer(Nstates, Ndet,Nint,delta_ij_loc, i_gen, minil
   double precision               :: phase, phase2
   integer                        :: degree, exc(0:2,2,2)
   integer(8), save :: diamond = 0
+
+!  call dress_with_alpha_buffer_old(Nstates, Ndet,Nint,delta_ij_loc, i_gen, minilist, det_minilist, n_minilist, alpha, iproc)
+!  return
+
+
   if(n_minilist == 1) return
   !check if not linked to reference
   do i=1,n_minilist
@@ -200,7 +205,7 @@ subroutine dress_with_alpha_buffer_neu(Nstates,Ndet,Nint,delta_ij_loc, i_gen, mi
       end do
       !diamond found
       diamond += 1
-      if(mod(diamond,10000) == 1) print *, "diam", diamond
+!      if(mod(diamond,10000) == 1) print *, "diam", diamond
 
       call get_excitation(psi_ref(1,1,i_I),det_minilist(1,1,j),exc,degree,phase,Nint)
       call get_excitation(alpha,det_minilist(1,1,i),exc,degree,phase2,Nint)
@@ -361,13 +366,11 @@ subroutine dress_with_alpha_buffer_old(Nstates,Ndet,Nint,delta_ij_loc, i_gen, mi
       
            
       if (perturbative_triples.and. (degree2 == 1) ) then
-          if(sum(popcnt(tmp_det(:,1))) /= elec_alpha_num) stop "STOP 1"
+          if (ok) then
-  if(sum(popcnt(tmp_det(:,2))) /= elec_beta_num) stop "STOP 2"
+            call i_h_j(psi_ref(1,1,i_I),tmp_det,Nint,hka)
-  if(sum(popcnt(tmp_det(:,1))) /= elec_alpha_num) stop "STOP 3"
+          else
-  if(sum(popcnt(tmp_det(:,2))) /= elec_beta_num) stop "STOP 4"
+            hka = 0.d0
-  
+          endif
-
-          call i_h_j(psi_ref(1,1,i_I),tmp_det,Nint,hka)
           hka = hij_cache_(k_sd,iproc) - hka
           if (dabs(hka) > 1.d-12) then
             call get_delta_e_dyall_general_mp(psi_ref(1,1,i_I),alpha,Delta_E_inv)

plugins/mrcepa0/dressing.irp.f

@@ -1228,7 +1228,11 @@ subroutine get_cc_coef(tq,c_alpha)
       endif
       
       if (perturbative_triples.and. (degree2 == 1) ) then
-        call i_h_j(psi_ref(1,1,i_I),tmp_det,N_int,hka)
+        if (ok) then
+          call i_h_j(psi_ref(1,1,i_I),tmp_det,N_int,hka)
+        else
+          hka = 0.d0
+        endif
         call i_h_j(tq,psi_non_ref(1,1,k_sd),N_int,hla)
         hka = hla - hka
         if (dabs(hka) > 1.d-12) then

tests/bats/mrcepa0.bats

@@ -5,7 +5,7 @@ source $QP_ROOT/tests/bats/common.bats.sh
 #=== H2O
 @test "MRCC-lambda H2O cc-pVDZ" {
   INPUT=h2o.ezfio
-  EXE=mrcc_zmq
+  EXE=mrcc
   test_exe $EXE || skip
   qp_edit -c $INPUT  
   ezfio set_file $INPUT
@@ -17,7 +17,7 @@ source $QP_ROOT/tests/bats/common.bats.sh
   ezfio set mrcepa0 n_it_max_dressed_ci 3
   cp -r $INPUT TMP ; qp_run $EXE TMP 
   ezfio set_file TMP
-  energy="$(ezfio get mrcepa0 energy_pt2)"
+  energy="$(ezfio get mrcc energy)"
   rm -rf TMP
   eq $energy -76.2382975461183 1.e-4
 }
@@ -25,7 +25,7 @@ source $QP_ROOT/tests/bats/common.bats.sh
 
 @test "MRCC H2O cc-pVDZ" {
   INPUT=h2o.ezfio
-  EXE=mrcc_zmq
+  EXE=mrcc
   test_exe $EXE || skip
   qp_edit -c $INPUT  
   ezfio set_file $INPUT
@@ -37,14 +37,14 @@ source $QP_ROOT/tests/bats/common.bats.sh
   ezfio set mrcepa0 n_it_max_dressed_ci 3
   cp -r $INPUT TMP ; qp_run $EXE TMP 
   ezfio set_file TMP
-  energy="$(ezfio get mrcepa0 energy_pt2)"
+  energy="$(ezfio get mrcc energy)"
   rm -rf TMP
   eq $energy -76.2382468380776 1.e-4
 }
 
 #@test "MRCC-stoch H2O cc-pVDZ" {
 #  INPUT=h2o.ezfio
-#  EXE=mrcc_zmq
+#  EXE=mrcc
 #  test_exe $EXE || skip
 #  qp_edit -c $INPUT  
 #  ezfio set_file $INPUT
@@ -60,41 +60,3 @@ source $QP_ROOT/tests/bats/common.bats.sh
 #  eq $energy -76.2379517543157 1.e-4
 #}
 
-@test "MRSC2 H2O cc-pVDZ" {
-  INPUT=h2o.ezfio
-  EXE=mrsc2
-  test_exe $EXE || skip
-  qp_edit -c $INPUT  
-  ezfio set_file $INPUT
-  ezfio set determinants threshold_generators 1.
-  ezfio set determinants threshold_selectors  1.
-  ezfio set determinants read_wf True
-  ezfio set mrcepa0 lambda_type 1
-  ezfio set mrcepa0 n_it_max_dressed_ci 3
-  ezfio set mrcepa0 perturbative_triples 0
-  cp -r $INPUT TMP ; qp_run $EXE TMP 
-  ezfio set_file TMP
-  energy="$(ezfio get mrcepa0 energy_pt2)"
-  rm -rf TMP
-  eq $energy -76.2358860928235 3.e-4
-}
-
-@test "MRCEPA0 H2O cc-pVDZ" {
-  INPUT=h2o.ezfio
-  EXE=mrcepa0
-  test_exe $EXE || skip
-  qp_edit -c $INPUT  
-  ezfio set_file $INPUT
-  ezfio set determinants threshold_generators 1.
-  ezfio set determinants threshold_selectors  1.
-  ezfio set determinants read_wf True
-  ezfio set mrcepa0 perturbative_triples 0
-  ezfio set mrcepa0 lambda_type 1
-  ezfio set mrcepa0 n_it_max_dressed_ci 3
-  cp -r $INPUT TMP ; qp_run $EXE TMP
-  ezfio set_file TMP
-  energy="$(ezfio get mrcepa0 energy_pt2)"
-  rm -rf TMP
-  eq $energy -76.2412031502384 2.e-4
-}
-