Merge branch 'develop'

ocaml/Pseudo.ml

@@ -124,23 +124,27 @@ let to_string t =
 let find in_channel element =
   In_channel.seek in_channel 0L;
   
-  let element_read, old_pos =
-     ref Element.X, 
+  let loop, element_read, old_pos =
+     ref true,
+     ref None, 
      ref (In_channel.pos in_channel)
   in
-  while !element_read <> element
+
+  while !loop
   do
+    try
       let buffer = 
         old_pos := In_channel.pos in_channel;
         match In_channel.input_line in_channel with
         | Some line -> String.split ~on:' ' line
           |> List.hd_exn
-      | None -> ""
+        | None -> raise End_of_file
       in
-    try
-      element_read := Element.of_string buffer
+      element_read := Some (Element.of_string buffer);
+      loop := !element_read <> (Some element) 
     with
     | Element.ElementError _ -> ()
+    | End_of_file -> loop := false
   done ;
   In_channel.seek in_channel !old_pos;
   !element_read
@@ -148,8 +152,9 @@ let find in_channel element =
 
 (** Read the Pseudopotential in GAMESS format *)
 let read_element in_channel element =
-  ignore (find in_channel element);
-
+  match find in_channel element with
+  | Some e when e = element -> 
+    begin
       let rec read result =
         match In_channel.input_line in_channel with
         | None -> result
@@ -264,7 +269,8 @@ let read_element in_channel element =
       decode_first_line data
       |> decode_local
       |> decode_non_local 
-
+    end
+  | _ -> empty element
 
  
 

plugins/MRPT_Utils/ezfio_interface.irp.f

@@ -0,0 +1,23 @@
+! DO NOT MODIFY BY HAND
+! Created by $QP_ROOT/scripts/ezfio_interface/ei_handler.py
+! from file /home/scemama/quantum_package/src/MRPT_Utils/EZFIO.cfg
+
+
+BEGIN_PROVIDER [ logical, do_third_order_1h1p  ]
+  implicit none
+  BEGIN_DOC
+! If true, compute the third order contribution for the 1h1p
+  END_DOC
+
+  logical                        :: has
+  PROVIDE ezfio_filename
+  
+  call ezfio_has_mrpt_utils_do_third_order_1h1p(has)
+  if (has) then
+    call ezfio_get_mrpt_utils_do_third_order_1h1p(do_third_order_1h1p)
+  else
+    print *, 'mrpt_utils/do_third_order_1h1p not found in EZFIO file'
+    stop 1
+  endif
+
+END_PROVIDER

plugins/Psiref_threshold/psi_ref.irp.f

@@ -6,19 +6,22 @@ use bitmasks
 &BEGIN_PROVIDER [ integer, N_det_ref ]
   implicit none
   BEGIN_DOC
-  ! Reference wave function, defined as determinants with coefficients > 0.05
+  ! Reference wave function, defined as determinants with amplitudes > 0.05
   ! idx_ref gives the indice of the ref determinant in psi_det.
   END_DOC
   integer                        :: i, k, l
   logical                        :: good
   double precision, parameter    :: threshold=0.05d0
+  double precision               :: t(N_states)
   N_det_ref = 0
-  t = threshold * abs_psi_coef_max
+  do l = 1, N_states
+    t(l) = threshold * abs_psi_coef_max(l)
+  enddo
   do i=1,N_det
     good = .False.
-    do l = 1, N_states
+    do l=1, N_states
      psi_ref_coef(i,l) = 0.d0
-     good = good.or.(dabs(psi_coef(i,l)) > t)
+     good = good.or.(dabs(psi_coef(i,l)) > t(l))
     enddo
     if (good) then
       N_det_ref = N_det_ref+1

plugins/QmcChem/e_curve_qmc.irp.f

@@ -1,7 +1,7 @@
 program e_curve
   use bitmasks
   implicit none
-  integer :: i,j,k, nab, m, l
+  integer :: i,j,k, kk, nab, m, l
   double precision :: norm, E, hij, num, ci, cj
   integer, allocatable :: iorder(:)
   double precision , allocatable :: norm_sort(:)
@@ -60,7 +60,7 @@ program e_curve
     num = 0.d0
     norm = 0.d0
     m = 0
-    !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(k,l,det_i,det_j,ci,cj,hij) REDUCTION(+:norm,m,num)
+    !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(k,kk,l,det_i,det_j,ci,cj,hij) REDUCTION(+:norm,m,num)
     allocate( det_i(N_int,2), det_j(N_int,2))
     !$OMP DO SCHEDULE(guided)
     do k=1,n_det
@@ -68,15 +68,19 @@ program e_curve
         cycle
       endif
       ci = psi_bilinear_matrix_values(k,1)
-      det_i(:,1) = psi_det_alpha_unique(:,psi_bilinear_matrix_rows(k))
-      det_i(:,2) = psi_det_beta_unique(:,psi_bilinear_matrix_columns(k))
+      do kk=1,N_int
+        det_i(kk,1) = psi_det_alpha_unique(kk,psi_bilinear_matrix_rows(k))
+        det_i(kk,2) = psi_det_beta_unique(kk,psi_bilinear_matrix_columns(k))
+      enddo
       do l=1,n_det
         if (psi_bilinear_matrix_values(l,1) == 0.d0) then
           cycle
         endif
         cj = psi_bilinear_matrix_values(l,1)
-        det_j(:,1) = psi_det_alpha_unique(:,psi_bilinear_matrix_rows(l))
-        det_j(:,2) = psi_det_beta_unique(:,psi_bilinear_matrix_columns(l))
+        do kk=1,N_int
+          det_j(kk,1) = psi_det_alpha_unique(kk,psi_bilinear_matrix_rows(l))
+          det_j(kk,2) = psi_det_beta_unique(kk,psi_bilinear_matrix_columns(l))
+        enddo
         call i_h_j(det_i, det_j, N_int, hij)
         num = num + ci*cj*hij
       enddo

plugins/mrcc_selected/dressing_slave.irp.f

@@ -0,0 +1,593 @@
+subroutine mrsc2_dressing_slave_tcp(i)
+  implicit none
+  integer, intent(in)            :: i
+  BEGIN_DOC
+! Task for parallel MR-SC2
+  END_DOC
+  call mrsc2_dressing_slave(0,i)
+end
+
+
+subroutine mrsc2_dressing_slave_inproc(i)
+  implicit none
+  integer, intent(in)            :: i
+  BEGIN_DOC
+! Task for parallel MR-SC2
+  END_DOC
+  call mrsc2_dressing_slave(1,i)
+end
+
+subroutine mrsc2_dressing_slave(thread,iproc)
+  use f77_zmq
+
+  implicit none
+  BEGIN_DOC
+! Task for parallel MR-SC2
+  END_DOC
+  integer,  intent(in)            :: thread, iproc
+!   integer                        :: j,l
+  integer                        :: rc 
+
+  integer                        :: worker_id, task_id 
+  character*(512)                :: task 
+ 
+  integer(ZMQ_PTR),external      :: new_zmq_to_qp_run_socket 
+  integer(ZMQ_PTR)               :: zmq_to_qp_run_socket 
+ 
+  integer(ZMQ_PTR), external     :: new_zmq_push_socket 
+  integer(ZMQ_PTR)               :: zmq_socket_push 
+
+  double precision, allocatable  :: delta(:,:,:) 
+  
+
+
+  integer                         :: i_state, i, i_I, J, k, k2, k1, kk, ll, degree, degree2, m, l, deg, ni, m2
+  integer                         :: n(2)
+  integer                         :: p1,p2,h1,h2,s1,s2, blok, I_s, J_s, kn
+  logical                         :: ok
+  double precision                :: phase_iI, phase_Ik, phase_Jl, phase_Ji, phase_al
+  double precision                :: diI, hIi, hJi, delta_JI, dkI, HkI, ci_inv(N_states), cj_inv(N_states)
+  double precision                :: contrib, wall, iwall
+  double precision, allocatable   :: dleat(:,:,:)
+  integer, dimension(0:2,2,2)     :: exc_iI, exc_Ik, exc_IJ
+  integer(bit_kind)               :: det_tmp(N_int, 2), det_tmp2(N_int, 2), inac, virt
+  integer, external               :: get_index_in_psi_det_sorted_bit, searchDet, detCmp
+  logical, external               :: is_in_wavefunction, isInCassd, detEq
+  integer,allocatable :: komon(:)
+  logical :: komoned
+  !double precision, external :: get_dij
+     
+  zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
+  zmq_socket_push      = new_zmq_push_socket(thread)
+
+  call connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread)
+
+  allocate (dleat(N_states, N_det_non_ref, 2), delta(N_states,0:N_det_non_ref, 2))
+  allocate(komon(0:N_det_non_ref))
+
+  do 
+    call get_task_from_taskserver(zmq_to_qp_run_socket,worker_id, task_id, task)
+    if (task_id == 0) exit
+    read (task,*) i_I, J, k1, k2
+    do i_state=1, N_states
+      ci_inv(i_state) = 1.d0 / psi_ref_coef(i_I,i_state)
+      cj_inv(i_state) = 1.d0 / psi_ref_coef(J,i_state)
+    end do
+    !delta = 0.d0
+    n = 0
+    delta(:,0,:) = 0d0
+    delta(:,:nlink(J),1) = 0d0
+    delta(:,:nlink(i_I),2) = 0d0
+    komon(0) = 0
+    komoned = .false.
+    
+    
+    
+    
+    do kk = k1, k2
+      k = det_cepa0_idx(linked(kk, i_I))
+      blok = blokMwen(kk, i_I)
+      
+      call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,k),exc_Ik,degree,phase_Ik,N_int)
+    
+      if(J /= i_I) then
+        call apply_excitation(psi_ref(1,1,J),exc_Ik,det_tmp2,ok,N_int)
+        if(.not. ok) cycle
+          
+        l = searchDet(det_cepa0(1,1,cepa0_shortcut(blok)), det_tmp2, cepa0_shortcut(blok+1)-cepa0_shortcut(blok), N_int)
+        if(l == -1) cycle
+        ll = cepa0_shortcut(blok)-1+l
+        l = det_cepa0_idx(ll)
+        ll = child_num(ll, J)
+      else
+        l = k
+        ll = kk
+      end if
+      
+        
+      if(.not. komoned) then
+        m = 0
+        m2 = 0
+          
+        do while(m < nlink(i_I) .and. m2 < nlink(J))
+          m += 1
+          m2 += 1
+          if(linked(m, i_I) < linked(m2, J)) then
+            m2 -= 1
+            cycle
+          else if(linked(m, i_I) > linked(m2, J)) then
+            m -= 1
+            cycle
+          end if
+          i = det_cepa0_idx(linked(m, i_I))
+          
+          if(h_(J,i) == 0.d0) cycle
+          if(h_(i_I,i) == 0.d0) cycle
+          
+          !ok = .false.
+          !do i_state=1, N_states
+          !  if(lambda_mrcc(i_state, i) /= 0d0) then
+          !    ok = .true.
+          !    exit
+          !  end if
+          !end do
+          !if(.not. ok) cycle
+!         
+          
+          komon(0) += 1
+          kn = komon(0)
+          komon(kn) = i
+          
+          
+!           call get_excitation(psi_ref(1,1,J),psi_non_ref(1,1,i),exc_IJ,degree2,phase_Ji,N_int)
+!           if(I_i /= J) call get_excitation(psi_ref(1,1,I_i),psi_non_ref(1,1,i),exc_IJ,degree2,phase_Ii,N_int)
+!           if(I_i == J) phase_Ii = phase_Ji
+          
+          do i_state = 1,N_states
+            dkI = h_(J,i) * dij(i_I, i, i_state)!get_dij(psi_ref(1,1,i_I), psi_non_ref(1,1,i), N_int)
+            !dkI = h_(J,i) * h_(i_I,i) * lambda_mrcc(i_state, i)
+            dleat(i_state, kn, 1) = dkI
+            dleat(i_state, kn, 2) = dkI
+          end do
+
+        end do
+          
+        komoned = .true.
+      end if
+      
+      
+      do m = 1, komon(0)
+        
+        i = komon(m)
+        
+        call apply_excitation(psi_non_ref(1,1,i),exc_Ik,det_tmp,ok,N_int)
+        if(.not. ok) cycle
+        if(HP(1,i) + HP(1,k) <= 2 .and. HP(2,i) + HP(2,k) <= 2) then
+!           if(is_in_wavefunction(det_tmp, N_int)) cycle
+          cycle
+        end if
+        
+        !if(isInCassd(det_tmp, N_int)) cycle
+          
+        do i_state = 1, N_states 
+          !if(lambda_mrcc(i_state, i) == 0d0) cycle
+          
+
+          !contrib = h_(i_I,k) * lambda_mrcc(i_state, k) * dleat(i_state, m, 2)! * phase_al
+          contrib =  dij(i_I, k, i_state) * dleat(i_state, m, 2)
+          delta(i_state,ll,1) += contrib
+          if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5) then
+            delta(i_state,0,1) -= contrib * ci_inv(i_state) * psi_non_ref_coef(l,i_state)
+          endif
+          
+          if(I_i == J) cycle
+          !contrib = h_(J,l) * lambda_mrcc(i_state, l) * dleat(i_state, m, 1)! * phase_al
+          contrib =  dij(J, l, i_state) * dleat(i_state, m, 1)
+          delta(i_state,kk,2) += contrib
+          if(dabs(psi_ref_coef(J,i_state)).ge.5.d-5) then
+            delta(i_state,0,2) -= contrib * cj_inv(i_state) * psi_non_ref_coef(k,i_state)
+          end if
+        enddo !i_state
+      end do ! while
+    end do ! kk
+
+      
+      call push_mrsc2_results(zmq_socket_push, I_i, J, delta, task_id) 
+      call task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id)
+ 
+!     end if
+    
+  enddo
+
+  deallocate(delta)
+
+  call disconnect_from_taskserver(zmq_to_qp_run_socket,zmq_socket_push,worker_id)
+  call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
+  call end_zmq_push_socket(zmq_socket_push,thread)
+
+end
+
+
+subroutine push_mrsc2_results(zmq_socket_push, I_i, J, delta, task_id) 
+  use f77_zmq 
+  implicit none 
+  BEGIN_DOC 
+! Push integrals in the push socket 
+  END_DOC 
+
+  integer, intent(in)            :: i_I, J
+  integer(ZMQ_PTR), intent(in)   :: zmq_socket_push 
+  double precision,intent(inout) :: delta(N_states, 0:N_det_non_ref, 2)
+  integer, intent(in)            :: task_id 
+  integer                        :: rc , i_state,  i, kk, li
+  integer,allocatable            :: idx(:,:)
+  integer                        :: n(2)
+  logical :: ok
+  
+  allocate(idx(N_det_non_ref,2))
+  rc = f77_zmq_send( zmq_socket_push, i_I, 4, ZMQ_SNDMORE) 
+  if (rc /= 4) then 
+    print *, irp_here,  'f77_zmq_send( zmq_socket_push, i_I, 4, ZMQ_SNDMORE)' 
+    stop 'error' 
+  endif 
+
+  rc = f77_zmq_send( zmq_socket_push, J, 4, ZMQ_SNDMORE) 
+  if (rc /= 4) then 
+    print *, irp_here,  'f77_zmq_send( zmq_socket_push, J, 4, ZMQ_SNDMORE)' 
+    stop 'error' 
+  endif 
+  
+  
+  do kk=1,2
+    n(kk)=0
+    if(kk == 1) li = nlink(j)
+    if(kk == 2) li = nlink(i_I)
+    do i=1, li
+      ok = .false.
+      do i_state=1,N_states
+        if(delta(i_state, i, kk) /= 0d0) then
+          ok = .true.
+          exit
+        end if
+      end do
+      
+      if(ok) then
+        n(kk) += 1
+!         idx(n,kk) = i
+        if(kk == 1) then
+          idx(n(1),1) = det_cepa0_idx(linked(i, J))
+        else
+          idx(n(2),2) = det_cepa0_idx(linked(i, i_I))
+        end if
+        
+        do i_state=1, N_states
+            delta(i_state, n(kk), kk) = delta(i_state, i, kk)
+        end do
+      end if
+    end do
+    
+    rc = f77_zmq_send( zmq_socket_push, n(kk), 4, ZMQ_SNDMORE) 
+    if (rc /= 4) then 
+      print *, irp_here,  'f77_zmq_send( zmq_socket_push, n, 4, ZMQ_SNDMORE)' 
+      stop 'error' 
+    endif 
+    
+    if(n(kk) /= 0) then
+      rc = f77_zmq_send( zmq_socket_push, delta(1,0,kk), (n(kk)+1)*8*N_states, ZMQ_SNDMORE)  ! delta(1,0,1) = delta_I   delta(1,0,2) = delta_J 
+      if (rc /=  (n(kk)+1)*8*N_states) then 
+        print *, irp_here,  'f77_zmq_send( zmq_socket_push, delta, (n(kk)+1)*8*N_states, ZMQ_SNDMORE)'
+        stop 'error' 
+      endif 
+      
+      rc = f77_zmq_send( zmq_socket_push, idx(1,kk), n(kk)*4, ZMQ_SNDMORE) 
+      if (rc /=  n(kk)*4) then 
+        print *, irp_here,  'f77_zmq_send( zmq_socket_push, delta, 8*n(kk), ZMQ_SNDMORE)'
+        stop 'error' 
+      endif 
+    end if
+  end do
+  
+  
+  rc = f77_zmq_send( zmq_socket_push, task_id, 4, 0) 
+  if (rc /= 4) then 
+    print *, irp_here,  'f77_zmq_send( zmq_socket_push, task_id, 4, 0)' 
+    stop 'error' 
+  endif 
+ 
+! ! Activate is zmq_socket_push is a REQ 
+!   integer :: idummy 
+!   rc = f77_zmq_recv( zmq_socket_push, idummy, 4, 0) 
+!   if (rc /= 4) then 
+!     print *, irp_here, 'f77_zmq_send( zmq_socket_push, idummy, 4, 0)' 
+!     stop 'error' 
+!   endif 
+end
+
+
+
+subroutine pull_mrsc2_results(zmq_socket_pull, I_i, J, n, idx, delta, task_id) 
+  use f77_zmq 
+  implicit none 
+  BEGIN_DOC 
+! Push integrals in the push socket 
+  END_DOC 
+
+  integer(ZMQ_PTR), intent(in)   :: zmq_socket_pull 
+  integer, intent(out)           :: i_I, J, n(2)
+  double precision, intent(inout) :: delta(N_states, 0:N_det_non_ref, 2)
+  integer, intent(out)           :: task_id 
+  integer                        :: rc , i, kk
+  integer,intent(inout) :: idx(N_det_non_ref,2)
+  logical :: ok
+  
+  rc = f77_zmq_recv( zmq_socket_pull, i_I, 4, ZMQ_SNDMORE) 
+  if (rc /= 4) then 
+    print *, irp_here,  'f77_zmq_recv( zmq_socket_pull, i_I, 4, ZMQ_SNDMORE)' 
+    stop 'error' 
+  endif 
+  
+  rc = f77_zmq_recv( zmq_socket_pull, J, 4, ZMQ_SNDMORE) 
+  if (rc /= 4) then 
+    print *, irp_here,  'f77_zmq_recv( zmq_socket_pull, J, 4, ZMQ_SNDMORE)' 
+    stop 'error' 
+  endif 
+  
+  do kk = 1, 2
+    rc = f77_zmq_recv( zmq_socket_pull, n(kk), 4, ZMQ_SNDMORE) 
+    if (rc /= 4) then 
+      print *, irp_here,  'f77_zmq_recv( zmq_socket_pull, n, 4, ZMQ_SNDMORE)' 
+      stop 'error' 
+    endif 
+    
+    if(n(kk) /= 0) then
+      rc = f77_zmq_recv( zmq_socket_pull, delta(1,0,kk), (n(kk)+1)*8*N_states, ZMQ_SNDMORE) 
+      if (rc /= (n(kk)+1)*8*N_states) then 
+        print *, irp_here,  'f77_zmq_recv( zmq_socket_pull, delta, (n(kk)+1)*8*N_states, ZMQ_SNDMORE)'
+        stop 'error' 
+      endif 
+      
+      rc = f77_zmq_recv( zmq_socket_pull, idx(1,kk), n(kk)*4, ZMQ_SNDMORE) 
+      if (rc /= n(kk)*4) then 
+        print *, irp_here,  'f77_zmq_recv( zmq_socket_pull, delta, n(kk)*4, ZMQ_SNDMORE)'
+        stop 'error' 
+      endif 
+    end if
+  end do
+  
+  rc = f77_zmq_recv( zmq_socket_pull, task_id, 4, 0) 
+  if (rc /= 4) then 
+    print *, irp_here,  'f77_zmq_recv( zmq_socket_pull, task_id, 4, 0)' 
+    stop 'error' 
+  endif 
+
+ 
+! ! Activate is zmq_socket_pull is a REP 
+!   integer :: idummy 
+!   rc = f77_zmq_send( zmq_socket_pull, idummy, 4, 0) 
+!   if (rc /= 4) then 
+!     print *, irp_here, 'f77_zmq_send( zmq_socket_pull, idummy, 4, 0)' 
+!     stop 'error' 
+!   endif 
+end
+
+
+
+subroutine mrsc2_dressing_collector(delta_ii_,delta_ij_)
+  use f77_zmq
+  implicit none
+  BEGIN_DOC 
+! Collects results from the AO integral calculation 
+  END_DOC 
+ 
+  double precision,intent(inout)               :: delta_ij_(N_states,N_det_non_ref,N_det_ref) 
+  double precision,intent(inout)               :: delta_ii_(N_states,N_det_ref)
+
+!   integer                        :: j,l
+  integer                        :: rc 
+   
+  double precision, allocatable  :: delta(:,:,:) 
+   
+  integer(ZMQ_PTR),external      :: new_zmq_to_qp_run_socket 
+  integer(ZMQ_PTR)               :: zmq_to_qp_run_socket 
+   
+  integer(ZMQ_PTR), external     :: new_zmq_pull_socket 
+  integer(ZMQ_PTR)               :: zmq_socket_pull 
+   
+  integer*8                      :: control, accu 
+  integer                        :: task_id, more 
+ 
+  integer                        :: I_i, J, l, i_state, n(2), kk
+  integer,allocatable :: idx(:,:)
+  
+  delta_ii_(:,:) = 0d0
+  delta_ij_(:,:,:) = 0d0
+
+  zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() 
+  zmq_socket_pull = new_zmq_pull_socket() 
+ 
+  allocate ( delta(N_states,0:N_det_non_ref,2) ) 
+  
+  allocate(idx(N_det_non_ref,2))
+  more = 1 
+  do while (more == 1) 
+         
+    call pull_mrsc2_results(zmq_socket_pull, I_i, J, n, idx, delta, task_id)
+    
+
+      do l=1, n(1)
+        do i_state=1,N_states
+          delta_ij_(i_state,idx(l,1),i_I) += delta(i_state,l,1)
+        end do
+      end do
+      
+      do l=1, n(2)
+        do i_state=1,N_states
+          delta_ij_(i_state,idx(l,2),J) += delta(i_state,l,2)
+        end do
+      end do
+
+    
+!     
+!       do l=1,nlink(J)
+!         do i_state=1,N_states
+!           delta_ij_(i_state,det_cepa0_idx(linked(l,J)),i_I) += delta(i_state,l,1)
+!           delta_ij_(i_state,det_cepa0_idx(linked(l,i_I)),j) += delta(i_state,l,2)
+!         end do
+!       end do
+!      
+       if(n(1) /= 0) then 
+       do i_state=1,N_states
+         delta_ii_(i_state,i_I) += delta(i_state,0,1)
+       end do
+       end if
+
+      if(n(2) /= 0) then 
+       do i_state=1,N_states
+         delta_ii_(i_state,J) += delta(i_state,0,2)
+       end do
+       end if
+
+
+    if (task_id /= 0) then 
+      call zmq_delete_task(zmq_to_qp_run_socket,zmq_socket_pull,task_id,more) 
+    endif 
+
+ 
+  enddo 
+  deallocate( delta )
+ 
+  call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
+  call end_zmq_pull_socket(zmq_socket_pull)
+
+end
+
+
+
+
+ BEGIN_PROVIDER [ double precision, delta_ij_old, (N_states,N_det_non_ref,N_det_ref) ]
+&BEGIN_PROVIDER [ double precision, delta_ii_old, (N_states,N_det_ref) ]
+  implicit none
+  
+  integer                         :: i_state, i, i_I, J, k, kk, degree, degree2, m, l, deg, ni, m2
+  integer                         :: p1,p2,h1,h2,s1,s2, blok, I_s, J_s, nex, nzer, ntot
+!   integer, allocatable            :: linked(:,:), blokMwen(:, :), nlink(:)
+  logical                         :: ok
+  double precision                :: phase_iI, phase_Ik, phase_Jl, phase_Ji, phase_al, diI, hIi, hJi, delta_JI, dkI(N_states), HkI, ci_inv(N_states), dia_hla(N_states)
+  double precision                :: contrib, wall, iwall ! , searchance(N_det_ref)
+  integer, dimension(0:2,2,2)     :: exc_iI, exc_Ik, exc_IJ
+  integer(bit_kind)               :: det_tmp(N_int, 2), det_tmp2(N_int, 2), inac, virt
+  integer, external               :: get_index_in_psi_det_sorted_bit, searchDet, detCmp
+  logical, external               :: is_in_wavefunction, isInCassd, detEq
+  character*(512)                :: task 
+  integer(ZMQ_PTR)               :: zmq_to_qp_run_socket 
+  
+  integer :: KKsize = 1000000
+  
+  
+  call new_parallel_job(zmq_to_qp_run_socket,'mrsc2')
+
+
+  call wall_time(iwall)
+!   allocate(linked(N_det_non_ref, N_det_ref), blokMwen(N_det_non_ref, N_det_ref), nlink(N_det_ref))
+  
+
+!   searchance = 0d0
+!   do J = 1, N_det_ref
+!     nlink(J) = 0
+!     do blok=1,cepa0_shortcut(0)
+!     do k=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1
+!       call get_excitation_degree(psi_ref(1,1,J),det_cepa0(1,1,k),degree,N_int)
+!       if(degree <= 2) then
+!         nlink(J) += 1
+!         linked(nlink(J),J) = k
+!         blokMwen(nlink(J),J) = blok
+!         searchance(J) += 1d0 + log(dfloat(cepa0_shortcut(blok+1) - cepa0_shortcut(blok)))
+!       end if
+!     end do
+!     end do
+!   end do
+
+     
+     
+!   stop
+  nzer = 0
+  ntot = 0
+  do nex = 3, 0, -1
+    print *, "los ",nex
+    do I_s = N_det_ref, 1, -1
+!         if(mod(I_s,1) == 0) then
+!           call wall_time(wall)
+!           wall = wall-iwall
+!           print *, I_s, "/", N_det_ref, wall * (dfloat(N_det_ref) / dfloat(I_s)), wall, wall * (dfloat(N_det_ref) / dfloat(I_s))-wall
+!         end if
+
+
+      do J_s = 1, I_s
+      
+        call get_excitation_degree(psi_ref(1,1,J_s), psi_ref(1,1,I_s), degree, N_int)
+        if(degree /= nex) cycle
+        if(nex == 3) nzer = nzer + 1
+        ntot += 1
+!           if(degree > 3) then
+!             deg += 1
+!             cycle
+!           else if(degree == -10) then
+!             KKsize = 100000
+!           else
+!             KKsize = 1000000
+!           end if
+        
+        
+        
+        if(searchance(I_s) < searchance(J_s)) then
+          i_I = I_s
+          J = J_s
+        else
+          i_I = J_s
+          J = I_s
+        end if
+        
+        KKsize = nlink(1)
+        if(nex == 0) KKsize = int(float(nlink(1)) / float(nlink(i_I)) * (float(nlink(1)) / 64d0))
+        
+        !if(KKsize == 0) stop "ZZEO"
+        
+        do kk = 1 , nlink(i_I), KKsize
+          write(task,*) I_i, J, kk, int(min(kk+KKsize-1, nlink(i_I)))
+          call add_task_to_taskserver(zmq_to_qp_run_socket,task)
+        end do
+        
+  !       do kk = 1 , nlink(i_I)
+  !         k = linked(kk,i_I)
+  !         blok = blokMwen(kk,i_I)
+  !         write(task,*) I_i, J, k, blok
+  !         call add_task_to_taskserver(zmq_to_qp_run_socket,task)
+  !         
+  !       enddo !kk
+      enddo !J
+          
+    enddo !I
+  end do ! nex
+  print *, "tasked"
+!  integer(ZMQ_PTR)               ∷ collector_thread
+!  external                       ∷ ao_bielec_integrals_in_map_collector
+!  rc = pthread_create(collector_thread, mrsc2_dressing_collector)
+  print *, nzer, ntot, float(nzer) / float(ntot)
+  provide nproc
+  !$OMP PARALLEL DEFAULT(none)  SHARED(delta_ii_old,delta_ij_old)  PRIVATE(i) NUM_THREADS(nproc+1)
+      i = omp_get_thread_num()
+      if (i==0) then
+        call mrsc2_dressing_collector(delta_ii_old,delta_ij_old)
+      else
+        call mrsc2_dressing_slave_inproc(i)
+      endif
+  !$OMP END PARALLEL
+
+!  rc = pthread_join(collector_thread)
+  call end_parallel_job(zmq_to_qp_run_socket, 'mrsc2')
+    
+
+END_PROVIDER
+
+
+

plugins/mrcc_selected/ezfio_interface.irp.f

@@ -0,0 +1,61 @@
+! DO NOT MODIFY BY HAND
+! Created by $QP_ROOT/scripts/ezfio_interface/ei_handler.py
+! from file /home/scemama/quantum_package/src/mrcc_selected/EZFIO.cfg
+
+
+BEGIN_PROVIDER [ double precision, thresh_dressed_ci  ]
+  implicit none
+  BEGIN_DOC
+! Threshold on the convergence of the dressed CI energy
+  END_DOC
+
+  logical                        :: has
+  PROVIDE ezfio_filename
+  
+  call ezfio_has_mrcc_selected_thresh_dressed_ci(has)
+  if (has) then
+    call ezfio_get_mrcc_selected_thresh_dressed_ci(thresh_dressed_ci)
+  else
+    print *, 'mrcc_selected/thresh_dressed_ci not found in EZFIO file'
+    stop 1
+  endif
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, n_it_max_dressed_ci  ]
+  implicit none
+  BEGIN_DOC
+! Maximum number of dressed CI iterations
+  END_DOC
+
+  logical                        :: has
+  PROVIDE ezfio_filename
+  
+  call ezfio_has_mrcc_selected_n_it_max_dressed_ci(has)
+  if (has) then
+    call ezfio_get_mrcc_selected_n_it_max_dressed_ci(n_it_max_dressed_ci)
+  else
+    print *, 'mrcc_selected/n_it_max_dressed_ci not found in EZFIO file'
+    stop 1
+  endif
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, lambda_type  ]
+  implicit none
+  BEGIN_DOC
+! lambda type
+  END_DOC
+
+  logical                        :: has
+  PROVIDE ezfio_filename
+  
+  call ezfio_has_mrcc_selected_lambda_type(has)
+  if (has) then
+    call ezfio_get_mrcc_selected_lambda_type(lambda_type)
+  else
+    print *, 'mrcc_selected/lambda_type not found in EZFIO file'
+    stop 1
+  endif
+
+END_PROVIDER

plugins/mrcc_selected/mrcc_selected.irp.f

@@ -0,0 +1,19 @@
+program mrsc2sub
+  implicit none
+  double precision, allocatable  :: energy(:)
+  allocate (energy(N_states))
+  
+  !mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrcc
+  mrmode = 3
+  
+  read_wf = .True.
+  SOFT_TOUCH read_wf
+  call print_cas_coefs
+  call set_generators_bitmasks_as_holes_and_particles
+  call run(N_states,energy)
+  if(do_pt2_end)then
+    call run_pt2(N_states,energy)
+  endif
+  deallocate(energy)
+end
+

plugins/mrcc_selected/mrcepa0_general.irp.f

@@ -0,0 +1,245 @@
+
+
+subroutine run(N_st,energy)
+  implicit none
+  
+  integer, intent(in) :: N_st 
+  double precision, intent(out) :: energy(N_st) 
+
+  integer :: i,j
+
+  double precision :: E_new, E_old, delta_e
+  integer :: iteration
+  double precision :: E_past(4)
+  
+  integer :: n_it_mrcc_max
+  double precision :: thresh_mrcc
+  double precision, allocatable :: lambda(:)
+  allocate (lambda(N_states))
+
+  
+  thresh_mrcc = thresh_dressed_ci
+  n_it_mrcc_max = n_it_max_dressed_ci
+
+  if(n_it_mrcc_max == 1) then
+    do j=1,N_states_diag
+      do i=1,N_det
+        psi_coef(i,j) = CI_eigenvectors_dressed(i,j)
+      enddo
+    enddo
+    SOFT_TOUCH psi_coef ci_energy_dressed
+    call write_double(6,ci_energy_dressed(1),"Final MRCC energy")
+    call ezfio_set_mrcepa0_energy(ci_energy_dressed(1))
+    call save_wavefunction
+    energy(:) = ci_energy_dressed(:)
+  else
+    E_new = 0.d0
+    delta_E = 1.d0
+    iteration = 0
+    lambda = 1.d0
+    do while (delta_E > thresh_mrcc)
+      iteration += 1
+      print *,  '===========================' 
+      print *,  'MRCEPA0 Iteration', iteration
+      print *,  '===========================' 
+      print *,  ''
+      E_old = sum(ci_energy_dressed)
+      call write_double(6,ci_energy_dressed(1),"MRCEPA0 energy")
+      call diagonalize_ci_dressed(lambda)
+      E_new = sum(ci_energy_dressed)
+      delta_E = dabs(E_new - E_old)
+      call save_wavefunction
+      call ezfio_set_mrcepa0_energy(ci_energy_dressed(1))
+      if (iteration >= n_it_mrcc_max) then
+        exit
+      endif
+    enddo
+    call write_double(6,ci_energy_dressed(1),"Final MRCEPA0 energy")
+    energy(:) = ci_energy_dressed(:)
+  endif
+end
+
+
+subroutine print_cas_coefs
+  implicit none
+
+  integer :: i,j
+  print *,  'CAS'
+  print *,  '==='
+  do i=1,N_det_cas
+    print *,  (psi_cas_coef(i,j), j=1,N_states)
+    call debug_det(psi_cas(1,1,i),N_int)
+  enddo
+  call write_double(6,ci_energy(1),"Initial CI energy")
+
+end
+
+
+
+
+subroutine run_pt2_old(N_st,energy) 
+  implicit none 
+  integer :: i,j,k 
+  integer, intent(in)          :: N_st 
+  double precision, intent(in) :: energy(N_st) 
+  double precision :: pt2_redundant(N_st), pt2(N_st)
+  double precision :: norm_pert(N_st),H_pert_diag(N_st)
+  
+  pt2_redundant = 0.d0 
+  pt2 = 0d0
+  !if(lambda_mrcc_pt2(0)  == 0) return
+  
+  print*,'Last iteration only to compute the PT2' 
+  
+  print * ,'Computing the redundant PT2 contribution'
+
+  if (mrmode == 1) then
+
+    N_det_generators = lambda_mrcc_kept(0) 
+    N_det_selectors = lambda_mrcc_kept(0) 
+
+    do i=1,N_det_generators
+      j = lambda_mrcc_kept(i)
+      do k=1,N_int
+        psi_det_generators(k,1,i) = psi_non_ref(k,1,j)
+        psi_det_generators(k,2,i) = psi_non_ref(k,2,j)
+        psi_selectors(k,1,i) = psi_non_ref(k,1,j)
+        psi_selectors(k,2,i) = psi_non_ref(k,2,j)
+      enddo
+      do k=1,N_st
+        psi_coef_generators(i,k) = psi_non_ref_coef(j,k)
+        psi_selectors_coef(i,k) = psi_non_ref_coef(j,k)
+      enddo
+    enddo
+
+  else
+
+    N_det_generators = N_det_non_ref
+    N_det_selectors = N_det_non_ref
+
+    do i=1,N_det_generators
+      j = i
+      do k=1,N_int
+        psi_det_generators(k,1,i) = psi_non_ref(k,1,j)
+        psi_det_generators(k,2,i) = psi_non_ref(k,2,j)
+        psi_selectors(k,1,i) = psi_non_ref(k,1,j)
+        psi_selectors(k,2,i) = psi_non_ref(k,2,j)
+      enddo
+      do k=1,N_st
+        psi_coef_generators(i,k) = psi_non_ref_coef(j,k)
+        psi_selectors_coef(i,k) = psi_non_ref_coef(j,k)
+      enddo
+    enddo
+
+  endif
+  
+  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_redundant, norm_pert, H_pert_diag,  N_st) 
+  
+  print * ,'Computing the remaining contribution'
+
+  threshold_selectors = max(threshold_selectors,threshold_selectors_pt2)
+  threshold_generators = max(threshold_generators,threshold_generators_pt2)
+
+  N_det_generators = N_det_non_ref + N_det_ref
+  N_det_selectors = N_det_non_ref + N_det_ref
+  
+  psi_det_generators(:,:,:N_det_ref) = psi_ref(:,:,:N_det_ref)
+  psi_selectors(:,:,:N_det_ref) = psi_ref(:,:,:N_det_ref)
+  psi_coef_generators(:N_det_ref,:) = psi_ref_coef(:N_det_ref,:)
+  psi_selectors_coef(:N_det_ref,:) = psi_ref_coef(:N_det_ref,:)
+  
+  do i=N_det_ref+1,N_det_generators
+    j = i-N_det_ref
+    do k=1,N_int
+      psi_det_generators(k,1,i) = psi_non_ref(k,1,j)
+      psi_det_generators(k,2,i) = psi_non_ref(k,2,j)
+      psi_selectors(k,1,i) = psi_non_ref(k,1,j)
+      psi_selectors(k,2,i) = psi_non_ref(k,2,j)
+    enddo
+    do k=1,N_st
+      psi_coef_generators(i,k) = psi_non_ref_coef(j,k)
+      psi_selectors_coef(i,k) = psi_non_ref_coef(j,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) 
+ 
+
+  print *, "Redundant PT2 :",pt2_redundant
+  print *, "Full      PT2 :",pt2
+  print *, lambda_mrcc_kept(0), N_det, N_det_ref, psi_coef(1,1), psi_ref_coef(1,1)
+  pt2 = pt2 - pt2_redundant
+  
+  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_mrcepa0_energy_pt2(energy(1)+pt2(1))
+
+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_mrcepa0_energy_pt2(energy(1)+pt2(1))
+
+end 
+

src/Determinants/determinants.irp.f

@@ -438,8 +438,12 @@ end
  do i=1,N_states
    psi_coef_min(i) = minval(psi_coef(:,i))
    psi_coef_max(i) = maxval(psi_coef(:,i))
-   abs_psi_coef_min(i) = dabs(psi_coef_min(i))
-   abs_psi_coef_max(i) = dabs(psi_coef_max(i))
+   abs_psi_coef_min(i) = minval( dabs(psi_coef(:,i)) )
+   abs_psi_coef_max(i) = maxval( dabs(psi_coef(:,i)) )
+   call write_double(6,psi_coef_max(i), 'Max coef')
+   call write_double(6,psi_coef_min(i), 'Min coef')
+   call write_double(6,abs_psi_coef_max(i), 'Max abs coef')
+   call write_double(6,abs_psi_coef_min(i), 'Min abs coef')
  enddo
 
 END_PROVIDER