From 5067e6948d369a4164111133a469f72071b2a800 Mon Sep 17 00:00:00 2001
From: BOULANGEOT Nathan <nathan.boulangeot@univ-lorraine.fr>
Date: Mon, 29 Jul 2024 07:08:17 +0000
Subject: [PATCH] Add AEM_image.py
---
Minimum Working Exemple/AEM_image.py | 137 +++++++++++++++++++++++++++
1 file changed, 137 insertions(+)
create mode 100644 Minimum Working Exemple/AEM_image.py
diff --git a/Minimum Working Exemple/AEM_image.py b/Minimum Working Exemple/AEM_image.py
new file mode 100644
index 0000000..60e2abb
--- /dev/null
+++ b/Minimum Working Exemple/AEM_image.py
@@ -0,0 +1,137 @@
+import numpy as np
+import ase
+import os
+import ase.io
+from ase.io import Trajectory
+import matplotlib.pyplot as plt
+from ase.data import covalent_radii
+plt.rcParams['font.size'] = 25
+plt.rcParams['axes.linewidth'] = 1.5
+
+
+def read_regu_all(result_path):
+ #Eslab and mu are the energy of the slab and the adsorbate alone.
+ #It needs to be manually changed to get the adsorption energy Eads=Etot-Eslab-mu
+ Eslab=-33.14301679*4
+ mu = -6.762788/2
+ at_train=Trajectory(result_path+'/train.traj')
+ y_mean=np.mean([at.get_potential_energy() for at in at_train])
+ ind=[]
+ Epred=[]
+ with open(result_path+'/resultat.txt','r') as f:
+ a=f.read().split('indice:')
+ for i,line in enumerate(a[1:]):
+ b=line.split()
+ ind.append(int(b[0]))
+ Epred.append(float(b[9]))
+ f.close()
+ ind,Epred = zip(*sorted(zip(ind,Epred)))
+ ind = list(ind)
+ Epred = np.array(Epred)+y_mean-Eslab-mu
+ return Epred
+
+def atomplot(ax,aX=0,aY=0,lw=1.3,POSCAR='POSCAR'):
+ c_radii=[]
+ el_colors=[]
+ at=ase.io.read(POSCAR,format='vasp')
+ at_number=at.get_atomic_numbers()
+ cell=at.cell
+ shift=[aX,aY,0]@cell
+ pos=at.get_positions()
+ zpos=pos[:,2]
+ zind=np.argsort(zpos)[-8:]
+ sub_surf_lim=13
+ for i in range(len(at)):
+ c_radii.append(covalent_radii[at_number[i]])
+ colors=['red','orange','purple']
+ el_colors=['k']*len(at_number)
+ for ian,an in enumerate(np.unique(at_number)):
+ for i in np.where(at_number==an)[0]:
+ el_colors[i]=colors[ian]
+ for i in zind:
+ if zpos[i]<sub_surf_lim:
+ ls=':'
+ else:
+ ls='-'
+ circle=plt.Circle((pos[i,0]+shift[0],pos[i,1]+shift[1]),c_radii[i],color=el_colors[i], fill=False,linestyle=ls,lw=lw)
+ ax.add_artist(circle)
+
+def replicate9(X,Y,E,cell):
+ X_tot=[]
+ Y_tot=[]
+ rep=[-1,0,1]
+ pos=np.array([[i,j,0] for i in rep for j in rep])
+ pos=pos@cell
+ xx=pos[:,0]
+ yy=pos[:,1]
+ for i in range(len(xx)):
+ X_tot=np.concatenate((X_tot,X+xx[i]))
+ Y_tot=np.concatenate((Y_tot,Y+yy[i]))
+ E=np.tile(E,9)
+ return X_tot,Y_tot,E
+
+
+def ads_mps(X,Y,E,ax,POSCAR='POSCAR',minE=None,maxE=None):
+ cmap='Blues_r'
+ if minE is None:
+ print('no minmaxE')
+ maxE=max(E)
+ minE=min(E)
+ ax.tricontourf(X,Y,E,100,cmap=cmap,vmax=maxE,vmin=minE)
+ ax.scatter(X,Y,s=1,c='k')
+ #-----------------
+ plotcell=(np.array([[0,1,1,0,0],[0,0,1,1,0],[0,0,0,0,0]]).T@cell).T[:2]
+ ax.plot(plotcell[0],plotcell[1],'k-')
+ ax.set_xlim(-cell[0,0]*0.15,1.55*cell[0,0])
+ ax.set_ylim(-cell[1,1]*0.35,1.35*cell[1,1])
+ ax.set_xlabel(r'$x$ ($\mathregular{\AA}$)',labelpad=-15)#,rotation=180)
+ ax.set_ylabel(r'$y$ ($\mathregular{\AA}$)',labelpad=-15)#,rotation=180)
+
+ ax.set_aspect('equal')
+ ax.minorticks_on()
+ ax.tick_params(which='major', width=1.5)
+ ax.tick_params(which='minor', width=1)
+ ax.tick_params(which='major', length=7)
+ ax.tick_params(which='minor', length=4)
+ ax.xaxis.set_tick_params(which='both', direction='in',top=True,labelbottom=True)
+ ax.yaxis.set_tick_params(which='both', direction='in',top=True,labelbottom=True)
+ for i in [-1,0,1]:
+ for j in [-1,0,1]:
+ atomplot(ax,i,j,POSCAR=POSCAR)
+
+cmap='Blues_r'
+postrain_props={'s':58,'c':'yellow','marker':'^','edgecolor':'k','linewidths':.9,'zorder':100}
+path_props={'s':58,'c':'green','marker':'o','edgecolor':'k','linewidths':.9,'zorder':100}
+path_ini='./'
+traj=Trajectory(path_ini+'/train.traj')
+E_tot=[t.get_potential_energy() for t in traj]
+Eslab=-33.14301679*4
+mu = -6.762788/2
+minE=min(E_tot)-Eslab-mu
+maxE=max(E_tot)-Eslab-mu
+print(minE,maxE)
+
+os.chdir(path_ini)
+POSCAR=path_ini+'/POSCAR'
+Xi,Yi,_=pos_to_relax(POSCAR)
+poscar=ase.io.read(POSCAR)
+cell=poscar.get_cell()
+fig,ax=plt.subplots(1,1,figsize=(8,6))
+result_path=path_ini
+Ei=read_regu_all(result_path)
+minE=min(Ei)
+maxE=max(Ei)
+X,Y,E=replicate9(Xi,Yi,np.array(Ei),cell)
+ads_mps(X,Y,E,ax,POSCAR=POSCAR,minE=minE,maxE=maxE)
+cax = plt.axes([0.2,-0.1,0.6, 0.03])
+m = plt.cm.ScalarMappable(cmap=cmap)
+cbar = plt.colorbar(m,boundaries=np.arange(minE,maxE,0.001),cax=cax,label=r'E$_{ads}$ H / Ag (eV)',orientation='horizontal')
+cbar.ax.tick_params(labelsize=20)
+cbar.formatter.set_powerlimits((0, 0))
+cax.xaxis.set_ticks_position("top")
+ticks=np.around(np.linspace(minE+.002,maxE-.002,5),2)
+cbar.set_ticks(ticks)
+
+ax.scatter(np.array(postrain)[:,0],np.array(postrain)[:,1],**postrain_props)
+fig.savefig(path_ini+'Ag111_AEM_'+str(i+1)+'pos.pdf', bbox_inches='tight',pad_inches=0.1, dpi = 300)
+
--
GitLab