Merge branch 'master' of github.com:scemama/quantum_package

ocaml/Address.ml

@@ -42,7 +42,7 @@ end = struct
     assert (String.is_prefix ~prefix:"inproc://" x);
     x
   let create name =
-    Printf.sprintf "ipc://%s" name
+    Printf.sprintf "inproc://%s" name
   let to_string x = x
 end
 

ocaml/TaskServer.ml

@@ -678,9 +678,9 @@ let run ~port =
 
               (** Debug input *)
               Printf.sprintf "q:%d  r:%d  n:%d  : %s\n%!"
-              (Queuing_system.number_of_queued  program_state.queue)
+              (Queuing_system.number_of_queued program_state.queue)
               (Queuing_system.number_of_running program_state.queue)
-              (Queuing_system.number_of_tasks   program_state.queue)
+              (Queuing_system.number_of_tasks program_state.queue)
               (Message.to_string message)
               |> debug;
 

plugins/Full_CI_ZMQ/fci_zmq.irp.f

@@ -134,7 +134,7 @@ subroutine ZMQ_selection(N_in, pt2)
 
   step = int(5000000.d0 / dble(N_int * N_states * elec_num * elec_num * mo_tot_num * mo_tot_num ))
   step = max(1,step)
-  do i= 1,N_det_generators, step
+  do i= N_det_generators, 1, -step
     i_generator_start = max(i-step+1,1)
     i_generator_max = i
     write(task,*) i_generator_start, i_generator_max, 1, N

plugins/MRCC_Utils/mrcc_dress.irp.f

@@ -271,7 +271,7 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge
           !delta_ii_(i_state,i_I)  = 0.d0
           do l_sd=1,idx_alpha(0)
             k_sd = idx_alpha(l_sd)
-            delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + dIa_hla(i_state,k_sd)
+            delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + 0.5d0 * dIa_hla(i_state,k_sd)
           enddo
         endif
       enddo

plugins/MRCC_Utils/mrcc_utils.irp.f

@@ -685,7 +685,7 @@ END_PROVIDER
         
 
   
-  do s = 1, N_states
+  do s=1, N_states
     
     A_val = 0d0
     A_ind = 0
@@ -698,61 +698,61 @@ END_PROVIDER
     
     !$OMP PARALLEL default(none) shared(psi_non_ref, hh_exists, pp_exists, N_int, A_val, A_ind)&
         !$OMP shared(s, hh_shortcut, psi_ref_coef, N_det_non_ref, psi_non_ref_sorted, psi_non_ref_sorted_idx, psi_ref, N_det_ref)&
-        !$OMP shared(active, active_hh_idx, active_pp_idx, nactive)&
+        !$OMP shared(active, active_hh_idx, active_pp_idx, nactive)  &
         !$OMP private(lref, pp, II, ok, myMask, myDet, ind, phase, wk, ppp, hh)
-     allocate(lref(N_det_non_ref))
+    allocate(lref(N_det_non_ref))
-    !$OMP DO schedule(static,10) 
+    !$OMP DO schedule(static,10)
-      do ppp=1,nactive
+    do ppp=1,nactive
-        pp = active_pp_idx(ppp)
+      pp = active_pp_idx(ppp)
-        hh = active_hh_idx(ppp)
+      hh = active_hh_idx(ppp)
-        lref = 0
+      lref = 0
-        do II = 1, N_det_ref
+      do II = 1, N_det_ref
-          call apply_hole_local(psi_ref(1,1,II), hh_exists(1, hh), myMask, ok, N_int)
+        call apply_hole_local(psi_ref(1,1,II), hh_exists(1, hh), myMask, ok, N_int)
-          if(.not. ok) cycle
+        if(.not. ok) cycle
-          call apply_particle_local(myMask, pp_exists(1, pp), myDet, ok, N_int)
+        call apply_particle_local(myMask, pp_exists(1, pp), myDet, ok, N_int)
-          if(.not. ok) cycle
+        if(.not. ok) cycle
-          ind = searchDet(psi_non_ref_sorted(1,1,1), myDet(1,1), N_det_non_ref, N_int)
+        ind = searchDet(psi_non_ref_sorted(1,1,1), myDet(1,1), N_det_non_ref, N_int)
-          if(ind /= -1) then
+        if(ind /= -1) then
-            call get_phase(myDet(1,1), psi_ref(1,1,II), phase, N_int)
+          call get_phase(myDet(1,1), psi_ref(1,1,II), phase, N_int)
-            if (phase  > 0.d0) then
+          if (phase  > 0.d0) then
-              lref(psi_non_ref_sorted_idx(ind)) =  II
+            lref(psi_non_ref_sorted_idx(ind)) =  II
-            else
+          else
-              lref(psi_non_ref_sorted_idx(ind)) = -II
+            lref(psi_non_ref_sorted_idx(ind)) = -II
-            endif
+          endif
-          end if
+        end if
-        end do
-        wk = 0
-        do i=1, N_det_non_ref
-          if(lref(i) > 0) then
-            wk += 1
-            A_val(wk, ppp) = psi_ref_coef(lref(i), s)
-            A_ind(wk, ppp) = i
-          else if(lref(i) < 0) then
-            wk += 1
-            A_val(wk, ppp) = -psi_ref_coef(-lref(i), s)
-            A_ind(wk, ppp) = i
-          end if
-        end do
-        A_ind(0,ppp) = wk
       end do
+      wk = 0
+      do i=1, N_det_non_ref
+        if(lref(i) > 0) then
+          wk += 1
+          A_val(wk, ppp) = psi_ref_coef(lref(i), s)
+          A_ind(wk, ppp) = i
+        else if(lref(i) < 0) then
+          wk += 1
+          A_val(wk, ppp) = -psi_ref_coef(-lref(i), s)
+          A_ind(wk, ppp) = i
+        end if
+      end do
+      A_ind(0,ppp) = wk
+    end do
     !$OMP END DO
     deallocate(lref)
-    !$OMP END PARALLEL   
+    !$OMP END PARALLEL
-
+    
-
+    
     print *,  'Done building A_val, A_ind'
     
     AtA_size = 0
     col_shortcut = 0
     N_col = 0
-    integer :: a_coll, at_roww
+    integer                        :: a_coll, at_roww
     
     
     !$OMP PARALLEL default(none) shared(k, psi_non_ref_coef, A_ind, A_val, x, N_det_ref, nex, N_det_non_ref)&
         !$OMP private(at_row, a_col, t, i, j, r1, r2, wk, A_ind_mwen, A_val_mwen, a_coll, at_roww)&
         !$OMP shared(col_shortcut, N_col, AtB, AtA_size, AtA_val, AtA_ind, s, nactive, active_pp_idx)
     allocate(A_val_mwen(nex), A_ind_mwen(nex))
-     
+    
     !$OMP DO schedule(dynamic, 100)
     do at_roww = 1, nactive ! nex
       at_row = active_pp_idx(at_roww)
@@ -762,8 +762,8 @@ END_PROVIDER
         j = active_pp_idx(i)
         AtB(at_row) = AtB(at_row) + psi_non_ref_coef(A_ind(i, at_roww), s) * A_val(i, at_roww)
       end do
-
+      
-      do a_coll = 1, nactive 
+      do a_coll = 1, nactive
         a_col = active_pp_idx(a_coll)
         t = 0d0
         r1 = 1
@@ -795,12 +795,12 @@ END_PROVIDER
         col_shortcut(at_roww) = AtA_size+1
         N_col(at_roww) = wk
         if (AtA_size+wk > size(AtA_ind,1)) then
-           print *, AtA_size+wk , size(AtA_ind,1)
+          print *, AtA_size+wk , size(AtA_ind,1)
-           stop 'too small'
+          stop 'too small'
         endif
         do i=1,wk
-         AtA_ind(AtA_size+i) = A_ind_mwen(i)
+          AtA_ind(AtA_size+i) = A_ind_mwen(i)
-         AtA_val(AtA_size+i) = A_val_mwen(i)
+          AtA_val(AtA_size+i) = A_val_mwen(i)
         enddo
         AtA_size += wk
         !$OMP END CRITICAL
@@ -822,41 +822,41 @@ END_PROVIDER
     rho_mrcc_init = 0d0
     
     allocate(lref(N_det_ref))
-    !$OMP PARALLEL DO default(shared) schedule(static, 1) &
+    !$OMP PARALLEL DO default(shared) schedule(static, 1)            &
         !$OMP private(lref, hh, pp, II, myMask, myDet, ok, ind, phase)
     do hh = 1, hh_shortcut(0)
-     do pp = hh_shortcut(hh), hh_shortcut(hh+1)-1
+      do pp = hh_shortcut(hh), hh_shortcut(hh+1)-1
-      if(active(pp)) cycle
+        if(active(pp)) cycle
-      lref = 0
+        lref = 0
-      do II=1,N_det_ref
+        do II=1,N_det_ref
-        call apply_hole_local(psi_ref(1,1,II), hh_exists(1, hh), myMask, ok, N_int)
+          call apply_hole_local(psi_ref(1,1,II), hh_exists(1, hh), myMask, ok, N_int)
-        if(.not. ok) cycle
+          if(.not. ok) cycle
-        call apply_particle_local(myMask, pp_exists(1, pp), myDet, ok, N_int)
+          call apply_particle_local(myMask, pp_exists(1, pp), myDet, ok, N_int)
-        if(.not. ok) cycle
+          if(.not. ok) cycle
-        ind = searchDet(psi_non_ref_sorted(1,1,1), myDet(1,1), N_det_non_ref, N_int)
+          ind = searchDet(psi_non_ref_sorted(1,1,1), myDet(1,1), N_det_non_ref, N_int)
-        if(ind == -1) cycle
+          if(ind == -1) cycle
-        ind = psi_non_ref_sorted_idx(ind)
+          ind = psi_non_ref_sorted_idx(ind)
-        call get_phase(myDet(1,1), psi_ref(1,1,II), phase, N_int)
+          call get_phase(myDet(1,1), psi_ref(1,1,II), phase, N_int)
-        X(pp) += psi_ref_coef(II,s)**2
+          X(pp) += psi_ref_coef(II,s)**2
-        AtB(pp) += psi_non_ref_coef(ind, s) * psi_ref_coef(II, s) * phase
+          AtB(pp) += psi_non_ref_coef(ind, s) * psi_ref_coef(II, s) * phase
-        lref(II) = ind
+          lref(II) = ind
-        if(phase < 0d0) lref(II) = -ind
+          if(phase < 0d0) lref(II) = -ind
+        end do
+        X(pp) =  AtB(pp) / X(pp)
+        do II=1,N_det_ref
+          if(lref(II) > 0) then
+            rho_mrcc_init(lref(II),s) = psi_ref_coef(II,s) * X(pp)
+          else if(lref(II) < 0) then
+            rho_mrcc_init(-lref(II),s) = -psi_ref_coef(II,s) * X(pp)
+          end if
+        end do
       end do
-      X(pp) =  AtB(pp) / X(pp)
-      do II=1,N_det_ref
-        if(lref(II) > 0) then
-          rho_mrcc_init(lref(II),s) = psi_ref_coef(II,s) * X(pp)
-        else if(lref(II) < 0) then
-          rho_mrcc_init(-lref(II),s) = -psi_ref_coef(II,s) * X(pp)
-        end if
-      end do
-     end do
     end do
     !$OMP END PARALLEL DO
     
     x_new = x
-     
+    
-    double precision :: factor, resold
+    double precision               :: factor, resold
     factor = 1.d0
     resold = huge(1.d0)
     do k=0,100000
@@ -882,10 +882,10 @@ END_PROVIDER
       !$OMP END PARALLEL
       
       res = 0.d0
-     
+      
       
       if (res < resold) then
-         do a_coll=1,nactive ! nex
+        do a_coll=1,nactive ! nex
           a_col = active_pp_idx(a_coll)
           do j=1,N_det_non_ref
             i = A_ind(j,a_coll)
@@ -894,39 +894,151 @@ END_PROVIDER
           enddo
           res = res + (X_new(a_col) - X(a_col))*(X_new(a_col) - X(a_col))
           X(a_col) = X_new(a_col)
-         end do
+        end do
-         factor = 1.d0
+        factor = 1.d0
       else
         factor = -factor * 0.5d0
       endif
       resold = res
-
+      
-      if(mod(k, 5) == 0) then
+      if(mod(k, 100) == 0) then
         print *, "res ", k, res
       end if
       
-      if(res < 1d-12) exit
+      if(res < 1d-9) exit
     end do
     
     
     
     norm = 0.d0
-     do i=1,N_det_non_ref
+    do i=1,N_det_non_ref
-       norm = norm + rho_mrcc(i,s)*rho_mrcc(i,s)
+      norm = norm + rho_mrcc(i,s)*rho_mrcc(i,s)
-     enddo
+    enddo
-     ! Norm now contains the norm of A.X
+    ! Norm now contains the norm of A.X
-
+    
-     do i=1,N_det_ref
+    do i=1,N_det_ref
-       norm = norm + psi_ref_coef(i,s)*psi_ref_coef(i,s)
+      norm = norm + psi_ref_coef(i,s)*psi_ref_coef(i,s)
-     enddo
+    enddo
-     ! Norm now contains the norm of Psi + A.X
+    ! Norm now contains the norm of Psi + A.X
-     
+    
-     print *, k, "res : ", res, "norm : ", sqrt(norm)
+    print *, k, "res : ", res, "norm : ", sqrt(norm)
-     
+        
-     !dIj_unique(:size(X), s) = X(:)
+!---------------
+! double precision               :: e_0, overlap
+! double precision, allocatable  :: u_0(:)
+! integer(bit_kind), allocatable :: keys_tmp(:,:,:)
+! allocate (u_0(N_det), keys_tmp(N_int,2,N_det) )
+! k=0
+! overlap = 0.d0
+! do i=1,N_det_ref
+!   k = k+1
+!   u_0(k) = psi_ref_coef(i,1)
+!   keys_tmp(:,:,k) = psi_ref(:,:,i)
+!   overlap += u_0(k)*psi_ref_coef(i,1)
+! enddo
+! norm = 0.d0
+! do i=1,N_det_non_ref
+!   k = k+1
+!   u_0(k) = psi_non_ref_coef(i,1)
+!   keys_tmp(:,:,k) = psi_non_ref(:,:,i)
+!   overlap += u_0(k)*psi_non_ref_coef(i,1)
+! enddo
+! 
+! call u_0_H_u_0(e_0,u_0,N_det,keys_tmp,N_int,1,N_det)
+! print *,  'Energy of |Psi_CASSD> : ', e_0 + nuclear_repulsion, overlap
+!
+! k=0
+! overlap = 0.d0
+! do i=1,N_det_ref
+!   k = k+1
+!   u_0(k) = psi_ref_coef(i,1)
+!   keys_tmp(:,:,k) = psi_ref(:,:,i)
+!   overlap += u_0(k)*psi_ref_coef(i,1)
+! enddo
+! norm = 0.d0
+! do i=1,N_det_non_ref
+!   k = k+1
+!   ! f is such that f.\tilde{c_i} = c_i
+!   f = psi_non_ref_coef(i,1) / rho_mrcc(i,1)
+!   
+!   ! Avoid numerical instabilities
+!   f = min(f,2.d0)
+!   f = max(f,-2.d0)
+!
+!   f = 1.d0
+!
+!   u_0(k) = rho_mrcc(i,1)*f
+!   keys_tmp(:,:,k) = psi_non_ref(:,:,i)
+!   norm += u_0(k)**2
+!   overlap += u_0(k)*psi_non_ref_coef(i,1)
+! enddo
+! 
+! call u_0_H_u_0(e_0,u_0,N_det,keys_tmp,N_int,1,N_det)
+! print *,  'Energy of |(1+T)Psi_0> : ', e_0 + nuclear_repulsion, overlap
+!
+! f = 1.d0/norm
+! norm = 1.d0
+! do i=1,N_det_ref
+!  norm = norm - psi_ref_coef(i,s)*psi_ref_coef(i,s)
+! enddo
+! f = dsqrt(f*norm)
+! overlap = norm
+! do i=1,N_det_non_ref
+!   u_0(k) = rho_mrcc(i,1)*f
+!   overlap += u_0(k)*psi_non_ref_coef(i,1)
+! enddo
+!
+! call u_0_H_u_0(e_0,u_0,N_det,keys_tmp,N_int,1,N_det)
+! print *,  'Energy of |(1+T)Psi_0> (normalized) : ', e_0 + nuclear_repulsion,  overlap
+!
+! k=0
+! overlap = 0.d0
+! do i=1,N_det_ref
+!   k = k+1
+!   u_0(k) = psi_ref_coef(i,1)
+!   keys_tmp(:,:,k) = psi_ref(:,:,i)
+!   overlap += u_0(k)*psi_ref_coef(i,1)
+! enddo
+! norm = 0.d0
+! do i=1,N_det_non_ref
+!   k = k+1
+!   ! f is such that f.\tilde{c_i} = c_i
+!   f = psi_non_ref_coef(i,1) / rho_mrcc(i,1)
+!   
+!   ! Avoid numerical instabilities
+!   f = min(f,2.d0)
+!   f = max(f,-2.d0)
+!
+!   u_0(k) = rho_mrcc(i,1)*f
+!   keys_tmp(:,:,k) = psi_non_ref(:,:,i)
+!   norm += u_0(k)**2
+!   overlap += u_0(k)*psi_non_ref_coef(i,1)
+! enddo
+! 
+! call u_0_H_u_0(e_0,u_0,N_det,keys_tmp,N_int,1,N_det)
+! print *,  'Energy of |(1+T)Psi_0> (mu_i): ', e_0 + nuclear_repulsion, overlap
+!
+! f = 1.d0/norm
+! norm = 1.d0
+! do i=1,N_det_ref
+!  norm = norm - psi_ref_coef(i,s)*psi_ref_coef(i,s)
+! enddo
+! overlap = norm
+! f = dsqrt(f*norm)
+! do i=1,N_det_non_ref
+!   u_0(k) = rho_mrcc(i,1)*f
+!   overlap += u_0(k)*psi_non_ref_coef(i,1)
+! enddo
+!
+! call u_0_H_u_0(e_0,u_0,N_det,keys_tmp,N_int,1,N_det)
+! print *,  'Energy of |(1+T)Psi_0> (normalized mu_i) : ', e_0 + nuclear_repulsion, overlap
+!
+! deallocate(u_0, keys_tmp)
+!
+!---------------
      
      norm = 0.d0
-     double precision :: f
+     double precision               :: f
      do i=1,N_det_non_ref
        if (rho_mrcc(i,s) == 0.d0) then
          rho_mrcc(i,s) = 1.d-32
@@ -969,8 +1081,9 @@ END_PROVIDER
      ! rho_mrcc now contains the product of the scaling factors and the
      ! normalization constant
     
-     dIj_unique(:size(X), s) = X(:)
+    dIj_unique(:size(X), s) = X(:)
   end do
+
 END_PROVIDER
 
 

plugins/loc_cele/loc_cele.irp.f

@@ -451,7 +451,7 @@
 
       enddo !big loop over symmetry
 
-      10 format (4E18.12)
+      10 format (4E19.12)
 
 
 !  Now we copyt the newcmo into the mo_coef

plugins/mrcepa0/EZFIO.cfg

@@ -23,7 +23,7 @@ interface: ezfio
 type: Threshold
 doc: Threshold on the convergence of the dressed CI energy
 interface: ezfio,provider,ocaml
-default: 1.e-4
+default: 5.e-5
 
 [n_it_max_dressed_ci]
 type: Strictly_positive_int

plugins/mrcepa0/dressing.irp.f

@@ -299,7 +299,7 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,i_generator,n_selected,det_buffe
           delta_ii_(i_state,i_I)  = 0.d0
           do l_sd=1,idx_alpha(0)
             k_sd = idx_alpha(l_sd)
-            delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + dIa_hla(i_state,k_sd)
+            delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + 0.5d0*dIa_hla(i_state,k_sd)
           enddo
         endif
       enddo
@@ -554,7 +554,6 @@ END_PROVIDER
         do k=1,N_det_non_ref
 
           call i_h_j(psi_ref(1,1,j), psi_non_ref(1,1,k),N_int,Hjk)
-          call i_h_j(psi_non_ref(1,1,k),psi_ref(1,1,i), N_int,Hki)
           
           delta_cas(i,j,i_state) += Hjk * dij(i, k, i_state) ! * Hki * lambda_mrcc(i_state, k)
           !print *, Hjk * get_dij(psi_ref(1,1,i), psi_non_ref(1,1,k), N_int), Hki * get_dij(psi_ref(1,1,j), psi_non_ref(1,1,k), N_int)
@@ -647,7 +646,7 @@ end function
   integer                         :: p1,p2,h1,h2,s1,s2, p1_,p2_,h1_,h2_,s1_,s2_, sortRefIdx(N_det_ref)
   logical                         :: ok
   double precision                :: phase_iI, phase_Ik, phase_Jl, phase_IJ, phase_al, diI, hIi, hJi, delta_JI, dkI(1), HkI, ci_inv(1), dia_hla(1)
-  double precision                :: contrib, HIIi, HJk, wall
+  double precision                :: contrib, contrib2,  HIIi, HJk, wall
   integer, dimension(0:2,2,2)     :: exc_iI, exc_Ik, exc_IJ
   integer(bit_kind)               :: det_tmp(N_int, 2), made_hole(N_int,2), made_particle(N_int,2), myActive(N_int,2)
   integer(bit_kind),allocatable   :: sortRef(:,:,:)
@@ -677,7 +676,7 @@ end function
     delta_mrcepa0_ij(:,:,:) = 0d0
 
     !$OMP PARALLEL DO default(none) schedule(dynamic) shared(delta_mrcepa0_ij, delta_mrcepa0_ii)       &
-    !$OMP private(m,i,II,J,k,degree,myActive,made_hole,made_particle,hjk,contrib)       &
+    !$OMP private(m,i,II,J,k,degree,myActive,made_hole,made_particle,hjk,contrib,contrib2)       &
     !$OMP shared(active_sorb, psi_non_ref, psi_non_ref_coef, psi_ref, psi_ref_coef, cepa0_shortcut, det_cepa0_active)     &
     !$OMP shared(N_det_ref, N_det_non_ref,N_int,det_cepa0_idx,lambda_mrcc,det_ref_active, delta_cas) &
     !$OMP shared(notf,i_state, sortRef, sortRefIdx, dij)
@@ -720,16 +719,18 @@ end function
           !$OMP ATOMIC
           notf = notf+1
 
-          call i_h_j(psi_non_ref(1,1,det_cepa0_idx(k)),psi_ref(1,1,J),N_int,HJk)
+!          call i_h_j(psi_non_ref(1,1,det_cepa0_idx(k)),psi_ref(1,1,J),N_int,HJk)
-          !contrib = delta_cas(II, J, i_state) * HJk * lambda_mrcc(i_state, det_cepa0_idx(k))
           contrib = delta_cas(II, J, i_state) * dij(J, det_cepa0_idx(k), i_state)
-          !$OMP ATOMIC
-          delta_mrcepa0_ij(J, det_cepa0_idx(i), i_state) += contrib
           
           if(dabs(psi_ref_coef(J,i_state)).ge.5.d-5) then
+            contrib2 = contrib / psi_ref_coef(J, i_state) * psi_non_ref_coef(det_cepa0_idx(i),i_state)
             !$OMP ATOMIC
-            delta_mrcepa0_ii(J,i_state) -= contrib / psi_ref_coef(J, i_state) * psi_non_ref_coef(det_cepa0_idx(i),i_state)
+            delta_mrcepa0_ii(J,i_state) -= contrib2 
+          else
+            contrib = contrib * 0.5d0
           end if
+          !$OMP ATOMIC
+          delta_mrcepa0_ij(J, det_cepa0_idx(i), i_state) += contrib
 
         end do kloop
       end do
@@ -753,7 +754,7 @@ END_PROVIDER
   integer                         :: p1,p2,h1,h2,s1,s2, p1_,p2_,h1_,h2_,s1_,s2_
   logical                         :: ok
   double precision                :: phase_Ji, phase_Ik, phase_Ii
-  double precision                :: contrib, delta_IJk, HJk, HIk, HIl
+  double precision                :: contrib, contrib2, delta_IJk, HJk, HIk, HIl
   integer, dimension(0:2,2,2)     :: exc_Ik, exc_Ji, exc_Ii
   integer(bit_kind)               :: det_tmp(N_int, 2), det_tmp2(N_int, 2)
   integer, allocatable            :: idx_sorted_bit(:)
@@ -778,7 +779,7 @@ END_PROVIDER
     !$OMP PARALLEL DO default(none) schedule(dynamic,10) shared(delta_sub_ij, delta_sub_ii)       &
     !$OMP private(i, J, k, degree, degree2, l, deg, ni)       &
     !$OMP private(p1,p2,h1,h2,s1,s2, p1_,p2_,h1_,h2_,s1_,s2_)     &
-    !$OMP private(ok, phase_Ji, phase_Ik, phase_Ii, contrib, delta_IJk, HJk, HIk, HIl, exc_Ik, exc_Ji, exc_Ii) &
+    !$OMP private(ok, phase_Ji, phase_Ik, phase_Ii, contrib2, contrib, delta_IJk, HJk, HIk, HIl, exc_Ik, exc_Ji, exc_Ii) &
     !$OMP private(det_tmp, det_tmp2, II, blok)    &
     !$OMP shared(idx_sorted_bit, N_det_non_ref, N_det_ref, N_int, psi_non_ref, psi_non_ref_coef, psi_ref, psi_ref_coef)   &
     !$OMP shared(i_state,lambda_mrcc, hf_bitmask, active_sorb)
@@ -827,13 +828,16 @@ END_PROVIDER
             delta_IJk = HJk * HIk * lambda_mrcc(i_state, k)
             call apply_excitation(psi_non_ref(1,1,i),exc_Ik,det_tmp,ok,N_int)
             if(ok) cycle
-              contrib = delta_IJk * HIl * lambda_mrcc(i_state,l)   
+            contrib = delta_IJk * HIl * lambda_mrcc(i_state,l)   
+            if(dabs(psi_ref_coef(II,i_state)).ge.5.d-5) then
+              contrib2 = contrib / psi_ref_coef(II, i_state) * psi_non_ref_coef(l,i_state)
+              !$OMP ATOMIC
+              delta_sub_ii(II,i_state) -= contrib2
+            else
+              contrib = contrib * 0.5d0
+            endif
             !$OMP ATOMIC
             delta_sub_ij(II, i, i_state) += contrib
-            if(dabs(psi_ref_coef(II,i_state)).ge.5.d-5) then
-              !$OMP ATOMIC
-              delta_sub_ii(II,i_state) -= contrib / psi_ref_coef(II, i_state) * psi_non_ref_coef(l,i_state)
-            endif
           end do
         end do
       end do

src/Davidson/davidson_parallel.irp.f

@@ -501,7 +501,7 @@ subroutine davidson_miniserver_end()
   integer                 rc
   character*(64)          buf
   
-  address = trim(qp_run_address_tcp)//':11223'
+  address = trim(qp_run_address)//':11223'
   requester = f77_zmq_socket(zmq_context, ZMQ_REQ)
   rc        = f77_zmq_connect(requester,address)
 
@@ -520,7 +520,7 @@ subroutine davidson_miniserver_get()
   character*(20)          buffer
   integer                 rc
   
-  address = trim(qp_run_address_tcp)//':11223'
+  address = trim(qp_run_address)//':11223'
   
   requester = f77_zmq_socket(zmq_context, ZMQ_REQ)
   rc        = f77_zmq_connect(requester,address)

src/Determinants/connected_to_ref.irp.f

@@ -109,8 +109,6 @@ integer function get_index_in_psi_det_sorted_bit(key,Nint)
       continue
     else
       in_wavefunction = .True.
-      !DIR$ IVDEP
-      !DIR$ LOOP COUNT MIN(3)
       do l=2,Nint
         if ( (key(l,1) /= psi_det_sorted_bit(l,1,i)).or.                           &
               (key(l,2) /= psi_det_sorted_bit(l,2,i)) ) then
@@ -175,7 +173,6 @@ logical function is_connected_to(key,keys,Nint,Ndet)
   do i=1,Ndet
     degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) +              &
         popcnt(xor( key(1,2), keys(1,2,i)))
-    !DEC$ LOOP COUNT MIN(3)
     do l=2,Nint
       degree_x2 = degree_x2 + popcnt(xor( key(l,1), keys(l,1,i))) +&
           popcnt(xor( key(l,2), keys(l,2,i)))
@@ -231,7 +228,6 @@ logical function is_connected_to_by_mono(key,keys,Nint,Ndet)
   do i=1,Ndet
     degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) +              &
         popcnt(xor( key(1,2), keys(1,2,i)))
-    !DEC$ LOOP COUNT MIN(3)
     do l=2,Nint
       degree_x2 = degree_x2 + popcnt(xor( key(l,1), keys(l,1,i))) +&
           popcnt(xor( key(l,2), keys(l,2,i)))
@@ -325,10 +321,12 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet)
     do i=N_past-1,1,-1
       degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) +              &
           popcnt(xor( key(1,2), keys(1,2,i)))
-      !DEC$ LOOP COUNT MIN(3)
       do l=2,Nint
         degree_x2 = degree_x2 + popcnt(xor( key(l,1), keys(l,1,i))) +&
             popcnt(xor( key(l,2), keys(l,2,i)))
+        if (degree_x2 > 4) then
+          exit
+        endif
       enddo
       if (degree_x2 > 4) then
         cycle
@@ -429,7 +427,6 @@ integer function connected_to_ref_by_mono(key,keys,Nint,N_past_in,Ndet)
     do i=N_past-1,1,-1
       degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) +              &
           popcnt(xor( key(1,2), keys(1,2,i)))
-      !DEC$ LOOP COUNT MIN(3)
       do l=2,Nint
         degree_x2 = degree_x2 + popcnt(xor( key(l,1), keys(l,1,i))) +&
             popcnt(xor( key(l,2), keys(l,2,i)))

src/Determinants/filter_connected.irp.f

@@ -300,7 +300,6 @@ subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx)
   else
     
 
-    integer, save :: icount(4) = (/0,0,0,0/)
     !DIR$ LOOP COUNT (1000)
     outer: do i=1,sze
       degree_x2 = 0
@@ -318,7 +317,6 @@ subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx)
       enddo
       idx(l) = i
       l = l+1
-      icount(3) = icount(3) + 1_8
     enddo outer
     
   endif

src/Determinants/slater_rules.irp.f

@@ -1320,7 +1320,6 @@ subroutine get_excitation_degree_vector(key1,key2,degree,Nint,sze,idx)
   l=1
   if (Nint==1) then
     
-    !DIR$ LOOP COUNT (1000)
     do i=1,sze
       d = popcnt(xor( key1(1,1,i), key2(1,1))) +       &
           popcnt(xor( key1(1,2,i), key2(1,2)))
@@ -1335,7 +1334,6 @@ subroutine get_excitation_degree_vector(key1,key2,degree,Nint,sze,idx)
     
   else if (Nint==2) then
     
-    !DIR$ LOOP COUNT (1000)
     do i=1,sze
       d = popcnt(xor( key1(1,1,i), key2(1,1))) +                     &
           popcnt(xor( key1(1,2,i), key2(1,2))) +                     &
@@ -1352,7 +1350,6 @@ subroutine get_excitation_degree_vector(key1,key2,degree,Nint,sze,idx)
     
   else if (Nint==3) then
     
-    !DIR$ LOOP COUNT (1000)
     do i=1,sze
       d = popcnt(xor( key1(1,1,i), key2(1,1))) +                     &
           popcnt(xor( key1(1,2,i), key2(1,2))) +                     &
@@ -1371,10 +1368,8 @@ subroutine get_excitation_degree_vector(key1,key2,degree,Nint,sze,idx)
     
   else
     
-    !DIR$ LOOP COUNT (1000)
     do i=1,sze
       d = 0
-      !DIR$ LOOP COUNT MIN(4)
       do m=1,Nint
         d = d + popcnt(xor( key1(m,1,i), key2(m,1)))                 &
               + popcnt(xor( key1(m,2,i), key2(m,2)))

src/ZMQ/utils.irp.f

@@ -17,8 +17,6 @@ END_PROVIDER
 
 
  BEGIN_PROVIDER [ character*(128), qp_run_address ]
-&BEGIN_PROVIDER [ character*(128), qp_run_address_ipc ]
-&BEGIN_PROVIDER [ character*(128), qp_run_address_tcp ]
 &BEGIN_PROVIDER [ integer, zmq_port_start ]
   use f77_zmq
   implicit none
@@ -36,22 +34,19 @@ END_PROVIDER
   integer                        :: i
   do i=len(buffer),1,-1
     if ( buffer(i:i) == ':') then
-      qp_run_address_tcp = trim(buffer(1:i-1))
+      qp_run_address = trim(buffer(1:i-1))
       read(buffer(i+1:), *) zmq_port_start
       exit
     endif
   enddo
-  qp_run_address_ipc = 'ipc:///tmp/qp_run'
-  qp_run_address = qp_run_address_ipc
 END_PROVIDER
 
-
  BEGIN_PROVIDER [ character*(128), zmq_socket_pull_tcp_address    ]
-&BEGIN_PROVIDER [ character*(128), zmq_socket_pull_inproc_address ]
 &BEGIN_PROVIDER [ character*(128), zmq_socket_pair_inproc_address ]
 &BEGIN_PROVIDER [ character*(128), zmq_socket_push_tcp_address    ]
+&BEGIN_PROVIDER [ character*(128), zmq_socket_pull_inproc_address ]
 &BEGIN_PROVIDER [ character*(128), zmq_socket_push_inproc_address ]
-&BEGIN_PROVIDER [ character*(128), zmq_socket_sub_address     ]
+&BEGIN_PROVIDER [ character*(128), zmq_socket_sub_tcp_address ]
   use f77_zmq
   implicit none
   BEGIN_DOC
@@ -59,12 +54,12 @@ END_PROVIDER
   END_DOC
   character*(8), external        :: zmq_port
 
+  zmq_socket_sub_tcp_address     = trim(qp_run_address)//':'//zmq_port(1)//' '
   zmq_socket_pull_tcp_address    = 'tcp://*:'//zmq_port(2)//' '
+  zmq_socket_push_tcp_address    = trim(qp_run_address)//':'//zmq_port(2)//' '
   zmq_socket_pull_inproc_address = 'inproc://'//zmq_port(2)//' '
-  zmq_socket_pair_inproc_address = 'inproc://'//zmq_port(3)//' '
-  zmq_socket_push_tcp_address    = trim(qp_run_address_tcp)//':'//zmq_port(2)//' '
   zmq_socket_push_inproc_address = zmq_socket_pull_inproc_address
-  zmq_socket_sub_address         = trim(qp_run_address)//':'//zmq_port(1)//' '
+  zmq_socket_pair_inproc_address = 'inproc://'//zmq_port(3)//' '
 
   ! /!\ Don't forget to change subroutine reset_zmq_addresses
 END_PROVIDER
@@ -77,13 +72,12 @@ subroutine reset_zmq_addresses
   END_DOC
   character*(8), external        :: zmq_port
 
+  zmq_socket_sub_tcp_address     = trim(qp_run_address)//':'//zmq_port(1)//' '
   zmq_socket_pull_tcp_address    = 'tcp://*:'//zmq_port(2)//' '
+  zmq_socket_push_tcp_address    = trim(qp_run_address)//':'//zmq_port(2)//' '
   zmq_socket_pull_inproc_address = 'inproc://'//zmq_port(2)//' '
-  zmq_socket_pair_inproc_address = 'inproc://'//zmq_port(3)//' '
-  zmq_socket_push_tcp_address    = trim(qp_run_address_tcp)//':'//zmq_port(2)//' '
   zmq_socket_push_inproc_address = zmq_socket_pull_inproc_address
-  zmq_socket_sub_address         = trim(qp_run_address)//':'//zmq_port(1)//' '
+  zmq_socket_pair_inproc_address = 'inproc://'//zmq_port(3)//' '
-
 end
 
 
@@ -111,7 +105,6 @@ subroutine switch_qp_run_to_master
       exit
     endif
   enddo
-  qp_run_address_tcp = qp_run_address
   call reset_zmq_addresses
 
 end
@@ -374,7 +367,7 @@ function new_zmq_sub_socket()
     stop 'Unable to subscribe new_zmq_sub_socket'
   endif
 
-  rc = f77_zmq_connect(new_zmq_sub_socket, zmq_socket_sub_address)
+  rc = f77_zmq_connect(new_zmq_sub_socket, zmq_socket_sub_tcp_address)
   if (rc /= 0) then
     stop 'Unable to connect new_zmq_sub_socket'
   endif