Homec4science
Diffusion BosonicCodes (master) Browse

R10106/hm.py
9cf2d06018ebmaster

Last Change
Raw File

hm.py
View Options

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Mar 5 15:53:55 2020
@author: david
"""
import numpy as np
import qutip
def hamiltonian(delta, U, K, J, F, Nmodes, Nmax, verbose=True,
trunc = 1e-12):
"""
Construct the Hamiltonian using the Qutip Python Package given indicated
parameters. The model consists of Nmodes modes each with Nmax Fockstates.
Coupling, Kerr nonlinearity cross-Kerr nonlinearity and pumping are
specified below.
For the connectivity of the modes, see :func:`~adjacency_mat`
Parameters
----------
delta :ArithmeticError float array_like of shape (Nmodes)
On-site detuning for each mode
If float: constant for every site
U : float or array_like of shape (Nmodes)
On-site Kerr nonlinearity for each mode
If float: constant for every site
K : array_like of shape (Nmodes, Nmodes) or 0
Cross-Kerr nonlinearity between modes
Since K MUST be symmetric, it is only necessary to specify the
upper triangle of the cross-Kerr nonlinearity matrix
If 0 : No Cross-Kerr nonlinearty
J : array_like of shape (Nmodes, Nmodes)
hopping amplitudes between modes
Since J MUST be symmetric, it is only necessary to specify the
upper triangle of the cross-Kerr nonlinearity matrix
If 0 : No hopping
F : float or array_like of shape (Nmodes)
pumping amplitudes for single-photon pumping
If float: constant for every site
Nmodes : int
number of modes
Nmax : int
number of Fock states for each mode (including vaccuum state)
verbose : bool (optional, default: True)
If True: verbose mode on (for logging and debugging)
trunc : float or None (optional, default: 1e-12)
truncation of the matrix elements with absolute value smaller than
specified by trunc. This reduces the number of stored matrix elements
in the sparse matrix representation.
Returns
-------
Hamiltonian : qutip.Qobj
Returns Hamiltonian as Qutip Object with the properties
dims = [[Nmax]*Nmodes, [Nmax]*Nmodes]
shape = (Nmax*Nmodes, Nmax*Nmodes)
The matrix is represented as a complex csr sparse matrix
"""
vprint = print if verbose else lambda *a, **k: None
# handle cases in which parameters are simply one-dimensional
if isinstance(delta, (int, float, complex)):
delta = np.ones(Nmodes) * delta
if isinstance(U, (int, float, complex)):
U = np.ones(Nmodes) * U
if J == 0:
J = np.zeros((Nmodes, Nmodes))
if K == 0:
K = np.zeros((Nmodes, Nmodes))
if isinstance(F, (int, float, complex)):
F = np.ones(Nmodes) * F
vprint("constructing operators...")
a = qutip.destroy(Nmax) # define destruction operator for a single mode
n = qutip.num(Nmax) # define number operator for a single mode
# Define annihilation and creation operator for entire Hilbert space
# initialize empty lists for creat. and annih. & numb ops for each mode
aa, aa_dag, nn = [], [], []
for i in np.arange(Nmodes):
a_tensor = qutip.tensor([qutip.qeye(Nmax)] * i
+ [a]
+ [qutip.qeye(Nmax)] * (Nmodes - i - 1))
nn.append(qutip.tensor([qutip.qeye(Nmax)] * i
+ [n]
+ [qutip.qeye(Nmax)] * (Nmodes - i - 1)))
aa.append(a_tensor)
aa_dag.append(a_tensor.dag())
del a_tensor # free memory
vprint("done.")
# initialize empty Hamilton operatornp.linspace(0.1, 10, 100)]
h = qutip.Qobj(dims=[[Nmax] * Nmodes, [Nmax] * Nmodes])
for i in range(Nmodes):
# all local terms here
h += delta[i] * nn[i] # onsite-energies
h += F[i] * aa_dag[i] + np.conj(F[i]) * aa[i] # pump
h += U[i] * (aa_dag[i] * aa[i] * aa_dag[i] * aa[i])
for j in np.arange(i+1, Nmodes): # non-local terms, ensure that j > i
# all non-local terms here
if J[i,j] != 0: # avoid unneccesary cases with zero element
h += J[i,j] * ( aa_dag[i] * aa[j] + aa_dag[j] * aa[i])
if K[i,j] !=0: # avoid unneccesary cases with zero element
h += K[i,j] * nn[i] * nn[j]
return h.tidyup(atol = trunc)
def adjacency_mat(Lx, Ly, pb=False):
"""
return upper triangle matrix of 2d adjacencies to specify the connectivity
of the 2d lattice with or without periodic boundary conditions
index counting:
[[0, 1, 2, ... , Lx-1], [Lx, Lx+1,...,2*Lx-1], [...], [... ,...,Ly*Lx-1]
Parameters
----------
Lx, Ly : int
number of sites inx x and y directions, respectively
pb : bool (optional, default: False)
If True: use periodic boundary conditions
Returns
-------
Adjacency Matrix : numpy nd_array size (Lx*Ly, Lx*Ly)
Returns upper triangle of adjacencies
"""
I_a, I_b = np.eye(Lx), np.eye(Ly)
offdi_a = np.diag(np.ones(Lx-1), k=1)
offdi_b = np.diag(np.ones(Ly-1), k=1)
if pb:
offdi_a[0,-1], offdi_b[0,-1] = 1, 1
return np.kron(offdi_a, I_b) + np.kron(I_a, offdi_b)
c4science · Help