# imports
import matplotlib.pyplot as plt
import numpy as np
from IPython.display import display
import qca
import coma
import configure_matplotlib
print(qca.__version__)
print(coma.__version__)
8.2.0 2.0.0-1.g4ad5ac5
# increases the resolution of the pngs displayed inline,
# so they are big enough even for a small figsize
%matplotlib inline
c = %config InlineBackend.rc
c['savefig.dpi'] = 150
%config InlineBackend.figure_format='png'
%config InlineBackend.rc = c
%matplotlib inline
configure_matplotlib.configure_matplotlib()
e1 = coma.Experiment('../experiments/experiment.000041/')
rs1 = e1.retrieve_results(
(
('Es','results/energies'),
),
(
('model', 'parameters/model'),
('V0', 'parameters/V0'),
('mu', 'parameters/mu'),
('N', 'parameters/layout/N'),
('V1', 'parameters/layout/V1'),
('boa', 'parameters/layout/boa'),
('q', 'parameters/q')
))
e2 = coma.Experiment('../experiments/experiment.000042/')
rs2 = e2.retrieve_results(
(
('T', 'parameters/T'),
('P', 'results/P'),
('N', 'results/N')
),
(
('model', 'parameters/model'),
('N', 'parameters/layout/N'),
('V1', 'parameters/layout/V1'),
('V0', 'parameters/V0')
))
rs1.sort(key=lambda i: (i.parameters.model,-i.parameters.N,-i.parameters.V1),reverse=True)
rs2.sort(key=lambda i: i.parameters.model, reverse=True)
#rs1,rs2
# function definitions
def dos(Es_, Emin=0, Emax=10, offset=0, deltaE=1E-4, Eoffset=None,gamma=None):
Es = Es_ - Es_[offset]
if Eoffset is not None:
Es += Eoffset
Es = Es[(Es >= Emin) & (Es <= Emax)]
# construct a matrix with zero everywhere in the range Emin to Emax,
# and only a few delta peaks where the energy eigenvalues are located
xs = np.arange(Emin,Emax,deltaE)
m = np.row_stack((
np.column_stack((xs,np.zeros(xs.shape))),
np.column_stack((Es,np.ones(Es.shape)))
))
m = m[m[:,0].argsort()]
# simply convolve with a Lorentzian: this gives us the density of states
# the resolution deltaE must be much better than the peak width gamma
lorentz = lambda x,x0,gamma: gamma*gamma / ((x-x0)*(x-x0) + gamma*gamma)
x0 = Emin+(Emax-Emin)/2.0
if not gamma:
gamma = (Emax-Emin)/1000.0
convolved = np.convolve(m[:,1], lorentz(m[:,0],x0,gamma),'same')
m = np.column_stack((m[:,0],convolved))
m = m[m[:,1] != 0]
return m
def dosplot(p, dos1, dos2, label1, label2, xlim=None, ylim=None, color1='blue', color2='green'):
p.plot(dos1[:,0],dos1[:,1],label=label1,color=color1)
p.plot(dos2[:,0],-dos2[:,1],label=label2,color=color2)
p.axhline(color='#666666')
p.set_xlabel('energy')
p.set_ylabel('state count')
if xlim: p.set_xlim(xlim[:2])
if ylim:
p.set_ylim(ylim[:2])
ylim = list(ylim)
ylim[1] += 1
p.set_yticks(range(*ylim))
p.set_yticklabels([abs(i) for i in range(*ylim)])
models = {'QcaBond': 'Bond','QcaFixedCharge': 'Fixed','QcaGrandCanonical': 'Grand','QcaIsing':'Ising'}
def PTplot(p,rs):
for r in rs:
d = r.table
p.plot(d[:,0], d[:,2], label=models[r.parameters.model])
p.set_xlabel('temperature')
p.set_ylabel('polarization $P_2$')
def NTplot(p,rs):
d = rs['model','QcaGrandCanonical'][0].table
p.plot(d[:,0], d[:,7],label='Grand. Cell 1.',color='blue')
p.plot(d[:,0], d[:,12],label='Grand. Cell 2.',color='red')
d = rs['model','QcaFixedCharge'][0].table
p.plot(d[:,0], d[:,7].round(2),label='Fixed. Both cells.',color='green')
p.set_xlabel('temperature')
p.set_ylabel('particle number')
def text(p,t,loc='upper right',dx=0.035,dy=0.035):
if loc=='upper right': x,y,v,h = 1-dx,1-dy,'top','right'
elif loc=='upper left': x,y,v,h = dx,1-dy,'top','left'
elif loc=='lower right': x,y,v,h = 1-dx,dy,'bottom','right'
elif loc=='lower left': x,y,v,h = dx,dy,'bottom','left'
p.text(x,y,t,verticalalignment=v,horizontalalignment=h,transform=p.transAxes)
def grand_fixed_figure():
T = 0.05
Emin,Emax,deltaE = (-2,50,1E-3)
E_f = rs1['mu',250]['model','QcaFixedCharge'][0].table[0,0]
E_g = rs1['mu',250]['model','QcaGrandCanonical'][0].table[0,0]
dos_f = dos(E_f,Emin,Emax,deltaE=deltaE,gamma=0.5*T)
dos_g = dos(E_g,Emin,Emax,deltaE=deltaE,gamma=0.5*T)
fig,ps = plt.subplots(1,2,figsize=(5,2.5))
dosplot(ps[0], dos_g,dos_f,'Grand','Fixed',(Emin,Emax),(-20,30,5))
NTplot(ps[1], rs2['N',2]['V0',1E3])
ps[0].set_yticks(np.arange(-20,31,10))
ps[0].set_yticklabels(abs(np.arange(-20,31,10)))
ps[0].legend(loc='upper left')
ps[1].set_yticks(np.arange(1.5,2.51,0.25))
ps[1].set_ylim((1.5,2.5))
ps[1].legend(loc='lower right')
for p,t in zip(ps,['(a)','(b)']):
text(p,t,'lower left')
fig.tight_layout()
fig.savefig('../plots/chapter02/fixed_charge_approximation.pdf')
plt.close()
return fig
def fixed_bond_spectrum_figure():
fig,ps = plt.subplots(1,2,figsize=(5,2.5))
Ts = [0.05,0.5]
ls = ['(a)','(b)']
for p,T,l in zip(ps,Ts,ls):
Emin,Emax,deltaE = (-1,14,1E-3)
E_f = rs1['mu',0]['N',1]['V1',20]['model','QcaFixedCharge'][0].table[0,0]
E_b = rs1['mu',0]['N',1]['V1',20]['model','QcaBond'][0].table[0,0]
dos_f = dos(E_f,Emin,Emax,deltaE=deltaE,gamma=0.5*T)
dos_b = dos(E_b,Emin,Emax,deltaE=deltaE,Eoffset=E_f[1]-E_f[0],gamma=0.5*T)
dosplot(p, dos_f,dos_b,'Fixed','Bond',(Emin,Emax),(-2,5))
text(p,'$T = {}$'.format(T),'upper right')
text(p,l,'lower left')
ps[0].legend(loc='upper left')
fig.tight_layout()
fig.savefig('../plots/chapter02/bond_approximation1.pdf')
plt.close()
return fig
def fixed_bond_spectrum_and_polarization_figure():
fig,ps = plt.subplots(2,2,figsize=(5,5))
V1s = [20,100]
for p,V1 in zip(ps[:,0],[20,100]):
T = 0.005
Emin,Emax,deltaE = (-0.5,4,1E-4)
E_f = rs1['mu',0]['N',2]['V1',V1]['model','QcaFixedCharge'][0].table[0,0]
E_b = rs1['mu',0]['N',2]['V1',V1]['model','QcaBond'][0].table[0,0]
dos_f = dos(E_f,Emin,Emax,deltaE=deltaE,gamma=0.5*T)
dos_b = dos(E_b,Emin,Emax,deltaE=deltaE,Eoffset=E_b[0]-E_f[0],gamma=0.5*T)
dosplot(p,dos_f,dos_b,'Fixed','Bond',(Emin,Emax),(-2,10,2))
p.set_xticks(range(5))
for p,V1 in zip(ps[:,1],V1s):
PTplot(p, rs2['N',2]['V1',V1]['V0',1E6])
p.set_xticks(range(5))
ps[1,1].set_yticks([0.0,0.2,0.4,0.6,0.8])
ps[0,1].legend(loc='upper right')
ls = ['(a)','(b)','(c)','(d)']
for p,l,V1,dy in zip(ps.reshape(-1).tolist(),ls,
[20,20,100,100],[0.035,0.27,0.035,0.035]):
text(p,l,'lower left')
text(p,'$V_1 = {}$'.format(V1),'upper right',dy=dy)
fig.tight_layout()
fig.savefig('../plots/chapter02/bond_approximation2.pdf')
plt.close()
return fig
def fixed_bond_ising_spectrum_figure():
fig,ps = plt.subplots(1,2,figsize=(5,2.5))
ls = ['(a)','(b)']
ms = [('QcaFixedCharge','QcaBond'),('QcaFixedCharge','QcaIsing')]
cs2 = ['green','red']
for p,(m1,m2),l,c2 in zip(ps,ms,ls,cs2):
T = 0.002
Emin,Emax,deltaE = (-0.05,0.4,1E-4)
E_f = rs1['mu',0]['q',0.5]['N',1]['V1',100]['model',m1][0].table[0,0]
E_b = rs1['mu',0]['q',0.5]['N',1]['V1',100]['model',m2][0].table[0,0]
dos_f = dos(E_f,Emin,Emax,deltaE=deltaE,gamma=0.5*T)
Eoffset = E_f[1]-E_f[0] if m2=='QcaBond' else 0
dos_b = dos(E_b,Emin,Emax,deltaE=deltaE,Eoffset=Eoffset,gamma=0.5*T)
dosplot(p, dos_f,dos_b,models[m1],models[m2],(Emin,Emax),(-2,5),'blue',c2)
text(p,l,'lower left')
p.legend(loc='upper right')
p.set_xticklabels(['',0,'',0.1,'',0.2,'',0.3,'',0.4])
fig.tight_layout()
#fig.savefig('../plots/chapter02/ising_approximation.pdf')
plt.close()
return fig
def fixed_ising_spectrum_figure():
fig,ps = plt.subplots(10,2,figsize=(5,25))
ms = [('QcaFixedCharge','QcaBond'),('QcaFixedCharge','QcaIsing')]
cs2 = ['green','red']
params = [(boa,V1) for V1 in [100,200] for boa in [1.2,2,3,4,5]]
for (boa,V1),ps_ in zip(params,ps):
for p,(m1,m2),c2 in zip(ps_,ms,cs2):
#T = 0.2 # for boa=1.2
#Emin,Emax,deltaE = (-0.05,50,1E-2)
T = 0.02
Emin,Emax,deltaE = (-0.05,10,1E-3)
E_f = rs1['mu',0]['q',0.5]['N',2]['boa',boa]['V1',V1]['model',m1][0].table[0,0]
E_b = rs1['mu',0]['q',0.5]['N',2]['boa',boa]['V1',V1]['model',m2][0].table[0,0]
dos_f = dos(E_f,Emin,Emax,deltaE=deltaE,gamma=0.5*T)
Eoffset = E_f[7]-E_f[0] if m2=='QcaBond' else 0
dos_b = dos(E_b,Emin,Emax,deltaE=deltaE,Eoffset=Eoffset,gamma=0.5*T)
dosplot(p, dos_f,dos_b,models[m1],models[m2],(Emin,Emax),(-2,5),'blue',c2)
text(p,'$b/a = {}$\n$V_1 = {}$'.format(boa,V1),'lower right')
p.legend(loc='upper right')
fig.tight_layout()
plt.close()
return fig
grand_fixed_figure()
# with temperature-broadened peaks to illustrate singlet-triplet splitting
fixed_bond_spectrum_figure()
fixed_bond_spectrum_and_polarization_figure()
fixed_bond_ising_spectrum_figure()
# Huge number of spectrum plots for various parameters
# TODO: the problem is that with these parameters (and P_D = 1) everything will be fully polarized
fixed_ising_spectrum_figure()
def P_gs(N,V1,boa,P_D=1):
ss = [qca.QcaFixedCharge(),qca.QcaBond(),qca.QcaIsing()]
Ps = []
for s in ss:
s.q = 0.5
s.t = 1
s.V0 = 1E6
s.mu = 0
s.T = 1E-9
s.l = qca.Wire(N,V1,boa,P_D)
s.init()
s.run()
Ps.append(s.results['P'][-1])
return Ps
def groundstate_polarization_over_boa():
fig,ps = plt.subplots(2,3,figsize=(7.5,4))
Ps = [0.1,1]
V1s = [20,100,200]
for ps_,P in zip(ps,Ps):
for p,V1 in zip(ps_,V1s):
boas = np.linspace(2,10,20)
Ps = np.array([P_gs(2,V1,boa,P) for boa in boas])
m = np.column_stack((boas,Ps))
p.plot(m[:,0],m[:,1],label='Fixed')
p.plot(m[:,0],m[:,2],label='Bond')
p.plot(m[:,0],m[:,3],label='Ising')
text(p,'$P_D = {}$\n$V_1 = {}$'.format(P,V1),'lower right')
#p.set_ylim((0,1.1))
p.set_ylabel('$P_2$')
p.set_xlabel('$b/a$')
ps[0,2].legend()
fig.tight_layout()
def groundstate_polarization_over_V1():
fig,ps = plt.subplots(2,3,figsize=(7.5,4))
Ps = [0.1,1]
boas = [2,3,4]
for ps_,P in zip(ps,Ps):
for p,boa in zip(ps_,boas):
V1s = np.linspace(20,200,20)
Ps = np.array([P_gs(2,V1,boa,P) for V1 in V1s])
m = np.column_stack((V1s,Ps))
p.plot(m[:,0],m[:,1],label='Fixed')
p.plot(m[:,0],m[:,2],label='Bond')
p.plot(m[:,0],m[:,3],label='Ising')
text(p,'$P_D = {}$\n$b/a = {}$'.format(P,boa),'lower right')
#p.set_ylim((0,1.1))
p.set_ylabel('$P_2$')
p.set_xlabel('$V_1$')
ps[0,2].legend()
fig.tight_layout()
groundstate_polarization_over_boa()
groundstate_polarization_over_V1()
I would say that the bottom line from these two plots is: Bond and Ising just do not reproduce the exact ground state at all, they are not even close. (Except when everything is fully polarized.)
def Es(T,N,V1,boa,P_D=1,V0=1E6):
ss = [qca.QcaFixedCharge(),qca.QcaBond(),qca.QcaIsing()]
Es = []
for s in ss:
s.q = 0.5
s.t = 1
s.V0 = V0
s.mu = 0.5
s.T = T
s.l = qca.Wire(N,V1,boa,P_D)
s.init()
s.run()
E = np.array(s.energies())
E.sort()
E = E[:10000]
Es.append(E)
return np.array(Es)
def P_all(T,N,V1,boa,P_D=1,V0=1E6,q=0.5):
ss = [qca.QcaFixedCharge(),qca.QcaBond(),qca.QcaIsing()]
Ps = []
for s in ss:
s.q = q
s.t = 1
s.V0 = V0
s.mu = 0.5
s.T = T
s.l = qca.Wire(N,V1,boa,P_D)
s.init()
s.run()
Ps.extend(s.results['P'])
return Ps
def relative_error(Ps):
N = Ps.shape[1] / 3
m = np.column_stack([
abs(Ps[:,0:N] - Ps[:,i*N:(i+1)*N]) / Ps[:,0:N]
for i in [1,2] ])
return m
def PETplot(T_,V1,boa,P_D):
fig,ps = plt.subplots(2,2,figsize=(2*2.5,2*2.5))
V0 = 1E3
Ts = np.logspace(-3,3,20)
Ps = np.array([P_all(T,2,V1,boa,P_D,V0) for T in Ts])
es = relative_error(Ps)
m_P = np.column_stack((Ts,Ps))
m_e = np.column_stack((Ts,es))
ls = ['Bond','Ising']
for p,o,l in zip(ps[0],[3,5],ls):
p.plot(m_P[:,0],m_P[:,1],label='Fixed 1')
p.plot(m_P[:,0],m_P[:,2],label='Fixed 2')
p.plot(m_P[:,0],m_P[:,o],label=l + ' 1')
p.plot(m_P[:,0],m_P[:,o+1],label=l + ' 2')
p.set_xscale('log')
p.set_xlabel('$T$')
p.set_ylabel('$P_2$')
p.set_ylim((0,1))
p.legend()
for p,o in zip(ps[1],[1,3]):
p.plot(m_e[:,0],m_e[:,o],label='Cell 1')
p.plot(m_e[:,0],m_e[:,o+1],label='Cell 2')
p.set_xscale('log')
p.set_xlabel('$T$')
p.set_ylabel('relative error')
p.set_ylim((0,1))
p.legend()
fig.tight_layout()
plt.close()
return fig
def PEV1plot(T,V1_,boa,P_D):
fig,ps = plt.subplots(2,2,figsize=(2*2.5,2*2.5))
V1s = np.logspace(1,3,20)
Ps = np.array([P_all(T,2,V1,boa,P_D) for V1 in V1s])
es = relative_error(Ps)
m_P = np.column_stack((V1s,Ps))
m_e = np.column_stack((V1s,es))
ls = ['Bond','Ising']
for p,o,l in zip(ps[0],[3,5],ls):
p.plot(m_P[:,0],m_P[:,1],label='Fixed 1')
p.plot(m_P[:,0],m_P[:,2],label='Fixed 2')
p.plot(m_P[:,0],m_P[:,o],label=l + ' 1')
p.plot(m_P[:,0],m_P[:,o+1],label=l + ' 2')
p.set_xscale('log')
p.set_xlabel('$V_1$')
p.set_ylabel('$P_2$')
p.set_ylim((0,1))
p.legend()
for p,o in zip(ps[1],[1,3]):
p.plot(m_e[:,0],m_e[:,o],label='Cell 1')
p.plot(m_e[:,0],m_e[:,o+1],label='Cell 2')
p.set_xscale('log')
p.set_yscale('log')
p.set_xlabel('$V_1$')
p.set_ylabel('relative error')
p.set_ylim((0,1))
p.legend()
fig.tight_layout()
plt.close()
return fig
def PEBOAplot(T,V1,boa_,P_D):
fig,ps = plt.subplots(2,2,figsize=(2*2.5,2*2.5))
boas = np.linspace(1.2,10,20)
Ps = np.array([P_all(T,2,V1,boa,P_D) for boa in boas])
es = relative_error(Ps)
m_P = np.column_stack((boas,Ps))
m_e = np.column_stack((boas,es))
ls = ['Bond','Ising']
for p,o,l in zip(ps[0],[3,5],ls):
p.plot(m_P[:,0],m_P[:,1],label='Fixed 1')
p.plot(m_P[:,0],m_P[:,2],label='Fixed 2')
p.plot(m_P[:,0],m_P[:,o],label=l + ' 1')
p.plot(m_P[:,0],m_P[:,o+1],label=l + ' 2')
p.set_xlabel('$b/a$')
p.set_ylabel('$P_2$')
p.set_ylim((0,1))
p.legend()
for p,o in zip(ps[1],[1,3]):
p.plot(m_e[:,0],m_e[:,o],label='Cell 1')
p.plot(m_e[:,0],m_e[:,o+1],label='Cell 2')
p.set_xlabel('$b/a$')
p.set_ylabel('relative error')
#p.set_ylim((0,1))
p.legend()
fig.tight_layout()
plt.close()
return fig
def PEPDplot(T,V1,boa,P_D_):
fig,ps = plt.subplots(2,2,figsize=(2*2.5,2*2.5))
PDs = np.linspace(0.01,1,10)
Ps = np.array([P_all(T,2,V1,boa,PD) for PD in PDs])
es = relative_error(Ps)
m_P = np.column_stack((PDs,Ps))
m_e = np.column_stack((PDs,es))
ls = ['Bond','Ising']
for p,o,l in zip(ps[0],[3,5],ls):
p.plot(m_P[:,0],m_P[:,1],label='Fixed 1')
p.plot(m_P[:,0],m_P[:,2],label='Fixed 2')
p.plot(m_P[:,0],m_P[:,o],label=l + ' 1')
p.plot(m_P[:,0],m_P[:,o+1],label=l + ' 2')
p.set_xlabel('$P_D$')
p.set_ylabel('$P_2$')
p.set_ylim((0,1))
p.legend()
for p,o in zip(ps[1],[1,3]):
p.plot(m_e[:,0],m_e[:,o],label='Cell 1')
p.plot(m_e[:,0],m_e[:,o+1],label='Cell 2')
p.set_xlabel('$P_D$')
p.set_ylabel('relative error')
#p.set_ylim((0,1))
p.legend()
fig.tight_layout()
plt.close()
return fig
def spectrumplot(T,V1,boa,P_D):
fig,ps = plt.subplots(1,2,figsize=(5,2.5))
es = Es(T,2,V1,boa,P_D)
cs2 = ['green','red']
T_s = 0.02
Emin,Emax,deltaE = (-0.5,4,1E-3)
for p,(E1,E2),m2,c2 in zip(ps,[(es[0],es[1]),(es[0],es[2])],['Bond','Ising'],cs2):
dos1 = dos(E1,Emin,Emax,deltaE=deltaE,gamma=0.5*T_s)
Eoffset = E1[7]-E1[0] if m2=='Bond' else 0
dos2 = dos(E2,Emin,Emax,deltaE=deltaE,Eoffset=Eoffset,gamma=0.5*T_s)
dosplot(p, dos1,dos2,'Fixed',m2,(Emin,Emax),(-2,5),'blue',c2)
#text(p,'$b/a = {}$\n$V_1 = {}$'.format(boa,V1),'lower right')
p.legend(loc='upper right')
fig.tight_layout()
plt.close()
return fig
T,V1,boa,P_D = 1,100,3.0,0.1
display(spectrumplot(T,V1,boa,P_D))
display(PETplot(T,V1,boa,P_D))
display(PEV1plot(T,V1,boa,P_D))
display(PEBOAplot(T,V1,boa,P_D))
display(PEPDplot(T,V1,boa,P_D))
T,V1,boa,P_D = 1,100,2.0,0.1
# display(spectrumplot(T,V1,boa,P_D))
# display(PETplot(T,V1,boa,P_D))
# display(PEV1plot(T,V1,boa,P_D))
# display(PEBOAplot(T,V1,boa,P_D))
# display(PEPDplot(T,V1,boa,P_D))
T,V1,boa,P_D = 1,200,3.0,0.1
display(spectrumplot(T,V1,boa,P_D))
display(PETplot(T,V1,boa,P_D))
display(PEV1plot(T,V1,boa,P_D))
display(PEBOAplot(T,V1,boa,P_D))
display(PEPDplot(T,V1,boa,P_D))
T,V1,boa,P_D = 1,40,2.0,0.5
display(spectrumplot(T,V1,boa,P_D))
display(PETplot(T,V1,boa,P_D))
#display(PEV1plot(T,V1,boa,P_D))
#display(PEBOAplot(T,V1,boa,P_D))
#display(PEPDplot(T,V1,boa,P_D))
# Have to fiddle with the resolution, but you can see, in the spectrum,
# the singlet ground state and then two singlet-triplet triplets and
# then the triplet-triplet states
T,V1,boa,P_D = 1,60,2.0,0.5
display(spectrumplot(T,V1,boa,P_D))
display(PETplot(T,V1,boa,P_D))
#display(PEV1plot(T,V1,boa,P_D))
#display(PEBOAplot(T,V1,boa,P_D))
#display(PEPDplot(T,V1,boa,P_D))
def P_q_plot(T_,V1,boa,P_D,q):
fig,ps = plt.subplots(2,2,figsize=(2*2.5,2*2.5))
V0 = 1E3
Ts = np.logspace(-3,3,20)
Ps = np.array([P_all(T,2,V1,boa,P_D,V0,q=q) for T in Ts])
es = relative_error(Ps)
m_P = np.column_stack((Ts,Ps))
m_e = np.column_stack((Ts,es))
ls = ['Bond','Ising']
for p,o,l in zip(ps[0],[3,5],ls):
p.plot(m_P[:,0],m_P[:,1],label='Fixed 1')
p.plot(m_P[:,0],m_P[:,2],label='Fixed 2')
p.plot(m_P[:,0],m_P[:,o],label=l + ' 1')
p.plot(m_P[:,0],m_P[:,o+1],label=l + ' 2')
p.set_xscale('log')
p.set_xlabel('$T$')
p.set_ylabel('$P_2$')
p.set_ylim((0,1))
p.legend()
for p,o in zip(ps[1],[1,3]):
p.plot(m_e[:,0],m_e[:,o],label='Cell 1')
p.plot(m_e[:,0],m_e[:,o+1],label='Cell 2')
p.set_xscale('log')
p.set_xlabel('$T$')
p.set_ylabel('relative error')
p.set_ylim((0,1))
p.legend()
fig.tight_layout()
plt.close()
return fig
P_q_plot(1,100,2,0.1,0.5)
P_q_plot(1,100,2,0.1,0.0)
I just wanted to briefly look at the Ising model with $q=0$ as I never did that before. Unsurprisingly, the Ising model does not work that well for non-charge neutral cells: Especially, the error of the last cell is much larger. This is because the Ising model only has two states per cell and in particular cannot represent "edge states" which however become more important when charge buildup is present. Likewise, for larger inter-cell distances the error of the Ising model (for $q = 0$) becomes smaller, as charge buildup plays a smaller role.
For the straight line of cells ($\theta = 0$), the Ising model should yield the same polarization, regardless of the compensation charge $q$ (because $J_{kl}$ does not depend on $q$ and $J^{\prime}_{kl}$ vanishes for $\theta = 0$). We briefly check that below.
T,V1,boa,P_D,V0 = 1,100,2,0.1,1E3
Ts = np.logspace(-3,3,10)
Ps_1 = np.array([P_all(T,2,V1,boa,P_D,V0,q=0) for T in Ts])
Ps_2 = np.array([P_all(T,2,V1,boa,P_D,V0,q=0.5) for T in Ts])
delta_f = abs(Ps_1[:,0:2] - Ps_2[:,0:2])
delta_b = abs(Ps_1[:,2:4] - Ps_2[:,2:4])
delta_i = abs(Ps_1[:,-2:] - Ps_2[:,-2:])
print("delta for fixed charge model: {}".format(sum(delta_f)))
print("delta for bond model: {}".format(sum(delta_b)))
print("delta for Ising model: {}".format(sum(delta_i)))
delta for fixed charge model: [ 0.1985545 0.91361819] delta for bond model: [ 0.00364888 0.68695917] delta for Ising model: [ 3.18664013e-13 3.21963196e-13]
def ising_one_cell_figure():
fig,ps = plt.subplots(1,2, figsize=(5,2.5))
T,N,V1,boa,PD = 1,1,100,3,0.1
# left plot
T_ = 0.005
es = Es(T,N,V1,boa)
ds = [dos(es_,-0.1,0.8,deltaE=1E-4,Eoffset=o,gamma=0.5*T_) for es_,o in zip(es,[0,es[1][0]-es[0][0],0])]
p = ps[0]
p.plot(ds[0][:,0],ds[0][:,1],label='Fixed')
p.plot(ds[1][:,0],-ds[1][:,1],label='Bond')
p.plot(ds[2][:,0],-ds[2][:,1],label='Ising')
p.axhline(color='#666666')
p.set_xlabel('energy')
p.set_ylabel('state count')
p.set_ylim(-2,4)
p.set_yticks([-2,-1,0,1,2,3,4])
p.set_yticklabels([2,1,0,1,2,3,4])
p.set_xticks(np.arange(0,0.81,0.2))
text(p,'(a)','lower left')
# right plot
Ts = np.logspace(-3,3,100)
Ps = np.array([P_all(T,N,V1,boa,PD) for T in Ts])
m = np.column_stack((Ts,Ps))
p = ps[1]
p.plot(m[:,0],m[:,1],label='Fixed')
p.plot(m[:,0],m[:,2],label='Bond')
p.plot(m[:,0],m[:,3],label='Ising')
p.set_xscale('log')
p.set_xlabel('temperature')
p.set_ylabel('polarization')
text(p,'(b)','lower left')
p.legend()
fig.tight_layout()
fig.savefig('../plots/chapter02/ising_approximation1.pdf')
plt.close()
return fig
display(ising_one_cell_figure())
def ising_two_cell_figure():
fig,ps = plt.subplots(2,2,figsize=(2*2.5,2*2.5))
for ps_,V1_ in zip(ps,[100,200]):
T,N,V1,boa,PD,V0 = 1,2,V1_,3,0.1,1E3
# left plot
T_ = 0.01
Emin,Emax = -0.25,1.6
es = Es(T,N,V1,boa,PD,V0)
ds = [dos(es_,Emin,Emax,deltaE=1E-4,Eoffset=o,gamma=0.5*T_)
for es_,o in zip(es,[0,es[0][7]-es[0][0],{100:-0.15,200:-0.08}[V1_]])]
p = ps_[0]
p.plot(ds[0][:,0],ds[0][:,1],label='Fixed')
p.plot(ds[1][:,0],-ds[1][:,1],label='Bond')
p.plot(ds[2][:,0],-ds[2][:,1],label='Ising')
p.axhline(color='#666666')
p.set_xlabel('energy')
p.set_ylabel('state count')
p.set_xlim((Emin,Emax))
p.set_ylim((-3,15))
ys = np.arange(-3,15.1,2)
p.set_yticks(ys)
p.set_yticklabels(abs(ys))
p.set_xticks(np.arange(0,1.61,0.5))
# right plot
Ts = np.logspace(-3,3,80)
Ps = np.array([P_all(T,N,V1,boa,PD,V0) for T in Ts])
es = relative_error(Ps)
m_e = np.column_stack((Ts,es))
p = ps_[1]
p.plot(m_e[:,0],m_e[:,2],label='Bond',color='green')
p.plot(m_e[:,0],m_e[:,4],label='Ising',color='red')
p.set_xscale('log')
p.set_xlabel('temperature')
p.set_ylabel('relative error for $P_2$')
p.set_ylim((0,1))
#p.legend()
handles,labels = ps[0,0].get_legend_handles_labels()
ps[1,1].legend(handles, labels, loc='upper left')
ls = ['(a)','(b)','(c)','(d)']
for p,l,V1 in zip(ps.reshape(-1).tolist(),ls,[100,100,200,200]):
text(p,l,'lower left')
text(p,'$V_1 = {}$'.format(V1),'upper right')
fig.tight_layout()
fig.savefig('../plots/chapter02/ising_approximation2.pdf')
plt.close()
return fig
ising_two_cell_figure()
# See notebook 7
e = coma.Experiment('../experiments/experiment.000031/')
rs = e.retrieve_results(
(
('P','results/P'),
),
(
('model','parameters/model'),
('N','parameters/layout/N'),
('V1','parameters/layout/V1'),
('boa','parameters/layout/boa'),
('P_D','parameters/layout/P'),
('T','parameters/T')
)
)
rs.sort(key=lambda r: r.parameters.T)
from copy import copy
def relative_error_over_wire_length(p,rs,V1,boa,P_D):
for T in [1.0,2.0]:
rows = [
np.concatenate((np.array([r1.parameters.N]),
r1.table[0][0:3],
r2.table[0][0:3],
abs(r1.table[0][0:3]-r2.table[0][0:3]) / r1.table[0][0:3]))
for r1,r2 in zip(
rs['V1',V1]['boa',boa]['T',T]['P_D',P_D]['model','QcaBond'],
rs['V1',V1]['boa',boa]['T',T]['P_D',P_D]['model','QcaIsing']
)
]
d = np.vstack(rows)
for i,c in zip(range(3),['blue','green','red']):
p.plot(d[:,0],d[:,i+7], '+-', label='Cell {}'.format(i+1), color=c)
p.set_xlabel('number of cells in wire')
p.set_ylabel('relative error')
p.set_xlim((2.9,5.1))
p.set_ylim((0.0,0.2))
p.text(4.5,0.123,'$T = 1.0$')
p.text(4.5,0.03,'$T = 2.0$')
def relative_error_over_cells(p,rs,V1,boa,P_D):
ls = sorted(rs['N',5], key=lambda r: (r.parameters.V1,r.parameters.boa,r.parameters.P_D,r.parameters.T))
ls = [ls[i:i+2] for i in range(0,len(ls)-2,2)]
rs_ = []
for l1,l2 in ls:
r = copy(l1)
r.table = abs(l1.table - l2.table) / l1.table
rs_.append(r)
rs_ = coma.ResultList(rs_)
for T in [1.0,2.0]:
r = rs_['V1',V1]['boa',boa]['T',T]['P_D',P_D][0]
p.plot(np.arange(1,6),r.table[0],'+-',label='T = {}'.format(T))
p.set_xlim((0,6))
p.set_xlabel('cell number')
p.set_ylabel('relative error')
p.set_ylim((0.0,0.2))
def ising_wire_figure():
V1,boa,P_D = 100,3.0,0.1
fig, ps = plt.subplots(1,2,figsize=(5,2.5))
p = ps[0]
relative_error_over_wire_length(p,rs,V1,boa,P_D)
p.set_xticks([3,4,5])
text(p,'(a)','lower left')
hs,ls = p.get_legend_handles_labels()
p.legend(hs[0:3],ls[0:3],loc='upper left')
p = ps[1]
relative_error_over_cells(p,rs,V1,boa,P_D)
text(p,'(b)','lower left')
p.legend(loc='upper left')
fig.tight_layout()
fig.savefig('../plots/chapter02/ising_approximation3.pdf')
plt.close()
return fig
display(ising_wire_figure())