ASE使用笔记(1)

ASE使用笔记(1

写在前面关于安装ASEASE的简介官网的第一个demo关于ASE调用其他第一原理计算软件如VASP对于材料结构上的操控(感觉还是和吸附有点像)旋转操作CIF和XYZ,POSCAR的转换总结

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吧。