假设一个正方晶格的二维光子晶体，晶格常数为0.5微米，有四个圆形介质柱，直径都为0.1微米，介质柱中心距离晶格中心的距离是0.18微米，用python计算晶体哈密顿量得到能带图

在二维光子晶体中，计算哈密顿量并得到能带图需要一些专业的光学和量子传输理论知识，以及使用专门针对这种问题的数值模拟库。在Python中，常用的是如photonic crystals库，尽管它不是直接提供哈密顿量计算功能，但我们可以通过建立适当的数学模型，然后使用数值求解方法来近似。

首先，我们需要了解光子晶体的基本原理，包括Bloch定理和传输矩阵方法。对于这个问题，我们可以假设每个圆柱形区域是一个“缺陷”或截止区域，周围的空气区域代表光子的自由传播空间。

下面是一个简化版的流程，实际计算会更复杂，需要考虑边界条件、周期性和离散化处理：

1. 定义物理参数：晶格常数a=0.5微米，圆柱直径d=0.1微米，中心间距p=0.18微米。
2. 离散化：将晶格分割成小单元格，每一步的大小会影响哈密顿ian的精确度。
3. 创建模式：在每个空隙（即圆柱之间的部分）建立自由传播模式，在圆柱内部则可能是零或截止模式。
4. 构建哈密顿ian：基于传递矩阵（Transfer Matrix Method），构造一个二维的Hermitean矩阵，反映光子在晶格中的传播情况。
5. 解线性方程组：求解这个哈密顿ian的特征值问题，得到每一列波矢下的能带。

以下是Python的一个简单示例框架（请注意这只是一个非常粗略的演示，实际代码会更复杂）：

import numpy as np
from scipy.sparse.linalg import eigsh

# 基本参数
a = 5e-6  # 晶格常数
d = 1e-6  # 圆柱直径
p = 9e-6  # 中心间距

# 离散化，这里以步长为1纳米
cell_size = 1e-9
n_cells = int((a - p - d) / cell_size + 1)

# 初始化哈密顿ian矩阵
hamiltonian = np.zeros((n_cells, n_cells))

# 填充矩阵元素（这里仅作示例，需根据具体模型填充）
for i in range(n_cells):
    for j in range(i, n_cells):
        if abs(i-j) <= (a-p-d)/cell_size:
            hamiltonian[i,j] = ...  # 对应于圆柱区域的阻抗差异
            hamiltonian[j,i] = hamiltonian[i,j].conjugate()

# 能带计算
eigenvalues, eigenvectors = eigsh(hamiltonian, k=n_cells//2+1)  # 获取前半部分能带
dispersion关系 = np.sqrt(eigenvalues) * cell_size  # 将频率转换为波矢

# 可视化能带图
plt.plot(dispersion关系, 'r-')
plt.xlabel('Wave vector along the crystal [1/nm]')
plt.ylabel('Energy [THz]')
plt.show()