ASE使用笔记(1)

ASE(Atomic Simulation Environment) is a set of tools and Python modules for setting up, manipulating, running, visualizing and analyzing atomistic simulations. The code is freely available under the GNU LGPL license.

写在前面

作为一个刚刚由机器学习转入第一原理的人，对一些物理图像还是不怎么清晰的，但是，工欲善其事必先利其器，既然要做DFT相关的事情了，把和DFT相关的这些Python lib琢磨一下，应该会提高后期工作的效率。

关于安装ASE

pip install 即可，中间遇到的报错如VC++14之类的问题可以网上搜索解决办法进行规避or下载提示需要的VC++14,很多与和第一原理做接口的Python库都需要这个（如Pymatgen也遇到了这样的问题）
对于VC++14的下载：复制这段内容后打开百度网盘手机App，操作更方便哦
链接：https://pan.baidu.com/s/1I0OTvHfeNhaqrnY_nCqfyQ 提取码：8n85
官网链接：https://wiki.fysik.dtu.dk/ase/index.html

ASE的简介

官网解释的很清楚了，不再赘述。

官网的第一个demo

计算吸附能，氮气分子在铜板板上的
本笔记是在Jupyter Notebook的IPython环境下运行的。

#### 如何得到一个氮气分子
from ase import Atoms
d = 1.10
molecule = Atoms('2N', positions=[(0., 0., 0.), (0., 0., d)])
#### 如何得到铜板板
from ase.build import fcc111
slab = fcc111('Cu', size=(4,4,2), vacuum=10.0)

#### 使用effective medium theory (EMT) calculator（有效介质理论）
from ase.calculators.emt import EMT
slab.calc = EMT()
molecule.calc = EMT()
#### 得到各自体系的孤立能量
e_slab = slab.get_potential_energy()
e_N2 = molecule.get_potential_energy()
#### 结构优化，这个例子是在on-top位置的
h = 1.85
add_adsorbate(slab, molecule, h, 'ontop')

#### 固定（fix）铜板板是不能动的，只有N可以调整
from ase.constraints import FixAtoms
constraint = FixAtoms(mask=[a.symbol != 'N' for a in slab])
slab.set_constraint(constraint)
#### 设置受力收敛标准为0.005
from ase.optimize import QuasiNewton
dyn = QuasiNewton(slab, trajectory='N2Cu.traj')
dyn.run(fmax=0.00005)

                Step[ FC]     Time          Energy          fmax
*Force-consistent energies used in optimization.
BFGSLineSearch:    0[  0] 14:33:13       11.689927*       1.0797
BFGSLineSearch:    1[  2] 14:33:13       11.670814*       0.4090
BFGSLineSearch:    2[  4] 14:33:13       11.625880*       0.0409
BFGSLineSearch:    3[  6] 14:33:13       11.625228*       0.0309
BFGSLineSearch:    4[  7] 14:33:13       11.625212*       0.0039
BFGSLineSearch:    5[  8] 14:33:13       11.625212*       0.0000

## 把结构信息存储和读取
from ase.io import write
write('slab.xyz', slab)
from ase.io import read
slab_from_file = read('slab.xyz')
## 结构的可视化
from ase.visualize import view
view(slab)

## 分子动力学，但是他的温度设置是如何的没有仔细看
from ase.md.verlet import VelocityVerlet
from ase import units
dyn = VelocityVerlet(molecule, dt=1.0 * units.fs)
for i in range(100):
     pot = molecule.get_potential_energy()
     kin = molecule.get_kinetic_energy()
     print('%2d: %.5f eV, %.5f eV, %.5f eV' % (i, pot + kin, pot, kin))
     dyn.run(steps=20)
0: 0.45841 eV, 0.32004 eV, 0.13837 eV
 1: 0.45945 eV, 0.41613 eV, 0.04332 eV
 2: 0.45790 eV, 0.31736 eV, 0.14054 eV
 3: 0.45831 eV, 0.34529 eV, 0.11302 eV
 4: 0.45924 eV, 0.39462 eV, 0.06462 eV
 5: 0.45870 eV, 0.34393 eV, 0.11477 eV
 6: 0.45933 eV, 0.40648 eV, 0.05285 eV
 7: 0.45741 eV, 0.27477 eV, 0.18264 eV
 8: 0.45975 eV, 0.44898 eV, 0.01077 eV
 9: 0.45782 eV, 0.27953 eV, 0.17828 eV
10: 0.46031 eV, 0.45943 eV, 0.00088 eV
11: 0.45762 eV, 0.26915 eV, 0.18847 eV
12: 0.45981 eV, 0.45633 eV, 0.00348 eV
13: 0.45737 eV, 0.26596 eV, 0.19141 eV
14: 0.45971 eV, 0.42801 eV, 0.03171 eV
15: 0.45845 eV, 0.32330 eV, 0.13515 eV
16: 0.45943 eV, 0.41328 eV, 0.04615 eV
17: 0.45795 eV, 0.32109 eV, 0.13686 eV
18: 0.45825 eV, 0.34118 eV, 0.11707 eV
19: 0.45928 eV, 0.39782 eV, 0.06146 eV
20: 0.45866 eV, 0.34048 eV, 0.11818 eV
21: 0.45939 eV, 0.41026 eV, 0.04913 eV
22: 0.45740 eV, 0.27294 eV, 0.18446 eV
23: 0.45976 eV, 0.45045 eV, 0.00931 eV
24: 0.45778 eV, 0.27753 eV, 0.18025 eV
25: 0.46032 eV, 0.45992 eV, 0.00040 eV
26: 0.45765 eV, 0.27054 eV, 0.18711 eV
27: 0.45981 eV, 0.45540 eV, 0.00441 eV
28: 0.45737 eV, 0.26702 eV, 0.19035 eV
29: 0.45965 eV, 0.42475 eV, 0.03491 eV
30: 0.45849 eV, 0.32661 eV, 0.13188 eV
31: 0.45940 eV, 0.41036 eV, 0.04903 eV
32: 0.45801 eV, 0.32490 eV, 0.13311 eV
33: 0.45819 eV, 0.33711 eV, 0.12107 eV
34: 0.45931 eV, 0.40097 eV, 0.05834 eV
35: 0.45862 eV, 0.33705 eV, 0.12157 eV
36: 0.45946 eV, 0.41395 eV, 0.04551 eV
37: 0.45739 eV, 0.27126 eV, 0.18612 eV
38: 0.45978 eV, 0.45182 eV, 0.00795 eV
39: 0.45775 eV, 0.27564 eV, 0.18211 eV
40: 0.46032 eV, 0.46022 eV, 0.00010 eV
41: 0.45768 eV, 0.27206 eV, 0.18562 eV
42: 0.45980 eV, 0.45435 eV, 0.00544 eV
43: 0.45737 eV, 0.26823 eV, 0.18914 eV
44: 0.45959 eV, 0.42137 eV, 0.03822 eV
45: 0.45853 eV, 0.32995 eV, 0.12858 eV
46: 0.45937 eV, 0.40739 eV, 0.05198 eV
47: 0.45806 eV, 0.32878 eV, 0.12928 eV
48: 0.45813 eV, 0.33310 eV, 0.12503 eV
49: 0.45934 eV, 0.40406 eV, 0.05527 eV
50: 0.45858 eV, 0.33364 eV, 0.12494 eV
51: 0.45952 eV, 0.41755 eV, 0.04197 eV
52: 0.45738 eV, 0.26973 eV, 0.18765 eV
53: 0.45979 eV, 0.45309 eV, 0.00670 eV
54: 0.45772 eV, 0.27386 eV, 0.18385 eV
55: 0.46032 eV, 0.46032 eV, 0.00000 eV
56: 0.45771 eV, 0.27370 eV, 0.18401 eV
57: 0.45979 eV, 0.45320 eV, 0.00658 eV
58: 0.45738 eV, 0.26960 eV, 0.18778 eV
59: 0.45953 eV, 0.41788 eV, 0.04164 eV
60: 0.45857 eV, 0.33332 eV, 0.12525 eV
61: 0.45934 eV, 0.40435 eV, 0.05499 eV
62: 0.45812 eV, 0.33273 eV, 0.12539 eV
63: 0.45807 eV, 0.32915 eV, 0.12892 eV
64: 0.45937 eV, 0.40710 eV, 0.05226 eV
65: 0.45854 eV, 0.33026 eV, 0.12827 eV
66: 0.45959 eV, 0.42105 eV, 0.03854 eV
67: 0.45737 eV, 0.26835 eV, 0.18902 eV
68: 0.45980 eV, 0.45425 eV, 0.00555 eV
69: 0.45768 eV, 0.27221 eV, 0.18548 eV
70: 0.46032 eV, 0.46024 eV, 0.00009 eV
71: 0.45775 eV, 0.27547 eV, 0.18228 eV
72: 0.45978 eV, 0.45195 eV, 0.00783 eV
73: 0.45739 eV, 0.27111 eV, 0.18627 eV
74: 0.45946 eV, 0.41429 eV, 0.04517 eV
75: 0.45861 eV, 0.33673 eV, 0.12189 eV
76: 0.45931 eV, 0.40126 eV, 0.05805 eV
77: 0.45818 eV, 0.33673 eV, 0.12145 eV
78: 0.45801 eV, 0.32526 eV, 0.13275 eV
79: 0.45939 eV, 0.41009 eV, 0.04931 eV
80: 0.45849 eV, 0.32692 eV, 0.13158 eV
81: 0.45965 eV, 0.42444 eV, 0.03521 eV
82: 0.45737 eV, 0.26713 eV, 0.19024 eV
83: 0.45981 eV, 0.45530 eV, 0.00450 eV
84: 0.45765 eV, 0.27068 eV, 0.18698 eV
85: 0.46032 eV, 0.45996 eV, 0.00036 eV
86: 0.45778 eV, 0.27735 eV, 0.18043 eV
87: 0.45976 eV, 0.45058 eV, 0.00918 eV
88: 0.45740 eV, 0.27278 eV, 0.18462 eV
89: 0.45940 eV, 0.41061 eV, 0.04879 eV
90: 0.45865 eV, 0.34016 eV, 0.11850 eV
91: 0.45928 eV, 0.39811 eV, 0.06116 eV
92: 0.45824 eV, 0.34079 eV, 0.11745 eV
93: 0.45796 eV, 0.32145 eV, 0.13651 eV
94: 0.45942 eV, 0.41301 eV, 0.04641 eV
95: 0.45845 eV, 0.32361 eV, 0.13484 eV
96: 0.45971 eV, 0.42771 eV, 0.03200 eV
97: 0.45737 eV, 0.26605 eV, 0.19132 eV
98: 0.45981 eV, 0.45625 eV, 0.00356 eV
99: 0.45762 eV, 0.26927 eV, 0.18835 eV

关于分子的一些优化和计算，本人不研究非周期体系，有时间了可以看一下

关于ASE调用其他第一原理计算软件如VASP

$ export ASE_VASP_COMMAND="mpiexec $HOME/vasp/bin/vasp_std"   ##可执行路径
$ export VASP_PP_PATH=$HOME/vasp/mypps   ##yanshiPP
$ export ASE_VASP_VDW=$HOME/<path-to-vdw_kernel.bindat-folder>  ##vdw文件路径

from ase.build import molecule
atoms = molecule('N2')    ## 这里还是拿个氮气作为例子
atoms.center(vacuum=5)
from ase.calculators.vasp import Vasp2

calc = Vasp2(xc='pbe',  # Select exchange-correlation functional
             encut=520, # Plane-wave cutoff
             kpts=(1, 1, 1)) # k-points

atoms.calc = calc
en = atoms.get_potential_energy()  # This call will start the calculation
print('Potential energy: {:.2f} eV'.format(en))


Potential energy: -16.59 eV

对于材料结构上的操控(感觉还是和吸附有点像)

from math import sqrt
from ase import Atoms
a = 3.55
atoms = Atoms('Ni4',
              cell=[sqrt(2) * a, sqrt(2) * a, 1.0, 90, 90, 120],
              pbc=(1, 1, 0),
              scaled_positions=[(0, 0, 0),
                                (0.5, 0, 0),
                                (0, 0.5, 0),
                                (0.5, 0.5, 0)])
atoms.center(vacuum=5.0, axis=2)

atoms.cell

Cell([[5.020458146424487, 0.0, 0.0], [-2.5102290732122423, 4.347844293440141, 0.0], [0.0, 0.0, 10.0]])

atoms.positions

array([[ 0.        ,  0.        ,  5.        ],
       [ 2.51022907,  0.        ,  5.        ],
       [-1.25511454,  2.17392215,  5.        ],
       [ 1.25511454,  2.17392215,  5.        ]])

行吧，和POSCAR基本吻合

from ase.visualize import view
atoms.write('slab.xyz')
view(atoms)

还可以看看扩胞后的，在命令行上使用的话需要输入

ase gui -r 3,3,2 slab.xyz

当然，我这个环境需要在命令行前加入一个感叹号

!ase gui -r 3,3,2 slab.xyz   ## 扩胞查看

## 如果加点东西进去（先加个原子进去）
h = 1.9
relative = (1 / 6, 1 / 6, 0.5)
absolute = np.dot(relative, atoms.cell) + (0, 0, h)
atoms.append('Ag')
atoms.positions[-1] = absolute
view(atoms)

## 创建了一个水分子进来
import numpy as np
from ase import Atoms
p = np.array(
    [[0.27802511, -0.07732213, 13.46649107],
     [0.91833251, -1.02565868, 13.41456626],
     [0.91865997, 0.87076761, 13.41228287],
     [1.85572027, 2.37336781, 13.56440907],
     [3.13987926, 2.3633134, 13.4327577],
     [1.77566079, 2.37150862, 14.66528237],
     [4.52240322, 2.35264513, 13.37435864],
     [5.16892729, 1.40357034, 13.42661052],
     [5.15567324, 3.30068395, 13.4305779],
     [6.10183518, -0.0738656, 13.27945071],
     [7.3856151, -0.07438536, 13.40814585],
     [6.01881192, -0.08627583, 12.1789428]])
c = np.array([[8.490373, 0., 0.],
              [0., 4.901919, 0.],
              [0., 0., 26.93236]])
W = Atoms('4(OH2)', positions=p, cell=c, pbc=[1, 1, 0])
W.write('WL.traj')

旋转操作

W.rotate(90, 'z', center=(0, 0, 0))
view(W)
W.wrap()
###########The wrap() method only works if periodic boundary conditions are enabled. We have a 2 percent lattice mismatch between Ni(111) and the water, so we scale the water in the plane to match the cell of the slab. The argument scale_atoms=True indicates that the atomic positions should be scaled with the unit cell. The default is scale_atoms=False indicating that the cartesian coordinates remain the same when the cell is changed
W.set_cell(slab.cell, scale_atoms=True)
zmin = W.positions[:, 2].min()
zmax = slab.positions[:, 2].max()
W.positions += (0, 0, zmax - zmin + 1.5)
interface = slab + W
interface.center(vacuum=6, axis=2)
interface.write('NiH2O.traj')
view(interface)

CIF和XYZ,POSCAR的转换

ASE的这个功能已经被Pymatgen整合了，两者用起来都很方便。
拿出你的一个cif文件，我这个从MP里直接下载了个钙钛矿。

f = open('TiPbO3.cif', 'r')
print(f.read())

# generated using pymatgen
data_TiPbO3
_symmetry_space_group_name_H-M   'P 1'
_cell_length_a   3.86121500
_cell_length_b   3.86121500
_cell_length_c   4.62706100
_cell_angle_alpha   90.00000000
_cell_angle_beta   90.00000000
_cell_angle_gamma   90.00000000
_symmetry_Int_Tables_number   1
_chemical_formula_structural   TiPbO3
_chemical_formula_sum   'Ti1 Pb1 O3'
_cell_volume   68.98476581
_cell_formula_units_Z   1
loop_
 _symmetry_equiv_pos_site_id
 _symmetry_equiv_pos_as_xyz
  1  'x, y, z'
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
  Ti  Ti0  1  0.50000000  0.50000000  0.52092300  1
  Pb  Pb1  1  0.00000000  0.00000000  0.96921200  1
  O  O2  1  0.50000000  0.50000000  0.14244900  1
  O  O3  1  0.00000000  0.50000000  0.62665800  1
  O  O4  1  0.50000000  0.00000000  0.62665800  1

import warnings
warnings.filterwarnings("ignore")
!ase gui TiPbO3.cif -o POSCAR  ##这个参数 -o表示out 即输出为什么
f = open('POSCAR', 'r')
print(f.read())

Ti Pb  O 
 1.0000000000000000
     3.8612150000000001    0.0000000000000000    0.0000000000000000
     0.0000000000000000    3.8612150000000001    0.0000000000000000
     0.0000000000000000    0.0000000000000000    4.6270610000000003
   1   1   3
Cartesian
  1.9306075000000000  1.9306075000000000  2.4103424973030001
  0.0000000000000000  0.0000000000000000  4.4846030459320003
  1.9306075000000000  1.9306075000000000  0.6591202123890000
  0.0000000000000000  1.9306075000000000  2.8995847921380005
  1.9306075000000000  0.0000000000000000  2.8995847921380005

发现这个POSCAR和我要的版本不一样，我的vasp是5.4.4，后来发现大师兄的一个博客提到了这个东西是需要去修改ASE的源代码的，大师兄的博客连接https://www.bigbrosci.com/2020/08/31/A17/
1） 找到ASE的安装目录，打开 ase/io 目录下的vasp5.py 文件；

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

3） 保存vasp.py 即可。

大师兄牛皮，解决问题。

总结

结合前面这堆操作，做吸附储氢，电池什么的伙伴们可以通过ASE来搭建输入的POSCAR文件了，下一篇主要来讨论下怎么使用ASE来扩胞和搭建自己需要的掺杂体系的POSCAR吧。