开边界格点系统对角化程序

最新推荐文章于 2022-12-23 20:30:53 发布

wwxy261

最新推荐文章于 2022-12-23 20:30:53 发布

阅读量187

点赞数

分类专栏：算法

本文链接：https://blog.csdn.net/wwxy1995/article/details/117600652

版权

算法专栏收录该内容

3633 篇文章 119 订阅

订阅专栏

首先，对于周期性系统，我们最好先解析推导出Hk的形式

其中Hk往往可以写成分块对角形式，因此我们可以采用如下代码，直接构建出稠密矩阵的哈密顿量

石墨烯的例子

Eigen::MatrixXcd Hamiltonian_Graphene(double k, int N, double t = -1.0) {
    Eigen::Matrix2cd T = Eigen::Matrix2cd::Zero(2,2);
    Eigen::Matrix2cd Dup = Eigen::Matrix2cd::Zero(2,2);
    T(0,1) = 1.0 + std::exp(-1i*k);
    T(1,0) = 1.0 + std::exp(-1i*k);
    Dup(1,0) = 1.0 + 0.0 * 1i;
    Eigen::Matrix2cd Ddown = Dup.transpose();
    Eigen::MatrixXcd H = Eigen::MatrixXcd::Zero(2*N,2*N);
    H.block<2,2>(0,0) = T;
    for(int i=2;i<2*N;i+=2){
        H.block<2,2>(i,i) = T;
        H.block<2,2>(i-2,i) = Dup;
        H.block<2,2>(i,i-2) = Ddown;
    }
    return t*H;
}

当然这一段代码有可能可以通过向量化矩阵操作写的更优雅，这里简单先这样写。

对于一维问题，给定路径对角化是十分容易的，通用的代码逻辑这样写：

void solve_H0_band(std::vector<double>& ks, int N) {
    std::ofstream out("H0_band.txt");
    for(auto k : ks) {
        Eigen::MatrixXcd H0 = Hamiltonian_Graphene(k, N);
        // std::cout << k << std::endl;
        // std::cout << H0 << std::endl;
        Eigen::SelfAdjointEigenSolver<Eigen::MatrixXcd> eigensolver(H0);
        out << eigensolver.eigenvalues().transpose() << std::endl;
    }
}

更复杂的例子：

其中C++的复数要熟悉std::complex，然后多用浮点数，Eigen多看文档。

对于复杂的缺陷，处理的方法是加入大的on_site，对于简单的缺陷，直接去掉对应的行和列即可。

Eigen::MatrixXcd Hamiltonian_H0_SU4(double k, int N, double t = 1){
    Eigen::Matrix4cd A;
    A << rt3,     0,       1,  std::exp(1i*k),
	      0,     rt3,      1,    1,
	      1,     1,       rt3,   0,
	     std::exp(-1i*k), 1, 0, rt3;
    Eigen::Matrix4cd B = Eigen::Matrix4cd::Zero();
    B(0, 2) = 1;
    B(1, 3) = -1;
    Eigen::MatrixXcd H = Eigen::MatrixXcd::Zero(4*N,4*N);
    Eigen::Matrix4cd Bdiag = B.transpose();
    H.block<4,4>(0,0) = A;
    for(int i = 4; i < 4 * N; i += 4) {
        H.block<4,4>(i,i) = A;
        H.block<4,4>(i - 4,i) = Bdiag;
        H.block<4,4>(i, i - 4) = B;
    }
    return H.block(1, 1, 4 * N - 1, 4 * N - 1);
}

顺便配合python numpy以及matplot画图

import matplotlib.pyplot as plt
import numpy as np
from math import pi


band = np.loadtxt("H0_band.txt")
nk, N = band.shape
k = np.linspace(0, 2*pi, nk)
for i in range(N):
    plt.plot(k, band[:, i])
plt.ylim(-1, 1)
plt.show()