diff --git a/codes/zopt.py b/codes/zopt.py deleted file mode 100644 index 175bf42ffb95afc3b49357d88b44390deaf74856..0000000000000000000000000000000000000000 --- a/codes/zopt.py +++ /dev/null @@ -1,86 +0,0 @@ -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]) - - -