diff --git a/devel/cc/denominators.irp.f b/devel/cc/denominators.irp.f
index d605086..3af8e4e 100644
--- a/devel/cc/denominators.irp.f
+++ b/devel/cc/denominators.irp.f
@@ -21,6 +21,7 @@ BEGIN_PROVIDER [ double precision, delta_OV, (spin_occ_num, spin_vir_num) ]
    do i=1,spin_occ_num
      delta_OV(i,a) = eO(i) - eV(a)
 !     delta_OV(i,a) = min(eO(i) - eV(a), -1.d-3)
+     if (delta_OV(i,a) == 0.d0) delta_OV(i,a) = 1.d-1
    enddo
  enddo
 
@@ -39,6 +40,7 @@ BEGIN_PROVIDER [ double precision, delta_OOVV, (spin_occ_num,spin_occ_num,spin_v
    do j=1,spin_occ_num
     do i=1,spin_occ_num
      delta_OOVV(i,j,a,b) = delta_OV(i,a) + delta_OV(j,b)
+     if (delta_OOVV(i,j,a,b) == 0.d0) delta_OOVV(i,j,a,b) = 1.d-1
     enddo
    enddo
   enddo
@@ -60,6 +62,7 @@ BEGIN_PROVIDER [ double precision, delta_OOOVVV, (spin_occ_num,spin_occ_num,spin
      do j=1,spin_occ_num
       do i=1,spin_occ_num
         delta_OOOVVV(i,j,k,a,b,c) = delta_OV(i,a) + delta_OV(j,b) + delta_OV(k,c)
+        if (delta_OOOVVV(i,j,k,a,b,c) == 0.d0) delta_OOOVVV(i,j,k,a,b,c) = 1.d-1
       enddo
      enddo
     enddo
diff --git a/devel/dmc_dress/EZFIO.cfg b/devel/dmc_dress/EZFIO.cfg
index 88bc4ff..af763fe 100644
--- a/devel/dmc_dress/EZFIO.cfg
+++ b/devel/dmc_dress/EZFIO.cfg
@@ -4,3 +4,9 @@ doc: Dressing matrix obtained from DMC
 size: (determinants.n_det)
 interface: ezfio, provider
 
+[dmc_delta_htc]
+type: double precision
+doc: Dressing matrix obtained from H_TC
+size: (determinants.n_det)
+interface: ezfio
+
diff --git a/devel/dmc_dress/dress.py b/devel/dmc_dress/dress.py
index 53f585b..4a5464a 100755
--- a/devel/dmc_dress/dress.py
+++ b/devel/dmc_dress/dress.py
@@ -12,21 +12,26 @@ sys.path.insert(0,QP_PATH)
 from ezfio import ezfio
 
 def read_hamiltonian(inp):
-    text = subprocess.run(["qmcchem", "result", inp], capture_output=True).stdout
+#    text = subprocess.run(["qmcchem", "result", inp], capture_output=True).stdout
+    text = subprocess.run(["cat", "OUT"], capture_output=True).stdout
     inside = None
     h = []
     s = []
     norm = None
     for line in text.splitlines():
-       line = str(line)
-       print (line)
+       line = line.decode("utf-8")
+       if line == "":
+         continue
        if "Psi_norm :" in line:
-           norm = float(line.split()[3])
+           norm = float(line.split()[2])
        if "]" in line:
            inside = None
        if inside == "H":
           data = line.split()
-          h.append(float(data[3]))
+#          h.append(float(data[2]))                            #ave
+          h.append(max(float(data[2]) - float(data[4]),0.))   #diff
+#          if float(data[2]) - float(data[4])>0.: h.append(float(data[2]))  
+#          else: h.append(0.)
        elif "Ci_dress" in line:
            inside = "H"
     h = np.array(h)/norm
@@ -38,3 +43,4 @@ h = read_hamiltonian(sys.argv[1])
 ezfio.set_file(sys.argv[1])
 ezfio.set_dmc_dress_dmc_delta_h(h)
 print(h)
+
diff --git a/devel/dmc_dress/dressing_vector.irp.f b/devel/dmc_dress/dressing_vector.irp.f
index 71c59c8..3498b90 100644
--- a/devel/dmc_dress/dressing_vector.irp.f
+++ b/devel/dmc_dress/dressing_vector.irp.f
@@ -11,7 +11,7 @@
  double precision, allocatable :: delta(:)
 
  allocate(delta(N_det))
- delta(1:N_det) = dmc_delta_h(1:N_det)
+ delta (1:N_det) = dmc_delta_h(1:N_det) + dmc_delta_htc(1:N_det) 
 
  call dset_order(delta,psi_bilinear_matrix_order_reverse,N_det)
 
@@ -31,3 +31,38 @@ END_PROVIDER
 
 
 
+BEGIN_PROVIDER [ double precision, dmc_delta_htc , (n_det) ]
+  implicit none
+  BEGIN_DOC
+! Dressing matrix obtained from H_TC
+  END_DOC
+
+  logical                        :: has
+  PROVIDE ezfio_filename
+  if (mpi_master) then
+    if (size(dmc_delta_htc) == 0) return
+
+    call ezfio_has_dmc_dress_dmc_delta_htc(has)
+    if (has) then
+      write(6,'(A)') '.. >>>>> [ IO READ: dmc_delta_htc ] <<<<< ..'
+      call ezfio_get_dmc_dress_dmc_delta_htc(dmc_delta_htc)
+    else
+      dmc_delta_htc(:) = 0.d0
+    endif
+  endif
+  IRP_IF MPI_DEBUG
+    print *,  irp_here, mpi_rank
+    call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+  IRP_ENDIF
+  IRP_IF MPI
+    include 'mpif.h'
+    integer :: ierr
+    call MPI_BCAST( dmc_delta_htc, (n_det), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+    if (ierr /= MPI_SUCCESS) then
+      stop 'Unable to read dmc_delta_htc with MPI'
+    endif
+  IRP_ENDIF
+
+  call write_time(6)
+
+END_PROVIDER
diff --git a/devel/trexio/import_trexio_integrals.irp.f b/devel/trexio/import_trexio_integrals.irp.f
index 4d8e3da..22ac61e 100644
--- a/devel/trexio/import_trexio_integrals.irp.f
+++ b/devel/trexio/import_trexio_integrals.irp.f
@@ -111,7 +111,6 @@ subroutine run
 
   integer*8 :: offset, icount
 
-!  open(unit=104,file='tmp')
   offset = 0_8
   icount = BUFSIZE
   rc = TREXIO_SUCCESS
@@ -125,7 +124,6 @@ subroutine run
       integral = V(m)
       call two_e_integrals_index(i, j, k, l, buffer_i(m) )
       buffer_values(m) = integral
-!  write(104,'(4(I5,X),2D22.15)') i,j,k,l, integral
     enddo
     call insert_into_ao_integrals_map(int(icount,4),buffer_i,buffer_values)
     offset = offset + icount
@@ -134,7 +132,6 @@ subroutine run
     endif
   end do
   n_integrals = offset
-  close(104)
 
   call map_sort(ao_integrals_map)
   call map_unique(ao_integrals_map)
diff --git a/fnmf/EZFIO.cfg b/fnmf/EZFIO.cfg
new file mode 100644
index 0000000..033d0f1
--- /dev/null
+++ b/fnmf/EZFIO.cfg
@@ -0,0 +1,17 @@
+[E_dmc]
+type: double precision
+doc: DMC energy
+interface: ezfio, provider
+
+[h_dmc]
+type: double precision
+doc: Dressing matrix obtained from DMC
+size: (determinants.n_det)
+interface: ezfio, provider
+
+[s_dmc]
+type: double precision
+doc: Dressing matrix obtained from H_TC
+size: (determinants.n_det)
+interface: ezfio, provider
+
diff --git a/fnmf/NEED b/fnmf/NEED
new file mode 100644
index 0000000..8d89a45
--- /dev/null
+++ b/fnmf/NEED
@@ -0,0 +1,2 @@
+determinants
+davidson_undressed
diff --git a/fnmf/README.rst b/fnmf/README.rst
new file mode 100644
index 0000000..71cd4ad
--- /dev/null
+++ b/fnmf/README.rst
@@ -0,0 +1,4 @@
+====
+fnmf
+====
+
diff --git a/fnmf/dmc_data.irp.f b/fnmf/dmc_data.irp.f
new file mode 100644
index 0000000..84257ce
--- /dev/null
+++ b/fnmf/dmc_data.irp.f
@@ -0,0 +1,70 @@
+ BEGIN_PROVIDER [ double precision, h_dmc_row , (N_det) ]
+&BEGIN_PROVIDER [ double precision, s_dmc_row , (N_det) ]
+ implicit none
+ BEGIN_DOC
+ ! Data sampled with QMC=Chem
+ END_DOC
+
+ h_dmc_row(:) = h_dmc(:)
+ s_dmc_row(:) = s_dmc(:)
+ call dset_order(h_dmc_row,psi_bilinear_matrix_order_reverse,N_det)
+ call dset_order(s_dmc_row,psi_bilinear_matrix_order_reverse,N_det)
+
+!  integer :: i
+!  do i=1,N_det
+!    s_dmc_row(i) = psi_coef(i,1)
+!    call i_h_psi(psi_det(1,1,i), psi_det, psi_coef, N_int, N_det, &
+!         N_det, 1, h_dmc_row(i) )
+!  enddo
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, mat_size ]
+ implicit none
+ BEGIN_DOC
+! Size of the matrices
+ END_DOC
+ mat_size = N_det+1
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, H_dmc_mat, (mat_size, mat_size) ]
+ implicit none
+ BEGIN_DOC
+! Hamiltonian extended with DMC data
+ END_DOC
+ integer :: i,j
+ do j=1,N_det
+   do i=1,N_det
+      call i_h_j(psi_det(1,1,i), psi_det(1,1,j), N_int, H_dmc_mat(i,j))
+   enddo
+ enddo
+
+ do i=1,N_det
+   call i_h_psi(psi_det(1,1,i), psi_det, psi_coef, N_int, N_det, &
+        N_det, 1, H_dmc_mat(i,N_det+1) )
+   H_dmc_mat(N_det+1,i) = h_dmc_row(i)
+ enddo
+ H_dmc_mat(mat_size,mat_size) = E_dmc - nuclear_repulsion
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, S_dmc_mat, (mat_size, mat_size) ]
+ implicit none
+ BEGIN_DOC
+! Overlap matrix extended with DMC data
+ END_DOC
+ integer :: i,j
+ S_dmc_mat = 0.d0
+ do i=1,mat_size
+   S_dmc_mat(i,i) = 1.d0
+ enddo
+
+ do i=1,N_det
+   S_dmc_mat(i,N_det+1) = psi_coef(i,1)
+   S_dmc_mat(N_det+1,i) = S_dmc_row(i)
+ enddo
+
+END_PROVIDER
+
+
+
diff --git a/fnmf/fnmf.irp.f b/fnmf/fnmf.irp.f
new file mode 100644
index 0000000..fff7a6a
--- /dev/null
+++ b/fnmf/fnmf.irp.f
@@ -0,0 +1,76 @@
+program fnmf
+  implicit none
+  BEGIN_DOC
+! TODO : Put the documentation of the program here
+  END_DOC
+  read_wf = .True.
+  TOUCH read_wf
+  call run
+end
+
+
+subroutine run
+ implicit none
+ integer :: i, n_real
+ double precision, allocatable :: beta(:), vr(:,:), vl(:,:)
+ double precision, allocatable :: htmp(:,:), stmp(:,:)
+
+ allocate(Htmp(N_det,N_det), Stmp(N_det,N_det))
+
+ htmp(:, :) = H_dmc_mat(1:N_det, 1:N_det)
+ stmp(:, :) = S_dmc_mat(1:N_det, 1:N_det)
+
+ htmp(1, 1) = H_dmc_mat(N_det+1, N_det+1)
+ stmp(1, 1) = S_dmc_mat(N_det+1, N_det+1)
+
+ htmp(2:N_det, 1) = H_dmc_mat(2:N_det, N_det+1)
+ stmp(2:N_det, 1) = S_dmc_mat(2:N_det, N_det+1)
+
+ htmp(1, 2:N_det) = H_dmc_mat(N_det+1, 2:N_det)
+ stmp(1, 2:N_det) = S_dmc_mat(N_det+1, 2:N_det)
+
+
+ allocate(beta(N_det), vr(N_det,N_det), vl(N_det, N_det))
+
+ call lapack_g_non_sym_real(N_det, Htmp, size(Htmp,1), &
+  Stmp, size(Stmp,1), beta, n_real, vl, size(vl,1), vr, size(vr,1))
+
+
+  vl(:,1) = psi_coef(:,1)
+
+  psi_coef(:,1) *= vr(1,1)
+  do i=2,N_det
+    psi_coef(i,1) += vr(i,1)
+  end do
+  double precision, external :: dnrm2, ddot
+  double precision :: norm
+
+  norm = dnrm2(N_det, psi_coef, 1)
+  psi_coef(:,1) *= 1.d0/norm
+  SOFT_TOUCH psi_coef
+
+  print *, '-------------------------------------------------'
+  print *, ' psi_old  psi_new'
+  print *, '-------------------------------------------------'
+  do i=1,N_det
+    print '(2(E16.8,X))', vl(i,1), psi_coef(i,1) 
+  enddo
+  print *, '-------------------------------------------------'
+
+  print *,  ''
+  print *, 'Overlap        : ', ddot(N_det, psi_coef, 1, vl, 1)
+  print *, 'DMC energy     : ', htmp(1,1) + nuclear_repulsion
+  print *, 'Updated energy : ', beta(1) + nuclear_repulsion
+
+  print *,  'H'
+  do i=1,N_det
+    print *,  htmp(i, 1:N_det)
+  enddo
+  print *,  'S'
+  do i=1,N_det
+    print *,  stmp(i, 1:N_det)
+  enddo
+  call save_wavefunction
+ 
+end
+
diff --git a/fnmf/non_hermit.irp.f b/fnmf/non_hermit.irp.f
new file mode 100644
index 0000000..058754c
--- /dev/null
+++ b/fnmf/non_hermit.irp.f
@@ -0,0 +1,74 @@
+subroutine lapack_g_non_sym_real(n, H, LDH, S, LDS, beta, &
+      n_real, vl, LDVL, vr, LDVR)
+  implicit none
+  integer, intent(in) :: n, LDH, LDS, LDVL, LDVR
+  double precision, intent(in) :: H(LDH,n), S(LDS,n)
+  double precision, intent(out) :: VL(LDVL,n), VR(LDVR,n), beta(n)
+  integer, intent(out) :: n_real
+
+  integer :: lwork, info, i,j 
+  double precision, allocatable :: work(:), Htmp(:,:), Stmp(:,:)
+  double precision, allocatable :: alphar(:), alphai(:), vltmp(:,:), vrtmp(:,:)
+  integer, allocatable :: iorder(:), list_good(:)
+
+
+  lwork = -1
+  allocate(work(1), alphar(n), alphai(n), vltmp(n,n), vrtmp(n,n), &
+  Htmp(n,n), Stmp(n,n), list_good(n))
+
+  Htmp(1:n,1:n) = H(1:n,1:n)
+  Stmp(1:n,1:n) = S(1:n,1:n)
+  call dggev('V', 'V', n, Htmp, size(Htmp,1), Stmp, size(Stmp,1), alphar, alphai, beta, &
+             vltmp, size(vltmp,1), vrtmp, size(vrtmp,1), work, lwork, info)
+
+  lwork = int(work(1))
+  deallocate(work)
+  allocate(work(lwork))
+  call dggev('V', 'V', n, Htmp, size(Htmp,1), Stmp, size(Stmp,1), alphar, alphai, beta, &
+             vltmp, size(vltmp,1), vrtmp, size(vrtmp,1), work, lwork, info)
+
+  deallocate(work)
+  if (info /= 0) then
+    stop 'DGGEV Diagonalization failed'
+  endif
+
+  allocate(iorder(n))
+  n_real = 0
+  do i=1,n
+    if (dabs(alphai(i)) < 1.d-10) then 
+       n_real += 1
+       list_good(n_real) = i
+    else
+       alphar(i) = dble(huge(1.0))
+    endif
+  enddo
+
+  do i=1,n
+    if (dabs(beta(i)/(alphar(i)+1.d-12)) < 1.d-10) then
+      beta(i) = 0.d0
+    else
+      beta(i) = alphar(i)/beta(i)
+    endif
+  enddo
+
+  do i=1,n_real
+    iorder(i) = i
+    do j=1,n
+      vrtmp(j,i) = vrtmp(j,list_good(i))
+      vltmp(j,i) = vltmp(j,list_good(i))
+    end do
+  end do
+
+  call dsort(beta, iorder, n_real)
+
+  do i=1,n_real
+    do j=1,n
+      vr(j,i) = vrtmp(j,iorder(i))
+      vl(j,i) = vltmp(j,iorder(i))
+    end do
+  end do
+
+  deallocate(vrtmp, vltmp, iorder, Htmp, Stmp)
+
+end
+