diff --git a/codes/position_selection.py b/codes/position_selection.py
new file mode 100644
index 0000000000000000000000000000000000000000..892cb5e910c2f3e31e124c451b849a49caef17e4
--- /dev/null
+++ b/codes/position_selection.py
@@ -0,0 +1,86 @@
+import os
+import numpy as np
+from sklearn.metrics import pairwise_kernels
+from ase.io import Trajectory
+from scipy.spatial import distance_matrix
+from sklearn.preprocessing import StandardScaler
+from create_descriptor import *
+
+
+def matrice_farthest_search(points, k,mode='euclidian',gamma=1e-5,n=3):
+         if mode is 'euclidian':
+             def distance_m(a):
+                return distance_matrix(a,a)
+         if mode is 'rbf':
+              def distance_m(a):
+                DM=pairwise_kernels(a,a,metric='rbf',gamma=gamma)
+                return 1-DM
+         DM=distance_m(points)
+         remaining_points = list(points)
+         solution_set = []
+         indices=[]
+         indices_cluster=[]
+         indiceleft=list(range(len(remaining_points)))
+         #init=indiceleft.pop(random.randint(0, len(indiceleft) - 1))
+         init=indiceleft.pop(np.argmin(np.mean(DM,axis=0)))         
+         indices.append(init)
+         for _ in range(k-1):
+             a=[min([DM[p][pl] for p in indices])for pl in indiceleft]
+             ind_rem=a.index(max(a))             
+             rem=indiceleft.pop(ind_rem)
+             indices.append(rem)   
+             solution_set.append(max(a))
+         solution_set=np.array(solution_set)
+         return solution_set, indices   
+
+
+
+
+if __name__ == '__main__':
+    #debut de l initialisation
+    path=os.getcwd()
+    atomseul='H' #atome depose
+    traj='final.traj' #bdd 1
+    data_folder='data'
+    gamma=1e-5
+    
+    sys=Trajectory(traj)  
+    species=[]
+    elem=sys[0].get_chemical_symbols() 
+    #identification de l atome seul dans la structure
+    for el in range(len(elem)): 
+        if elem[el] not in species:
+            species.append(elem[el])
+            if elem[el]== atomseul:
+                ats=el
+
+    # bases entieres
+    params={'species':species,'l_max':2,'n_max':2,'r_cut':7}
+    desc=create_descriptor(method='soap',params=params,ats=ats)
+    X_ini=desc.create(traj,load=True,save_file=data_folder)
+     
+    # bases poscars
+    index_pos_100=np.load(data_folder+'/ind_pos.npy')
+    X_poscar=[X_ini[i[0]] for i in index_pos_100]
+    len100=len(X_poscar)
+    print(len100)
+         
+    #scaling
+    scaler=StandardScaler()
+    scaler.fit(X_poscar)
+    X_poscar_T=scaler.transform(X_poscar)
+    ncomp=144        
+    #initialisation
+    distances,indices=matrice_farthest_search(X_poscar_T, ncomp,)
+    np.save("distances.npy",distances)
+    for i in range(1,145):           
+        if i in [4,9,16,25,36,49,64,81,100,]:
+            #resultpath    
+            resultpath=path+'/'+str(i)      
+            try:
+                os.mkdir(resultpath)
+            except:
+                pass
+            #save
+            print(len(indices[:i]))
+            np.save(resultpath+'/indices.npy',indices[:i])