在量子力学中，电子的s、p、d、f轨道最初是通过求解氢原子的薛定谔方程来定义的，氢原子的总波函数可以写为径向部分和角向部分的乘积，即：

其中是径向波函数，是球谐函数，并且n，l，m分别为主量子数、角量子数和磁量子数，不同的轨道对应了不同的角量子数，而磁量子数的个数又由角量子数确定，即-l≤m≤l，因此

s轨道：l=0，m=0，只有一个轨道；

p轨道：l=1，m=-1，0，1，有px，py，pz三个轨道；

d轨道：l=2；m=-2，-1，0，1，2，有dxy，dyz，dxz，dx2，dz2五个轨道。

而能级都是由主量子数决定，因此对于同一能级，s、p、d轨道的能量是简并的。但是这个结论只有在球对称的场中才会成立。

那么对于一般的晶体场该如何分析呢？以八面体晶体场为例，配体所具有的孤对电子被抽象为了点负电荷并与金属原子的d轨道互相排斥，但是由于轨道的方向不同，所受到的排斥不同，最终导致能量不再简并，如下图所示：

根据配位场理论，在局域八面体晶体场作用下， Cr的d轨道分裂成能量较高的两重简并eg(dz2和dx2－y2)和能量较低的三重简并t2g(dxz、dyz和dxy)两组轨道［如上图c所示][1]。

如果这个八面体被拉长，导致eg轨道受到的排斥力不同，就会导致eg轨道发生劈裂，如下图所示[2]：

因此，通过晶体场，我们很容易分析轨道成分的变化。

但是，我们用vasp计算的时候，常常遇到这种问题，例如我下面计算的CoCl2结构，这是一个典型的八面体晶体场的结构：

我们使用vasp计算能带，并画出轨道成分，会得到如下图所示的能带结构：

可以看到，轨道成分似乎没有像我们预想的那样分布，t2g的轨道和eg的耦合到了一起。

其实，这是因为vasp计算的轨道成分是基于实际的空间坐标系xyz进行投影得到的，因此我们需要转换一下POSCAR，将八面体晶体的轴心对准z轴，这样计算得到的轨道成分才严格按照晶体场理论进行劈裂，如下图所示：

此时，重新计算能带结构并画出轨道成分，如下图所示：

可以看到，相比于前面的能带，在费米面附近，eg和t2g明显干净了许多。

我们可以再试试其他例子，下面是二维材料MoTe2的结构：

选择这种坐标系，画出的能带如下图所示：

同样，修改z轴方向：

能带如下图所示：

可以看到，此时eg的能量高于t2g，和前面理论的情况相同。

那么，如何修改坐标系，让z轴指向正八面体的轴呢？

我们可以通过下面的小脚本实现：

#!/usr/bin/env python

import numpy as np, sys, itertools

if sys.version_info[0] == 2: input = raw_input

else: input = input

# **************************

# * parse vasp poscar file *

# **************************

filename = sys.argv[1] if len(sys.argv) > 1 else 'POSCAR'

vasp = open(filename, 'r').readlines()

lattice = np.array([line.split() for line in vasp[2:5]], dtype='float')

atom_type, atom_num = vasp[5].split(), np.array(vasp[6].split(), dtype='int')

atoms = [t for t, n in zip(atom_type, atom_num) for _ in range(n)]

natom = atom_num.sum()

idx = 8 if vasp[7][0].lower() == 'd' else 9

position = np.array([i.split() for i in vasp[idx:idx + natom]], dtype='float')

position_xyz = np.dot(position, lattice)

# **********************

# * get final position *

# **********************

#

# Step 1. get Vec(oz) normalize vector

#

oz_string = input('Please enter one or two vector [fractional coordinates]\n'

                 '\tto define Vector of oz: \n')

temp = oz_string.split()

if len(temp) == 3:

    oz = np.dot(np.array(temp, dtype='float'), lattice)

elif len(temp) == 6:

    vec1 = np.array(temp[:3], dtype='float')

    vec2 = np.array(temp[3:], dtype='float')

    oz = np.dot(vec2 - vec1, lattice)

oz /= np.linalg.norm(oz)

#

# Step 2. project Vec(ox) to plane(oxy) and normlize

#         Vec_new(ox) = Vec(ox) - [Vec(ox) * Vec(oz)] * Vec[oz]

#

ox_string = input('Please enter one or two vector [fractional coordinates]\n'

                 '\tto define Vector of ox: \n')

temp = ox_string.split()

if len(temp) == 3:

    ox = np.dot(np.array(temp, dtype='float'), lattice)

elif len(temp) == 6:

    vec1 = np.array(temp[:3], dtype='float')

    vec2 = np.array(temp[3:], dtype='float')

    ox = np.dot(vec2 - vec1, lattice)

ox = ox - np.dot(ox, oz) * oz

ox /= np.linalg.norm(ox)

#

# Step 3. cross production of oz, ox ot get oy

#

oy = np.cross(oz, ox)

mat1 = np.concatenate((ox, oy, oz)).reshape((3, 3))

rot = np.linalg.solve(mat1, np.identity(3))

lattice_new = np.dot(lattice, rot)

# *********************************************

# * output to new position and check xyz file *

# *********************************************

# output to rotated.vasp

vasp[2:5] = ['%20.12F %20.12F %20.12F\n' % tuple(line) for line in lattice_new]

open('rotated.vasp', 'w').writelines(vasp)

# supercell and output to rotated.xyz

position_new = np.dot(position_xyz, rot)

with open('rotated.xyz', 'w') as fb:

    fb.write(' %d\nmolecule\n' % (natom * 3**3))

    for i, j, k in itertools.product(*map(range, [3, 3, 3])):

        position_new2 = position_new + np.dot([i, j, k], lattice_new)

        for _type, (x, y, z) in zip(atoms, position_new2):

            fb.write('%-2s %14.6F %14.6F %14.6F\n' % (_type, x, y, z))

使用方法也很简单，直接运行代码，根据提示输入z轴和x轴的向量的起始和终点坐标即可，而z轴和x轴确定之后，y轴会自动确定，因此只需要定义两个方向。

例如我下面的POSCAR，注意需要分数坐标下的：

Cl2 Co

1.0

     -2.158964588602     -0.196442203497     6.524749936631

      1.249606232921     -1.771497077862     5.396893118714

      0.909358355680      1.967939281359     5.509477221611

 Co   Cl

 1    2

Direct

     0.000000000     0.000000000     0.000000000

    0.250000000     0.250000000     0.250000000

    0.750000000     0.750000000     0.750000000

把上面的脚本命名为rot.py，依次在集群中输入命令：

python rot.py

Enter  

0.000000000 0.000000000 0.000000000 0.750000000 0.750000000 0.750000000  

Enter  

0.000000000 0.000000000 0.000000000 0.250000000 0.250000000 0.250000000  

Enter

运行结束后就会出现一个新的POSCAR，用这个进行后面的计算即可。

参考文献：

[1]李树宗,司君山,吴绪才,等.巧用Wannier函数分析局域坐标系下的晶体场劈裂[J].大学物理,2022,41(12):31-35.DOI:10.16854/j.cnki.1000-0712.220192.

[2]Jiang C, Lu K, Wang R, et al. Magnetic Anisotropy and Jahn–Teller Effect in Ferromagnetic Two-Dimensional CrGa2Te4[J]. ACS Applied Electronic Materials, 2022, 4(7): 3220-3225.