In [33]:
from __future__ import print_function
%matplotlib inline
import collections
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
import matplotlib as mpl
from ase.io import read
from pyqstem.util import atoms_plot
from pyqstem import PyQSTEM
from pyqstem.potentials import poisson_solver,create_potential_slices
mpl.rc('font',**{'size' : 13})
In [34]:
epsilon0=0.07957747154594767 # permitivity of vacuum in atomic units
a0=0.52917721067 # bohr radius in angstrom
factor = 2/a0/epsilon0 # conversion factor between ASE units and QSTEM units
defocus = 0 # defocus [Angstrom]
v0 = 60 # acceleration voltage [keV]
alpha = 18 # convergence angle [mrad]
In [35]:
atoms=read('graphene.cif',index=0) # atomic configuration
tiled_atoms=atoms*(3,3,1)
cell=atoms.get_cell()
Lx,Ly,Lz=np.diag(cell)
qstem=PyQSTEM('STEM')
qstem.set_atoms(tiled_atoms)
scan_range=[[cell[0,0]*1,2*cell[0,0],30],
[cell[1,1]*1,2*cell[1,1],30]]
In [37]:
probes=[]
images=[]
spacings=[0.08,0.04,0.02]
for h in spacings:
print(h)
V=np.load('graphene_ae_potential_{0}.npy'.format(h)) # all-electron potential from GPAW using LDA
Nx,Ny,Nz=V.shape
res_x,res_y=Lx/Nx,Ly/Ny
V_dft=create_potential_slices(V,10,(Lx,Ly,Lz)) # create projected potential slices
V_dft.array=np.tile(V_dft.array,(3,3,1))
V_dft.array/=-factor
resolution = (res_x,res_y) # the resolution shoudl match the potential
samples = (int(300*0.02/h),int(300*0.02/h)) # samples in x and y-direction
qstem.build_probe(v0,alpha,samples,resolution=resolution,defocus=defocus)
probes.append(qstem.get_wave())
qstem.set_potential(V_dft,scan_range)
qstem.add_detector('det1',(70,200))
qstem.add_detector('det2',(0,70))
qstem.run()
images.append(qstem.read_detector('det1'))
images.append(qstem.read_detector('det2'))
In [43]:
fig,axarr=plt.subplots(1,len(probes),figsize=(10,6))
if not isinstance(axarr,collections.Iterable):
axarr=(axarr,)
for probe,ax,h in zip(probes,axarr,spacings):
probe.view(method='real',ax=ax,title='probe real (h={0})'.format(h),
cmap='seismic',vmin=-8,vmax=8)
plt.tight_layout()
In [39]:
fig,axarr=plt.subplots(len(images)//2,2,figsize=(8,9))
extent=[0,scan_range[0][1]-scan_range[0][0],0,scan_range[1][1]-scan_range[1][0]]
axarr=axarr.flatten()
for i,(img,ax) in enumerate(zip(images,axarr)):
ims=ax.imshow(img.T,extent=extent,interpolation='nearest',cmap='gray')
ax.set_title('h={0}'.format(spacings[i//2]))
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
plt.colorbar(ims, cax=cax)
plt.tight_layout()
In [7]:
images_charge=[]
spacings=[0.1,0.05]
for h in spacings:
print(h)
rho=np.load('graphene_ae_density_{0}.npy'.format(h)) # all-electron density from GPAW using LDA
V=poisson_solver(rho,atoms,smooth=0,units='QSTEM')
Nx,Ny,Nz=V.shape
res_x,res_y=Lx/Nx,Ly/Ny
V_dft=create_potential_slices(V,8,(Lx,Ly,Lz)) # create projected potential slices
V_dft.array=np.tile(V_dft.array,(3,3,1))
resolution = (res_x,res_y) # the resolution shoudl match the potential
samples = (int(400*0.02/h),int(400*0.02/h)) # samples in x and y-direction
qstem.build_probe(v0,alpha,samples,resolution=resolution,defocus=defocus)
probes.append(qstem.get_wave())
qstem.set_potential(V_dft,scan_range)
qstem.add_detector('det1',(70,200))
qstem.add_detector('det2',(0,70))
qstem.run()
images_charge.append(qstem.read_detector('det1'))
images_charge.append(qstem.read_detector('det2'))
In [8]:
fig,axarr=plt.subplots(len(images_charge)//2,2,figsize=(8,9))
extent=[0,scan_range[0][1]-scan_range[0][0],0,scan_range[1][1]-scan_range[1][0]]
axarr=axarr.flatten()
for i,(img,ax) in enumerate(zip(images_charge,axarr)):
ims=ax.imshow(img.T,extent=extent,interpolation='nearest',cmap='gray')
ax.set_title('h={0}'.format(spacings[i//2]))
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
plt.colorbar(ims, cax=cax)
plt.tight_layout()
In [32]:
fig,((ax1,ax2),(ax3,ax4))=plt.subplots(2,2,figsize=(8,6))
ims=ax1.imshow(images[2].T,extent=extent,cmap='gray')
divider = make_axes_locatable(ax1)
cax = divider.append_axes("right", size="5%", pad=0.05)
plt.colorbar(ims, cax=cax)
ax1.set_title('from ps2ae \n h=0.04')
ims=ax2.imshow(images[3].T,extent=extent,cmap='gray')
divider = make_axes_locatable(ax2)
cax = divider.append_axes("right", size="5%", pad=0.05)
plt.colorbar(ims, cax=cax)
ax2.set_title('from ps2ae \n h=0.04')
ims=ax3.imshow(images_charge[2].T,extent=extent,cmap='gray')
divider = make_axes_locatable(ax3)
cax = divider.append_axes("right", size="5%", pad=0.05)
plt.colorbar(ims, cax=cax)
ax3.set_title('from charge \n h=0.05')
ims=ax4.imshow(images_charge[3].T,extent=extent,cmap='gray')
divider = make_axes_locatable(ax4)
cax = divider.append_axes("right", size="5%", pad=0.05)
plt.colorbar(ims, cax=cax)
ax4.set_title('from charge \n h=0.05')
plt.tight_layout()