Experimental structure Refernce for Bi2Te3 (starting point):
['Atuchin, V.V.; Gavrilova, T.A.; Kokh, K.A.; Kuratieva, N.V.; Pervukhina, '
'N.V.; Surovtsev, N.V.',
'Structural and vibrational properties of PVT grown Bi2 Te3 microcrystals',
'Solid State Communications (2012) 152, (13) Article ID * (p. 1119-1122)',
'10.1016/j.ssc.2012.04.007']
%load_ext autoreload
%autoreload 2
from aiida import load_profile
from aiida import orm
import numpy as np
import matplotlib.pyplot as plt
from aiida_kkr.tools import plot_kkr
profile = load_profile()
def only_Pd(struc):
"""
Create auxiliary structure keeping only the Pd sites
"""
s = orm.StructureData(cell=struc.cell)
for site in struc.sites:
# print('Pd' in site.kind_name, site.kind_name)
if 'Pd' in site.kind_name:
s.append_atom(position=site.position, symbols='Pd')
return s
Relax_workchain = orm.load_node('fe892256-1575-479f-8bb7-3c6a54f7dc18')
struc_inp = Relax_workchain.called[0].called[1].called[0].outputs.fleurinpData.get_structuredata()
struc_out = Relax_workchain.outputs.optimized_structure
plot_kkr(struc_out, silent=True, repeat_uc=(1,1,2))
p1_all = np.array([i.position for i in struc_inp.sites])
isort = (p1_all[:,0] + p1_all[:,1]*1000 + p1_all[:,2]*1000**2).argsort()
p1_all = np.round(p1_all[isort], 2)
sites1 = np.array(struc_inp.sites)[isort]
sites2 = np.array(struc_out.sites)[isort]
for i, s1 in enumerate(sites1):
# if i<28:
s2 = sites2[i]
p1, p2 = np.round(np.array(s1.position), 2), np.round(np.array(s2.position), 2)
if i==0:
print(' atom dist pos input diff vec rel diff')
print('---------------------------------------------------------------------')
print(f'{i:3} {s1.kind_name.split("(")[1].split(")")[0]:3} {np.linalg.norm(p1-p2):.2f} {p1} {p1-p2} {np.round(abs((p1-p2)/p1), 2)}')
scale = 0.1
p1 = np.array([i.position for i in sites1[p1_all[:,2]==-4.43]])
p2 = np.array([i.position for i in sites2[p1_all[:,2]==-4.43]])
plt.plot(p1[:,0], p1[:,1], 'ro', ms=10, label='Bi')
plt.plot(p2[:,0], p2[:,1], 'bs', ms=5, label='Bi, relax')
plt.quiver(p1[:,0], p1[:,1], (p2-p1)[:,0], (p2-p1)[:,1], color='r', scale=scale, units='xy', angles='xy', scale_units='xy' )
p1 = np.array([i.position for i in sites1[p1_all[:,2]==-3.41]])
p2 = np.array([i.position for i in sites2[p1_all[:,2]==-3.41]])
plt.plot(p1[:,0], p1[:,1], 'o', ms=10, color='C0', label='Pd')
plt.plot(p2[:,0], p2[:,1], 's', ms=5, color='C1', label='Pd, relax')
plt.quiver(p1[:,0], p1[:,1], (p2-p1)[:,0], (p2-p1)[:,1], color='C0', scale=scale, units='xy', angles='xy', scale_units='xy' )
p1 = np.array([i.position for i in sites1[p1_all[:,2]==-2.76]])
p2 = np.array([i.position for i in sites2[p1_all[:,2]==-2.76]])
plt.plot(p1[:,0], p1[:,1], 'o', ms=10, color='C2', label='Te')
plt.plot(p2[:,0], p2[:,1], 's', ms=5, color='C3', label='Te, relax')
plt.quiver(p1[:,0], p1[:,1], (p2-p1)[:,0], (p2-p1)[:,1], color='C2', scale=scale, units='xy', angles='xy', scale_units='xy' )
plt.legend(ncol=3, loc=3)
plt.axis('equal')
plt.xlabel('x ($\AA$)')
plt.ylabel('y ($\AA$)')
# cell = np.array(struc_inp.cell)
# box = np.array([
# [0, 0],
# (cell[0])[:2],
# (cell[0]+cell[1])[:2],
# (cell[1])[:2],
# [0, 0],
# ])
# box[:,0] += p1_all[:,0].min() - 1.0
# box[:,1] += p1_all[:,1].min() + 1.5
# plt.plot(box[:,0], box[:,1], 'k-', lw=1)
# plt.savefig('relax_shift_xy.pdf')
plt.show()
scale = 0.3
dz = np.array(struc_inp.cell[2])
p1 = np.array([i.position for i in sites1[p1_all[:,2]==-4.43]]) + dz
p2 = np.array([i.position for i in sites2[p1_all[:,2]==-4.43]]) + dz
p1 = p1[p1[:,0]>-4]; p2 = p2[p2[:,0]>-4]; p1 = p1[p1[:,0]<7]; p2 = p2[p2[:,0]<7]
plt.plot(p1[:,0], p1[:,2], 'ro', ms=10, label='Bi')
plt.plot(p2[:,0], p2[:,2], 'bs', ms=5, label='Bi, relax')
plt.quiver(p1[:,0], p1[:,2], (p2-p1)[:,0], (p2-p1)[:,2], color='r', scale=scale, units='xy', angles='xy', scale_units='xy' )
p1 = np.array([i.position for i in sites1[p1_all[:,2]==-3.41]]) + dz
p2 = np.array([i.position for i in sites2[p1_all[:,2]==-3.41]]) + dz
p1 = p1[p1[:,0]>-4]; p2 = p2[p2[:,0]>-4]; p1 = p1[p1[:,0]<7]; p2 = p2[p2[:,0]<7]
plt.plot(p1[:,0], p1[:,2], 'o', ms=10, color='C0', label='Pd')
plt.plot(p2[:,0], p2[:,2], 's', ms=5, color='C1', zorder=50, alpha=1, label='Pd, relax')
plt.quiver(p1[:,0], p1[:,2], (p2-p1)[:,0], (p2-p1)[:,2], color='C0', scale=scale, units='xy', angles='xy', scale_units='xy' )
p1 = np.array([i.position for i in sites1[p1_all[:,2]==-2.76]]) + dz
p2 = np.array([i.position for i in sites2[p1_all[:,2]==-2.76]]) + dz
p1 = p1[p1[:,0]>-4]; p2 = p2[p2[:,0]>-4]; p1 = p1[p1[:,0]<7]; p2 = p2[p2[:,0]<7]
plt.plot(p1[:,0], p1[:,2], 'o', ms=10, color='C2', label='Te')
plt.plot(p2[:,0], p2[:,2], 's', ms=5, color='C3', label='Te, relax')
plt.quiver(p1[:,0], p1[:,2], (p2-p1)[:,0], (p2-p1)[:,2], color='C2', scale=scale, units='xy', angles='xy', scale_units='xy' )
p1 = np.array([i.position for i in sites1[p1_all[:,2]==0]]) + dz
p2 = np.array([i.position for i in sites2[p1_all[:,2]==0]]) + dz
p1 = p1[p1[:,0]>-4]; p2 = p2[p2[:,0]>-4]; p1 = p1[p1[:,0]<7]; p2 = p2[p2[:,0]<7]
plt.plot(p1[:,0], p1[:,2], 'o', ms=10, color='C4', label='Te')
plt.plot(p2[:,0], p2[:,2], 's', ms=5, color='C5', label='Te, relax')
plt.quiver(p1[:,0], p1[:,2], (p2-p1)[:,0], (p2-p1)[:,2], color='C4', scale=scale, units='xy', angles='xy', scale_units='xy' )
p1 = np.array([i.position for i in sites1[p1_all[:,2]==0]])
p2 = np.array([i.position for i in sites2[p1_all[:,2]==0]])
p1 = p1[p1[:,0]>-4]; p2 = p2[p2[:,0]>-4]; p1 = p1[p1[:,0]<7]; p2 = p2[p2[:,0]<7]
plt.plot(p1[:,0], p1[:,2], 'o', ms=10, color='C4')
plt.plot(p2[:,0], p2[:,2], 's', ms=5, color='C5')
plt.quiver(p1[:,0], p1[:,2], (p2-p1)[:,0], (p2-p1)[:,2], color='C4', scale=scale, units='xy', angles='xy', scale_units='xy' )
p1 = np.array([i.position for i in sites1[p1_all[:,2]==1.67]])
p2 = np.array([i.position for i in sites2[p1_all[:,2]==1.67]])
p1 = p1[p1[:,0]>-4]; p2 = p2[p2[:,0]>-4]; p1 = p1[p1[:,0]<7]; p2 = p2[p2[:,0]<7]
plt.plot(p1[:,0], p1[:,2], 'o', ms=10, color='C6', label='Bi')
plt.plot(p2[:,0], p2[:,2], 's', ms=5, color='C7', label='Bi, relax')
plt.quiver(p1[:,0], p1[:,2], (p2-p1)[:,0], (p2-p1)[:,2], color='C6', scale=scale, units='xy', angles='xy', scale_units='xy' )
p1 = np.array([i.position for i in sites1[p1_all[:,2]==3.7]])
p2 = np.array([i.position for i in sites2[p1_all[:,2]==3.7]])
p1 = p1[p1[:,0]>-4]; p2 = p2[p2[:,0]>-4]; p1 = p1[p1[:,0]<7]; p2 = p2[p2[:,0]<7]
plt.plot(p1[:,0], p1[:,2], 'o', ms=10, color='C8', label='Te')
plt.plot(p2[:,0], p2[:,2], 's', ms=5, color='C9', label='Te, relax')
plt.quiver(p1[:,0], p1[:,2], (p2-p1)[:,0], (p2-p1)[:,2], color='C8', scale=scale, units='xy', angles='xy', scale_units='xy' )
plt.axhline(8.5, color='grey', zorder=-100, ls='--')
plt.legend(ncol=1)
plt.ylabel('z ($\AA$)')
plt.xlabel('x ($\AA$)')
plt.axis('equal')
plt.xlim(-4,12)
# plt.savefig('relax_shift_xz.pdf')
plt.show()
# 1 Pd in 3x3 cell
Relax_workchain = orm.load_node(uuid = '024e7d89-473e-46dd-bee4-ba364439fc81')
# struc_inp = Relax_workchain.called[0].called[1].called[0].outputs.fleurinpData.get_structuredata()
struc_out = Relax_workchain.outputs.optimized_structure
plot_kkr(struc_out, silent=True, repeat_uc=(1,1,2), canvas_size=(800,600))
plot_kkr(only_Pd(struc_out), silent=True, repeat_uc=(1,1,2), canvas_size=(800,600))
plot_kkr(only_Pd(struc_out), silent=True, repeat_uc=(1,1,2), canvas_size=(800,600), rotations=('-10x,0y,180z'))
# second calculation with 2Pd atoms in 3x3 unit cell
Relax_workchain2 = orm.load_node(uuid = 'dcc785ae-fee5-40a2-a568-da301b4f3dd4')
struc_out = Relax_workchain2.outputs.optimized_structure
plot_kkr(struc_out, silent=True, repeat_uc=(1,1,2), canvas_size=(800,600))
plot_kkr(only_Pd(struc_out), silent=True, repeat_uc=(1,1,2), canvas_size=(800,600))
plot_kkr(only_Pd(struc_out), silent=True, repeat_uc=(1,1,2), canvas_size=(800,600), rotations=('-10x,0y,180z'))
# third calculation with 3Pd atoms in 3x3 unit cell
Relax_workchain3 = orm.load_node(uuid = '6551418a-39fc-4006-bf2d-9f41b1860d76')
struc_out = Relax_workchain3.outputs.optimized_structure
plot_kkr(struc_out, silent=True, repeat_uc=(1,1,2), canvas_size=(800,600))
plot_kkr(only_Pd(struc_out), silent=True, repeat_uc=(1,1,2), canvas_size=(800,600))
plot_kkr(only_Pd(struc_out), silent=True, repeat_uc=(1,1,2), canvas_size=(800,600), rotations=('-10x,0y,180z'))
# fourth calculation with 2Pd atoms in 3x3 unit cell, vertically above each other
Relax_workchain4 = orm.load_node(uuid = '214e2068-4aea-4d87-84d3-8596ff27bfc8')
struc_out = Relax_workchain4.outputs.optimized_structure
plot_kkr(struc_out, silent=True, repeat_uc=(1,1,2), canvas_size=(800,600))
plot_kkr(only_Pd(struc_out), silent=True, repeat_uc=(1,1,2), canvas_size=(800,600))
plot_kkr(only_Pd(struc_out), silent=True, repeat_uc=(1,1,2), canvas_size=(800,600), rotations=('-10x,0y,180z'))
# plot_kkr(Relax_workchain.inputs.scf.structure, silent=True, repeat_uc=(2,2,2))
# relax = Relax_workchain; zPd, zPd2 = -1.48, None # 1 Pd in vdW gap
# relax = Relax_workchain2; zPd, zPd2 = -1.58, -1.38 # 2 Pd in vdW gap
# relax = Relax_workchain3; zPd, zPd2 = -1.58, -1.38 # 3 Pd in vdW gap
relax = Relax_workchain4; zPd, zPd2 = -2.08, -0.88 # 2 Pd in vdW gap vertically above each other
struc_inp = relax.called[0].called[1].called[0].outputs.fleurinpData.get_structuredata()
# plot_kkr(struc_inp, silent=True, repeat_uc=(2,2,2))
# plot_kkr(struc_inp, silent=True, repeat_uc=(1,1,2))
struc_out = relax.outputs.optimized_structure
# plot_kkr(struc_out, silent=True, repeat_uc=(2,2,2))
plot_kkr(struc_out, silent=True, repeat_uc=(1,1,2))
p1_all = np.array([i.position for i in struc_inp.sites])
isort = (p1_all[:,0] + p1_all[:,1]*1000 + p1_all[:,2]*1000**2).argsort()
p1_all = np.round(p1_all[isort], 2)
sites1 = np.array(struc_inp.sites)[isort]
sites2 = np.array(struc_out.sites)[isort]
for i, s1 in enumerate(sites1):
# if i<28:
s2 = sites2[i]
p1, p2 = np.round(np.array(s1.position), 2), np.round(np.array(s2.position), 2)
if i==0:
print(' atom dist pos input diff vec rel diff')
print('---------------------------------------------------------------------')
print(f'{i:3} {s1.kind_name.split("(")[1].split(")")[0]:3} {np.linalg.norm(p1-p2):.2f} {p1} {p1-p2} {np.round(abs((p1-p2)/p1), 2)}')
scale = 0.1
p1 = np.array([i.position for i in sites1[p1_all[:,2]==-4.43]])
p2 = np.array([i.position for i in sites2[p1_all[:,2]==-4.43]])
plt.plot(p1[:,0], p1[:,1], 'ro', ms=10, label='Bi')
plt.plot(p2[:,0], p2[:,1], 'bs', ms=5, label='Bi, relax')
plt.quiver(p1[:,0], p1[:,1], (p2-p1)[:,0], (p2-p1)[:,1], color='r', scale=scale, units='xy', angles='xy', scale_units='xy' )
p1 = np.array([i.position for i in sites1[p1_all[:,2]==zPd]])
p2 = np.array([i.position for i in sites2[p1_all[:,2]==zPd]])
plt.plot(p1[:,0], p1[:,1], 'o', ms=10, color='C0', label='Pd')
plt.plot(p2[:,0], p2[:,1], 's', ms=5, color='C1', label='Pd, relax')
plt.quiver(p1[:,0], p1[:,1], (p2-p1)[:,0], (p2-p1)[:,1], color='C0', scale=scale, units='xy', angles='xy', scale_units='xy' )
if zPd2 is not None:
p1 = np.array([i.position for i in sites1[p1_all[:,2]==zPd2]])
p2 = np.array([i.position for i in sites2[p1_all[:,2]==zPd2]])
# p1 = np.array([i.position for i in sites1[p1_all[:,2]==zPd2]])+cell[0]+cell[1]
# p2 = np.array([i.position for i in sites2[p1_all[:,2]==zPd2]])+cell[0]+cell[1]
plt.plot(p1[:,0], p1[:,1], 'o', ms=10, color='C0', label='Pd')
plt.plot(p2[:,0], p2[:,1], 's', ms=5, color='C1', label='Pd, relax')
plt.quiver(p1[:,0], p1[:,1], (p2-p1)[:,0], (p2-p1)[:,1], color='C0', scale=scale, units='xy', angles='xy', scale_units='xy' )
p1 = np.array([i.position for i in sites1[p1_all[:,2]==-2.76]])
p2 = np.array([i.position for i in sites2[p1_all[:,2]==-2.76]])
plt.plot(p1[:,0], p1[:,1], 'o', ms=10, color='C2', label='Te')
plt.plot(p2[:,0], p2[:,1], 's', ms=5, color='C3', label='Te, relax')
plt.quiver(p1[:,0], p1[:,1], (p2-p1)[:,0], (p2-p1)[:,1], color='C2', scale=scale, units='xy', angles='xy', scale_units='xy' )
# plt.legend(ncol=3, loc=3)
plt.axis('equal')
plt.xlabel('x ($\AA$)')
plt.ylabel('y ($\AA$)')
# cell = np.array(struc_inp.cell)
# box = np.array([
# [0, 0],
# (cell[0])[:2],
# (cell[0]+cell[1])[:2],
# (cell[1])[:2],
# [0, 0],
# ])
# box[:,0] += p1_all[:,0].min() - 1.0
# box[:,1] += p1_all[:,1].min() + 1.5
# plt.plot(box[:,0], box[:,1], 'k-', lw=1)
# plt.savefig('relax_shift_xy.pdf')
plt.show()
scale = 0.3
dz = np.array(struc_inp.cell[2])
p1 = np.array([i.position for i in sites1[p1_all[:,2]==-4.43]]) + dz
p2 = np.array([i.position for i in sites2[p1_all[:,2]==-4.43]]) + dz
p1 = p1[p1[:,0]>-4]; p2 = p2[p2[:,0]>-4]; p1 = p1[p1[:,0]<7]; p2 = p2[p2[:,0]<7]
plt.plot(p1[:,0], p1[:,2], 'ro', ms=10, label='Bi')
plt.plot(p2[:,0], p2[:,2], 'bs', ms=5, label='Bi, relax')
plt.quiver(p1[:,0], p1[:,2], (p2-p1)[:,0], (p2-p1)[:,2], color='r', scale=scale, units='xy', angles='xy', scale_units='xy' )
p1 = np.array([i.position for i in sites1[p1_all[:,2]==zPd]]) + dz
p2 = np.array([i.position for i in sites2[p1_all[:,2]==zPd]]) + dz
p1 = p1[p1[:,0]>-4]; p2 = p2[p2[:,0]>-4]; p1 = p1[p1[:,0]<7]; p2 = p2[p2[:,0]<7]
plt.plot(p1[:,0], p1[:,2], 'o', ms=10, color='C0', label='Pd')
plt.plot(p2[:,0], p2[:,2], 's', ms=5, color='C1', zorder=50, alpha=1, label='Pd, relax')
plt.quiver(p1[:,0], p1[:,2], (p2-p1)[:,0], (p2-p1)[:,2], color='C0', scale=scale, units='xy', angles='xy', scale_units='xy' )
if zPd2 is not None:
# p1 = np.array([i.position for i in sites1[p1_all[:,2]==zPd2]])+cell[0]+cell[1]+dz
# p2 = np.array([i.position for i in sites2[p1_all[:,2]==zPd2]])+cell[0]+cell[1]+dz
p1 = np.array([i.position for i in sites1[p1_all[:,2]==zPd2]])+dz
p2 = np.array([i.position for i in sites2[p1_all[:,2]==zPd2]])+dz
plt.plot(p1[:,0], p1[:,2], 'o', ms=10, color='C0', label='Pd')
plt.plot(p2[:,0], p2[:,2], 's', ms=5, color='C1', label='Pd, relax')
plt.quiver(p1[:,0], p1[:,2], (p2-p1)[:,0], (p2-p1)[:,2], color='C0', scale=scale, units='xy', angles='xy', scale_units='xy' )
p1 = np.array([i.position for i in sites1[p1_all[:,2]==-2.76]]) + dz
p2 = np.array([i.position for i in sites2[p1_all[:,2]==-2.76]]) + dz
p1 = p1[p1[:,0]>-4]; p2 = p2[p2[:,0]>-4]; p1 = p1[p1[:,0]<7]; p2 = p2[p2[:,0]<7]
plt.plot(p1[:,0], p1[:,2], 'o', ms=10, color='C2', label='Te')
plt.plot(p2[:,0], p2[:,2], 's', ms=5, color='C3', label='Te, relax')
plt.quiver(p1[:,0], p1[:,2], (p2-p1)[:,0], (p2-p1)[:,2], color='C2', scale=scale, units='xy', angles='xy', scale_units='xy' )
p1 = np.array([i.position for i in sites1[p1_all[:,2]==0]]) + dz
p2 = np.array([i.position for i in sites2[p1_all[:,2]==0]]) + dz
p1 = p1[p1[:,0]>-4]; p2 = p2[p2[:,0]>-4]; p1 = p1[p1[:,0]<7]; p2 = p2[p2[:,0]<7]
plt.plot(p1[:,0], p1[:,2], 'o', ms=10, color='C4', label='Te')
plt.plot(p2[:,0], p2[:,2], 's', ms=5, color='C5', label='Te, relax')
plt.quiver(p1[:,0], p1[:,2], (p2-p1)[:,0], (p2-p1)[:,2], color='C4', scale=scale, units='xy', angles='xy', scale_units='xy' )
p1 = np.array([i.position for i in sites1[p1_all[:,2]==0]])
p2 = np.array([i.position for i in sites2[p1_all[:,2]==0]])
p1 = p1[p1[:,0]>-4]; p2 = p2[p2[:,0]>-4]; p1 = p1[p1[:,0]<7]; p2 = p2[p2[:,0]<7]
plt.plot(p1[:,0], p1[:,2], 'o', ms=10, color='C4')
plt.plot(p2[:,0], p2[:,2], 's', ms=5, color='C5')
plt.quiver(p1[:,0], p1[:,2], (p2-p1)[:,0], (p2-p1)[:,2], color='C4', scale=scale, units='xy', angles='xy', scale_units='xy' )
p1 = np.array([i.position for i in sites1[p1_all[:,2]==1.67]])
p2 = np.array([i.position for i in sites2[p1_all[:,2]==1.67]])
p1 = p1[p1[:,0]>-4]; p2 = p2[p2[:,0]>-4]; p1 = p1[p1[:,0]<7]; p2 = p2[p2[:,0]<7]
plt.plot(p1[:,0], p1[:,2], 'o', ms=10, color='C6', label='Bi')
plt.plot(p2[:,0], p2[:,2], 's', ms=5, color='C7', label='Bi, relax')
plt.quiver(p1[:,0], p1[:,2], (p2-p1)[:,0], (p2-p1)[:,2], color='C6', scale=scale, units='xy', angles='xy', scale_units='xy' )
p1 = np.array([i.position for i in sites1[p1_all[:,2]==3.7]])
p2 = np.array([i.position for i in sites2[p1_all[:,2]==3.7]])
p1 = p1[p1[:,0]>-4]; p2 = p2[p2[:,0]>-4]; p1 = p1[p1[:,0]<7]; p2 = p2[p2[:,0]<7]
plt.plot(p1[:,0], p1[:,2], 'o', ms=10, color='C8', label='Te')
plt.plot(p2[:,0], p2[:,2], 's', ms=5, color='C9', label='Te, relax')
plt.quiver(p1[:,0], p1[:,2], (p2-p1)[:,0], (p2-p1)[:,2], color='C8', scale=scale, units='xy', angles='xy', scale_units='xy' )
plt.axhline(dz[2]-2.76/2, color='grey', zorder=-100, ls='--')
plt.legend(ncol=1)
plt.ylabel('z ($\AA$)')
plt.xlabel('x ($\AA$)')
plt.axis('equal')
plt.xlim(-4,12)
# plt.savefig('relax_shift_xz.pdf')
plt.show()
scale = 0.3
p1 = np.array([i.position for i in sites1[p1_all[:,2]==-4.43]])
p2 = np.array([i.position for i in sites2[p1_all[:,2]==-4.43]])
plt.plot(p1[:,1], p1[:,2], 'ro', ms=10)
plt.plot(p2[:,1], p2[:,2], 'bs', ms=5)
plt.quiver(p1[:,1], p1[:,2], (p2-p1)[:,1], (p2-p1)[:,2], color='r', scale=scale, units='xy', angles='xy', scale_units='xy' )
p1 = np.array([i.position for i in sites1[p1_all[:,2]==zPd]])
p2 = np.array([i.position for i in sites2[p1_all[:,2]==zPd]])
plt.plot(p1[:,1], p1[:,2], 'o', ms=10, color='C0')
plt.plot(p2[:,1], p2[:,2], 's', ms=5, color='C1')
plt.quiver(p1[:,1], p1[:,2], (p2-p1)[:,1], (p2-p1)[:,2], color='C0', scale=scale, units='xy', angles='xy', scale_units='xy' )
if zPd2 is not None:
p1 = np.array([i.position for i in sites1[p1_all[:,2]==zPd2]])
p2 = np.array([i.position for i in sites2[p1_all[:,2]==zPd2]])
# p1 = np.array([i.position for i in sites1[p1_all[:,2]==zPd2]])+cell[0]+cell[1]
# p2 = np.array([i.position for i in sites2[p1_all[:,2]==zPd2]])+cell[0]+cell[1]
plt.plot(p1[:,1], p1[:,2], 'o', ms=10, color='C0', label='Pd')
plt.plot(p2[:,1], p2[:,2], 's', ms=5, color='C1', label='Pd, relax')
plt.quiver(p1[:,1], p1[:,2], (p2-p1)[:,1], (p2-p1)[:,2], color='C0', scale=scale, units='xy', angles='xy', scale_units='xy' )
p1 = np.array([i.position for i in sites1[p1_all[:,2]==-2.76]])
p2 = np.array([i.position for i in sites2[p1_all[:,2]==-2.76]])
plt.plot(p1[:,1], p1[:,2], 'o', ms=10, color='C2')
plt.plot(p2[:,1], p2[:,2], 's', ms=5, color='C3')
plt.quiver(p1[:,1], p1[:,2], (p2-p1)[:,1], (p2-p1)[:,2], color='C2', scale=scale, units='xy', angles='xy', scale_units='xy' )
plt.axhline(-2.76/2, color='grey', ls='--')
p1 = np.array([i.position for i in sites1[p1_all[:,2]==0]])
p2 = np.array([i.position for i in sites2[p1_all[:,2]==0]])
p1 = p1[p1[:,0]>-4]; p2 = p2[p2[:,0]>-4]; p1 = p1[p1[:,0]<7]; p2 = p2[p2[:,0]<7]
plt.plot(p1[:,1], p1[:,2], 'o', ms=10, color='C4', label='Te')
plt.plot(p2[:,1], p2[:,2], 's', ms=5, color='C5', label='Te, relax')
plt.quiver(p1[:,1], p1[:,2], (p2-p1)[:,1], (p2-p1)[:,2], color='C4', scale=scale, units='xy', angles='xy', scale_units='xy' )
plt.axis('equal')
plt.show()