diff --git a/codes/zopt.py b/codes/zopt.py
new file mode 100644
index 0000000000000000000000000000000000000000..175bf42ffb95afc3b49357d88b44390deaf74856
--- /dev/null
+++ b/codes/zopt.py
@@ -0,0 +1,86 @@
+import os
+import sys
+import ase.io
+import numpy as np
+from ase import Atoms, Atom
+from ase.io import Trajectory
+from sklearn.metrics import mean_squared_error
+from dscribe.descriptors import SOAP
+from scipy.interpolate import UnivariateSpline
+from sklearn.gaussian_process import GaussianProcessRegressor
+from sklearn.gaussian_process.kernels import RBF
+from sklearn.model_selection import train_test_split, learning_curve, validation_curve, GridSearchCV
+from joblib import Parallel, delayed
+from ase import neighborlist 
+from tool import *
+from sklearn.preprocessing import StandardScaler
+from create_descriptor import *
+
+def searchEminzgrid(desc,p,ind=0,atomseul='H',poscar='POSCAR'):   
+            pos=p.copy()        
+            zetap=[]
+            traj_tot=[]
+            NB=False
+            traj_ini=ase.io.read(poscar,format='vasp')
+            cell=traj_ini.get_cell()
+            while NB is False:
+                traj=traj_ini+Atoms(atomseul,[pos],pbc=True,cell=cell)             
+                traj_tot.append(traj)
+                zetap.append(pos[2])
+                cutoffs=neighborlist.natural_cutoffs(traj)
+                nl = neighborlist.NeighborList(cutoffs,skin=0,bothways=True,self_interaction=False)
+                nl.update(traj)
+                indices, offsets = nl.get_neighbors(-1)
+                pos[2]-=0.1
+                if len(indices)>0:
+                    NB=True
+            zlim=pos[2]+0.1    
+            for _ in range(5):
+                traj=traj_ini+Atoms(atomseul,[pos],pbc=True,cell=cell)
+                traj_tot.append(traj)
+                zetap.append(pos[2])
+                pos[2]-=0.1
+            #print(traj,traj[0])   
+            print(len(zetap)) 
+            x_pred=desc.create(traj_tot)
+            return x_pred,zetap,ind,zlim
+        
+        
+    
+    
+def pos_reduc(zetap,Eetap,Estd,ind,zlim,xpos,ypos): 
+        y_spl = UnivariateSpline(np.flip(zetap),np.flip(Eetap),s=0,k=3)
+        y_spl_2d = y_spl.derivative(n=2)
+        y_spl_1d = y_spl.derivative(n=1)
+        minis=[]
+        yspld2m1=y_spl_2d(zetap[2])
+        yspld1m1=y_spl_1d(zetap[2])
+        for it in range(200):
+            zz=zetap[0]+it/200*(zetap[-1]-zetap[0])           
+            if y_spl_2d(zz)>0:
+                         if sign(y_spl_1d(zz))!=sign(yspld1m1):
+                                 j=np.where(abs(zetap-zz)==min(abs(zetap-zz)))[0][0]
+                                 minis.append(j)
+            yspld2m1=y_spl_2d(zz)
+            yspld1m1=y_spl_1d(zz)
+        try:
+            minE=Eetap[minis[np.argsort(Eetap[minis]+Estd[minis])[0]]]
+        except:
+            minE=Eetap[np.argsort(Eetap+Estd)[0]]
+            print('no min for i=',ind)
+        E=minE
+        z=zetap[np.where(Eetap==minE)[0][0]]
+        Estd=Estd[np.where(Eetap==minE)[0][0]]                
+        
+        with open('resultat.txt', 'a') as f:
+                  f.write('indice: '+str(ind)+' pos: '+str(xpos)+' '+str(ypos)+' zlim: '+str(zlim)+' zpred: '+str(z)+' Epred: '+str(E)+' Estd: '+str(Estd)+'\n')
+        f.close()
+
+def run_para(p,ind,desc,scaler,g_search):
+        x_pred,zetap,ind,zlim=searchEminzgrid(desc,p,ind=ind,atomseul='H')  
+        x_pred=scaler.transform(x_pred)          
+        Eprd,Estd=g_search.best_estimator_.predict(x_pred,return_std=True)
+        pos_reduc(np.array(zetap),Eprd,Estd,ind,zlim,p[0],p[1])
+
+   
+