Quantum mechanism study

最近学习下《朗道量子力学教程》,并想通过这个blog记录下学习过程中的疑问和感想。

§2 叠加态

位形空间: 经典力学中,表示包含系统可能所有状态的一个空间。 比如描述系统的坐标集为q(这里粗体是表明这是多个坐标参量的集合而不仅仅是单个坐标,与下后面括号中的q区分开来), q可能包含多个坐标参量 $(q, q’, q’’, q’’’…)$,那么系统的任一状态就可以表示为$(q_0,q_0’,q_0’’,q_0’’’ …)$。 假如只考虑一个单粒子的状态,不考虑其速度,只考虑其坐标作为其状态的表示,那么它的位形空间就是普通空间,该单粒子的任意一个状态,即该单粒子处于空间的某个位置$(x_0,y_0,z_0)$,$d\bf{q}$就是普通空间的一个体积元dV, 即$dxdydz$.

假如一个体系的状态可以通过坐标函数在$\psi(\textbf{q})$表示,也就是说在系统在$d\textbf{q}_0$的体积元内,表示系统状态的值为$\psi(\textbf{q}_0)$. 通过一个经典的例子理解:假如在一个封闭盒子充满了NaCl溶液,且$Na^+$的浓度是不均匀的。 如果用$Na^+$的状态表示系统的状态,表示为$\rho(\bf{r})$,则在$x_0,y_0,z_0$位置的小体积元内,$Na^+$的浓度为$\rho(x_0,y_0,z_0)$, $Na^+$的数量为$\rho(x_0,y_0,z_0)dxdydz$.

疑问

  1. §11矩阵中,常常用到定态的本征函数$\psi_n$, $n$是对应本征能量$E_n$,即$\psi_n$是该体系中本征能量为$E_n$的电子在空间的分布。那么本征能量是一个值吗? 能带图中一个能带(能带图中的一条线)是具有一定的宽度的。 能带图中的一个条能带对应的电子的空间分布$\psi_n$是同一个吗?

答:本书所说的典型的定态结构一般是指如球方势阱等非周期结构,定态中的本征函数$\psi_n$对应的就是一个能量值,而本征函数$\psi_n$就是电子的一种空间分布。

而能带的概念是适用于周期性的结构,需要结合bloch理论进行分析,而§11所讲的内容是基于定态(封闭系统或处于恒定外场中的系统),当然所涵盖的范围应该包括封闭系统下的周期结构。周期性结构中,我们计算的过程中,先指定K点,再解KS方程,得到本征能量和本征函数,讲所有选取的K点(K mesh)都计算完成后,每个K点都有对应的一系列的本征值和本征函数。 所以周期结构中,能量的表示方法一般为 $E_n,_k$, 而本征函数 $\psi_n,_k(x,y,z)$。 $E_n,_k$沿着高对称点所取的路径选择K点计算得到的一系列本征值所连的曲线就是能带结构中的一条能带。 那么一个能带所对应的电子的空间分别不是将所有的$\psi_n$$_0,_k(x,y,z)$进行叠加吗?

人工智能学习

《Artificial Intelligence A modern approach》 学习笔记

2 Intelligence Agent

Agent: 是指接受外界信息(percept),进行处理后,做出对应反应(action)的事物。
Intelligence Agent: 主要强调rationality。

本章思考

agent 的概念实际上是适用范围极广,即接收外界信息做出反应的anything。 那么人就是一个极为复杂的agent,一个人无时无刻不在接受外界信息并做出反应,简单的叫应激反应,复杂的需要通过长时间的思考再行动。
agent的设计由简单到复杂有:

  • simple reflex agents
  • model-based reflex agents
  • goal-based agents
  • utility-based agent

Appendix A.1.1

Asymptotic analysis

主要讲的是数据结构中的复杂度分析,主要用来评估算法的效率。 具体见数据结构之算法时间复杂度数据结构之算法时间复杂度

Matlab application in my research

NIXSW Simulation

NIXSW的原理

垂直入射X射线驻波技术的原理简单可以概括为: 单色X射线垂直入射晶体表面时,会形成周期与晶面间距相同的驻波,且该驻波从晶体内部延伸至真空,故表面吸附的原子分子亦处在驻波场中,可以通过调节驻波的波长观察原子分子的光电子发射变化,获得原子分子相对驻波场的具体位置,又因驻波的波节(波峰?)的位置与晶面位置一致,故可以确定原子分子与最表层原子的垂直距离,即原子在表面的高度的定量确定。

等式右侧为晶面间距,z为吸附原子相对于驻波波节的距离,即原子在表面的高度,这两个参数可从优化后得模型获得(z在模型中为分立的值,因此实际计算的时候右侧积分会转变为求和的计算),因此可以通过理论计算的模型获得右侧的值,而左侧D和也可通过二元一次方程解得。

matlab读取DFT优化后结构参数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
function [Direct2Cartesian, AllElements, Element_Num,Num_Element_PositionD,Num_Element_PositionC_Height] = Read_CONTCAR(NumOftopCu,filename)
%This function is used to read the information of optimized structure from
%VASP calculation.

%Direct2Cartesian is the matrix to transform the Direct coordinates to
%Cartesian coordiantes, i.e. the lattice parameters of the supercell.

%ALLElements contains all elemetns in the structure, for F4TCNQ/Cu(111) Cu
%C N F.

%Element_Num contains the Number of atoms of each element according to the
%order in AllElements.

%Num_Element_PositionD contains the Direct coordinates of all atoms with
%the same order in CONTCAR

%,Num_Element_PositionC_Height contains the heights of all atoms with
%the same order in CONTCAR

%NumOftopCu give the serial number of top-most atoms .

%filename the output from VASP, i.e. CONTCAR.

fid=fopen(filename);
Direct2Cartesian = cell2mat(textscan(fid, '%f %f %f',3, 'Headerlines', 2,'CollectOutput',1));
A=textscan(fid, '%s %s %s %s',1 , 'CollectOutput', 1);
%***** the number of %s should be changed to be the same as number of elements

AllElements=A{1,1};
B=textscan(fid, '%d %d %d %d',1, 'CollectOutput', 1);
Element_Num=B{1,1}; Total_atoms = sum(Element_Num);
C=textscan(fid, '%f %f %f %*s %*s %*s',Total_atoms, 'Headerlines', 3, 'CollectOutput',1);
% %*s %*s %*s is to find and ignore "T T T" or "F F F"
% The 'Headerlines' has to be 3 not 2.

AtomsPositionD = C{1,1};
Num_Element_PositionD=[[1:Total_atoms]' string([repelem(AllElements, Element_Num)]') AtomsPositionD];
AtomsPositionC= AtomsPositionD*Direct2Cartesian;
TopLayerSubHeight=mean(AtomsPositionC(NumOftopCu,3));
Atom_Height= AtomsPositionC(:,3)-TopLayerSubHeight;
disp(Atom_Height)
Num_Element_PositionC_Height=[[1:Total_atoms]' string([repelem(AllElements, Element_Num)]') AtomsPositionC Atom_Height];


fclose(fid);

#

DFT calculation (VASP) in surface science

Simulation of dI/dV spectrum and mapping

dI/dV 对应于LDOS

STM 表征中dI/dV直接与理论计算中局域电子态密度LDOS相对应,LDOS为空间坐标和能量的函数$\rho(r,E)$,dI/dV谱对应特定位置r0的LDOS分布$\rho(r_0,E)$,dI/dV mapping 对应特定能量的DOS的空间分布$\rho(r,E_0)$。
LDOS定义为:

即能量为E的电子态的空间概率分布。回想VASP中模拟STM图的过程正是将$E_F$至$E_F-eV$范围内的态密度进行叠加,我们也可以通过类似的过程进行di/dV模拟。
整个计算过程包括

  • 计算体系的DOSCAR文件,确定DOS中态密度范围
  • 计算DOS中态密度下限至不同能量区间$E_i$的积分LDOS对应得PARCHG文件
  • 对PARCHG文件进行分析,获得特定位置积分LDOS,并微分得LDOS

    计算体系的DOSCAR文件,确定DOS中态密度范围

    DOS 计算得时候参数和静态计算一致,只不过得加上LORBIT=10或11。
    INCAR
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    ISTART = 0  
    ICHARG = 2
    PREC = Accurate
    ENCUT = 480
    EDIFF = 1.0e-06
    ALGO = Fast
    ISMEAR = 0; SIGMA = 0.01
    LREAL = .FALSE.
    LCHARG = .TRUE.
    LAECHG = .TRUE.
    LORBIT = 10

计算得到的DOSCAR包含了DOS信息,其中DOS是通过积分DOS微分而来。

用p4v软件将DOSCAR作图获得DOS图。(直接打开DOSCAR获得的DOS图费米能级未移到横坐标0的位置,而打开vasprun.xml获得的DOS图费米能级位于0)

由图可知,DOS中态密度不为零的能量范围为(-20,6)。(此处费米能级位于0的位置)

计算积分LDOS

复制上述DOS计算后的文件夹,并进入该文件夹,用下面的脚本进行积分LDOS的计算。

Shell script for PARCHG calculation

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
#!/bin/bash
for i in {1..260}
do
Vi=-20
Vf=$(echo "scale=2; -20+26/260*$i"|bc)
echo $i $Vi $Vf
cat >INCAR<<!
ISTART = 1
ICHARG = 0
ENMAX = 480
ISMEAR = 0 ;SIGMA = 0.01
ALGO= Fast
PREC = Single
IVDW = 11
EDIFF = 1E-4

partial charge densities:
LPARD = .TRUE.
LSEPK = .FALSE.
LSEPB = .FALSE.
NBMOD = -3 % Vi and Vf are vs. fermi energy
EINT = $Vi $Vf

mechine:
LPLANE = .TRUE
NCORE = 8
NSIM = 4
!
mpirun -np 32 vasp_std
mv PARCHG $i-PARCHG_$Vf
done

该脚本计算得到-20 到-19.9, 19.8, 19.7…5.8, 5.9, 6.0的260个能量范围的空间态密度积分,对应得文件为1-PARCHG-19.9, 2-PARCHG-19.8 … 259-PARCHG_5.9, 260-PARCHG_6.0。
PARCHG 记录了Vi至Vf区间态密度的积分。

这其实是个空间坐标的函数,在PARCHG中以一个三维矩阵记录态密度信息。

该文件中将Unit cell 分成84x84x140的空间格子,并将各个格子内的态密度与unite cell体积的乘积在PARCHG中记录,记录的顺序为(1,0,0)(2,0,0)…(84,0,0)(0,1,0)…(84,84,0)(0,0,1)…(84,84,140)

计算LDOS谱和LDOS mapping

LDOS谱的计算需要选定空间位置r,此处我们选择距离分子中心高度为$2\overset{\circ}{A}$的位置,对应的分数坐标为:0.5508 0.6530 0.3203,对应的空间网格为46 55 46,对应PARCHG中第(46-1)x84x84+(55-1)x84 +46=322102个数据,PARCHG中密度数据记录是按照一行十列,因此位于第32210行,第2列。按照下面的脚本可以从计算的i-PARCHG-Vf 文件中读取该点的数据,获得该位置积分LDOS曲线。

dIdV_dataAnalyse

1
2
3
4
5
6
7
8
9
10
11
12
13
#!/bin/bash
rm LDOS_data
e=0.01
for i in {1..260}
do
Vi=-20
Vf=$(echo "scale=2; -20+26/260*$i"|bc)
LDOS=$(awk -v line=$(awk '/84 84 140/{print NR}' $i-PARCHG_$Vf) '{if(NR==line+32210){print $2}}' $i-PARCHG_$Vf)
echo $Vf $LDOS >>LDOS_data
echo $i $Vf $LD0S
done
gnuplot
plot LDOS_data

在origin或者matlab中对LDOS_data中的数据进行作图获得积分LDOS曲线,作图的时候数据除去了Unit cell 的体积。

进行微分可获得LDOS曲线

LDOS maping

dI/dV mapping 对应特定能量的LDOS的空间分布$\rho(r,E_0)$可以通过下面的公式求得:

通过下面的脚本可以计算上式积分部分,计算时需要给定$E_0$和极小值$\varepsilon$。

计算$E0±\varepsilon$范围电子态密度空间分布

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
#!/bin/bash
E0=-2
for e in 0.005 0.01 0.025 0.05 0.075
do
Vi=$(echo "scale=6; $E0-$e"|bc)
Vf=$(echo "scale=6; $E0+$e"|bc)
echo $Vi $Vf
cat >INCAR<<!
ISTART = 1
ICHARG = 0
ENMAX = 480
ISMEAR = 0 ;SIGMA = 0.01
ALGO= Fast
PREC = Single
IVDW = 11
EDIFF = 1E-4

partial charge densities:
LPARD = .TRUE.
LSEPK = .FALSE.
LSEPB = .FALSE.
NBMOD = -3
EINT = $Vi $Vf

mechine:
LPLANE = .TRUE
NCORE = 8
NSIM = 4
!
echo "EINT = $i "; mpirun -np 48 vasp_std
mv PARCHG PARCHG_$E0_$e
done

由于分母$2\varepsilon$为常数,所以PARCHG中的数据直接对应$\rho(r,E_0)$,因此直接通过PARCHG可以模拟dI/dV mapping。

Phonopy study

Phonopy的学习主要参照官网的内容: https://atztogo.github.io/phonopy

Phonopy 安装

Step 1: 安装Phonopy 所需的phonon环境包。
sudo apt-get install python-dev python-numpy python-matplotlib python-yaml python-h5py
Step 2: 在 https://pypi.org/project/phonopy/#files 中下载phonopy安装包phonopy-1.13.2.32.tar.gz。然后解压该软件,进入解压后的目录打开终端,输入下面的代码进行phonopy安装即可。
sudo python setup.py install
Step 3:查看phonopy位置,如输出phonopy的路径即安装成功。
which phonopy

Phonopy 使用

依据官网的教程,使用VASP和Phonopy进行结构热力学性质和声子谱等一系列性质的计算,并同官网的结果相检验。

有限位移方法计算结构力常数

计算力常数是计算结构声子性质的第一步,计算力常数的方法有两种:有限位移方法和密度泛函微扰理论。有限位移方法通过计算结构中各个原子微小位移前后能量的变化获得结构力常数。
Step 1:准备需要计算到结构的VASP输入文件,此处命名为POSCAR-unitcell。(一般为已经优化完成的结构)

POSCAR-unitcell

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
NaCl
1.00000000000000
5.6903014761756712 0.0000000000000000 0.0000000000000000
0.0000000000000000 5.6903014761756712 0.0000000000000000
0.0000000000000000 0.0000000000000000 5.6903014761756712
4 4
Direct
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 0.5000000000000000 0.5000000000000000
0.5000000000000000 0.0000000000000000 0.5000000000000000
0.5000000000000000 0.5000000000000000 0.0000000000000000
0.5000000000000000 0.5000000000000000 0.5000000000000000
0.5000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 0.5000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 0.5000000000000000

使用下面命令生成2x2x2的超胞SPOSCAR文件(扩胞是消除虚频常用的方法),包含位移信息的disp.yaml文件,以及对应各个原子各个方向位移后的POSCAR-001 … POSCAR-00N文件(如无对称性,这些文件的数量等于POSCAR中原子数量的6倍,对应于各个原子三个自由度上正负方向的位移)。
phonopy -d --dim="2 2 2" -c POSCAR-unitcell
step2:将各个POSCAR-001 … POSCAR-00N分别作为POSCAR用VASP进行静态计算(用shell编个循环语句的shell脚本计算即可),得到vasprun.xml-001 … vasprun.xml-00N 文件。

静态计算INCAR

1
2
3
4
5
6
7
8
9
PREC = Accurate
IBRION = -1
ENCUT = 500
EDIFF = 1.0e-08
ISMEAR = 0; SIGMA = 0.01
IALGO = 48
LREAL = .FALSE.
LWAVE = .FALSE.
LCHARG = .FALSE

step3:再用phonopy处理即可得到力常数FORCE_SETS文件

1
phonopy -f vasprun.xml-{001..00N}

密度泛函微扰理论计算结构力常数

用VASP-DFPT计算力常数更为简单,同时还能给出体系简正振动频率,伯恩有效电荷,压电常数等信息(LEPSILON = .TRUE.),所以我们计算的时候主要采取这种方法。
step1:将超胞SPOSCAR重命名为POSCAR作为VASP的输入文件计算体系力常数(用超胞进行计算能获得更准确到结果)。
mv SPOSCAR POSCAR
Step 2: 用IBRION = 8 和 NSW = 1计算体系力常数。为使受力计算更准确,提高电子步迭代收敛要求(EDIFF = 1E-8,并将是空间积分LREAL设为FALSE。具体INCAR如下。

INCAR

1
2
3
4
5
6
7
8
9
10
11
12
PREC = Accurate
ENCUT = 500
IBRION = 8
EDIFF = 1.0e-08
IALGO = 38
ISMEAR = 0; SIGMA = 0.01
LREAL = .FALSE.
ISYM = 0
ADDGRID = .TRUE.
LWAVE = .FALSE.
LCHARG = .FALSE.
LEPSILON = .TRUE.

准备好KPOINTS和POTCAR文件就可以进行VASP计算了,计算后用下面的命令检查是否存在虚頻(此处存在虚頻,是因为我们直接用的官网的POSCAR,该结构没有用我们计算用到参数和赝势进行优化)。
grep 'cm-1' OUTCAR
Step 3: 用下面到命令从输出文件vasprun.xml得到包含力常数的文件FORCE_CONSTANTS。
phonopy --fc vasprun.xml

计算声子谱

Step 1: 新建band.conf文件,设定计算声子谱的参数。由于我们关心的是POSCAR-unitcell的声子性质,而计算受力时是用超胞进行计算,所以在band.conf中需要注明DIM = 2 2 2,而在step2和step3计算时用参数 -c POSCAR-unitcell使计算的结果对应扩胞前的结构。

band.conf

1
2
3
4
5
6
ATOM_NAME = Na Cl
DIM = 2 2 2
PRIMITIVE_AXIS = 0 0.5 0.5 0.5 0 0.5 0.5 0.5 0
FORCE_CONSTANTS = READ
BAND = 0.00 0.00 0.00 0.50 0.00 0.00 0.50 0.50 0.00 0.00 0.00 0.00 0.50 0.50 0.50
BAND_POINTS = 200

Step 2: 输入下面命令计算声子谱。

Step 3: 通过下面两个命令使用phonopy-bandplot脚本可直接画出声子谱或得到文本格式的声子谱数据。
phonopy-bandplot
phonopy-bandplot --gnuplot >PhononBAND.dat

计算声子谱态密度和体系热力学性质

step 1:新建mesh.conf文件。

mesh.conf

1
2
3
4
ATOM_NAME = Na Cl
DIM = 2 2 2
MP = 8 8 8
FORCE_CONSTANTS = READ

step 2:使用下面的命令计算体系声子谱的态密度。
phonopy -p mesh.conf -c POSCAR-unitcell

step 3:使用下面的命令计算体系的热力学性质(自由能,零点能,熵,热容)。
phonopy -t -p mesh.conf -c POSCAR-unitcell >thermal_properties.dat

surface science technology

晶体表面吸附物LEED表征

使用LEEDpat通过吸附矩阵模拟LEED图案

LEEDpat软件可以直接从网上下载,直接安装即可。
使用方法:

Enjoy Life

这个博客算是基本搭好了,最近很忙,但总还是该写点东西。

Physics

我是个学物理的,这个blog的初衷是记录在学习过程中的问题与感悟,并希望以最直观的方式表达出来。

Music

学物理的人或多或少有些超然傲气甚至是偏执,因而也容易感觉没有归属感,而音乐就是一种表达自己,理解自己的方式。

Emotion

自我的情绪有的时候是很容易感觉到,但是难以控制。

END

由于智商有限,物理学的不好,也没有音乐细胞,以前情绪还不好。但一直希望在这三个方面有所提升,也断断续续地在学东西。

很多东西做不到没有关系,但至少得知道自己应该要做什么。


,