Delete zopt.py

codes/zopt.py

-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])
-
-   
-