diff --git a/bin/jast_opt.py b/bin/jast_opt.py
index 3816be3..1725ddf 100755
--- a/bin/jast_opt.py
+++ b/bin/jast_opt.py
@@ -29,6 +29,8 @@ def main():
     filename = sys.argv[1]
     ezfio.set_file(filename)
 
+    jast_type = ezfio.jastrow_jast_type
+    print (jast_type)
 
     def make_atom_map():
         labels = {}
@@ -47,18 +49,62 @@ def main():
     atom_map = make_atom_map()
 
 
-    def get_params_b():
-        return ezfio.jastrow_jast_b_up_dn
+    if jast_type == 'Simple':
 
-    def get_params_pen():
-        d = ezfio.jastrow_jast_pen
-        print(atom_map)
-        for m in atom_map:
-          print(m[0])
-          print (d[m[0]])
-        sys.stdout.flush()
-        return np.array([d[m[0]] for m in atom_map])
+      def get_params_b():
+          return ezfio.jastrow_jast_b_up_dn
 
+      def set_params_b(x):
+          x = np.abs(x)
+          ezfio.set_jastrow_jast_b_up_up(x)
+          ezfio.set_jastrow_jast_b_up_dn(x)
+
+    elif jast_type == 'Mu':
+
+      def get_params_b():
+          return ezfio.jastrow_mu_erf
+
+      def set_params_b(x):
+          x = np.abs(x)
+          ezfio.set_jastrow_mu_erf(x)
+
+    if jast_type == 'Simple':
+
+      def get_params_pen():
+          d = ezfio.jastrow_jast_pen
+          return np.array([d[m[0]] for m in atom_map])
+
+      def set_params_pen(x):
+          x = np.abs(x)
+          y=list(ezfio.jastrow_jast_pen)
+          for i,m in enumerate(atom_map):
+              for j in m:
+                  y[j] = x[i]
+          ezfio.set_jastrow_jast_pen(y)
+
+    elif jast_type == 'Mu':
+
+      def get_params_pen():
+          d = ezfio.jastrow_jast_1bgauss_pen
+          return np.array([d[m[0]] for m in atom_map])
+
+      def set_params_pen(x):
+          x = np.abs(x)
+          y=list(ezfio.jastrow_jast_1bgauss_pen)
+          for i,m in enumerate(atom_map):
+              for j in m:
+                  y[j] = x[i]
+          ezfio.set_jastrow_jast_1bgauss_pen(y)
+
+    def get_norm():
+        return 1.,0.
+#        buffer = subprocess.check_output(['qmcchem', 'result', '-e', 'psi_norm', filename],
+#                                         encoding='UTF-8')
+#        if buffer.strip() != "":
+#            buffer = buffer.splitlines()[-1]
+#            _, energy, error = [float(x) for x in buffer.split()]
+#        else:
+#            return None, None
 
     def get_energy():
         buffer = subprocess.check_output(['qmcchem', 'result', '-e', 'e_loc', filename],
@@ -82,19 +128,6 @@ def main():
             return None, None
 
 
-    def set_params_b(x):
-        x = np.abs(x)
-        ezfio.set_jastrow_jast_b_up_up(x)
-        ezfio.set_jastrow_jast_b_up_dn(x)
-
-    def set_params_pen(x):
-        x = np.abs(x)
-        y=list(ezfio.jastrow_jast_pen)
-        for i,m in enumerate(atom_map):
-            for j in m:
-                y[j] = x[i]
-        ezfio.set_jastrow_jast_pen(y)
-
 
     def run_qmc():
         return subprocess.check_output(['qmcchem', 'run', filename])
@@ -105,7 +138,7 @@ def main():
 
 
     def set_vmc_params():
-#        subprocess.check_output(['qmcchem', 'edit', '-c', '-j', 'Simple',
+#        subprocess.check_output(['qmcchem', 'edit', '-c', '-j', jast_type,
 #                                 '-m', 'VMC',
 #                                 '-l', str(block_time),
 #                                 '--time-step=0.3',
@@ -113,7 +146,7 @@ def main():
 #                                 '--norm=1.e-5',
 #                                 '-w', '10',
 #                                 filename])
-        subprocess.check_output(['qmcchem', 'edit', '-c', '-j', 'Simple',
+        subprocess.check_output(['qmcchem', 'edit', '-c', '-j', jast_type,
                                  '-l', str(block_time),
                                  filename])
 
@@ -142,9 +175,12 @@ def main():
             while err > local_thresh:
                 time.sleep(block_time)
                 e, e_err = get_energy()
-                variance, v_err = get_variance()
+                variance, v_err = get_variance() 
                 if e is None or variance is None:
                     continue
+#                norm, _ = get_norm() 
+#                e, e_err = e, e_err/norm
+#                variance, v_err =variance/norm, v_err/norm
                 energy = e #+ variance
                 err = e_err #sqrt(e_err*e_err+v_err*v_err)
                 print(" %f  %f    %f  %f  %f  %f"%(e, e_err, variance, v_err, energy, err))
@@ -173,7 +209,9 @@ def main():
         if sum(x) == 0.:
           jast_a_up_dn = ezfio.jastrow_jast_a_up_dn
           x += jast_a_up_dn
-        opt = sp.optimize.minimize(f,x,method="Powell",
+#        opt = sp.optimize.minimize(f,x,method="Powell",
+#                             options= {'disp':True, 'ftol':thresh,'xtol':0.02})
+        opt = sp.optimize.minimize(f,x,method="Nelder-Mead",
                              options= {'disp':True, 'ftol':thresh,'xtol':0.02})
         print("x = "+str(opt))
         sys.stdout.flush()
diff --git a/src/JASTROW/jastrow_1b.irp.f b/src/JASTROW/jastrow_1b.irp.f
index 6c74b82..9acbc84 100644
--- a/src/JASTROW/jastrow_1b.irp.f
+++ b/src/JASTROW/jastrow_1b.irp.f
@@ -49,8 +49,7 @@ BEGIN_PROVIDER [ double precision, jast_1b_value, (elec_num_8)  ]
 
     elseif( jast_1b_type .eq. 4 ) then ! add 1body-RSDFT  Jastrow
                                        ! J(i) = - \sum_A [ -z_A r_iA erfc(mu*r_iA) + z_A exp(-(mu*r_iA)^2)/(mu*sqt_pi) ]
-!      mu    = jast_mu_erf
-      mu    = mu_erf
+      mu    = jast_mu_erf
       mu_pi = 1.d0 / ( dsqpi * mu )
       !DIR$ LOOP COUNT (100)
       do j = 1, nucl_num
@@ -151,8 +150,7 @@ END_PROVIDER
 
     elseif( jast_1b_type .eq. 4 ) then ! add 1body-RSDFT  Jastrow
                                        ! J(i) = - \sum_A [ -z_A r_iA erfc(mu*r_iA) + z_A exp(-(mu*r_iA)^2)/(mu*sqt_pi) ]
-!      mu = jast_mu_erf
-      mu = mu_erf
+      mu = jast_mu_erf
       !DIR$ LOOP COUNT (100)
       do j = 1, nucl_num
         rij = nucl_elec_dist(j,i)
@@ -248,8 +246,7 @@ BEGIN_PROVIDER [ double precision, jast_1b_lapl, (elec_num_8) ]
 
     elseif( jast_1b_type .eq. 4 ) then ! add 1body-RSDFT  Jastrow
                                        ! J(i) = - \sum_A [ -z_A r_iA erfc(mu*r_iA) + z_A exp(-(mu*r_iA)^2)/(mu*sqt_pi) ]
-!      mu    = jast_mu_erf
-      mu    = mu_erf
+      mu    = jast_mu_erf
       mu_pi = mu / dsqpi
       !DIR$ LOOP COUNT (100)
       do j = 1, nucl_num
@@ -278,7 +275,7 @@ BEGIN_PROVIDER [ double precision, jast_1b_lapl, (elec_num_8) ]
         a   = jast_1bGauss_pen(j)
         rij = nucl_elec_dist(j,i)
         c   = a * rij * rij
-        tmp = 2.d0 * a * dexp(-c) * (3.d0-2.d0*a*c)
+        tmp = 2.d0 * a * dexp(-c) * (3.d0-2.d0*c)
         jast_1b_lapl(i) -= tmp
       enddo
 
diff --git a/src/JASTROW/jastrow_mu.irp.f b/src/JASTROW/jastrow_mu.irp.f
index b88f81a..56efefd 100644
--- a/src/JASTROW/jastrow_mu.irp.f
+++ b/src/JASTROW/jastrow_mu.irp.f
@@ -13,10 +13,10 @@ implicit none
  double precision :: a, b, rij, tmp
  include '../constants.F'
  double precision :: mu
- mu = mu_erf
+ mu = jast_mu_erf
 
  do i=1,elec_num
-  jast_elec_Mu_value(i) = 0.d0
+   jast_elec_Mu_value(i) =  jast_1b_value(i)
  enddo
 
  do j=1,elec_num 
@@ -25,10 +25,9 @@ implicit none
     if(j==i)cycle
     rij = elec_dist(i,j)
     tmp = 0.5d0 * rij * (1.d0 - derf(mu*rij)) - 0.5d0/(dsqpi*mu) * dexp(-mu*mu*rij*rij)
-    jast_elec_Mu_value(i) += tmp 
+    jast_elec_Mu_value(j) += 0.5d0*tmp ! symmetrization 
   enddo
  enddo
- jast_elec_Mu_value = jast_elec_Mu_value * 0.5d0 ! symmetrization 
 END_PROVIDER
 
  BEGIN_PROVIDER [ double precision , jast_elec_Mu_grad_x, (elec_num_8) ]
@@ -44,13 +43,12 @@ END_PROVIDER
  double precision :: a, b, rij, tmp, x, y, z
  include '../constants.F'
  double precision :: mu
- mu = mu_erf
+ mu = jast_mu_erf
 
  do i=1,elec_num
-  jast_elec_Mu_grad_x(i) = 0.d0
-  jast_elec_Mu_grad_y(i) = 0.d0
-  jast_elec_Mu_grad_z(i) = 0.d0
-  !DIR$ LOOP COUNT (100)
+  jast_elec_Mu_grad_x(i) = jast_1b_grad_x(i)
+  jast_elec_Mu_grad_y(i) = jast_1b_grad_y(i)
+  jast_elec_Mu_grad_z(i) = jast_1b_grad_z(i)
  enddo
 
  ! (grad of J(r12) with respect to xi, yi, zi)
@@ -76,10 +74,10 @@ BEGIN_PROVIDER [ double precision , jast_elec_Mu_lapl, (elec_num_8) ]
  double precision :: a, b, rij, tmp, x, y, z
  include '../constants.F'
  double precision :: mu, x_ij, y_ij, z_ij, rij_inv
- mu = mu_erf
+ mu = jast_mu_erf
 
  do i=1,elec_num
-  jast_elec_Mu_lapl(i) = 0.d0
+  jast_elec_Mu_lapl(i) = jast_1b_lapl(i)
  enddo
 
  do i=1, elec_num
@@ -96,10 +94,6 @@ BEGIN_PROVIDER [ double precision , jast_elec_Mu_lapl, (elec_num_8) ]
  enddo
 END_PROVIDER
 
-BEGIN_PROVIDER [double precision, mu_erf ]
- implicit none
- mu_erf = 0.5d0
-END_PROVIDER 
 
  BEGIN_PROVIDER [double precision, grad_j_mu_x,(elec_num, elec_num)]
 &BEGIN_PROVIDER [double precision, grad_j_mu_y,(elec_num, elec_num)]
@@ -110,7 +104,7 @@ END_PROVIDER
  END_DOC
  integer :: i,j
  double precision :: rij, mu,scal
- mu = mu_erf
+ mu = jast_mu_erf
  grad_j_mu_x = 0.d0
  grad_j_mu_y = 0.d0
  grad_j_mu_z = 0.d0
@@ -126,3 +120,4 @@ END_PROVIDER
  enddo
 
 END_PROVIDER 
+
diff --git a/src/JASTROW/jastrow_param.irp.f b/src/JASTROW/jastrow_param.irp.f
index f0dca09..eac69b6 100644
--- a/src/JASTROW/jastrow_param.irp.f
+++ b/src/JASTROW/jastrow_param.irp.f
@@ -156,16 +156,24 @@ END_PROVIDER
 
 END_PROVIDER
 
+BEGIN_PROVIDER [ real, jast_mu_erf ]
+  implicit none
+  include '../types.F'
+  jast_mu_erf = 0.5
+  call get_jastrow_mu_erf(jast_mu_erf)
+END_PROVIDER
+
+
 BEGIN_PROVIDER  [ integer, jast_1b_type ]
  implicit none
  include '../types.F'
- jast_1b_type = 0 ! no  1body Jastrow
+! jast_1b_type = 0 ! no  1body Jastrow
  !jast_1b_type = 2 ! add 1body-Tanh   Jastrow
  !jast_1b_type = 3 ! add 1body-Simple Jastrow
  !jast_1b_type = 4 ! add 1body-RSDFT  Jastrow
  !jast_1b_type = 5 ! add 1body-erf    Jastrow
- !jast_1b_type = 6 ! add 1body-Gauss  Jastrow
- call get_jastrow_jast_1b_type(jast_1b_type)
+ jast_1b_type = 6 ! add 1body-Gauss  Jastrow
+! call get_jastrow_jast_1b_type(jast_1b_type)
 END_PROVIDER
 
 ! useful if  jast_1b_type = 2
diff --git a/src/Makefile b/src/Makefile
index 68e3a23..be94d05 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -5,3 +5,8 @@ export
 
 irpf90.make: $(filter-out IRPF90_temp/%, $(wildcard */*.irp.f)) $(wildcard *.irp.f) $(wildcard *.inc.f) Makefile 
 	$(IRPF90)
+
+IRPF90_temp/irp_touches.irp.o: IRPF90_temp/irp_touches.irp.F90
+	$(FC) -O2 -c -g  IRPF90_temp/irp_touches.irp.F90 -o IRPF90_temp/irp_touches.irp.o
+
+
diff --git a/src/PROPERTIES/properties_ci.irp.f b/src/PROPERTIES/properties_ci.irp.f
index 033ee2c..d7f8f8c 100644
--- a/src/PROPERTIES/properties_ci.irp.f
+++ b/src/PROPERTIES/properties_ci.irp.f
@@ -331,3 +331,71 @@ return
   ci_dress_max = max(ci_dress_max,maxval(ci_dress))
   SOFT_TOUCH ci_dress_min ci_dress_max
 END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, ci_dress_mu_opt ]
+  BEGIN_DOC
+  ! Use for optimizing mu
+  END_DOC
+  implicit none
+  integer          :: i, j, k, l
+  double precision :: T, dij, f, E_noJ, dE
+!  E_noJ = -0.5d0*psidet_lapl*psidet_inv + E_pot + E_nucl 
+!  ci_dress_mu_opt = (E_loc - energy_mu) * psi_value_inv * jast_value_inv * &
+!                    det_alpha_value(1) * det_beta_value(1)
+  ! energy_mu = H_mu \Phi / \Phi
+  dE = (E_loc - energy_mu) * psi_value_inv * jast_value_inv
+!  dE = (E_loc - E_noJ) * psi_value_inv * jast_value_inv
+  k = 1
+  i = det_coef_matrix_rows(   k)
+  j = det_coef_matrix_columns(k)
+  f = det_alpha_value(i) * det_beta_value(j)
+  ci_dress_mu_opt = dE * f
+  ci_dress_mu_opt_min = min(ci_dress_mu_opt_min, ci_dress_mu_opt)
+  ci_dress_mu_opt_max = max(ci_dress_mu_opt_max, ci_dress_mu_opt)
+  SOFT_TOUCH ci_dress_mu_opt_min ci_dress_mu_opt_max
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, ci_dress_Htilde, (size_ci_dress) ]
+  implicit none
+  BEGIN_DOC
+  ! < det(i) e^{-J} |H| Psi >
+  !
+  ! Dimensions : det_num
+  END_DOC
+
+  integer          :: i, j, k, l
+  double precision :: T, h_psidet, dij, f, E_noJ, dE
+
+  h_psidet = -0.5d0*psidet_lapl*psidet_inv + E_pot + E_nucl 
+  E_noJ = h_psidet 
+  dE = E_loc - E_noJ
+  do k=1,det_num
+     i = det_coef_matrix_rows(k)
+     j = det_coef_matrix_columns(k)
+     f = det_alpha_value(i)*det_beta_value(j) * psi_value_inv * jast_value_inv
+     ci_dress(k) = dE * f
+  enddo
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, ci_dress_H, (size_ci_dress) ]
+  implicit none
+  BEGIN_DOC
+  ! < det(i) e^{-J} |H| Psi >
+  !
+  ! Dimensions : det_num
+  END_DOC
+
+  integer          :: i, j, k, l
+  double precision :: T, h_psidet, dij, f, E_noJ, dE
+
+  E_noJ= -0.5d0*psidet_lapl*psidet_inv + E_pot + E_nucl 
+  do k=1,det_num
+     i = det_coef_matrix_rows(k)
+     j = det_coef_matrix_columns(k)
+     f = det_alpha_value(i)*det_beta_value(j) * psi_value_inv * jast_value_inv
+     ci_dress(k) = E_noJ * f
+  enddo
+
+END_PROVIDER
+
diff --git a/src/PROPERTIES/properties_mu.irp.f b/src/PROPERTIES/properties_mu.irp.f
index 5d9ae6b..6125149 100644
--- a/src/PROPERTIES/properties_mu.irp.f
+++ b/src/PROPERTIES/properties_mu.irp.f
@@ -57,8 +57,7 @@ END_PROVIDER
   integer :: i,j
   double precision :: rij, mu
 
-  mu = mu_erf
-!  mu = jast_mu_erf
+  mu = jast_mu_erf
   Eff_pot_mu_elec = 0.d0
 
 
@@ -85,18 +84,21 @@ END_PROVIDER
   enddo
 
   ! 1-body Jastrow
-!  if( jast_1b_type .gt. 0 ) then
-!    do i = 1, elec_num
-!      Eff_pot_mu_elec(i) -= 0.5d0 * jast_1b_lapl(i) 
-!      Eff_pot_mu_elec(i) -= 0.5d0 * jast_1b_grad_sq(i) 
-!      do j = 1, elec_num
-!        if(i==j) cycle
-!        Eff_pot_mu_elec(i) -= ( jast_elec_Mu_grad_x(i) * ( jast_1b_grad_x(i) - jast_1b_grad_x(j) ) &
-!                              + jast_elec_Mu_grad_y(i) * ( jast_1b_grad_y(i) - jast_1b_grad_y(j) ) &
-!                              + jast_elec_Mu_grad_z(i) * ( jast_1b_grad_z(i) - jast_1b_grad_z(j) ) )
-!      enddo
-!    enddo
-!  endif
+  if( jast_1b_type .gt. 0 ) then
+    do i = 1, elec_num
+      Eff_pot_mu_elec(i) -= 0.5d0 * jast_1b_lapl(i) 
+      Eff_pot_mu_elec(i) -= 0.5d0 * jast_1b_grad_sq(i) 
+      do j = 1, elec_num
+        if(i==j) cycle
+        ! + sign for i <--> j
+        ! 0.5d0  for double counting
+        Eff_pot_mu_elec(i) += 0.5d0 *                                                      &
+                            ( grad_j_mu_x(j,i) * ( jast_1b_grad_x(i) - jast_1b_grad_x(j) ) &
+                            + grad_j_mu_y(j,i) * ( jast_1b_grad_y(i) - jast_1b_grad_y(j) ) &
+                            + grad_j_mu_z(j,i) * ( jast_1b_grad_z(i) - jast_1b_grad_z(j) ) )
+      enddo
+    enddo
+  endif
 
 END_PROVIDER 
 
@@ -143,8 +145,7 @@ BEGIN_PROVIDER [double precision, Eff_pot_deriv_mu_elec, (elec_num) ]
   integer          :: i, j
   double precision :: rij, mu
 
-!  mu = jast_mu_erf
-  mu = mu_erf
+  mu = jast_mu_erf
   Eff_pot_deriv_mu_elec = 0.d0
 
   ! 2body-Jastrow: (eq A4)
@@ -237,3 +238,4 @@ BEGIN_PROVIDER [ double precision, ci_dress_mu, (size_ci_dress_mu) ]
   SOFT_TOUCH ci_dress_mu_min ci_dress_mu_max
 END_PROVIDER
 
+
diff --git a/src/ezfio_interface.irp.f b/src/ezfio_interface.irp.f
index ee71245..ac038f4 100644
--- a/src/ezfio_interface.irp.f
+++ b/src/ezfio_interface.irp.f
@@ -18,7 +18,7 @@ data = [ \
 ("ao_basis_ao_power"               , "integer"       , "(ao_num,3)"            ),
 ("ao_basis_ao_expo"                , "real"          , "(ao_num,ao_prim_num_max)" ),
 ("ao_basis_ao_coef"                , "real"          , "(ao_num,ao_prim_num_max)" ),
-("jastrow_mu_erf"            , "real"          , ""                   ),
+("jastrow_mu_erf"                  , "real"          , ""                   ),
 ("jastrow_jast_a_up_up"            , "real"          , ""                   ),
 ("jastrow_jast_a_up_dn"            , "real"          , ""                   ),
 ("jastrow_jast_b_up_up"            , "real"          , ""                   ),