merge zopt+tool def

codes/tool.py

+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