pymatgen的transformations模塊提供了許多用于結(jié)構(gòu)轉(zhuǎn)換的方法(類)芽腾，使用這些類可以輕易的實(shí)現(xiàn)摻雜竹伸、建超胞奇唤、添加氧化態(tài)等功能辫樱，這里介紹一下用于有序無序轉(zhuǎn)換的兩個(gè)類：EnumerateStructureTransformation 和 OrderDisorderedStructureTransformation峭拘。這兩個(gè)類可以實(shí)現(xiàn)類似supercell軟件的功能，將無序結(jié)構(gòu)轉(zhuǎn)變?yōu)橛行蚪Y(jié)構(gòu)狮暑。
下面我們以立方體結(jié)構(gòu)鈣鈦礦為例鸡挠，學(xué)習(xí)一下如何使用該模塊建模。
的結(jié)構(gòu)文件可以去 materials project 網(wǎng)站下載(mp-5827)搬男。
關(guān)于有序無序轉(zhuǎn)變
首先來了解一下什么是有序(order)與無序(disorder)結(jié)構(gòu)拣展，打開我們在mp上下載的cif文件：
# generated using pymatgen
data_CaTiO3
_symmetry_space_group_name_H-M Pm-3m
_cell_length_a 3.88947141
_cell_length_b 3.88947141
_cell_length_c 3.88947141
_cell_angle_alpha 90.00000000
_cell_angle_beta 90.00000000
_cell_angle_gamma 90.00000000
_symmetry_Int_Tables_number 221
_chemical_formula_structural CaTiO3
_chemical_formula_sum 'Ca1 Ti1 O3'
_cell_volume 58.83987623
_cell_formula_units_Z 1
loop_
_symmetry_equiv_pos_site_id
_symmetry_equiv_pos_as_xyz
1 'x, y, z'
2 '-x, -y, -z'
3 '-y, x, z'
4 'y, -x, -z'
5 '-x, -y, z'
6 'x, y, -z'
7 'y, -x, z'
8 '-y, x, -z'
9 'x, -y, -z'
10 '-x, y, z'
11 '-y, -x, -z'
12 'y, x, z'
13 '-x, y, -z'
14 'x, -y, z'
15 'y, x, -z'
16 '-y, -x, z'
17 'z, x, y'
18 '-z, -x, -y'
19 'z, -y, x'
20 '-z, y, -x'
21 'z, -x, -y'
22 '-z, x, y'
23 'z, y, -x'
24 '-z, -y, x'
25 '-z, x, -y'
26 'z, -x, y'
27 '-z, -y, -x'
28 'z, y, x'
29 '-z, -x, y'
30 'z, x, -y'
31 '-z, y, x'
32 'z, -y, -x'
33 'y, z, x'
34 '-y, -z, -x'
35 'x, z, -y'
36 '-x, -z, y'
37 '-y, z, -x'
38 'y, -z, x'
39 '-x, z, y'
40 'x, -z, -y'
41 '-y, -z, x'
42 'y, z, -x'
43 '-x, -z, -y'
44 'x, z, y'
45 'y, -z, -x'
46 '-y, z, x'
47 'x, -z, y'
48 '-x, z, -y'
loop_
_atom_site_type_symbol
_atom_site_label
_atom_site_symmetry_multiplicity
_atom_site_fract_x
_atom_site_fract_y
_atom_site_fract_z
_atom_site_occupancy
Ca Ca0 1 0.500000 0.500000 0.500000 1
Ti Ti1 1 0.000000 0.000000 0.000000 1
O O2 3 0.000000 0.000000 0.500000 1
內(nèi)容很多，前面的內(nèi)容我們只需要關(guān)注這么幾個(gè)內(nèi)容：
- 空間群
Pm-3m - 化學(xué)式
CaTiO3 - 各原子數(shù)目
Ca1 Ti1 O3
下面我們來看一下最后幾行
_atom_site_type_symbol
_atom_site_label
_atom_site_symmetry_multiplicity
_atom_site_fract_x
_atom_site_fract_y
_atom_site_fract_z
_atom_site_occupancy
Ca Ca0 1 0.500000 0.500000 0.500000 1
Ti Ti1 1 0.000000 0.000000 0.000000 1
O O2 3 0.000000 0.000000 0.500000 1
上面七行對應(yīng)的是下面七列的信息缔逛，分別是：
| 元素符號 | 原子序號 | 原子數(shù) | x軸坐標(biāo) | y軸坐標(biāo) | z軸坐標(biāo) | 占據(jù)比 |
|---|---|---|---|---|---|---|
| Ca | Ca0 | 1 | 0.5 | 0.5 | 0.5 | 1 |
| Ti | Ti1 | 1 | 0 | 0 | 0 | 1 |
| O | O2 | 3 | 0 | 0 | 0.5 | 1 |
可以看到CaTiO3中每個(gè)原子的占據(jù)比都是1备埃，且其結(jié)構(gòu)符合 空間群的對稱性，因此它是一個(gè)有序結(jié)構(gòu)褐奴，也可以說是一個(gè)完美晶體按脚。而無序結(jié)構(gòu)顧名思義就是不是完美晶體的結(jié)構(gòu)，在這里一般指的是占據(jù)比 occupancy != 1的結(jié)構(gòu)敦冬。
在研究摻雜的過程中辅搬，我們會經(jīng)常遇到這種情況，某個(gè)位置發(fā)生取代摻雜脖旱，部分原子被取代后占據(jù)比就會小于1堪遂。還有一種情況就是某種材料具有鈣鈦礦介蛉、尖晶石等結(jié)構(gòu)，但是沒有這種材料的結(jié)構(gòu)文件溶褪，比如(LLTO)币旧，鋰鑭鈦氧是一種離子電導(dǎo)率很高的固態(tài)電解質(zhì)材料，通過解析XRD數(shù)據(jù)發(fā)現(xiàn)它具有鈣鈦礦結(jié)構(gòu)竿滨，其中Li和La占據(jù)A site佳恬，也就是Ca的位置，那么我們就可以以CaTiO3的結(jié)構(gòu)為基礎(chǔ)于游，通過取代Ca來得到 LLTO 的結(jié)構(gòu)文件毁葱。
但是問題就來了，我們?nèi)〈蟮奈募卣紦?jù)比小于1贰剥，變成了無序結(jié)構(gòu)(disordered structure)倾剿，而我們使用VASP的軟件進(jìn)行計(jì)算的時(shí)候，原子的占據(jù)比必須是1蚌成，那我們怎么辦呢前痘？
一種通用的做法是建一個(gè)超胞(supercell)，也就是把現(xiàn)有的晶胞在三個(gè)方向上擴(kuò)展若干倍担忧，將分?jǐn)?shù)占據(jù)比的元素按比例分配到每個(gè)位置上芹缔，得到一系列的結(jié)構(gòu)后計(jì)算它們的能量，取能量最低者即為我們需要的結(jié)構(gòu)瓶盛。
比如我們在三個(gè)方向上都擴(kuò)展三倍最欠，得到一個(gè) 3x3x3 的超胞，此時(shí)結(jié)構(gòu)化學(xué)式就變成了 惩猫，對于電導(dǎo)率最高的
：
- Li : 0.33 x 27 = 8.91
9
- La : 0.56 x 27 = 15.12
- Ti : 1 x 27 = 27
- O : 3x27 = 81
此時(shí)LLTO的化學(xué)式可以寫作 芝硬，也就是說，我們用9個(gè)Li轧房、15個(gè)La和3個(gè)空位取代了A site上的27個(gè)Ca拌阴，此時(shí)我們得到了一個(gè)有132個(gè)原子的超胞，我們可以把它作為 LLTO 的晶胞奶镶，這個(gè)過程稱之為 ordering迟赃，通過這個(gè)過程我們把一個(gè)分?jǐn)?shù)占據(jù)的disordered的結(jié)構(gòu)變成了一個(gè)ordered的結(jié)構(gòu)。
但是問題就來了厂镇，我們?nèi)〈罂梢缘玫蕉嗌俜N結(jié)構(gòu)呢捺氢？做個(gè)排列組合： = 17383860 x 220 = 3824449200，對于簡單地取代剪撬，我們可以通過枚舉將所有結(jié)構(gòu)列出來，但對這種有天文數(shù)字的組合悠反，顯然我們不可能全部列出來残黑。
對這種情況馍佑，一種可行的方法是隨機(jī)取代得到若干種結(jié)構(gòu)，然后通過靜電能(electrostatic energy) 排序取能量最低的幾個(gè)結(jié)構(gòu)梨水，再對這幾個(gè)結(jié)構(gòu)進(jìn)行DFT計(jì)算以得到能量最低者拭荤。
兩種方法
下面我們分別來介紹一下兩種方法，我們以LSTF為例疫诽，即 舅世，這是一種固態(tài)電解質(zhì)，具有立方體型鈣鈦礦結(jié)構(gòu)奇徒。為了方便雏亚，我們忽略F，只考慮Li摩钙、Sr罢低、Ta、Hf胖笛、O五種元素网持，化合價(jià)分別為+1、+2长踊、+5功舀、+4、-2身弊。細(xì)心的同學(xué)可能會發(fā)現(xiàn)電荷不平衡辟汰，不過沒有關(guān)系。
首先我們要設(shè)計(jì)一下結(jié)構(gòu)佑刷，為了對稱這里我選擇建一個(gè) 3x3x3 的超胞莉擒，則該結(jié)構(gòu)化學(xué)式為：
pymatgen提供了多種建超胞的方法，可以通過structure的make_supercell()方法直接建瘫絮，也可以通過SupercellTransformation類返回新的結(jié)構(gòu)而不改變原結(jié)構(gòu)涨冀。
OrderDisorderedStructureTransformation類
我們首先來看一下如何使用OrderDisorderedStructureTransformation類進(jìn)行ordering，該類位于pymatgen.transformations.standard_transformations模塊中：
class OrderDisorderedStructureTransformation(algo=0, symmetrized_structures=False, no_oxi_states=False)
可以看到麦萤，該類有三個(gè)可選參數(shù)：
- algo鹿鳖，即algorithm，排序所用算法壮莹，默認(rèn)為0翅帜，可以不用管
- 0：ALGO_FAST
- 1：ALGO_COMPLETE
- 2：ALGO_BEST_FIRST
- symmetrized_structures：是否為對稱性結(jié)構(gòu)，默認(rèn)False
- no_oxi_states：是否在ordering前去除氧化態(tài)命满，默認(rèn)False涝滴，即不去除
氧化態(tài)即是否帶電荷，我們的結(jié)構(gòu)文件中是原子，因此不帶電荷歼疮，而計(jì)算靜電能需要知道每種元素的電荷杂抽，因?yàn)殪o電能是整體的電荷相互作用，沒有電荷也就無法計(jì)算靜電能了韩脏。我們可以通過Structure類的add_oxidation_state_by_element方法給元素添加氧化態(tài)缩麸，也可以通過該模塊下的OxidationStateDecorationTransformation類進(jìn)行添加。
屬性與方法
屬性：
- lowest_energy_structure：具有最低能量的結(jié)構(gòu)
- _all_structures：存有ordering后結(jié)構(gòu)的列表赡矢，結(jié)構(gòu)以字典的形式儲存杭朱。
- inverse： 沒什么意義，不用管
- is_one_to_many： 是否返回多個(gè)結(jié)果吹散，不用管
方法：apply_transformation(structure, return_ranked_list=False)
ordering的方法弧械，structure為我們要排序的結(jié)構(gòu)，return_ranked_list控制返回值：
- False：默認(rèn)值送浊，返回Ewald能最低的結(jié)構(gòu)
- 整數(shù)：返回該數(shù)量的結(jié)構(gòu)梦谜，以字典列表的形式，每個(gè)結(jié)構(gòu)儲存在一個(gè)字典中袭景，字典還包括其他信息唁桩，可以通過
structure鍵來訪問。
代碼：
import os
from pymatgen.io.cif import CifParser, CifWriter
from pymatgen.transformations.standard_transformations import SubstitutionTransformation, OrderDisorderedStructureTransformation, SupercellTransformation
# read in the CaTiO3 structure from cif file
filename = 'CaTiO3_mp-5827_symmetrized.cif'
parser = CifParser(filename)
init_structure = parser.get_structures()[0]
# add oxidation states to the structure
data = {"Ca":2, "Ti":4, "O":-2}
init_structure.add_oxidation_state_by_element(data)
# substitute to get the partial occupancied structure
species_map = {"Ca2+":{"Li1+":0.38, "Sr2+":0.44},"Ti4+":{"Ta5+":0.7,"Hf4+":0.3}}
substitutuin = SubstitutionTransformation(species_map)
structure = substitutuin.apply_transformation(init_structure)
``
# make a 3x3x3 supercell
structure.make_supercell(3)
'''
# you can also use the SupercellTransformation class to achieve it
sc = SupercellTransformation().from_scaling_factors(3,3,3)
supercell = sc.apply_transformation(structure)
'''
# ordering耸棒，here we set return_ranked_list as 100 to get 100 structures
order = OrderDisorderedStructureTransformation()
standard_structures = order.apply_transformation(structure,return_ranked_list=100)
# save ordered structures to files
if !os.path.exists('ordering'):
os.mkdir('ordering')
i = 1
for s in ordered_structures:
cif = CifWriter(s['structure'])
cif.write_file('ordering/structure_{:0>3d}.cif'.format(i))
i+=1
當(dāng)你運(yùn)行這段代碼的時(shí)候荒澡，你會發(fā)現(xiàn)程序會報(bào)錯(cuò)：
ValueError: Occupancy fractions not consistent with size of unit cell
這是為什么呢？我們查看一下該類的源碼,找到出錯(cuò)的位置：
for k, v in total_occupancy.items():
if abs(v - round(v)) > 0.25:
raise ValueError("Occupancy fractions not consistent "
"with size of unit cell")
再來看看文檔中對該類的說明：
simple rounding of the occupancies are performed, with no attempt made to achieve a target composition. This is usually not a problem for most ordering problems, but there can be times where rounding errors may result in structures that do not have the desired composition.
龜龜与殃，文檔中說會簡單地取整单山，應(yīng)該可以理解為四舍五入，然鵝幅疼，代碼中設(shè)置了限制米奸，超過0.25就不給算了，而對Li：0.38 x 27 = 10.26 比我們?nèi)〉?0大了0.26爽篷，已經(jīng)超過了0.25悴晰，我們對代碼做一下修改，把0.38改為0.37逐工，此時(shí)：0.37 x 27 = 9.99 就沒得問題了：
# substitute to get the partial occupancied structure
species_map = {"Ca2+":{"Li1+":0.37, "Sr2+":0.44},"Ti4+":{"Ta5+":0.7,"Hf4+":0.3}}
substitutuin = SubstitutionTransformation(species_map)
structure = substitutuin.apply_transformation(init_structure)
再次運(yùn)行程序铡溪，我的電腦陷入了沉思...
看一下說明：
The algorithm can currently compute approximately 5,000,000 permutations per minute.
這個(gè)算法每分鐘能處理大約5百萬個(gè)結(jié)果，但是我們有上億種結(jié)構(gòu)泪喊，因此需要程序跑幾分鐘棕硫。當(dāng)然如果你用超算或者服務(wù)器的話應(yīng)該會快一點(diǎn)。
程序運(yùn)行完后袒啼，你會發(fā)現(xiàn)腳本目錄下多了一個(gè) ordering 文件夾哈扮，里面有我們order后的結(jié)構(gòu)纬纪，文件名形如 structure_***.cif 的形式，結(jié)構(gòu)是根據(jù)ewald energy / atom由低到高排列的灶泵。
