diff --git a/devel/dmc_dress/dmc_dress.irp.f b/devel/dmc_dress/dmc_dress.irp.f
index f79b98b..72a78b4 100644
--- a/devel/dmc_dress/dmc_dress.irp.f
+++ b/devel/dmc_dress/dmc_dress.irp.f
@@ -1,4 +1,4 @@
-program diagonalize_h
+program dmc_dress
  implicit none
  BEGIN_DOC
 ! Program that extracts the lowest states of the Hamiltonian dressed by the QMC
@@ -7,13 +7,30 @@ program diagonalize_h
  END_DOC
  read_wf = .True.
  touch read_wf
+ call pre
  call routine
  call save_wavefunction_general(N_det,N_states,psi_det_sorted,size(psi_coef_sorted,1),psi_coef_sorted)
 end
 
+subroutine pre
+ implicit none
+ double precision, allocatable :: left(:,:), right(:,:), tmp(:,:), res(:,:)
+ integer :: i
+ allocate (left(1,1:N_det), right(1:N_det,1), tmp(1:N_det,1), res(1,1))
+ left(1,1:N_det) = psi_coef(1:N_det,1)
+ right(1:N_det,1) = psi_coef(1:N_det,1)
+
+ tmp(1:N_det,1:1) = matmul(h_matrix_dressed(1:N_det,1:N_det), right(1:N_det,1:1))
+ res(1:1,1:1) = matmul(left(1:1,1:N_det), tmp(1:N_det,1:1))
+ print *,  'E_in = ', res(1,1)
+ do i=1,N_det
+ print *,  'HPsi/c0 = ', tmp(i,1)/psi_coef(i,1)  
+ enddo
+end
 subroutine routine
  implicit none
  psi_coef(1:N_det,1) = ci_eigenvectors_dressed(1:N_det,1)
  print*,'N_det = ',N_det
+ print *, 'E = ', ci_energy_dressed(1) + nuclear_repulsion
  SOFT_TOUCH psi_coef
 end
diff --git a/devel/dmc_dress/dress.py b/devel/dmc_dress/dress.py
new file mode 100755
index 0000000..53f585b
--- /dev/null
+++ b/devel/dmc_dress/dress.py
@@ -0,0 +1,40 @@
+#!/usr/bin/env python
+
+import numpy as np
+import subprocess
+import sys
+import os
+
+QP_PATH=os.environ["QP_ROOT"]+"/external/ezfio/Python/"
+
+sys.path.insert(0,QP_PATH)
+
+from ezfio import ezfio
+
+def read_hamiltonian(inp):
+    text = subprocess.run(["qmcchem", "result", inp], capture_output=True).stdout
+    inside = None
+    h = []
+    s = []
+    norm = None
+    for line in text.splitlines():
+       line = str(line)
+       print (line)
+       if "Psi_norm :" in line:
+           norm = float(line.split()[3])
+       if "]" in line:
+           inside = None
+       if inside == "H":
+          data = line.split()
+          h.append(float(data[3]))
+       elif "Ci_dress" in line:
+           inside = "H"
+    h = np.array(h)/norm
+
+    return h
+
+
+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 9f79996..71c59c8 100644
--- a/devel/dmc_dress/dressing_vector.irp.f
+++ b/devel/dmc_dress/dressing_vector.irp.f
@@ -3,21 +3,30 @@
  implicit none
  BEGIN_DOC
  ! \Delta_{state-specific}. \Psi
+ ! Diagonal element is divided by 2 because Delta = D + D^t
  END_DOC
 
  integer :: i,ii,k,j, l
  double precision :: f, tmp
- double precision, external :: u_dot_v
- logical, external :: detEq
+ double precision, allocatable :: delta(:)
+
+ allocate(delta(N_det))
+ delta(1:N_det) = dmc_delta_h(1:N_det)
+
+ call dset_order(delta,psi_bilinear_matrix_order_reverse,N_det)
 
  dressing_column_h(:,:) = 0.d0
  dressing_column_s(:,:) = 0.d0
 
  l = dressed_column_idx(1)
  do j = 1, n_det
-   dressing_column_h(j,1) = 0.5d0*dmc_delta_h(j)
-   dressing_column_h(l,1) -= 0.5d0 * psi_coef(j,1) * dmc_delta_h(j) /psi_coef(l,1)
+   if (j == l) cycle
+   dressing_column_h(j,1)  = delta(j) 
+   dressing_column_h(l,1) -= psi_coef(j,1) * delta(j) / psi_coef(l,1)
  enddo
+ dressing_column_h(l,1) += delta(l) 
+ dressing_column_h(l,1) *= 0.5d0
+
 END_PROVIDER