diff --git a/codes/tool.py b/codes/tool.py
new file mode 100644
index 0000000000000000000000000000000000000000..b536607c02dba0e88deea1611f1d8acf0fc4c730
--- /dev/null
+++ b/codes/tool.py
@@ -0,0 +1,144 @@
+import numpy as np
+import os
+import ase
+import ase.io
+from ase import Atoms, Atom
+from ase.io import Trajectory
+from ase.constraints import FixScaled,FixAtoms
+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 GridSearchCV
+from joblib import Parallel, delayed
+from ase import neighborlist
+from sklearn.preprocessing import StandardScaler
+from create_descriptor import *
+
+
+
+def sign(x):
+ if x<0:
+ return -1
+ else:
+ return 1
+
+def get_train_data(traj,data_folder,train_pos_indices):
+ sys=Trajectory(traj)
+ ind_slab=sys[0].constraints[1].index
+ ind_free=np.delete(np.arange(len(sys[0])),ind_slab)
+ species=[]
+ elem=sys[0].get_chemical_symbols()
+ #identification des espèces
+ for el in range(len(elem)):
+ if elem[el] not in species:
+ species.append(elem[el])
+ y_ini=np.load(data_folder+'/E.npy')
+ index_pos=np.load(data_folder+'/ind_pos.npy')
+ ind_train=[]
+ #select only 3 config per position
+ for i in range(len(index_pos)):
+ if i in train_pos_indices:
+ ind_train.append(index_pos[i,0])
+ ind_train.append((index_pos[i,1]-index_pos[i,0]+1)//2)
+ ind_train.append(index_pos[i,1])
+
+ y_train=np.array([y_ini[i] for i in ind_train])
+ return [sys[i] for i in ind_train],y_train,species
+
+def pos_to_relax(poscar_file,n=20):
+ poscar=ase.io.read(poscar_file)
+ #position à relaxer
+ cell=poscar.get_cell()
+ a=cell[0,0]
+ b=cell[1,1]
+ c=cell[2,2]
+ #maillage
+ n=20
+ Xpos=[i*a/n for i in range(n)]
+ Ypos=[i*b/n for i in range(n)]
+
+ pos=np.array([[i/n,j/n,0] for i in range(n) for j in range(n)])
+ pos=pos@cell
+ xx,yy=pos[:,0],pos[:,1]
+ zz=np.empty(len(xx))
+ #zpos
+ atom_top=poscar.get_positions()
+ for i in range(len(xx)):
+ dx = atom_top[:,0] - xx[i]
+ dy = atom_top[:,1] - yy[i]
+ dz = atom_top[:,2] - 18
+ dxy=np.linalg.norm([dx,dy,dz],axis=0)
+ indice=np.where(dxy==np.min(dxy))[0][0]
+ zz[i]=atom_top[indice,2]+1.6
+ return xx,yy,zz
+
+def searchEminzgrid(desc,p,ind=0,atomseul='H',poscar='POSCAR'):
+ pos=p.copy()
+ zetap=[]
+ traj_tot=[]
+ NB=False
+ try:
+ traj_ini=ase.io.read(poscar,format='vasp')
+ except:
+ traj_ini=poscar.copy()
+ 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,poscar='POSCAR'):
+ x_pred,zetap,ind,zlim=searchEminzgrid(desc,p,ind=ind,atomseul='H',poscar=poscar)
+ 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])
\ No newline at end of file