anchor.diff:    
     integer :: II, blok
     
  +  call wall_time(wall)
  +  print *, "cepa0", wall
     provide det_cepa0_active delta_cas lambda_mrcc
     provide mo_bielec_integrals_in_map
     allocate(idx_sorted_bit(N_det))
  
positive.code:   enddo
    
    
    
    
    do i_state = 1, N_states 
  
      delta_ii_old(:,:) = 0 
      delta_ij_old(:,:,:) = 0 
      
      !$OMP PARALLEL DO default(none) schedule(dynamic,10) shared(delta_ij_old, delta_ii_old)       &
      !$OMP private(i, J, k, degree, degree2, l, deg, ni, inac, virt)       &
      !$OMP private(ok,p1,p2,h1,h2,s1,s2, p1_,p2_,h1_,h2_,s1_,s2_)     &
      !$OMP private(phase_iI, phase_Ik, phase_Jl, phase_IJ, phase_al, diI, hIi, hJi, delta_JI, dkI, HkI, ci_inv, dia_hla) &
      !$OMP private(contrib, exc_iI, exc_Ik, exc_IJ, det_tmp, det_tmp2)    &
      !$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)
      do i = 1 , N_det_non_ref
        if(mod(i,1000) == 0) print *, i, "/", N_det_non_ref
        if(lambda_mrcc(i_state, i) == 0d0) cycle
        do i_I = 1 , N_det_ref
          call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,i),exc_iI,degree2,phase_iI,N_int)
          if(degree2 == -1) cycle
          ci_inv(i_state) = 1.d0 / psi_ref_coef(i_I,i_state)
  
          call decode_exc(exc_iI,degree2,h1,p1,h2,p2,s1,s2)
        
          call i_h_j(psi_non_ref(1,1,i), psi_ref(1,1,i_I),N_int,hIi)
          diI = hIi * lambda_mrcc(i_state, i)
          do J = 1 , i_I ! N_det_ref !!! 
            call get_excitation(psi_ref(1,1,i_I),psi_ref(1,1,J),exc_IJ,degree,phase_IJ,N_int)
            call i_h_j(psi_non_ref(1,1,i), psi_ref(1,1,J),N_int,hJi)
            if(hJi == 0) cycle
            delta_JI = hJi * diI 
            do k = 1 , N_det_non_ref
              if(lambda_mrcc(i_state, k) == 0d0) cycle
              
              call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,k),exc_Ik,degree,phase_Ik,N_int)
              if(degree == -1) cycle
  
              call decode_exc(exc_Ik,degree,h1_,p1_,h2_,p2_,s1_,s2_)
  
              
              det_tmp(:,:) = 0_bit_kind
              det_tmp2(:,:) = 0_bit_kind
              
              
        
              call apply_excitation(psi_non_ref(1,1,i),exc_Ik,det_tmp,ok,N_int)
              
              
              if(.not. ok) cycle
              !if(is_in_wavefunction(det_tmp, N_int)) cycle
              inac = iand(HF_bitmask(1,1), not(active_sorb(1,1)))
              virt = iand(not(HF_bitmask(1,1)), not(active_sorb(1,1)))
  
              deg = 0
              deg += popcnt(xor(iand(inac,det_tmp(1,1)), inac))
              deg += popcnt(xor(iand(inac,det_tmp(1,2)), inac))
              if(deg <= 2) then
                deg = 0
                deg += popcnt(iand(virt, det_tmp(1,1)))
                deg += popcnt(iand(virt, det_tmp(1,2)))
                if(deg <= 2) then
                  cycle
                end if
              end if
  
              
              !call get_excitation(psi_non_ref(1,1,i), det_tmp, exc_Ik, degree, phase_al, N_int)
  
        
              !if(is_in_wavefunction(det_tmp, N_int)) cycle
        
  
              call apply_excitation(psi_ref(1,1,J),exc_Ik,det_tmp,ok,N_int)
              if(.not. ok) cycle
        
              call get_excitation(psi_ref(1,1,J), det_tmp, exc_Ik, degree, phase_Jl, N_int)
        
              l = get_index_in_psi_det_sorted_bit(det_tmp, N_int)
              if(l == 0) cycle
              l = idx_sorted_bit(l)
              if(l ==-1) cycle
              
              
  
              call i_h_j(psi_non_ref(1,1,k), psi_ref(1,1,i_I),N_int,HkI)
              dkI(i_state) = HkI * lambda_mrcc(i_state, k)! * phase_Jl * phase_Ik
              !if( phase_Jl * phase_Ik < 0d0 ) stop "STOOOOOOP"
              contrib = dkI(i_state) * delta_JI
              !$OMP ATOMIC  
              delta_ij_old(i_I,l,i_state) += contrib
              if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5) then
                !$OMP ATOMIC  
                delta_ii_old(i_I,i_state) -= contrib * ci_inv(i_state) * psi_non_ref_coef(k,i_state)
              endif
  
            enddo
  
negative.code:     !$OMP END PARALLEL DO
    end do
    deallocate(idx_sorted_bit)
    call wall_time(wall)
    print *, "cepa0", wall
    stop
  END_PROVIDER
  
  
  