物理问题：空间中$n$个电子，质量为$m$，相距为$a$，放置在一维势阱($U=0.0$)中，只计最近邻相互作用，求其耦合能谱分布。$n$大于3时可以考虑周期性边界条件。

第一步，转化为矩阵方程。

因为求定态问题，列定态薛定谔方程如下：

$[-frac{hbar^2}{2m}frac{partial^2 }{partial x^2}+U(x)]psi=Epsi$

其中含有一个二阶导数，因为是数值求解，我们用二阶差分代替它，

$frac{partial^2psi(x) }{partial x^2}|_{x=x_n}rightarrow frac{psi(x_{n+1})-2psi(x_{n})+psi(x_{n-1})}{a^2}$

这里采用$a$作为微分步长，因为$a$是个常数，不妨把它与其它常数写在一起，定义

$t_0=frac{hbar^2}{2ma^2}$

同时在一维势阱中$U(x_n)$也是个常数，简化成$U_0$。稍做整理，薛定谔方程就变成了差分形式的递推式

$(U_0+2t_0)psi_n-t_0psi_{n-1}-t_0psi_{n+1}=Epsi_n$

由于波函数的正交性有$psi_m^*psi_n=delta_{m,n}$，因为我们在上式两边同乘以$psi_m^*$并对$m$求和

$sum_{m}[(U_n+2t_0)delta_{n,m}-t_0delta_{n,m+1}-t_0delta_{n,m-1}]psi_m=Esum_{m}delta_{n,m}$

可见只有当$|n-m|<=1$时矩阵元不为零：

$[H]=begin{bmatrix}U+2t_0 & -t_0 & 0 &cdots   &0 \ -t_0 & U+2t_0 & -t_0 &cdots &0 \ 0 & -t_0 & U+2t_0  & cdots &0 \ vdots & vdots  & vdots  & ddots  & vdots \0 & 0 & 0 & cdots & U+2t_0end{bmatrix}$

要找到向量$E$使得$[H]-[EI]=0$，这样转化成了求$[H]$本征值问题。

第二步，利用数学库求解本征值问题。

观察$[H]$矩阵可见，这是一个实对称矩阵，其本征值为一系列实数。于是我们在MLK数学库中寻找（实）对称矩阵本征值的子程序。打开MKL的REFERENCE MANUAL

http://software.intel.com/sites/products/documentation/doclib/mkl_sa/11/mklman/index.htm

找到LAPACK Routines: Least Squares and Eigenvalue Problems 下的

从说明可见这个函数正是我们所需要的。?syev是子程序的名字，其中前面的问号可以用不同字母替换来适用于不同精度。这个问题中我们将采用双精度，所以按理说应该使用dsyev。不过在FORTRAN 95中，我们不需要手动选择精度，程序会根据我们定义的矩阵数据类型来确定所使用的子程序，所以我们只要无视问号选用syev就可以了。

最简单的情况是体系中有两个电子，即求解矩阵

程序段很简单，有几个注意的地方：

1. 第二行 USE mkl95_LAPACK 是必要的，因为需要调用外部函数；

2. mtH是哈密顿矩阵（MatrixHamiltonian,2$times$2实矩阵），egW是存放本征值的列向量（EigenValue,2$times$1实向量），它们作为输入参数被SYEV调用；

3. syev是前述的MKL子程序，从手册中我们可以看到这个程序可以输入的参数很多，即使在FORTRAN95中也有多达5个参数

call syev(a, w [,jobz] [,uplo] [,info])

不过后面三个是可选的，也就是说如果我们只需要本征值的话，后面的三个可以完全不写。不过有时候我们还想知道本征向量，那么就要使用后面的选项，将jobz设置为'V'，这样在运行过后原先的输入矩阵a就被本征向量矩阵覆盖了。

第三步，编译运行。

程序完成后需要编译运行，不过因为调用了MKL库，我们需要告诉程序MKL库的位置。这一步稍微复杂一些。先找到计算机中MKL的安装位置，通常会跟INTEL FORTRAN 编译器放在一起。以下是以笔者的系统为例，在系统终端输入

which ifort

得到

/opt/intel/bin/ifort

那么在/opt/intel这个文件夹寻找MKL的安装路径，发现为

/opt/intel/mkl/lib/intel64

注意区分32位和64位的库。然后开始编译，假定刚才写的源程序文件名为a.f90，用命令

ifort a.f90  -I/opt/intel/mkl/include/intel64/lp64 -L/opt/intel/mkl/lib/intel64 -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -lmkl_blas95_lp64 -lmkl_lapack95_lp64 -liomp5

嗯，看起来比较长，不过一般不需要改动。大体说一下参数的意义：

首先 -L是标明后面要紧接着给出库文件所在的路径，这里我们加上MKL的lib文件夹。如果采用这种方法写的话，后面的库文件名就需要一个小的变动。比如我们在lib文件夹中可以发现名叫libmkl_core.a和libmkl_core.so的文件（其实是同样功能，只不过前面是静态库，后面是动态库），但是在编译中我们要把lib写成-l，因此就变成了-lmkl_core。这个库是MKL的核心库，一般要使用MKL时都需要调用。其它重要的库还有blas95和lapack95这两个，lapack95就是我们需要的FORTRAN95版本的LAPACK，它依赖于blas95，所以要一起调用。

另一块以-I开头的是include，这里我们在程序中用到了use mkl95_LAPACK，这个被使用的mkl95_LAPACK.mod文件就放在MKL的INCLUDE文件夹下，大家可以找找看。

其它的几个库都是起辅助作用的，比如线程管理之类，只要包含进来就可以了。

编译通过之后，产生了a.out文件，运行它就输出两行

这里打印出了体系的本征能量。还记得计算之前两个电子的能量吗？$U_0+2t_0=2.0$。所以在两个电子相互耦合之后，共同构成了两个新的能级，刚好是1.0和3.0。是不是很简单呢？

这样我们就可以求更大一些的体系，比如五电子链体系。这时候除了刚才的做法，还可以考虑周期性边界条件，具体来说就是在设置$[H]$矩阵时将两个角上的矩阵元$left langle 1|H |5 right rangle$和$left langle 5|H |1 right rangle$也设置为$-t_0$就可以了（不考虑的时候是零）。这时的矩阵依然是实对称矩阵，还是可以用上面的MKL子程序来求解。