Learn VASP The Hard Way

在ipad，手机上运行终端，连接服务器可以让我们随时随地查看，提交任务。在这里推荐一款ipad和iphone上的终端APP: SSHaking。如果你用的安卓系统，一些推荐的APP可以直接到文章底部去找。
先依照个人使用的经验，列出来以下SSHaking的优点：

1） 界面简洁，没广告；

2） 可以同时打开好几个terminal界面；

3）连接服务器设置也很简单；

4） 支持连接多个Hosts； （平时常用的terminus只支持一个，连多个的话要花钱，服务器如果只有一个，也推荐terminus）

5） 免费！！！！

6）APP的开发者很牛，可能也很帅！

7）希望大家多多下载这款APP，有什么bug及时给作者提出来，而且作者看到后修复的也很及时。

如果你用的是安卓系统，可以有以下几个选择：

1） juiceSSH (轻和浅吟🐳)
2） Admin hands (HeyTamas)

括号里是提供APP名字的好心人。如果你也有推荐的其他好用的APP(免费的)，可以在微信公众号留言。

https://mp.weixin.qq.com/s/jvghzIlV6NOsNFjUa445xQ

Citrine 也是一个材料学数据库，自带了机器学习工具，大家也可以研究琢磨一下。官网是 https://citrination.com/search/simple?searchMatchOption=fuzzyMatch。这里，我们将通过几个例子了解如何在这个数据库中获取数据。在Citrine中获取数据需要用到 matminer.data_retrieval.retrieve_Citrine.CitrineDataRetrieval 这个工具，其跟materials project类似，也学要一个API_key。大家可以在官网注册，并记录下API_Key。

1
2
3

#请先在终端中输入： pip install citrination_client ； 否则会报错
#然后我们回到代码，导入 CitrineDataRetrieval 这个模块
from matminer.data_retrieval.retrieve_Citrine import CitrineDataRetrieval

示例1：获取化学式为 Si 的所有材料的实验测得的能带

1
2

# 先实例化这个工具，需要提供你的API_Key
cdr = CitrineDataRetrieval(api_key='')

1
2
3
4

# 接下来，类似从materials project获取数据，我们需要用get_dataframe函数
df = cdr.get_dataframe(criteria={'formula':'Si', 'data_type': 'EXPERIMENTAL'}, 
                       properties=['Band gap'],
                       secondary_fields=True)

1
2
3
4
5
6
7
8
9
10
11

  0%|                                                                                            | 0/7 [00:00<?, ?it/s]d:\programs\anaconda3\envs\py36\lib\site-packages\matminer\data_retrieval\retrieve_Citrine.py:103: FutureWarning:

pandas.io.json.json_normalize is deprecated, use pandas.json_normalize instead

100%|████████████████████████████████████████████████████████████████████████████████████| 7/7 [00:00<00:00, 87.33it/s]

all available fields:
['Band gap', 'uid', 'Band gap-units', 'category', 'references', 'chemicalFormula', 'Crystallinity', 'Band gap-conditions', 'Band gap-methods', 'Band gap-dataType']

suggested common fields:
['references', 'chemicalFormula', 'Crystallinity', 'Band gap', 'Band gap-units', 'Band gap-conditions', 'Band gap-methods', 'Band gap-dataType']

1

df #查看表格

从上面结果可以看出，数据库里有化学式为Si的材料共7中实验测得的能带数据。这些数据来源的文献，带隙的值，实验方法都被记录了下来。

示例2：获取材料对 *OH 和 *O物种的吸附能

想不到citrine数据库还有各种材料对*OH 和 *O 物种的吸附能，一起来看看怎么获取把。要查询更多的数据库，请探索官网https://citrination.com/datasets

1
2
3

# 我们发现这里不需要指定筛选条件，只需要返回所需要的性质即可。
df_OH = cdr.get_dataframe(criteria={}, properties=['adsorption energy of OH'], secondary_fields=True)
df_O = cdr.get_dataframe(criteria={}, properties=['adsorption energy of O'], secondary_fields=True)

  0%|                                                                                            | 0/9 [00:00<?, ?it/s]d:\programs\anaconda3\envs\py36\lib\site-packages\matminer\data_retrieval\retrieve_Citrine.py:103: FutureWarning:

pandas.io.json.json_normalize is deprecated, use pandas.json_normalize instead

100%|███████████████████████████████████████████████████████████████████████████████████| 9/9 [00:00<00:00, 128.63it/s]


all available fields:
['uid', 'Adsorption energy of OH', 'Adsorption energy of OH-units', 'Morphology', 'Surface facet', 'category', 'references', 'chemicalFormula', 'Adsorption energy of OH-conditions', 'Adsorption energy of OH-dataType']

suggested common fields:
['references', 'chemicalFormula', 'Surface facet', 'Adsorption energy of OH', 'Adsorption energy of OH-units', 'Adsorption energy of OH-dataType', 'Morphology', 'Adsorption energy of OH-conditions']


100%|█████████████████████████████████████████████████████████████████████████████████| 21/21 [00:00<00:00, 140.70it/s]

all available fields:
['uid', 'Surface facet', 'Adsorption energy of O', 'category', 'references', 'chemicalFormula', 'Reconstruction', 'Adsorption energy of O-units', 'Adsorption energy of O-conditions']

suggested common fields:
['references', 'chemicalFormula', 'Surface facet', 'Adsorption energy of O', 'Adsorption energy of O-units', 'Adsorption energy of O-conditions', 'Reconstruction']

1

df_OH.head()

1

df_O.head()

大家可以看到，这些数据库里面的数据还不够丰富。大家也可以把自己计算得到的数据上传到citrine数据库中，促进科学知识的开放获取。

• 实际上 citrine 数据库包含了很多实验数据，下次的分享中，我们会把citrine数据库中的实验数据和materials project中的计算数据做个对比。
• 强烈推荐大家在学习的时候使用jupyter notebook。
• 附上jupyter notebook文件共大家参考：【链接：https://pan.baidu.com/s/1W1XKj-Fnjmx1yIe2QB7_MQ 提取码：jvoi 】

本次分享到此结束。

https://nbviewer.jupyter.org/github/hackingmaterials/matminer_examples/blob/master/matminer_examples/data_retrieval-nb/data_retrieval_basics.ipynb

从Materials Project 数据库获取

从Materials Project数据库获取数据得用到 matminer.data_retrieval.retrieve_MP.MPDataRetrieval 功能。

1
2
3
4
5

# 首先导入模块
from matminer.data_retrieval.retrieve_MP import MPDataRetrieval
# 实例化 MPDataRetrieval 这个类
mpdr = MPDataRetrieval(api_key='your_own_api_key') 
# 在实例化 MPDataRetrieval 时需要输入用户在 Material Project 网站的 API Key.

示例1：获取所有单元素材料的密度

上一次的介绍说到，matminer 处理的数据跟 pandas 库一样，是 dataframe，所以我们要用 get_dataframe 方法获取数据; dataframe 看起来就像excel 表格一样。

1
2
3

df = mpdr.get_dataframe(criteria={"nelements": 1}, properties=['density', 'pretty_formula'])
print("There are {} entries on MP with 1 element".format(df['density'].count())) # 计算有多少材料
#请自行关注 df['density'].count() 语句

criteria 是搜索条件，它是一个字典，字典里面可以有多个键值对; properties 是想要获取的材料性质，它是个列表，里面可以有多个性质。其实criteria 也是材料的性质。想知道所有可用的性质，请参考 https://github.com/materialsproject/pymatgen/blob/master/pymatgen/ext/matproj.py

上述代码意味着 我们想在数据库中搜索材料，这些材料元素种类为1，对于这些材料，我们想得到它们的密度和化学式信息。

结果如下：

让我们看一下这个 dataframe 表格长什么样

1
2

df.head()  # 这个语句会显示 df 这个表格的前几行，可以看到，我们获得了 716 个单元素材料。
           # 表格里记录了我们需要的密度 和 化学式的性质，以及它们的ID

结果如下：

示例2：获取所有带隙大于 4.0 eV的材料

方法挺简单的，就是在上个例子的 get_dataframe 方法中设置新的搜索条件: 带隙大于4.0 eV.

1
2
3
4
5

# 带隙是 band_gap，大于4怎么写呢？ 是 greater than, 缩写成 gt. 在 matminer中要写成 $gt 
# 大家注意体会 criteria 的字典形式。criteria=, properties= 可以省略不写。
df = mpdr.get_dataframe({"band_gap": {"$gt": 4.0}}, ['pretty_formula', 'band_gap'])

print("There are {} entries on MP with a band gap larger than 4.0".format(df['band_gap'].count()))

结果如下：

再看一下表格长什么样子

1

df.head()

结果如下：

如果想把表格保存成excel可读写的格式，可以用下面的语句：

1

df.to_csv('materials_bg_gt_4.csv')

示例3：获取 VRH 剪切模量和体积模量¶

首先得要求这些材料存在弹性的性质，其次，我们还想要求这些弹性数据没有警告信息

1
2
3
4
5
6

# 存在弹性常数信息用 "elasticity": {"$exists": True} 表示，没有警告信息用 一个空列表表示："elasticity.warnings": []

df = mpdr.get_dataframe({"elasticity": {"$exists": True}, "elasticity.warnings": []},
                        ['pretty_formula', 'elasticity.K_VRH', 'elasticity.G_VRH'])

print("There are {} elastic entries on MP with no warnings".format(df['elasticity.K_VRH'].count()))

搜索结果如下：

我们想统计一下这些数据

1

df.describe() # 该语句可以很方便地对表格中的每一列数据进行统计，给出数量，平均数，方差，最小、最大值等信息

得到：

接下来看一个更复杂的例子

1
2
3
4
5
6
7
8
9
10
11
12
13

'''
除了上次的搜索条件外，我们想搜索包含 Pb 和 Te 的材料："elements": {"$all": ["Pb", "Te"]} 
材料的稳定性也在考虑之中，energy above hull 必须在1e-6以下："e_above_hull": {"$lt": 1e-6}
'''
df = mpdr.get_dataframe(criteria={"elasticity": {"$exists": True}, 
                         "elasticity.warnings": [],
                         "elements": {"$all": ["Pb", "Te"]},
                         "e_above_hull": {"$lt": 1e-6}}, # to limit the number of hits for the sake of time
                        properties = ["elasticity.K_VRH", "elasticity.G_VRH", "pretty_formula", 
                                      "e_above_hull", "bandstructure", "dos"])

print("There are {} elastic entries on MP with no warnings that contain "
      "Pb and Te with energy above hull ~ 0.0 eV".format(df['elasticity.K_VRH'].count()))

结果是：

There are 3 elastic entries on MP with no warnings that contain Pb and Te with energy above hull ~ 0.0 eV

1

df.head()

下面来查看一下其中一个材料的能带和DOS图

1
2
3
4
5
6

from pymatgen.electronic_structure.plotter import BSDOSPlotter # 调用 Pymatgen 作图

mpid = 'mp-20740'
idx = df.index[df.index==mpid][0]  # 获取 mp-20740 这个材料所在的行数
plt = BSDOSPlotter().get_plot(bs=df.loc[idx, 'bandstructure'], dos=df.loc[idx, 'dos']);
plt.show()

结尾

• matminer获得的数据(本教程保存在 df 变量中)可以可视化，也可以直接用来做机器学习。有兴趣的同学请自学 pandas 这个 python 库，学习 dataframe 数据的处理方法。

• Pymatgen 做的能带和DOS图不够美观。大家可以留言，我有机会可以跟大家分享修改pymatgen代码来美化这些图。

• 官网的教程是基于 jupyter notebook的。jupyter notebook 的好处是既能运行代码，又能自由地给代码加注释，公式，图片。所得到的文档便于阅读和理解。希望大家也能学会这个python 库。

• 如果想调用 jupyter notebook来运行上述代码，建议大家在具有图形界面的终端上试验。以windows 系统为例，再按照上篇教程创建好 py36 虚拟环境后，需要

  1	pip install jupyter # 具体的使用方法请参考其它网络教程

• 本文还附带jupyter notebook 版本的文档供大家阅读和使用 【链接：https://pan.baidu.com/s/153PcfiM3vaAG71MCCTFuCQ
  提取码：w4us 】

本次的教程到这就结束了。

今天向大家隆重介绍一款材料类的机器学习开源软件：Matminer！机器学习预测材料性质，乃至逆向设计材料是材料领域的热门研究话题。很多同学有心做机器学习方面的研究，奈何不知道有哪些实用的代码，从而使项目一拖再拖。材料类的机器学习软件有很多，我比较喜欢Matminer这一款，因为它是由Pymatgen开发者开发的，非常易于使用。现在就向大家介绍一下这款软件（本教程根据Matminer官网https://hackingmaterials.lbl.gov/matminer/#installing-matminer撰写）

简介

Matminer是用于材料数据挖掘的基于Python的开源软件。它可以从各种数据库获取材料属性数据，将复杂材料属性(如成分、晶体结构、能带结构)表征为与物理相关的特征量，训练机器学习模型，并分析数据挖掘的结果。

Matminer使用pandas数据格式。它还有一个进阶的自动化版本叫Automatminer，可以自动化的训练机器学习模型并得到结果。

Matminer 可以生产可交互的图像。下面的流程图生动展示了matminer的工作内容。Matminer 可以访问Citrine，Materials Project，MDF等数据库，获得材料结构、能带、力学性能等多种性质。这些性质可以被转换成数值化的、可视化的特征量，该特征量可以被用于训练机器学习模型。

安装

Matminer的安装过程非常简单，大家只要安装好Anaconda3 (Python 3.6) 后，在终端输入 pip install matminer 就行了。建议大家一并安装Pymatgen。

下面以天河超算为例，展示matminer的安装过程：

1
2

# 创建基于python 3.6 的虚拟环境，名字叫 py36
(base) [test@ln0%tianhe2 ~]$ conda create -n py36 python=3.6

安装好环境后，激活该环境，并安装pymatgen 和 matminer

1
2
3
4
5
6
7
8
9

# 激活环境用conda activate 命令
(base) [test@ln0%tianhe2 ~]$ conda activate py36

# 括号里的 base 就会变成 py36
(py36) [test@ln0%tianhe2 ~]$ conda install -c conda-forge pymatgen

# 下面这个命令为安装 matminer
(py36) [test@ln0%tianhe2 ~]$ pip install matminer

下面测试安装结果，输入下列命令，如果不出错，就是安装成功了。

1
2
3
4
5
6
7

(py36) [test@ln0%tianhe2 ~]$ python
Python 3.6.12 |Anaconda, Inc.| (default, Sep  8 2020, 23:10:56) 
[GCC 7.3.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import pymatgen
>>> import matminer
>>> 

现在可以尽情地去测试啦！后续教程将陆续分享。

大师兄在9月8日的文章中讲了如何用 Pymatgen 计算离子的电导率（该方法可以用来研究固态电解质等问题），本次，大师兄介绍一下如何使用 Pymatgen 可视化离子的迁移概率密度。

先举个例子，

在“Design principles for solid-state lithium superionic conductors”一文中（Wang et al., Nature Materials 2015, 14 , 1026–1031. ），作者用Ab Initio Molecular Dynamic （AIMD）计算了Li 离子在Li$\mathrm{_1}$$\mathrm{_0}$GeP$\mathrm{_2}$S$\mathrm{_1}$$\mathrm{_2}$， Li$\mathrm{_7}$P$\mathrm{_3}$S$\mathrm{_1}$$\mathrm{_1}$，Li$\mathrm{_2}$S，和 Li$\mathrm{_4}$GeS$\mathrm{_4}$ 四种材料中的迁移概率密度（Probability Density），结果如下图所示：

• 从图中可以看出Li离子在图a所示材料中主要沿c轴方向的通道迁移，而且由于这个通道连通得比较好，Li离子的迁移势垒会比较低（0.22~0.25 eV）。
• Li离子在图b所示的材料的迁移路径形成了一个三维网格，而且由于这个概率密度比图b中的概率密度分布得更加均匀，Li离子的迁移势垒更低（0.18~0.19 eV）。
• 图b所示的材料就完全不行了，因为Li离子的概率密度仅分布在特定的位点附近，说明离子不能有效地移动。
• Li离子在图d所示材料中也存在迁移局域化的行为。
• 作者总结说 “A general principle for the design of Li-ion conductors with low activation energy can be distilled from the above findings: all of the sites within the diffusion network should be energetically close to equivalent, with large channels connecting them.”

那么我们如何在自己的计算中画出这样的图呢？Pymatgen 举手说，它可以帮忙！

但是在开始之前，我们要安装Pymatgen的插件：Pymatgen-diffusion（https://github.com/materialsvirtuallab/pymatgen-diffusion）。

安装 Pymatgen-diffusion

推荐大家使用最新版的Anaconda安装Pymatgen及其插件。点击上面的链接，进入官网后，点击最新版本链接，

我们可以下载.zip文件，

下载完成后，大家可以解压这个文件，得到 pymatgen-diffusion-2019.8.18文件夹。

我们把其中的 pymatgen_diffusion 文件夹放到 Anaconda的site-packages文件夹下，路径是 Windows 系统：……\Anaconda\Lib\site-packages；Linux系统：……/anaconda3/lib/pythonx.x/site-packages，就算安装好了。

接下来我们可以启动python，导入这个模块，如果不报错就没有问题了。

1
2
3
4
5
6

[test@ln0%tianhe2 li_sn_s]$ python
Python 3.8.3 (default, Jul  2 2020, 16:21:59) 
[GCC 7.3.0] :: Anaconda, Inc. on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import pymatgen_diffusion
>>> 

学习用法

我们可以在其github网站上通过例子学习这个模块的用法。

点击打开 probbility_analysis.ipynb 文件。

其内容如下（有所删减）：如果不想看的话可直接查看 开始作图 部分

1
2
3
4
5
6
7
8
9
10
11
12
13
14

from pymatgen.analysis.diffusion_analyzer import DiffusionAnalyzer
from pymatgen_diffusion.aimd.pathway import ProbabilityDensityAnalysis
import json

#ProbabilityDensityAnalysis object
filename="/Users/iekhengchu/repos/pymatgen-diffusion/pymatgen_diffusion/aimd/tests/cNa3PS4_pda.json"

data = json.load(open("../pymatgen_diffusion/aimd/tests/cNa3PS4_pda.json", "r"))
diff_analyzer = DiffusionAnalyzer.from_dict(data) # 初始化DiffusionAnalyzer类

pda = ProbabilityDensityAnalysis.from_diffusion_analyzer(diff_analyzer, interval=0.5, 
                                                         species=("Na", "Li")) #可以指定离子
#Save probability distribution to a CHGCAR-like file
pda.to_chgcar(filename="CHGCAR_new2.vasp") #保存概率密度文件

开始作图

代码（test.py）如下：

1
2
3
4
5
6
7
8
9

from pymatgen_diffusion.aimd.pathway import ProbabilityDensityAnalysis
from pymatgen.core.trajectory import Trajectory
from pymatgen.io.vasp.outputs import Xdatcar
from pymatgen.analysis.diffusion_analyzer import DiffusionAnalyzer

traj = Trajectory.from_file('XDATCAR')
diff = DiffusionAnalyzer.from_structures(traj,'Li',900,2,1)
pda = ProbabilityDensityAnalysis.from_diffusion_analyzer(diff,interval=0.5,species=("Li"))
pda.to_chgcar(filename="pda.vasp") #保存概率密度文件

此处理过程大概耗时8分钟，因机器而异。

在VESTA中可视化如下：

前面一节我们介绍了Ubuntu虚拟机的安装，今天就介绍一些Ubuntu常用软件的软装。大部分极其简单，敲个命令，输入下密码就能解决。少部分非常简单，设置下环境变量就可以解决。有些安装比较蛋疼的可以借助第三方平台比如Anaconda等来安装，也可以归纳到非常简单的行列中。所以，So easy！！！

1 第一件事：

打开终端（CONTROL + ALT + T），创建~/bin 文件夹。

1

qli@bigbrosci:~$ mkdir ~/bin

这个文件夹用于存储一些日常的脚本文件，可执行的小程序。 平时勤备份也不会占用太大空间。

2 命令安装程序

复制下面的sudo apt-get命令，输入到Terminal中，回车，输入密码即可。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

sudo apt-get update
sudo apt-get install python3-pip
sudo apt-get install curl
curl https://bootstrap.pypa.io/get-pip.py --output get-pip.py
sudo python get-pip.py

sudo apt-get install vim 
sudo apt-get install vim-gtk
sudo apt-get install sshfs
sudo apt-get install rename 
sudo apt-get install xclip 
sudo apt-get install git 
sudo apt-get install gnuplot
sudo apt-get install xcrysden
sudo apt-get tree
sudo apt-get install inkscape  

sudo apt-get install fish 
sudo apt-get install tmux 
sudo apt-get install texlive-full
sudo apt-get install kile  

pip 用来安装python程序；

vim就不介绍了；

sshfs 用来把服务器挂载到本地电脑，方便传输数据；

rename 可以很方便批量命名；

xclip 可以把命令的输出复制到剪切板里面，方便鼠标粘贴；

git 用来从git-hub下载，上传，管理一些脚本程序等，

gnuplot 画图用的

xcrysden 可视化软件

tree 方便显示文件目录

inkscape 一款非常好用的做图软件

fish，tmux 装逼效果明显，没对象的重点关注下，好好学习一番，当然功能没的说，很强大。

Latex 写文章的可以使用最后两个命令安装。

其他的暂时也想不起来，如果你有推荐的，可以发邮件给我，后面再添上去。（lqcata@gmail.com）

3 安装Anaconda

它可以帮助你节省很多时间和精力去安装软件。下载后，按照下面命令，然后一路狂摁回车或者输入yes即可，官网也有安装教程，这里就不累赘了。通过Anaconda安装其他软件，可以参考前面的教程：https://www.bigbrosci.com/2020/09/19/A20/

1
2
3
4
5
6
7
8
9
10

qli@bigbrosci:~/Downloads$ ls
Anaconda3-2020.07-Linux-x86_64.sh  p4vasp-0.3.30.tgz 
qli@bigbrosci:~/Downloads$ bash Anaconda3-2020.07-Linux-x86_64.sh 
Welcome to Anaconda3 2020.07

In order to continue the installation process, please review the license
agreement.
Please, press ENTER to continue

>>> 

Anaconda安装的软件可以通过https://anaconda.org 来查询，都会列出来安装的命令，比如RDkit：

1

conda install -c rdkit rdkit

另一个需要注意的是，Anaconda安装后，最喜爱的p4vasp不能正常使用了。解决办法如下：

1） 找到p4v的目录，一般在/usr/bin/

2）将p4v中的python 全部改为python2

3. 将p4v.py 第一行中的python改为python2。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

(base) qli@bigbrosci:~$ p4v 
  File "/usr/bin/p4v.py", line 9
    print """You need to get version 2.0 (or later) of PyGTK for this to work. You
             can get source code from http://www.pygtk.org """
          ^
SyntaxError: invalid syntax
(base) qli@bigbrosci:~$ which p4v 
/usr/bin/p4v
(base) qli@bigbrosci:~$ sudo vi /usr/bin/p4v
[sudo] password for qli: 
(base) qli@bigbrosci:~$ sudo vi /usr/bin/p4v.py 
(base) qli@bigbrosci:~$ grep python  /usr/bin/p4v
export PYTHONPATH=$PYTHONPATH:/usr/lib/python2.7/dist-packages
exec python2 /usr/bin/p4v.py "$@"
(base) qli@bigbrosci:~$ grep python  /usr/bin/p4v.py 
#!/usr/bin/env python2

4 其他可执行的程序

这里我们以VESTA为例，这种方法非常简单实用。

4.1） 下载解压VESTA程序包。（https://jp-minerals.org/vesta/en/） 直接通过官网下载即可，没必要去群里面求助，也没必要往群里面上传。

4.2） 本人习惯性地将其解压到/home/qli/Documents/VESTA-gtk3

4.3) 设置环境变量，将下面的export 两行复制到~/.bashrc 文件中， source 下 .bashrc 文件，然后直接敲命令VESTA即可打开可视化界面。

1
2
3
4
5
6
7
8
9
10
11
12
13
14

qli@bigbrosci:~/Documents/VESTA-gtk3$ ls
asfdc           Library_License  RIETAN.sh          style         VESTA-core
bvparm2016.cif  libVESTA.so      spgra.dat          style.ini     VESTA_def.ini
elements        LICENSE          spgro.dat          template.ins  VESTA-gui
elements.ini    mspgr.dat        STRUCTURE_TIDY     tmp           VESTA.ini
img             PowderPlot       STRUCTURE_TIDY.sh  VESTA         wyckoff.dat
qli@bigbrosci:~/Documents/VESTA-gtk3$ pwd
/home/qli/Documents/VESTA-gtk3
qli@bigbrosci:~/Documents/VESTA-gtk3$ cd 
qli@bigbrosci:~$ vi .bashrc
export VESTA='/home/qli/Documents/VESTA-gtk3'
export PATH=$PATH:$VESTA
qli@bigbrosci:~$ source .bashrc	 
qli@bigbrosci:~$ VESTA

注意 上面VESTA的路径记得换成你自己电脑里面的。

5 Ubuntu的备份

本人喜欢瞎折腾，也经常会把电脑搞崩溃，备份就显得尤其重要。对于Ubuntu来说，前面我们在虚拟机的安装中提到存储Ubuntu系统的vdi文件，如果你有足够大的硬盘，可以定期将一个稳定版本的vdi保存起来，万一系统崩溃，删掉坏的，导入最新的备份即可。

如果没有用虚拟机，直接装的Ubuntu系统，或者双系统，可以尝试定期保存~/bin 和~/.bashrc 文件，本人之前就是这样做的（现在懒得折腾，系统不坏了），它们不怎么占空间，打包压缩下，邮件基本就可以上传。如果系统坏掉了，非得安装一个新的才能解决的时候，只需要重复本教程，下载一些常用的软件，以及用备份的~/bin 和~/.bashrc 文件还原一下设置即可。数据一般来说都直接存储在服务器中，相关备份不在本教程讨论范围之内。

当然还有其他的办法备份和恢复虚拟机，比如创建当前Ubuntu版本的iso文件，定期备份当前整个系统以便引导修复等等，精力有限，不再累赘。有兴趣的可以尝试。

6 打赏

如果感觉本文对你的相关研究有帮助，欢迎打赏，支持作者的热心付出。如果你也有自己擅长的操作等，热烈欢迎无私分享，可以通过QQ（122103465）或者邮件（lqcata@gmail.com）联系。

为方便新手入坑，这里我们先介绍一下Ubuntu虚拟机的安装细节。VMWARE和VirtualBox是两个主要的选择，前者收费，后者免费。当然介绍免费的喽。 本节主要把安装的过程展示一遍，然后把需要注意的细节强调一下。安装虚拟机要量力而行，优先在台式机上安装虚拟机。如果你的电脑配置很普通的话，运行起来可能不会很流畅。

第一步： 准备工作

1 下载Ubuntu的安装包

https://ubuntu.com/download/desktop

2 下载Virtual Box，并安装

https://www.virtualbox.org/wiki/Downloads

3 下载Virtual Box的扩展包并安装（Extension Pack）

第二步：安装虚拟机(初始化)

2.1） 创建新的虚拟机

2.2）设置虚拟机的名字，名字里面只要你输入Ubuntu字样，后面的Type和Version会自动更新，Machine Folder 最好选择一个硬盘比较大的盘。

2.3） 设置虚拟机的内存，本人喜欢设置为host的一半。

2.4）创建虚拟机，按照图片选择即可。

2.5） 设置虚拟机的vdi文件的目录，如果前面选择好了，这个采用默认即可。虚拟机的大小：如果使用Ubuntu比较频繁，安装的软件、数据处理等比较多的话，根据自己的情况酌情增加下。但是尽量别太小，到时候存储不够容易导致死机崩溃。

另外一个需要注意的是：Ubuntu的系统以及后面你安装的软件，下载的东西，写的脚本等都会存在这个vdi文件里面，如果你换了一台电脑，还想继续使用虚拟机，大可不用再重复本教程的操作，直接导入vdi即可（前面2.4步骤的第一个图片）。

2.6） 经过前面步骤，虚拟机初始化基本完成了（界面上可以看到），按照箭头，点Settings。

2.7) System设置处理器的数目，同样也是host的一半。太小Ubuntu容易卡死，太大host容易卡死。

2.8）指定Ubuntu安装iso文件的目录，点Storage –> 光盘 Empty –> 右侧光盘图标

2.9） 在弹出的对话框，选择前面下载好的Ubuntu安装文件。

3 安装虚拟机 （Ubuntu的安装）

3.1） 点击下方图片的Start绿色图表，开始启动安装

3.2） 选择Install Ubuntu，后面的一路Continue

3.3） 设置键盘类型

3.4）按照下面的选择，Continue

3.5) 设置用户名，计算机名，以及密码。

3.6） Ubuntu 终于开始安装了，电脑配置好的别走开，配置不好的一边凉快等着去。（PS：配置不好不建议虚拟机，卡的你怀疑人生）

3.7） 安装完成后，会提示让你先移除安装的ISO文件，然后再点ENTER

3.8） 这里本人把安装包改了下名字，然后点ENTER。

3.9) Ubuntu 初步安装完成。

4 后续设置

前面只是顺利安装成功，但距离真正可以使用，还有几个设置需要继续操作下。

4.1）安装增强扩展包。（前面准备工作中，我们在host中安装了扩展包，Ubuntu里面也需要安装一下）

最明显的效果：

安装增强包前，不管你怎么拖放VirtualBox的界面，Ubuntu的界面总是那么一点。

安装后：Ubuntu的界面可以随意调整。（重启后才会更新效果）

4.2）设置虚拟机的网络。

Setting –> Network，然后按照下面的设置即可。

4.2) 设置剪切板共享：host复制的文字可以直接粘贴在Ubuntu里，反之亦然。(重启后生效)

4.3） 设置共享文件夹。

共享文件夹可以让你在host和Ubuntu之间及其方便地进行数据分享。

4.3.1 Virtual Box界面，点Settings –> Shared Folders，然后按照下面图中设置即可。

Folder Path为host的文件夹目录

Mount Point 为虚拟机中所挂在的位置。

点OK， 然后我们进入到虚拟机中，会发现一个类似硬盘的图标，名为：Shared，点一下，发现没有权限。

解决办法

在Ubuntu中，打开终端Terminal，输入命令。

1
2
3
4
5

qli@bigbro:~/Desktop$ sudo adduser $USER vboxsf 
[sudo] password for qli: 
Adding user 'qli' to group 'vboxsf' ...
Adding user qli to group vboxsf
Done.

完事后，重启虚拟机。如果不能启动，显示USB相关的报错信息。取消USB的选项，再次启动即可。

到这里，虚拟机的安装，设置基本就完成了，剩下的就是在Ubuntu里面安装一些常用的软件程序了。这个会在后续教程里面简单介绍。

打赏

如果感觉本文对你的相关研究有帮助，欢迎打赏，支持作者的热心付出。如果你遇到相关的问题，或者也有自己擅长的操作，热烈欢迎无私分享，可以通过QQ（122103465）或者邮件（lqcata@gmail.com）联系。

P4vasp是一款基于Python2的，用于可视化，搭建结构的软件。由于Python2已经被时代所淘汰，p4vasp官网也不见任何的针对python3的更新改进，这个软件估计要凉，跟我们再见了。但Ubuntu20安装也不是很难，能撑到博士毕业应该问题不大。经过多次的测试，下面是p4vasp的安装过程，非常简单。

下载并解压

下载链接：http://www.p4vasp.at/#/download

这里本人把安装包解压缩到了~/Documents/ 目录下。

安装步骤-1

1
2
3
4
5
6
7
8
9
10
11

qli@bigbrosci:~/Documents/p4vasp-0.3.30$ ls
BUGS           install           odpdom                setenvironment.sh
ChangeLog      install-local.sh  p4v                   src
data           install.sh        p4vasp-0.3.30-1.spec  test
diagnostic.py  lib               p4vasp.log            uninstall.sh
doc            LICENSE           p4v.py                utils
ext            Makefile          README                vinfo.py
FAQS           Makefile.MacOS    README.MacOS          vinfo.pyc
qli@bigbrosci:~/Documents/p4vasp-0.3.30$ sudo apt-get install gcc make
qli@bigbrosci:~/Documents/p4vasp-0.3.30$ sudo apt-get install make
qli@bigbrosci:~/Documents/p4vasp-0.3.30$ sudo apt-get install python2-dev

安装步骤-2：

http://archive.ubuntu.com/ubuntu/pool/universe/p/pygtk/

从这个网站下载：python-gtk2 和 python-glade2的deb文件，下载后，双击即可安装。

注意：先安装gtk2，再安装glade2。

安装步骤-3：

1
2
3
4

qli@bigbrosci:~/Documents/p4vasp-0.3.30$ sudo apt-get install libxft-dev
qli@bigbrosci:~/Documents/p4vasp-0.3.30$ sudo apt-get install libglfw3-dev libgl1-mesa-dev libglu1-mesa-dev
qli@bigbrosci:~/Documents/p4vasp-0.3.30$ bash install/install-ubuntu-dependencies.sh
qli@bigbrosci:~/Documents/p4vasp-0.3.30$ bash install.sh

注意：

bash install/install-ubuntu-dependencies.sh 这一步最后可能会出python-epydoc的报错，可以忽略。

反馈：

上面的命令是Ubuntu20里面安装用的，如果安装过程出现问题，或者有更好的安装方法，欢迎交流讨论：

QQ：122103465

E-mail：lqcata@gmail.com

至于其他的Linux系统，本人没有测试，如果你有比较好的安装方法，也欢迎分享讨论。

打赏

如果感觉本文对你的相关研究有帮助，欢迎打赏，支持作者的热心付出。如果你也有自己的骚操作，热烈欢迎无私分享，可以通过QQ（122103465）或者邮件（lqcata@gmail.com）联系。

今天我们介绍ASE的另一个骚操作：结合Openbabel自动搭分子结构。与其说是ASE的骚操作，不如把这个骚字放在Openbabel头上。Openbabel是一款功能极其强大的结构转化工具。如果你在用它，我就不在这里啰嗦了。如果你没听过，可以通过这个链接来初步了解一下。（http://openbabel.org/wiki/Main_Page）

鉴于Openbabel的安装比较蛋疼，先送一波福利（Mac，Linux用户）。下面是本人通过Anaconda（已经路转粉）创建虚拟环境，并且安装这些软件的一些命令。通过它们可以让你安全顺利，麻溜滴地安装Openbabel, RDkit, Ase，和pymatgen。

1
2
3
4
5
6

conda create --name qrobot
conda activate qrobot
conda install -c conda-forge openbabel
conda install -c conda-forge rdkit
conda install -c conda-forge ase
conda install --channel conda-forge pymatgen

RDkit也是一款功能及其强大地软件，尤其是化学信息学，生物化学相关的领域。如果不了解，点击这个官网学习。（https://www.rdkit.org/）顺便吐槽一下，百度里面你搜索一堆RDkit的使用，大部分都是把官网例子复制过来，然后贴个二维码要钱的。交不交智商税决定权在你手上，如果手痒，请跳到文末最后，把钱打赏给我。

废话不多说，转到我们今天的主题： SMILES to XYZ。又是一个新知识，什么是SMILES？

SMILES是simplified molecular-input line-entry system 的缩写，旨在用ASCII字符串来标识，描述分子的结构。如果你不了解，可以通过下面两个链接初步了解一下：

https://en.wikipedia.org/wiki/Simplified_molecular-input_line-entry_system

https://www.daylight.com/dayhtml/doc/theory/theory.smiles.html

通过SMILES我们可以将一个分子用字符串来标识出来，比如我们后面练习中的：

乙醇（CCO），甲烷（C），甲基（[CH3]）, 苯（c1ccccc1）。

所以，在操作之前，你需要知道所搭结构的SMILES字符串。通过上面两个链接，耐心钻研一下，应该不到1天就可以搞定大部分的分子。将SMILES转化为3D的结构，可以通过Openbabel，也可以通过RDkit。今天我们就介绍前者。下面是脚本内容（smiles_to_xyz.py）以及使用的效果。最后，我们再详细解释这个脚本是怎么工作的。

脚本内容

下载链接：https://github.com/BigBroSci-LVTHW/LVTHW/tree/master/source/_posts/A20

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

#!/usr/bin/env python3 
# -*- coding: utf-8 -*-
"""
convert SMILES to xyz
"""
import sys
from openbabel import pybel
import numpy as np
import ase
import ase.io
from ase import Atoms
from ase.constraints import FixAtoms
 
smiles, out_name = sys.argv[1:3]

'''
Introduction to SMILES:
https://en.wikipedia.org/wiki/Simplified_molecular-input_line-entry_system
https://www.daylight.com/dayhtml/doc/theory/theory.smiles.html
'''

## Use openbabel to convert SMILES to xyz. 
mol = pybel.readstring("smi", smiles)
mol.make3D(forcefield='mmff94', steps=100)
# mol.write("xyz", filename=out_name+'_pybel.xyz', overwrite=True)

# USE ase to make POSCAR
''' https://wiki.fysik.dtu.dk/ase/ase/atoms.html'''
geo = Atoms()
geo.set_cell(np.array([[16.0, 0.0, 0.0],[0.0, 17.0, 0.0],[0.0, 0.0, 18.0]]))
geo.set_pbc((True,True,True))
for atom in mol:
    atom_type = atom.atomicnum
    atom_position = np.array([float(i) for i in atom.coords])
    geo.append(atom_type)
    geo.positions[-1] = atom_position
geo.center()

'''https://wiki.fysik.dtu.dk/ase/ase/constraints.html'''
c = FixAtoms(indices=[atom.index for atom in geo if atom.symbol == 'XX'])
geo.set_constraint(c)

ase.io.write(out_name + '_ase.xyz', geo, format='xyz')
ase.io.write(out_name + '_POSCAR', geo, format='vasp', vasp5='True')

脚本运行

1
2
3
4
5
6

 (qrobot) qli@MacBook-Pro test_openbabel % python3 smiles_to_xyz.py '[CH3]' Methyl
 (qrobot) qli@MacBook-Pro test_openbabel % python3 smiles_to_xyz.py 'CCO' ET
 (qrobot) qli@MacBook-Pro test_openbabel % python3 smiles_to_xyz.py 'C' Methane
 (qrobot) qli@MacBook-Pro test_openbabel % python3 smiles_to_xyz.py 'c1ccccc1' BZ
 (qrobot) qli@MacBook-Pro test_openbabel % ls
BZ_POSCAR         BZ_ase.xyz        ET_POSCAR         ET_ase.xyz        Methane_POSCAR    Methane_ase.xyz   Methyl_POSCAR     Methyl_ase.xyz    smiles_to_xyz.py

结果展示

内容详解

1. Openbabel部分其实就三行代码，非常简单。

1
2
3

1 mol = pybel.readstring("smi", smiles)
2 mol.make3D(forcefield='mmff94', steps=100)
3 # mol.write("xyz", filename=out_name+'_pybel.xyz', overwrite=True)

1. 读取SMILES，并创建一个mol的对象（object）。

2. make3D 通过MMFF94力场对生成的3D结构优化了100步。

3. 保存结构，这里我把它注释掉了，因为后面我们要通过ASE来生成xyz和POSCAR。

2. ASE部分，代码多一些，可能本人水平有限吧。不管怎么样，这里主要告诉大家的一点就是：我们可以将Openbabel和ASE无缝连接起来。同样，你也可以把ASE与你自己写的一些东西关联，这样就可以调用ASE的功能来快速实现一些小目标了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

1 geo = Atoms()
2 geo.set_cell(np.array([[16.0, 0.0, 0.0],[0.0, 17.0, 0.0],[0.0, 0.0, 18.0]]))
3 geo.set_pbc((True,True,True))
4 for atom in mol:
5    atom_type = atom.atomicnum
6    atom_position = np.array([float(i) for i in atom.coords])
7    geo.append(atom_type)
8    geo.positions[-1] = atom_position
9 geo.center()

'''https://wiki.fysik.dtu.dk/ase/ase/constraints.html'''
10 c = FixAtoms(indices=[atom.index for atom in geo if atom.symbol == 'XX'])
11 geo.set_constraint(c)

12 ase.io.write(out_name + '_ase.xyz', geo, format='xyz')
13 ase.io.write(out_name + '_POSCAR', geo, format='vasp', vasp5='True')

1）创建一个ASE的Atom空对象；

2）设置放分子的格子大小，这里我们用的是16x17x18 $\AA^3$的。如果嫌小，自己改大一下就可以了。

3）设置周期性；

4-8）将Openbabel的mol对象中的原子添加到ASE的Atom对象里面。

9）将分子放到格子的中心。

10-11）固定或者放开原子，这里我们打算放开所有的原子，所以用到了一个XX 。如果原子符号是XX，那么就固定住。其实世界上根本就没有XX原子，因此所有的原子就会被放开了。

12-13） 保存xyz和VASP的POSCAR。

小结

通过SMILES来搭建结构非常方便，当然也可以通过数据库下载，可视化程序例如GaussianView等自己手动搭建。

值得注意的是，不管通过什么方式搭建的结构，提交任务之前，都要尽可能保证结构的合理性。所以不要搭建完结构就立刻提交任务，先认真检查一遍，没有任何问题之后再提交。

打赏

如果感觉本文对你的相关研究有帮助，欢迎打赏，支持作者的热心付出。如果你也有自己的骚操作，热烈欢迎无私分享，可以通过QQ（122103465）或者邮件（lqcata@gmail.com）联系。

本节介绍一下ASE的另一个骚操作—扩胞。说到扩胞，这个大家都不陌生，可以有很多种方式来实现，本教程中就介绍了p4vasp的一种操作，今天就看一下ASE的骚操作吧。废话不多说，直接上例子和脚本。

例子：

本例子中，主要把优化好的一个(1x1)-Ru(0001)的slab 扩胞为（4x4）的slab。为什么这样做呢？因为直接优化刚刚切好的（4x4)slab会有些费时间。先优化（1x1)的，在CONTCAR的基础上扩展为（4x4）的，再优化会相对来说快一些。但现在大家服务器一般都很给力，从刚切好的（1x1) slab直接扩展到（4x4），然后再优化也可以的。不管怎么样，扩胞这个操作都是必须的。如果你硬要抬扛说，我可以直接从bulk切出来（4x4)的，不用扩胞这么麻烦。OK，这也是可以的。但等你需要扩胞操作的时候，记得回来看这个骚操作就行。

先展示一下效果吧：

1
2
3
4
5
6
7

qli@bigbro:~/Desktop/A19$ ls
CONTCAR  expand.py 
qli@bigbro:~/Desktop/A19$ python3 expand.py  CONTCAR  4 4 1 
qli@bigbro:~/Desktop/A19$ ls
CONTCAR  expand.py  POSCAR
qli@bigbro:~/Desktop/A19$ ase gui "-R -90x" CONTCAR 
qli@bigbro:~/Desktop/A19$ ase gui "-R -90x" POSCAR 

扩胞前：（CONTCAR）

扩胞后：（POSCAR）

脚本：

脚本可以通过git-hub下载，链接：

https://github.com/BigBroSci-LVTHW/LVTHW/tree/master/source/_posts/A19

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Oct  2 10:54:36 2019

@author: qli
"""
import sys, os
import ase.io.vasp

file_read = sys.argv[1]
x,y,z     = [int(i) for i in sys.argv[2:5]]
try:
    cell = ase.io.vasp.read_vasp("CONTCAR")
    ase.io.vasp.write_vasp("POSCAR",cell*(x, y, z), direct=True,sort=True)
except:
    print(os.getcwd())

前面的扩胞命令：expand.py CONTCAR 4 4 1

CONTCAR 对应的为 sys.argv[1]

4 4 1 分别为在x,y,z 三个方向扩胞的倍数。

打赏

如果感觉本文对你的相关研究有帮助，欢迎打赏，支持作者的热心付出。如果你也有自己的骚操作，热烈欢迎无私分享，可以通过QQ（122103465）或者邮件（lqcata@gmail.com）联系。

ASE 在 DFT 计算中的骚操作令人印象深刻。类似 ASE 的科研工具还有很多，比如 VASPKIT 和 QVASP。熟练使用这些科研工具肯定会助力大家快速完成科研任务，发出更好更高质量的文章。

此次，本文向大家隆重介绍另一款出色的科研工具：Pymatgen。

Pymatgen 是 python materials genomics 的缩写，它是一款基于 python 的、开源的、强大的材料分析软件（https://pymatgen.org/）。

Pymatgen 包含一系列能够表示元素（Element）、位点（Site）、分子（Molecule）、和结构（Structure）的类（Class）。它具有为很多计算软件提供前处理和后处理的能力。这些计算软件包括VASP，ABINIT，exciting，FEFF，QCHEM，LAMMPS，ADF，AIIDA，ASE，Gaussian，Lobster，Phonopy，Shengbte，Pwscf，和Zeo++等等。它能实现科研狗的众多后处理需求，包括生成相图（Phase diagram）和布拜图（Pourbaix diagrams），分析态密度和能带等等。

Pymatgen 还提供了很多数据库（Materials Project REST API，Crystallography Open Database，and other external data sources）的接口，方便大家从数据库中查询结构和其他数据。

真是科研狗快乐科研之利器呀！

以下是Pymatgen官网提供的后处理的例子：

Top: (left) Phase and (right) Pourbaix diagram from the Materials API. Bottom left: Calculated bandstructure plot using pymatgen’s parsing and plotting utilities. Bottom right: Arrhenius plot using pymatgen’s DiffusionAnalyzer.

此次，本文就介绍一下如何使用 Pymatgen 的 DiffusionAnalyzer 类去计算锂离子固态电解质中锂离子电导率。

计算离子电导率的理论与公式

目前，比较准确的计算离子电导率的方法是先用NVT系综第一性原理分子动力学（AIMD，ab initio molecular dynamics）模拟材料中离子在不同温度下的运动，然后计算出离子的平均（average）均方位移（MSD，mean square displacement），再计算出自扩散系数（D$_s$，self-diffusion coefficient），最后求得离子在某温度下的电导率（$\sigma$，conductivity）。

如何进行AIMD计算

AIMD计算通常非常耗时，所以，为了减少计算成本，我们可以适当放宽计算精度。如果用 VASP 进行计算，具体的，大家可以

• 采用较小的截断能。氧化物用 400 eV，硫化物用 280 eV，硒化物用 270 eV
• 采用Gamma点作为K点设置，并使用gam版本的 VASP 进行计算
• 采用单胞计算，如果材料的单胞包含比较多的原子
• 采用合适的步长，比如2 fs，即 POTIM = 2

后处理的基本公式

一旦AIMD计算完成，大家就可以着手计算离子电导率了。本文首先先介绍以下计算过程中使用的公式，方便有兴趣的同学自己开发脚本。

平均均方位移（average MSD）可以通过以下公式计算：
$$
averageMSD=\langle[, \mathbf r(t)]^2,\rangle=\frac{1}{N}\sum_{i}^{N}\langle[, \mathbf r_i(t+t_0)],^2 - [, \mathbf r_i(t_0)],^2\rangle
$$
$\mathbf r_i(t)$ 是第 $i$ 个离子在 $t$ 时刻的位移。

自扩散系数（$D_s$）可以通过以下公式计算：
$$
D_s=\frac{averageMSD}{2dt}
$$
$d$ 是离子在材料中的扩散维度（一般地，$d=3$），$t$ 是离子扩散的时间。

最后，离子电导率（$\sigma$）可以这样计算：
$$
\sigma=\frac{ne^2Z^2}{k_BT}D_s
$$
$n$ 是材料中的离子密度，$e$ 是元电荷，$Z$ 是离子的价态，$k_B$ 是玻尔兹曼常数，$T$ 是温度。

电导率计算的例子

现在我们通过一个 Li_Sn_S 材料的例子来详细了解一下整个计算和处理的过程。该材料的结构显示如下：

本例中采用单胞做计算，INCAR 设置如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

[test@ln0%tianhe2 li_sn_s]$ vi INCAR

ISTART = 0
ICHARG = 2
IBRION = 0
ISIF = 2
NPAR = 8
NSW = 30000
TEBEG = 900 #还要设置成 1500K 等等
PREC = N
POTIM = 2
SMASS = 0.0
NELMIN = 4
LWAVE = F
LCHARG = F
IALGO = 48
LREAL = A

AIMD 计算结束之后会得到 XDATCAR 文件。很多时候，由于超算的时间限制，一个完整的AIMD计算需要提交两三次，从而产生两三个 XDATCAR 文件，这时，我们只要把它们按顺序通过 cat 命令合并在一起就行。例如我们有三个 XDATCAR 文件，分别命名成 XDATCAR01，XDATCAR02，和 XDATCAR03。

1

[test@ln0%tianhe2 li_sn_s]$ cat XDATCAR01 XDATCAR02 XDATCAR03 > XDATCAR 

新得到的XDATCAR文件，注意删掉重复的与晶格信息相关的行，一般续算的次数也不多，在使用上面命令的时候，手动把XDATCAR02, XDATCAR03 中的删除即可。

Pymatgen 大显身手

安装pymatgen

首先让我们安装 pyamtgen，推荐大家参考官网，使用 anaconda 安装，否则会出现问题。安装好了anaconda之后，不管是 linux 还是 windows, 安装 pyamtgen 的指令是一样的。下面以吕梁天河超算为例：

1

[test@ln0%tianhe2 li_sn_s]$ conda install --channel conda-forge pymatgen

安装完成后，我们可以试着运行 python，导入 Pyamtgen 模块，如果像下面一样没有出错，就是安装成功了。

1
2
3
4
5
6

[test@ln0%tianhe2 li_sn_s]$ python
Python 3.7.3 (default, Mar 27 2019, 22:11:17) 
[GCC 7.3.0] :: Anaconda, Inc. on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import pymatgen
>>> 

查看 DiffusionAnalyzer 的类

大家可以通过官方文档（https://pymatgen.org/pymatgen.analysis.diffusion_analyzer.html）查看接下来要使用的类，熟悉一下代码的用法。

1
2
3
4

class DiffusionAnalyzer(MSONable):
    def __init__(self, structure, displacements, specie, temperature,
                 time_step, step_skip, smoothed="max", min_obs=30,
                 avg_nsteps=1000, lattices=None):

这段代码显示，运行这个类需要一系列的输入信息，包括材料结构（structure），位移（displacements），要研究的离子（specie），温度（temperature）等等。

但是这个类提供了很多方法让大家可以通过读取 XDATCAR 或者 vasprun 文件的方式来实例化，例如

1
2
3
4
5
6
7
8
9
10
11
12
13

@classmethod
    def from_structures(cls, structures, specie, temperature,
                        time_step, step_skip, initial_disp=None,
                        initial_structure=None, **kwargs):
        """
        Convenient constructor that takes in a list of Structure objects to
        perform diffusion analysis.
        Args:
            structures ([Structure]): list of Structure objects (must be
                ordered in sequence of run). E.g., you may have performed
                sequential VASP runs to obtain sufficient statistics.
        ... ...
        """

好了，废话不多说，直接上代码，开始进行后处理。

代码示例

新建一个文件，名字为li_conductivity.py

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

'''
分析AIMD结果，计算MSD 和 conductivity
'''
import os
from pymatgen.core.trajectory import Trajectory
from pymatgen.io.vasp.outputs import Xdatcar
from pymatgen import Structure
from pymatgen.analysis.diffusion_analyzer import DiffusionAnalyzer
import numpy as np
import pickle

# 这一步是读取 XDATCAR，得到一系列结构信息
traj = Trajectory.from_file('XDATCAR')

# 这一步是实例化 DiffusionAnalyzer 的类
# 并用 from_structures 方法初始化这个类； 900 是温度，2 是POTIM 的值，1是间隔步数
# 间隔步数（step_skip）不太容易理解，但是根据官方教程:
# dt = timesteps * self.time_step * self.step_skip

diff = DiffusionAnalyzer.from_structures(traj,'Li',900,2,1)

# 可以用内置的 plot_msd 方法画出 MSD 图像
# 有些终端不能显示图像，这时候可以调用 export_msdt() 方法，得到数据后再自己作图
diff.plot_msd()

# 接下来直接得到 离子迁移率， 单位是 mS/cm
C = diff.conductivity

with open('result.dat','w') as f:
    f.write('# AIMD result for Li-ion\n')
    f.write('temp\tconductivity\n')
    f.write('%d\t%.2f\n' %(900,C))

在终端运行该文件

1

[test@ln0%tianhe2 li_sn_s]$ python li_conductivity.py

一段时间后就会得到MSD图像和离子电导率

1
2
3
4
5

[test@ln0%tianhe2 li_sn_s]$ vi result.dat

# AIMD result for Li-ion
temp	conductivity
900	884.05

可见，该材料在 900K 时的锂离子电导率为 884.05 mS/cm。

例子下载：

链接：https://pan.baidu.com/s/1WGzOVJBoe6Ym8mvR1uWanA
提取码：jhc5

思考

• 简短几行代码就可以计算出离子电导率，那么如何得出材料在300K下的电导率呢？
• 如何计算离子在材料中的迁移势垒？
• 如何可视化离子在材料中的扩散路径？

本文只是浅谈离子电导率的计算，欢迎大家指出计算过程中的不足之处。

打赏

如果感觉本文对你的相关研究有帮助，欢迎打赏，支持作者的热心付出。如果你也有自己的骚操作，热烈欢迎无私分享，可以通过QQ联系（122103465）

对科研狗们来说，ASE不是Automotive Service Excellence的缩写，而是Atomic Simulation Environment。
是一款基于Python语言的工具软件，可以方便地处理DFT计算中结构搭建，准备输入文件准备，读取输出文件，既可以可视化，又可以用来转换文件格式。本文主要介绍ASE在DFT计算中的骚操作。这个”骚”字让我想起来小学的时候，同桌读了屈原的《离骚》，然后转过头来用这个字来笑话我，说我真骚。我被憋的满脸通红却不知如何反驳。骚就骚吧，骚前面加个离可以千古流传，那么我相信，它后面加个操作也可以对大家日常的计算有所帮助。
首先声明：本教程具有自动过滤功能。每次过滤的时候，我都会提醒一下，以防你被滤纸卡住。但如果你读完教程却不能顺利进行操作，请不要找我，请分解自己让自己滤下去。

滤纸-1： ASE成功安装了没？ （https://wiki.fysik.dtu.dk/ase/install.html）被卡住就不要继续喽！

今天我们讲的骚操作就是用ASE来进行格式转换。这也是本人日常用的最多的一个，最稳定也最方便的一个操作。既然是操作，那么就举个栗子吧：

例子1）：将mol文件转化为xyz。

首先在Chemspider数据库中，下载一个乙烷3D构型的mol文件。

滤纸-2：Chemspider 是啥玩意，怎么下载mol文件？ 之前教程有讲过，不会的找百度。被卡住就不要继续喽！

操作过程：

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

qli@bigbro:~/Downloads$ ls
6084.mol
qli@bigbro:~/Downloads$ head -n 13 6084.mol 
6324
  Marvin  12300703363D          

  8  7  0  0  0  0            999 V2000
    0.7686    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.7686   -0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    1.1430    1.0253    0.0000 H   0  0  0  0  0  0  0  0  0  0  0  0
    1.1430   -0.5127    0.8880 H   0  0  0  0  0  0  0  0  0  0  0  0
    1.1430   -0.5127   -0.8880 H   0  0  0  0  0  0  0  0  0  0  0  0
   -1.1430   -1.0253    0.0000 H   0  0  0  0  0  0  0  0  0  0  0  0
   -1.1430    0.5127   -0.8880 H   0  0  0  0  0  0  0  0  0  0  0  0
   -1.1430    0.5127    0.8880 H   0  0  0  0  0  0  0  0  0  0  0  0
  1  2  1  0  0  0  0
qli@bigbro:~/Downloads$ ase gui 6084.mol  -o C2H6.xyz 
/home/qli/.local/lib/python3.8/site-packages/ase/gui/ag.py:85: UserWarning: You should be using "ase convert ..." instead!
  warnings.warn('You should be using "ase convert ..." instead!')
(qrobot) qli@bigbro:~/Downloads$ cat C2H6.xyz 
8
Properties=species:S:1:pos:R:3 pbc="F F F"
C        0.76860000       0.00000000       0.00000000
C       -0.76860000      -0.00000000       0.00000000
H        1.14300000       1.02530000       0.00000000
H        1.14300000      -0.51270000       0.88800000
H        1.14300000      -0.51270000      -0.88800000
H       -1.14300000      -1.02530000       0.00000000
H       -1.14300000       0.51270000      -0.88800000
H       -1.14300000       0.51270000       0.88800000
qli@bigbro:~/Downloads$ ase gui C2H6.xyz 
qli@bigbro:~/Downloads$ 

我们成功用ase gui 将mol转化为了xyz文件。

1） 每次运行都会出现那行警示信息，让你用ase convert转化。命令如下：

ase convert -i mol -o xyz 6084.mol C2H6.xyz

这个命令比较长，输入起来烦人，远远不如ase gui 干脆利落。警告信息就用你那高度近视镜自行滤过吧。

2. -o 中的字母o是out的缩写，意思是输出。

3） ase gui 后面跟一个文件，可以直接可视化。

学会了ase gui的命令操作，下面我们需要进一步分解自己，以便透过更多的滤纸。

滤纸-3：

DFT计算中有很多的文件格式，VASP有POSCAR(CONTCAR)，Gaussian有gjf或者com，Material Studio导出来可以有cif啥的。这些格式之间转换着实让一些新手头疼。以至于经常在群里看到有人询问类似的问题。如果你没有用文本编辑器打开这些格式的文件，并认真分析下它们的数据结构，那么你被滤纸卡住了。请把自己变小一些再继续阅读。

分析一下上面的操作：ase gui 6084.mol -o C2H6.xyz

6084.mol 是我们读取的对象，C2H6.xyz 是输出的文件。回到大家常用的VASP计算上来，cif文件文件都不陌生，那么我们就转化一下它吧。

例子-2： 将CIF转化为XYZ，POSCAR

首先从Materials Project 下载一个干冰的cif文件。

滤纸-4：如果你不知道啥是MP。打开这个网址：https://materialsproject.org，注册，随便下载个cif文件。

操作过程：

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
46
47
48
49
50
51
52
53
54
55

qli@bigbro:~/Downloads$ ls
6084.mol  C2H6.xyz  CO2_mp-20066_conventional_standard.cif
qli@bigbro:~/Downloads$ head -n 5 CO2_mp-20066_conventional_standard.cif 

# generated using pymatgen

data_CO2
_symmetry_space_group_name_H-M   'P 1'
_cell_length_a   5.80269900
_cell_length_b   5.80269900
qli@bigbro:~/Downloads$ ase gui CO2_mp-20066_conventional_standard.cif -o CO2_mp-20066.xyz 
/home/qli/.local/lib/python3.8/site-packages/ase/gui/ag.py:85: UserWarning: You should be using "ase convert ..." instead!
  warnings.warn('You should be using "ase convert ..." instead!')
qli@bigbro:~/Downloads$ ls
6084.mol  C2H6.xyz  CO2_mp-20066_conventional_standard.cif  CO2_mp-20066.xyz
qli@bigbro:~/Downloads$ cat CO2_mp-20066.xyz 
12
Lattice="5.802699 0.0 0.0 0.0 5.802699 0.0 0.0 0.0 5.802699" Properties=species:S:1:pos:R:3 spacegroup="P 1" unit_cell=conventional pbc="T T T"
C        0.00000000       0.00000000       0.00000000
C        2.90134950       0.00000000       2.90134950
C        2.90134950       2.90134950       0.00000000
C        0.00000000       2.90134950       2.90134950
O        0.67853280       0.67853280       0.67853280
O        2.22281670       5.12416620       3.57988230
O        3.57988230       2.22281670       5.12416620
O        5.12416620       3.57988230       2.22281670
O        5.12416620       5.12416620       5.12416620
O        3.57988230       0.67853280       2.22281670
O        2.22281670       3.57988230       0.67853280
O        0.67853280       2.22281670       3.57988230
qli@bigbro:~/Downloads$ ase gui CO2_mp-20066_conventional_standard.cif -o POSCAR 
/home/qli/.local/lib/python3.8/site-packages/ase/gui/ag.py:85: UserWarning: You should be using "ase convert ..." instead!
  warnings.warn('You should be using "ase convert ..." instead!')
qli@bigbro:~/Downloads$ cat POSCAR 
 C  O 
 1.0000000000000000
     5.8026989999999996    0.0000000000000000    0.0000000000000000
     0.0000000000000000    5.8026989999999996    0.0000000000000000
     0.0000000000000000    0.0000000000000000    5.8026989999999996
   4   8
Cartesian
  0.0000000000000000  0.0000000000000000  0.0000000000000000
  2.9013494999999998  0.0000000000000000  2.9013494999999998
  2.9013494999999998  2.9013494999999998  0.0000000000000000
  0.0000000000000000  2.9013494999999998  2.9013494999999998
  0.6785328048660000  0.6785328048660000  0.6785328048660000
  2.2228166951340000  5.1241661951339994  3.5798823048659996
  3.5798823048659996  2.2228166951340000  5.1241661951339994
  5.1241661951339994  3.5798823048659996  2.2228166951340000
  5.1241661951339994  5.1241661951339994  5.1241661951339994
  3.5798823048659996  0.6785328048660000  2.2228166951340000
  2.2228166951340000  3.5798823048659996  0.6785328048660000
  0.6785328048660000  2.2228166951340000  3.5798823048659996
qli@bigbro:~/Downloads$ ase gui POSCAR 
qli@bigbro:~/Downloads$ 

OK，大功告成。可以拖到服务器提交任务了。值得注意的是，ASE输出的POSCAR总是把元素行放在第一行的位置，类似于VASP4的POSCAR格式，如果不爽的话，

1） 自行修改为VASP5/6的POSCAR格式。

2） 也可以修改ASE的源代码，

2.1） 找到ASE的安装目录，打开 ase/io 目录下的vasp5.py 文件；

2.2） 将vasp.py中vasp5 = False 全部改为vasp5 = True；

2.3） 保存vasp.py 即可。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

qli@bigbro:~/anaconda3/lib/python3.8/site-packages/ase/io$ pwd
/home/qli/anaconda3/lib/python3.8/site-packages/ase/io
qli@bigbro:~/anaconda3/lib/python3.8/site-packages/ase/io$ grep vasp5 vasp.py
    vasp5 = True
        vasp5 = True
    if not vasp5:
               symbol_count=None, long_format=True, vasp5=True,
    if vasp5:
qli@bigbro:~/anaconda3/lib/python3.8/site-packages/ase/io$ cd ~/Downloads/
qli@bigbro:~/Downloads$ ls
6084.mol  C2H6.xyz  CO2_mp-20066_conventional_standard.cif  CO2_mp-20066.xyz  POSCAR
qli@bigbro:~/Downloads$ ase gui CO2_mp-20066_conventional_standard.cif  -o POSCAR
/home/qli/.local/lib/python3.8/site-packages/ase/gui/ag.py:85: UserWarning: You should be using "ase convert ..." instead!
  warnings.warn('You should be using "ase convert ..." instead!')
qli@bigbro:~/Downloads$ head -n 6 POSCAR
 C  O
 1.0000000000000000
     5.8026989999999996    0.0000000000000000    0.0000000000000000
     0.0000000000000000    5.8026989999999996    0.0000000000000000
     0.0000000000000000    0.0000000000000000    5.8026989999999996
   C   O

如果有其他的需求，可以继续研究ASE或者等后面的骚操作。再继续分解自己：怎么才能把其他格式的文件转化成我需要的格式呢？

上文中的例子已经给你足够的启发了。

妹妹你大胆地往前走
往前走 莫回呀头
通天的大路 九千九百
九千九百九呀
……

本节主要给大家介绍一下中间体计算的具体流程。

计算中间体前的准备：

1）催化剂表面模型的结构以及性质自己多多了解；表面原子排列，谁带正电荷，谁带负电荷，原子间距等等，这些基本的东西自己心里要有数；没事可以整理到表格或者ppt里面，方便以后写文章用；

2） 表面上可能发生的反应路径，各个基元反应，自己列举出来，整理归类方便以后分析。

3）反应路径中各种可能分子的吸附结构先优化完成，通过这些基本分子的吸附能，吸附结构，对体系大体了解一下，一方面参考文献中别人提到的反应路径，另一方面添加自己认为可能的情况。

关键点：

中间体的优化和前面我们讲到的分子或者原子在表面上的优化是一样的。意思就是计算分子吸附的时候，基本的操作同样适用于中间体的优化。可以分布算，先从粗糙的精度出发，然后慢慢提高精度。注意的一点：

对于计算速度的影响：

1） KPOINTS：先把slab固定住，用（1x1x1）的K点对刚刚搭建的中间体结构预优化一下（基本上跑个NSW = 60就差不多了）得到一个理想的初始结构，然后再增加K点密度，放开slab的表层原子，继续优化；

2）EDIFF，EDIFFG：这两个参数主要影响电子步和离子步的步数，当我们用低密度的K点进行计算的时候，电子步很快收敛，离子步收不收敛我们也不在乎，所以在这一步的时候，不用对这两个参数太较真；关键在于增加K点后，我们选取一个什么样的标准，既得到理想的结构和能量，也避免了服务器做无意义的计算；

3）ISPIN：体系需要考虑磁性的话，粗优化的时候，可以先关掉。后面提高精度后，记得加上就可以了。

4）NCORE等并行参数也不要忘记加上，详细参考：

熟练掌握不同精度下这几个参数的使用，可以极大地提高你的计算效率，节省机时。前面说的这些，基本上试用大部分的体系，但还是要靠大家自己多多测试，摸索一下最适合自己体系的方法。

本节主要取表面上：CH3OH –> CH3O –> CH2O 这两步来简单介绍，下图中圈出来的两步；其他的路径照着葫芦画瓢就行了。

第一步：CH3OH –> CH3O + H

1）反应产物为表面上的甲氧基（CH3O*）和氢原子（H*)， CH3O* 可以继续进行C–H键断裂，或者C–O键断裂，生成CH3*+O* 或者CH2O* + H*。

注意：这里你应该可以想象到反应网络的复杂性，从一个反应物出发可以有多个不同类型的基元反应，产物作为下一步的反应物，也会有不同类型的反应，导致整个反应网络变得越来越复杂，最后交叉合并，汇总到最基本的分解产物上来。这里我们用的是甲醇，当分子C链增长的时候，基元反应的数目就会急剧增加，增加到我们挨个算不能完成的情况，这也是我们为啥要研究BEP关系的原因。另一个原因就是，成千上万个反应，整天瞎算过渡态，有时候你并不会学到太多新鲜事物，只是重复搬砖的工作，对我们以后的成长不利。

2）H原子的吸附，我们就跳过吧，相信大家经过前面的学习，应该都会了。计算CH3O的结构，我们可以从CH3OH的POSCAR出发，进行下面几点思考：

A) 删除POSCAR中对应H原子的坐标，保存为新的POSCAR，即表面上只有CH3O，提交任务进行优化；

B) CH3OH的结构中，O是在一个top的位置上，但O–H 键断裂之后，O的成键能力增强了，还会继续呆在top位置上么？所以，我们要把OCH3挪到bridge, hcp, fcc 位点上，得到几个不同的结构，提交任务进行优化；

C) 前面我们的结构OCH3是斜着在表面上；会不会有直立的结构？ OCH3是通过O与表面结合，我们是不是可以参考O原子在表面的吸附结构；搭建OCH3直立的吸附结构？

D）前面说的这几点，我们都可以通过前面所说的关键点，用一个粗糙的模型扫描一遍。然后大致分析下结果，然后进行下一步的计算，最终得到稳定的OCH3在表面上的稳定结构。

第二步：CH3O –> CH2O + H

这一步，从甲氧基进行了一步C–H反应，得到了甲醛和H。H的话继续跳过；甲醛的话，这里我们应该也需要跳过；因为甲醛是一个分子，在计算中间体前的准备这一步中，我们应该已经完成。安全起见，我们也列举一下几点需要注意的结构搭建工作：

1） 甲醛分子以C，以O还是C=O双键与表面结合？

2） 以不同部分与表面结合的时候，会在哪个位点上？

3）垂直？斜着，还是平行吸附？

总之，尽可能多思考不同的表面结合方式，采用粗糙的模型快速进行筛选，也要实时思考，粗糙模型会不会对这样的快速筛选产生影响？会不会漏掉什么可能的结构？

后面的过程依次类推，直至甲醇分子分解为最基本的化学单元(H，C, O)。

下面介绍另外一种中间体模型，首先回顾一下上面搭建CH$_3$OH–> CH$_3$O + H这一个反应中间体的过程：（这里我们先定义这个中间体模型为：A）

1. 从CH$_3$OH的POSCAR出发，

i) 删除POSCAR中与O相连的H原子的坐标，此时表面上只有CH$_3$O的结构，且CH$_3$O位于top位上；

ii) 稍微修改一下OCH$_3$的结合位点，尝试一下附近的hcp，fcc，或者bridge位点上CH$_3$O的结构，并分别保存到不同文件夹里面；

iii) 提交任务计算，并对比不同OCH3结构的能量，取最低的那个。

iV) 第ii）步中我们暴力直接删掉H原子，基于的假设是一旦O–H键断裂后，产物之间（CH$_3$O*, H*）互不影响。

另外一种中间体模型（设为：B）长什么样子呢？

1） OCH$_3$的结构与前面A的模型一样，有top，fcc，hcp，bri各种尝试的结构；

2） 除此之外，CH$_3$OH 中与O相连的H也在OCH$_3$ 附近。

也就是说，在A 模型中，我们直接删掉H原子，而在B模型中，我们只是断开O–H 键，H还留在表面上。也即是B模型保持了与反应物分子相一致的原子种类以及数目。

这两个办法有什么区别呢？哪个更好些?

我们先分析下B)的模型：它比A模型多了一个H，可以认为H 吸附在A模型上；因此A中CH$_3$O的存在会对H的吸附能产生影响，一方面来源CH$_3$O对催化剂电子结构的影响，另一方面来自与空间的影响，类似于前面我们学习过的覆盖度对O原子吸附的影响。因此我们可以计算并对比H在纯净slab和A上面的吸附能。

$\Delta E_{ads}$ =$E_{ads}^{slab}$ -$E_{ads}^{A}$ = [$E_{H*}$ - $E_{Slab}$ - $E_{H_2}/2$] - [$E_{CH_3O*+H*}$ - $E_{CH_3O*}$ - $E_{H_2}/2$]

= ($E_{H*}$ + $E_{CH_3O*}$) - ($E_{Slab}$ + $E_{CH_3O*+H*}$)

A和B模型的区别为： 将$CH_{3}O$ 和 H 分离到无穷远时候的体系能量的变化。

再看一下A和B两个模型中反应热的计算：

A)模型：$\Delta$E(A) = $E_{CH_3O*}$ + $E_{H*}$ - $E_{CH_3OH*}$ - $E_{slab}$

B)模型：$\Delta$E(B) = $E_{CH_3O* + H*}$ - $E_{CH_3OH*}$

反应热的差：

$\Delta$E(A) - $\Delta$E(B) = $E_{CH_3O*}$ + $E_{H*}$ - ($E_{Slab}$ + $E_{CH_3O*+H*}$) = $\Delta E_{ads}$

通过上面的分析，我们得到同样的结果，即两个体系模型的区别主要在于产物之间的相互作用。如果产物之间相互作用弱，可以认为彼此相差甚远，互不影响，也就是A模型，这也是大部分人在计算过程中常常采用的方法。例子有很多，就不一一介绍了，随便取一个： ACS Catal. 2014, 4, 11, 4178–4188。有兴趣的可以去看看。

对于B模型，可以从表面覆盖度对反应的影响这一角度来分析，在键断裂后，如果产物与催化剂结合较强，来不及扩散至互不影响，那么局部的覆盖度就会增加，从而对后面的反应产生影响。有兴趣的可以参考一下这两篇文章：

ACS Catal. 2015, 5, 1, 104–112； ACS Catal. 2014, 4, 6, 1991–2005。

覆盖度对表面反应的影响很复杂，一般人也不愿意碰，但是在计算过渡态相关的操作过程中，B模型却发挥着重要的作用。可以从下面2个角度来分析：1）过渡态插点的技术角度 以及2）获取合理过渡态的角度。这些将会在下一节分析。

上一节我们简单介绍了一下搬砖的原则性问题。这一节，我们主要分析下搬什么砖的问题。继续采用甲醇在金属表面上分解的这个例子。物化书中我们学过表面反应的基本步骤： 1）反应物分子的吸附；2）表面反应； 3） 产物脱附。

甲醇分解反应第一步是先从气相吸附到催化剂金属表面。 对于这一步， 我们自然需要从甲醇的气相分子出发，研究甲醇分子在表面上的吸附结构。其中涉及到一些基本的计算细节：

1） 如何优化气相分子的结构，计算频率呢？

2）如何优化体相结构，切面，优化，计算表面相关性质：比如表面能、功函数等。

这些我们前面也已经讲过了，所以直接跳过。

这里还有一个问题在于，你去用什么表面去模拟催化剂？ 比如说Cu，可以有(111), (110), (100), (211)….等各种不同的指数面，你选哪个来算？为什么选这个来算？这涉及到：我们的研究对象是什么？通过研究这一对象，我们要解释什么？期望获得什么样的结果？对于催化反应来说，我们的研究对象有2个：催化剂和它所催化的反应。

如何获得催化剂的合理模型就是一个至关重要的问题了。一般来说，有下面几个供大家参考的办法：

1）选稳定的表面，也就是表面能低的那个，对于纳米颗粒来说，尺寸如果够大，那么表面能低的表面原子在所有暴露的表面原子中的比例会高，那么可以近似用这个面来模拟催化剂的情况；那么尺寸多大才算足够大呢？

我们可以尝试着从两个极端来进行思考：我们的催化剂如果是金属的话，如果采用一个单金属原子作为催化剂模型，可以想象吸附能肯定很大。当原子慢慢堆积成团簇，团簇再慢慢长到无穷大的时候，吸附能也随之发生变化。我们取无穷大时（最稳定的表面）的吸附能（E0）作为参考，看看金属团簇多大尺寸的时候，吸附能跟E0接近。那么该尺寸以及更大的纳米颗粒就可以用无穷大的催化剂模型来模拟，也就是最稳定的表面。一般来说，纳米颗粒的直径大于2nm（胆儿肥的）的时候，可以用最稳定的表面来模拟催化剂。保守些，直径大约3nm的时候可以放心地去用，不会有什么问题。

下面推荐一本书： Fundamental Concepts in Heterogeneous Catalysis 。书中第21页， 作者用O在Pt(111)表面的吸附能作为参考，测试了不同纳米颗粒尺寸时的吸附能，发现直径在2nm的时候，已经非常接近了。前面说的那个2nm，也是从这里来的。如果你想问我这本书从哪里可以下载，可以关注公众号：BigBroScience，然后回复：ex84 。

2）对于前面的一番说法， 也会有人说，催化剂的活性位不一定在平坦的表面上，边边角角，台阶上可能更容易发生反应。如果你单纯把工作聚集在最稳定的表面上，肯定会遗漏些什么？这也是很多审稿人喜欢折腾我们这些小白鼠的地方。一个金属有很多不同的指数面，单纯考虑最稳定的是一个保险折中的办法。因为不管你的体系搞的再复杂，跟真实情况相比，模型终究是简单的。我们只是尽自己可能去尝试着理解，解释实验的现象，描述下体系的性质，以及进行一些比较简单直接的预测。但是，如果你要深究这些，兴趣也恰好落在在这一块，那么可以对体系的各个相关结构进行一个全面的计算，分析不同表面结构，性质，反应活性有什么不同。在这个全面的计算体系中，

除了各种台阶面，最稳定的表面你还是要算的，这样才能有个对比。原则上说，只要你力气足够，资源充足，可以选任何你想要的表面来计算；但问题是你总会有吐的时候，如果想多算几个面的时候，要留意下面的几点：

注意-1：开放的表面在计算表面能，slab优化或者吸附的时候，层数的影响会很大，这个前期的测试工作要做好；

注意-2：生存法则：

A）灌水第一条：同一个反应， 这个指数面算一下发一篇文章，换个面再发篇文章，再换个面继续发文章….

B）灌水第二条：同一个面： 算一个反应发一篇文章，换个反应算一下再发一篇文章，再换反应；

C）灌水第三条：反应，表面和反应交替，排列组合式发文章。

3）还有一种情况，就是我们手上有实验数据，想模拟一下实验的结果。这个时候，就可以参考下催化剂的表征结果，有针对性地选取催化剂的模型，使其尽可能地去描述催化剂的结构。

注意的有以下几点：

A) 如果做实验的人找你帮忙算个东西，能推掉就推掉，别在上头浪费自己的功夫，顶死也就是个共同一作，国内基本不认可，做无用功；如果刚开始推不掉，就看看他的实验表征结果，一般般的表征以模型搭不出来为由直接拒掉；好的模型是成功的重要因素。前面我们说了，模型再复杂也是简单的。再考虑到实验表征的结果不确定性，你有很大的几率是跟实验上对不上的，到头还是白忙活。很多人会拿工作要冲子刊，杰克斯，俺狗娃为由诱惑你上船，要掂掂自己的分量，别头脑发热用屁股做决定。

B) 对于做实验的人来说，99%会低估计算的复杂性，总以为鼠标点点就OK的事情，其中也不乏一些所谓的大牛们。话虽然难听，却也是大实话：不要瞎几把评论自己不懂的领域。真想找人算一下，要找个靠谱的合作伙伴，找到不靠谱的瞎算，投文章更会给整个工作拉分。

4）最后一个及其重要的就是，认真进行文献调研。看看前人都是怎么选取的，为什么这样选？根据什么依据来选？最终确定自己的计算模型。

关于模型的一些就先到这儿， 本篇的工作主要是重复一遍前人已经算烂的反应，然后在此基础上，分析BEP和TSS线性关系。就直接拿各个金属最稳定的表面来计算了。那么， 什么是BEP线性关系？ 这里的BEP的全称是：Bell–Evans–Polanyi principle. 指的是对于同一类型的N个反应，这些反应的能垒和反应热之间存在的一个线性关系: $E_a$= $\alpha \Delta E$ + $\beta$ 。更详细点的介绍，参考维基百科搜索：Bell–Evans–Polanyi principle。起初这个关系是基于均相反应来说的。在2000年的时候，Neurock将这个关系引入多相催化体系，发现表面反应也存在这一关系。随后，经过了Norskov等人的继续研究，相关的文章就越来越多了。有兴趣的可以通过下面的几篇经典文献开始接触，然后再逐渐深入：J. Catal. 191, 301 (2000)， J. Catal. 197, 229 (2001)， J. Chem. Phys. 114, 8244 (2001)

BEP关系有啥用呢？主要有两点:

关键的一点就是用来节约机时：过渡态不是一个稳定的结构，它的优化需要的计算时间比直接优化反应物种的要多得多。插8个点，也就是8倍的单个优化。再加上计算还经常容易失败。因此我们可以通过BEP关系，跳过TS的优化，直接预测反应的能垒，从而达到节约机时和生命的作用。这个方法简单粗暴，与直接进行反应路径扫描相比，在时间成本上具有压倒性的优势，但精度和准确度上就会稍逊一筹。

另外一点是来预测反应的能垒，这个其实和上面的有些重复，但角度不同。比如一个复杂的催化体系，基元反应数目成千上万，我们挨个算过渡态显然不可行。但是中间体的数目明显少的多，那么我们通过中间体获取反应热之后，通过BEP关系得到反应能垒，也就有了动力学的一些基本数据，从此出发，可以做一些机器学习以及微观反应动力学的模拟工作。从而跨越理论计算的原子分子尺度来到达一个更为宏观的尺寸和时间的尺度，也就是所谓的多尺度模拟。

继续前面的分析，如果中间体的数目也很多，算都算不完，那该怎么办？这个也可以通过一些线性组合，数据分析的方法来进行预测。比如通过Group Additivity的方法，可以将对分子中不同group的能量的求和来获得。具体的内容在后面其他章节再慢慢进行介绍。由于篇幅关系，本节（ex84的上半部分）主要是让大家对自己的研究对象有一个明确的认识，为什么研究它？为什么不研究别的？下半部分介绍中间体计算的一些注意事项。

在线弹性理论中，体系被施加一个无穷小的应变，从而体系应变后的能量可以以小应变为自变量泰勒展开，并忽略二阶导以上的高阶项。通过对不同特定的独立的应变模式求解能量应变方程或者应力应变方程，我们可以获得材料所有的二阶弹性常数（SOECs）。SOECs反映了体系的简谐弹性特征。通过SOECs，一方面我们可以获得体系的基于Vogit-Reuss-Hill平均的多晶力学性能，包括杨氏模量，剪切模量，体积模量，泊松比和维氏硬度等，以及在此基础之上的各向异性系数等等特征；另一方面，我们也可以通过SOECs获得单晶力学性能在空间中的具体分布，并将其用图的模式直观的表示出来。

图中我们展示了某种金属的杨氏模量和剪切模量（每个面的平均值）的空间分布图，感兴趣的读者可以猜猜这是哪种布拉维晶型的金属。

在实际负载下，材料承受的往往是有限形变，这时候材料的非简谐弹性将称为不可忽视的因素。很多实验数据有去测量体系的三阶弹性常数（TOECs），甚至四阶弹性常数（FOECs）。注意到TOECs和FOECs相对比SOECs在实验中更加受到实验精度的影响，具体的测量非常困难。如果我们能在理论上对其进行预测，将变得非常有意义。实际上已经有若干文献进行了相关研究，笔者也抱着好奇的态度尝试进行了TOECs和FOECs的计算。原理上来说相对简单，我们只需要将能量对应变进行泰勒展开，将阶次保留到四次，即可同时获得SOECs，TOECs和FOECs。

将所有的应变η设为同一个数值，我们就可以用一个包含四次，三次和二次项的多项式来拟合计算的数据，获得包含各个弹性常数的多元一次方程组。求解这个这个方程组即可获得所有的弹性常数。

不过实际操作起来并没有这么简单。即使是最简单的立方体系，其TOECs也达到了6个，而FOECs更是达到了11个。对于更低对称性的体系，其计算量更是一个天文数字（开玩笑的，天文倒不至于，很大就是了）。一种避让的方法就是，我们可以用应力对应变进行展开，这样相对需要的独立应变数少一些。笔者使用Au作为测试体系，根据文献对体系施加了11种独立应变，如下图。

感兴趣的童鞋可以直接阅读文末给的参考文献。笔者通过测试发现，首先，体系的应变量必须足够大。笔者一开始犯上了计算SOECs的习惯，使用的很小的应变，计算结果简直南辕北辙。后来将应变加到15%后，获得的TOECs相对理想。细想也是，应变很小时，体系还处于线弹性阶段，这时候的非谐弹性非常微弱，极容易被数值误差干扰。将应变加大到非谐区后，数值就非常稳定了。同时，尽管我们使用了大应变条件，但是实际上高阶弹性常数的具体数值依旧会对数值精度的波动非常敏感。

同文献一样，笔者使用了非常高精度的k点半径和截断能。随着精度的提高，结果趋于较好的收敛。笔者在此展示了笔者自己算的TOECs与实验的比较，两者相对来说还是比较接近的（相对实验的误差与文献差不多）。总的来说，如果有小伙伴需要计算材料的TOECs或者FOECs，笔者建议一定要对应变，k点，截断能三者都要做仔细地收敛测试。

Wang,Hao, and Mo Li. “Ab initio calculations of second-, third-, andfourth-order elastic constants for single crystals.” PhysicalReview B 79.22 (2009): 224102.

本文经ponychen授权发布，版权属于ponychen。

在上文中，我们介绍了用于解释Hall Petch关系的位错阻塞模型和晶界位错源模型。两种理论都是较为早期出现的模型，原理上相对简单，因此文章中仔细介绍了相对应的推导过程。本文将继续介绍剩下的相关模型，由于推理过程相对复杂，读者如果对具体推导过程有兴趣可以阅读对应参考文献。

上文介绍的两种模型都未在公式中直接体现材料变形中塑性应变对Hall Petch系数的影响。因此，位错阻塞模型和晶界位错源模型仅仅只能描述多晶材料在屈服极限之下的情况。我们知道，当多晶材料所受应变超过其屈服极限，材料由于位错的交截和阻碍出现了加工硬化。加工硬化势必会影响Hall Petch系数的数值。在Meakin和Petch（文献1）提出的work haraening（加工硬化）模型中，Hall Petch系数被看成是一个与应变有关的因变量。

由于Hall Petch系数中包含了应变，该模型同时成功预测了材料中的加工硬化行为。如果我们将晶粒尺寸d固定，可以发现该模型描述了一种应力与应变成抛物线形状的加工硬化关系，这与很多实验是相符合的。

目前为止介绍的所有模型都默认认为位错来自于位错源的开动。然而实际上，当晶体被弯曲，或者在晶粒之间的交叠以及孔洞处，由于非均匀变形，晶体为了降低变形带来的巨大应力，会自发形成一些几何必须位错来释放应力。下图展示:

当晶体被弯曲时，晶体中自发在易滑移面上出现的几何必须位错林。在Ashby（文献2）等提出的几何必须位错模型中，作者认为晶界作为晶体中的应力协调区，其上的位错类型主要为几何必须位错主导，并且晶体中绝大部分的几何必须位错都位于晶界上，因此材料的强度由晶界上的几何必须位错密度决定。而平衡位错密度又由晶粒尺寸决定，通过结合加工硬化模型，作者获得了如下的Hall Petch关系式。

位错阻塞模型将材料的强度归功于晶粒内部位错源开动的位错受阻，而目前介绍的其他模型则都将材料的强度归功于晶界处的位错密度受阻。在这两种极端情况的基础上，Meyersm（文献3）等人认为（复合模型），一个多晶材料应该被分为两部分：连续的拥有较高流变强度的晶界部分以及被晶界包围的离散的拥有较低流变强度的块体部分（见下图4，t为晶界厚度）。

当材料承受外力时，晶界部分首先发生变形，形成大量几何必须位错；随着进一步变形，晶界首先发生局部屈服形成加工硬化层，并形核大量位错向晶粒内运动；最后，随着晶粒内位错密度的上升，晶粒部分发生屈服。至此，材料整体完成了整个塑性变形过程。基于这一过程假设，作者得出了如下关系。

可以发现，该关系式并非符合正常的Hall Petch关系，式子中同时出现了一次方和二次方项。作者进一步假设随着晶粒尺寸的增大，其容纳几何必须位错的数量也在上升，因此其厚度t也在上升。基于这一假设，作者得到了下图的强度与晶粒尺寸的关系。

可以发现，复合模型在晶粒尺寸较大的时候可以较好的匹配Hall Petch关系，而当晶粒尺寸小到纳米级的时候，复合模型成功预测了超细晶材料的强度反常降低现象的存在。这是之前所有的理论模型都无法给出的结果。

实际上，在复合模型之后材料学家依旧提出了大量的理论模型，其目的也不仅仅是在于解释Hall Petch关系。由于篇幅所限，本文仅仅介绍了一些具有代表性的模型。实际上，从实验的角度来说，很多材料也并非严格遵循Hall Petch关系中-1/2指数，该指数存在一个很大的范围。对于我们来说，-1/2只是代表材料强度与晶粒尺寸联系的一个优美符号。

1. Meakin, J., and N. J. Petch. “Report ASD-TDR-63-324.” Symposium on the Role of Substructure in the Mechanical Behavior of Metals,(ASD-7DR-63-324, Orlando, 1963) pp. 1963.
2. Ashby, M. F. “The deformation of plastically non-homogeneous materials.” The Philosophical Magazine: A Journal of Theoretical Experimental and Applied Physics21.170 (1970): 399-424.
3. Meyersm, Marc A., and E. Ashworth. “A model for the effect of grain size on the yield stress of metals.” Philosophical Magazine A 46.5 (1982): 737-759.
4. Kato, Masaharu. “Hall–Petch relationship and dislocation model for deformation of ultrafine-grained and nanocrystalline metals.” Materials Transactions 55.1 (2014): 19-24.

本文经ponychen授权发布，版权属于ponychen。

Hall Petch关系是材料学中描述多晶材料最为重要的一个相关关系，它反映了材料强度与晶粒的尺寸的开方成反比。同时材料在极小晶粒尺寸下表现出的反常Hall Petch 关系也一度成为材料学家的研究热点。关于反常Hall Petch关系的研究发现屡现于高水平期刊，本文将不做更多介绍。从Hall Petch关系的第一次报导开始，已经有大量的文章针对各种材料体系（纯金属，金属间化合物，多相合金等等）进行了相关研究。然而，从Hall Petch 关系提出至今，依旧没有出现一个完美简洁的理论可以解释这个现象。研究人员为了解释这个现象已经提出了大量的理论，在接下来的篇幅中，作者将从中提取了几个典型的理论做简要介绍

1951年，Hall（文献1）发表的论文中首次报导了钢中存在的Hall Petch关系:

$\sigma=\sigma_0 + k_{HP}d^{-1/2} $

其中$\sigma$表示流变应力，d表示晶粒尺寸，k是Hall Petch 系数，其对不同材料，温度和应变条件都非常敏感。

与此同时，他也在文章中首次提出了试图解释Hall Petch关系的pile up（位错阻塞）模型。

该理论认为材料的强度来自于外力作用下导致的晶粒内部位错源开动，位错源启动后形成的位错林被晶粒之间的界面（晶界）阻碍后引起了材料的应力强化。简单起见，我们不妨假设两个近似六边形的晶粒以晶界AB接触，左边晶粒的中心O点存在一个位错源，在外加应力的作用下位错源沿着某一滑移面释放出一系列位错。随着第一个发射的位错冲撞到晶界AB上，整个位错列队被禁锢在位错源与晶界之间。由于晶界墙的阻力和新产生位错的弹性应变，整个位错列最后会形成一个稳定分布，假设晶粒的平均尺寸为d，则位错列中位错的数目为：

$n=\frac{(1-v)d\tau_e}{2\mu{b}}$

其中$\tau_e$表示位错的有效剪切应力，它可以表示为：

$\tau_e=\tau-\tau_i$

显然，它实际上是在加载应力中扣除了弹性应力部分。试想如果你要产生一根位错，首先晶体要发生塑性变形，因此塑性变形前的弹性部分并不是我们应该要考虑的部分。既然我们已经知道位错列中位错的数目，如果我们粗略地认为，所有的位错作用到晶界上的应力都是等效的，那么晶界上受到的应力$\tau_p$为

$\tau_p = n\tau_e$

Hall认为，材料屈服发生于位错列冲破晶界的时候，也就是说 $\tau_p$大于临界值$\tau_c$。如果我们取等号，再将前面的n代入，我们就可以得到：

$\tau = \tau_i +\surd{\frac{2\tau_c\mu{b}}{1-v}} \times d^{-1/2}$

该公式看起来非常合理简洁，-1/2次方的关系也出来了。但是如果材料的屈服依靠的是位错像裂纹一般冲破晶界，我们应该在实验中大量发现，但实际上，只有少数几篇文献报导过相关现象。更重要的是，依据该理论计算的屈服强度跟实际数据还是有较大出入。

在pile up模型的基础上，Cotrell（文献2）认为，材料的屈服并非来自于位错林冲破晶界，而是由于位错林在晶界的终端的应力场诱发了邻近晶粒中的位错源(如下图的R处)的开动。这样子需要的 临界应力要小于pile-up模型。

晶粒内部出现位错源固然可行，但是另一方面，晶界本身就是天生的位错源，在晶粒内部没有经受剧烈变形的情况下，晶界上的位错形核位置要远远多于晶粒内部。Li （文献3）于1963年提出了boundary source（晶界位错源）模型。为了方便起见，我们不妨假设晶粒是边长为d的立方体。假设晶粒的每一个面上的台阶密度为m，假设每一个台阶都有一定几率发射位错。进一步简化，我们不妨认为每个台阶都会发射位错，则一个面的总位错数目为$md^2$。每一个晶粒与毗邻晶粒公用一个晶面，则一个晶粒可以认为有三个晶界面属于这个晶粒。所以:

我们可以认为一个晶粒有3$md^2$个位错。这些位错在外加应力作用下继续向晶内扩散，形成位错林。位错林的密度为:
$\rho = 3md^2 / d^3 = 3m/d$
另一方面，晶体材料的位错密度和强度之间可以由著名的Taylor公式来描述:
$\sigma = \sigma_0 + M\alpha\mu{b}\surd\rho$
代入前者即得:

$\sigma = \sigma_0 + M\alpha\mu{b}3^{1/2}m^{1/2}d^{-1/2}$

我们再一次获得了Hall Petch关系。在原文中也有说明这个方式计算的理论屈服强度与pile up模型基本一致。目前介绍的这两个模型，前者认为位错完全起源于晶粒内部，后者则认为位错完全起源于晶界，两者分别是实际情况下的两个极端。但是殊途同归，两个极端最后得到了相同的结论和近似的屈服强度。不过，这两个模型的局限性也是明显的，实验中Hall Petch曲线的斜率对晶粒结构和化学成分非常敏感，但是上述模型中并未考虑这些实际因素。

值得注意的是，在2002年的一篇文章4中，该文作者成功在MD模拟中发现了晶界处位错的形核和发射过程。

（未完待续）

1 Hall, E. O. The deformation and ageing of mild steel: III discussion of results. Proceedings of the Physical Society. Section B 64.9 (1951): 747.
2 A.H. Cottrell,: The Mechanical Properties of Matter, Wiley, New York (1964)
3 Li, James CM. Petch relation and grain boundary sources. Transactions of the Metallurgical Society of AIME 227.1 (1963): 239.
4 Van Swygenhoven, H., P. M. Derlet, and A. Hasnaoui. Atomic mechanism for dislocation emission from nanosized grain boundaries. Physical Review B 66.2 (2002): 024101.
本文经ponychen授权发布，版权属于ponychen。

*基于DFT+U研究应力调控CeO2(111)催化分解水制氢的活性*

研究背景

氢能清洁无污染，是最具潜力的能量载体。使用太阳能、风能和水能转换得到的清洁电能电解水制氢是人来未来最具发展潜力的、可持续的、绿色环保的制氢技术。传统低温碱性电解水制氢技术的效率受限于液体电解质较低的传质过程和氧电极较慢的反应动力。相比于液体传质技术，新兴发展的高温固体燃料电池 (SOEC)技术以固体作为传质媒介，充分利用工业废热提供反应所需的热能，具有极高的传质效率和水电解效率。然而人们对SOEC制氢反应机理尚缺乏深刻认识。计算模拟技术是揭示能源材料领域中多相催化反应过程及机理的重要手段，可以辅助了解SOEC电解水的关键反应步骤，以开发、低廉且稳定的催化剂用于增强氢电极的稳定性以及提高制氢反应在燃料电池中的效率。兼具优异的电子和离子传质能力的氧化铈（CeO2）是SOEC水分解重要的催化剂和稳定剂。CeO2优越的催化活性归因于Ce3+与Ce4+的氧化还原循环。基于密度泛函理论，该工作使用Hubbard-U（Ueff=4.5 eV）修正的方法来描述Ce 4f轨道电子的强关联作用，研究表面应力对CeO2(111)催化制氢的影响。

研究内容

水在CeO2表面分解制氢包括三个基本反应步骤：（1）CeO2(111)表面生成一个氧空位；（2）H2O吸附于氧空位附近最终形成表面OH-；（3）表面OH-分解形成H2并留下一个晶格氧，如Figure 1所示水分解过程中形成的反应中间体以及最终生成H2的反应过程。该工作主要研究了三种可能的反应路径，即不同OH- 覆盖率下形成H2的过程如Figure 1所示。通过研究应力对三种不同反应路径中反应中间体形成的影响，以及应力对各反应步骤反应能垒的调控，来探究应力如何增强CeO2催化制氢的性能。

首先，研究发现拉伸应变有效增强了氧空位和OH-等中间体的稳定性，随之增加了各反应步骤的反应能垒。与之相反的是压缩应变可以通过减弱OH-与表面的结合能，从而降低生成H2的能垒。通过计算不同反应温度下中间体的自由能可以获得沿着2H，4H和5H的反应路径能量图。研究发现不管处于何种应力条件下生成氧空位和氧空位的扩散，以及OH-的形成不需要很高的能量势垒，而OH-分解形成H2是水分解制氢反应路径中能垒最高的一个步骤。然后，通过对比不同反应温度下，研究发现在温度低于1000 K时，沿着5H路径可以最快速地水分解生成H2 如图2所示，且最佳反应路径不受应力的影响。然而，当温度升高到1200 K时，应力有效地调控了制氢的反应路径，即拉伸应力下最优路径沿着5H，而在压缩应力下最优反应路径沿着2H或者4H。此外，研究还发现反应速度随着压缩应力的增大而增快，且当压缩应力大于3%时，反应速度获得显著提升。因此，计算模拟理论研究显示压缩应力（>3%）有利于增强CeO2(111)催化水分解制氢的反应活性。

图1. 2H，4H和5H覆盖的CeO2(111)表面生成H2的反应路径示意图。黄色, 红色，蓝色，白色，黑色和灰色原子球分别代表Ce，表面O，次表面O，H，次表面氧空位和表面氧空位。

图2不同反应温度与应力条件下最优反应路径生成H2的反应速度对比图。

作者简介

本文第一作者，武甜甜为丹麦科技大学能源转化与储存系博士后研究员。

Google Scholar：

https://scholar.google.com/citations?user=rm84jfsAAAAJ&hl=en

指导老师为：Tejs Vegge教授。

课题组主页：https://www.dtu.dk/english/service/phonebook/person?id=5334&tab=1。

本文经作者授权发布，版权属于第一作者。

理论模拟导向的催化剂理性高效设计策略用于醇的选择性胺化

第一作者：王涛

通讯作者：Philippe Sautet

研究背景

胺是工业上生产药物、农用化学品、生物活性剂、聚合物和染料的重要中间体化合物，其中伯胺是衍生化反应的关键试剂。传统合成胺的方法包括氨与烷基卤的直接反应，或通过还原硝基络合物和腈，但是这些工艺的缺陷在于产生有害的副产物、废弃物并消耗大量氢气。而利用氨气通过借氢反应机理直接胺化脂肪类醇用于合成胺被认为是一种有潜力的绿色合成路线并且已经有数千吨规模的工业化应用。然而，合成胺工业仍然面临诸多挑战，尤其是选择性合成伯胺。在此背景下，探索催化剂的本征性质与合成胺活性和选择性的依赖关系变得尤为重要。近年来，理论模拟的飞速发展使得揭示催化过程中的构效关系成为可能。本文结合第一性原理微观动力学模拟和实验，揭示了决定醇的胺化活性和选择性的关键因素并构建了其与碳和氧原子吸附能的依赖关系，实现了该反应金属合金催化剂的理性高效筛选，为实验上醇胺化催化剂的设计和升级提供了重要理论参考。

内容和讨论

A. 醇的胺化反应机理: 首先甲醇的胺化被作为模型反应进行研究，该反应机理主要包括三部分(如图1)：醇的脱氢生成醛；醛与NH3迅速反应生成亚胺；亚胺的加氢生成胺。基于上述反应机理，DFT计算在九种过渡金属的密堆积表面系统展开。由于实际反应气氛中具有大量的NH3，因此催化剂的表面会覆盖一定浓度的NH3，早期的研究发现表面覆盖1/9 ML的NH3有助于准确描述胺化反应催化剂的活性顺序。其中，第二步醛与NH3反应生成亚胺被认为是气相里的快速平衡反应而且是理论模拟上的难点，该工作通过详细的反应器动力学模拟证明采用Eley–Rideal 反应机理能够准确描述该关键步骤。最终为微观动力学模拟提供了可靠的反应机理。

图1: 甲醇与氨气的胺化反应机理

B. 线性降维: 基于Norskov课题组提出的d带中心理论和scaling关系，反应过程中表面物种的吸附能并非相互独立的而是受制于scaling关系，并且基元反应的能垒也可以通过BEP关系与反应热关联，因此通过详细的描述符的筛选以及线性拟合可以实现降维和最优描述符的确认。由于该反应涉及含C、N、O的三类中间体，scaling关系可以最终将该反应中涉及的中间体的能量与C、N和O原子的吸附能进行关联。然而，包含三种不同描述符导致无法利用典型的二维火山曲线进行活性描述。深入研究发现：N原子的吸附能同样依赖于C和O原子的吸附能并且具有完美的线性关系（如图2）。而该依赖关系的物理化学根源一方面为N(3.0)原子的电负性恰好介于C(2.5)和O(3.5)原子之间，另一方面N原子的最高占据轨道能量(HOMO)可以线性分解为C和O原子的相应能量。基于上述分析，甲醇胺化反应的能量维度被降低为二(C和O原子的吸附能)。

**

图2: N原子的吸附能与C和O原子吸附能的线性关系

C. 醇的胺化活性和选择性与简单描述符的依赖关系: 基于图1的反应机理和图2的线性降维，通过构建微观动力学模型，最终建立了醇胺化的反应速率(TOF)与C和O原子吸附能的依赖关系。如图3a所示，由于对C和O原子吸附能力的区别，九种金属分布于火山图的不同位置而呈现出不同的反应活性。然而，针对醇的胺化反应，人们更加关心且关键的问题是选择性。因此，该工作进一步系统考察了导致伯胺选择性降低的不同可能性，最终发现CH3NH和CH2NH物种的耦合是导致选择性丢失的关键因素之一，并巧妙的将其与简单描述符进行了关联，最终构建了醇胺化生成伯胺的选择性与C和O原子吸附能的依赖关系如图3b。随着CH3NH/CH2NH耦合能垒(Ea)的增加，副反应变得困难，进而催化剂的选择性增加。至此，通过计算催化剂表面C和O原子的吸附能，结合图3中的火山图便可定性预测相应催化剂的醇胺化反应的活性和选择性，从而大大提高催化剂的筛选效率。

图3: 醇的胺化活性(a)和选择性(b)与C和O原子吸附能的依赖关系

D. 实验室论证预测: 基于上述活性和选择性图，通过对元素周期表中不同过渡金属进行两两排列组合，最终计算了三百余种双金属合金的C和O原子吸附能，并预测出系列有潜力的醇选择性胺化催化剂。基于上述理论预测，索尔维公司研发中心对不同合金进行了详细的实验室测试以论证图3的可靠性。具体的：本文合成了系列Co基合金催化剂，一方面由于Co金属的相对廉价，另一方面Co具有较高的活性和伯胺选择性。根据图3a的理论预测，提高Co的活性需要设法降低其O原子的吸附能，由于Pt和Pd具有弱的O吸附能力，因此将Co与Pt和Pd形成合金将可以实现反应活性提高的，然而基于图3b的信息，这将同时伴随着伯胺选择性的降低。事实上，图4和图5的实验结果也证实了CoPt和CoPd合金相比Co金属的高活性和低伯胺选择性。此外，图3表明：将Co与Ag和Ru组合将会维持甚至提高活性且提高选择性。该预测结果同样与图4和图5中的实验结果一致，并且成功设计出活性选择性都大大提高的CoAg合金催化剂。本文预测的其他有潜力的合金催化剂也将进一步推动实验上对醇胺化反应研究。

图4: 实验上不同催化剂醇的胺化生成伯胺的TOF

图5: 实验上不同催化剂醇的胺化生成伯胺的选择性

总结和展望

通过精确的第一性原理计算和微观动力学方法，准确的模拟了醇的胺化这一复杂反应的详细机理，揭示了决定金属合金催化剂醇的胺化活性和选择性的关键因素，并与简单易得的C和O原子吸附能进行关联，系统的实验测试论证了理论预测的可靠性。最终大大提高了醇胺化反应金属合金催化剂的筛选效率，为工业上醇胺化催化剂的设计和升级提供了参考。此外，该工作采用的改进的动力学模拟方法能够引入反应器环境对复杂反应机理的影响，一方面实现对实验条件更加精准的描述，另一方面能够准确定性预测金属合金催化剂的活性和选择性的变化趋势，为后续高通量催化剂筛选提供了基础。

作者简介

文章第一作者王涛博士目前为美国斯坦福大学SUNCAT中心博士后研究员，通讯作者Philippe Sautet为美国加州大学洛杉矶分校教授、ACS Catalysis副主编、法国科学院院士、美国能源局IMASC能源前沿研究中心副主任，长期致力于理论催化方向的基础研究。此外，法国里昂高等师范学院、法国里尔大学、比利时根特大学以及跨国公司索尔维研发中心E2P2实验室也为该工作提供了支持。本文经大师兄涛本人授权发布。

文章链接：https://www.nature.com/articles/s41929-019-0327-2

关注公众号：大师兄科研网，回复：胺化，可以获取下载链接。

工作简介

上海交通大学密西根学院朱虹老师课题组诚聘博士后 1~2 名，负责或协助团队进行金属合金体系的腐蚀性能的理论预测和腐蚀机理的研究，将与腐蚀实验团队以及世界各地的材料基因组计算研究团队开展紧密合作 (UC Berkeley, University of Michigan, UC San Diego, Georgia Tech, University of Maryland, etc)，并在现有的高通量计算平台上开发针对腐蚀相关基础参数的高通量计算工作流。诚挚邀请具有理论计算、电化学背景同学的加入。

导师简介

朱虹 (hong.zhu@sjtu.edu.cn) 助理教授，上海交通大学材料基因组联合研究中心成员； 材料科学与工程学院，上海交通大学密西根联合学院双聘。

Google Scholar: https://scholar.google.com/citations?user=x1BGIfEAAAAJ

课题组主页: http://umji.sjtu.edu.cn/~hzhu。

工作待遇

1. 年收入 22w+， 享受相应上海交通大学福利待遇；
2. 根据上海市博士后管理政策办理有关落户事宜；
3. 优秀者可转入专职科研序列。

招聘要求

1. 拥有材料、化学或相关专业的博士学位
2. 具有理论模拟计算经验或具有丰富的机器学习特别是材料研究领域的经验
3. 能够独立开展研究，在知名期刊上以第一作者发表论文
4. 在同等条件下，优先考虑具有较高编程能力的同学

申请方式

发邮件至 hong.zhu@sjtu.edu.cn，并附上简历以及可以到站的时间。