Plots and analysis for chapter 2: Approximations

In [1]:
# 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
In [2]:
# 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()

Reproduce old plots (from candidacy)

In [3]:
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
In [18]:
# 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')
In [21]:
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
In [20]:
grand_fixed_figure()
Out[20]:
In [7]:
# with temperature-broadened peaks to illustrate singlet-triplet splitting
fixed_bond_spectrum_figure()
Out[7]:
In [22]:
fixed_bond_spectrum_and_polarization_figure()
Out[22]:
In [9]:
fixed_bond_ising_spectrum_figure()
Out[9]:
In [10]:
# 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()
Out[10]:

Ground state polarization for bond and Ising

In [6]:
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
In [7]:
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()
In [13]:
groundstate_polarization_over_boa()
In [14]:
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.)

Explore interesting plots for various parameters (Ising, bond)

In [6]:
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
In [7]:
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
In [17]:
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))