Add AEM_image.py

Minimum Working Exemple/AEM_image.py

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