ABACUS ASE接口使用技巧（持更

jiyy

ASE 是一个用来构建模型、进行可视化以及后处理分析的Python工具包。 ABACUS 提供了相应的接口（https://gitlab.com/1041176461/ase-abacus），使 ASE 可以对 ABACUS 的输入输出进行分析，并调用 ABACUS 为Calculator进行计算。

本贴主要用于记录我用过的一些方法 （欢迎补充！！！）

1. 文件处理

之前 HongriTianqi 已经介绍了如何把cif文件转换为STRU文件 （详见：怎样把cif文件转换为STRU文件？）。我们也可以采用相同的方式读取ABACUS的running_*log文件，只需要修改format参数为"abacus-out"即可，例如：

1. from ase.io import read
   
2. atoms = read("running_scf.log", format='abacus-out')

复制代码

这个方法对于分析MD结果有比较大的帮助（详见：https://gitlab.com/1041176461/ase-abacus/-/issues/3），接口还支持读取MD续算后的running_md.log

2. 能带计算

（1）调用ABACUS calculator做SCF计算，输出电荷密度：

1. from ase.calculators.abacus import Abacus, AbacusProfile
   
2. from ase.build import bulk
   
3. 
   
4. Si = bulk('Si')
   
5. pp = {'Si':'Si_ONCV_PBE-1.0.upf'}
   
6. pseudo_dir = r'/home/jiyy/potential/SG15/SG15_upf'
   
7. basis_dir = r'/home/jiyy/GD_orb/dpsi-T-S'
   
8. basis = {'Si':'dpsi_Si.dat'}
   
9. kpts = {'size': [6, 6, 6], 'gamma': True}
   
10. parameters = {
    
11.    'xc': 'pbe',
    
12.    'ecutwfc': 100,
    
13.    'smearing_method': 'gaussian',
    
14.    'smearing_sigma': 0.01,
    
15.    'basis_type': 'lcao',
    
16.    'ks_solver': 'genelpa',
    
17.    'mixing_type': 'pulay',
    
18.    'scf_thr': 1e-8,
    
19.    'out_chg': 1,
    
20.    'calculation': 'scf',
    
21. }
    
22. 
    
23. profile = AbacusProfile(argv=['mpirun', '-n', '8', 'abacus'])
    
24. calc = Abacus(profile=profile, pp=pp, pseudo_dir=pseudo_dir, basis=basis, basis_dir=basis_dir, kpts=kpts, **parameters)
    
25. Si.calc = calc
    
26. Si.get_potential_energy()
    

复制代码

（2）获取SCF计算得到的Fermi能级，并生成高对称路径

1. fermi_level = calc.get_fermi_level()  # 读取Fermi能级
   
2. kpath = Si.cell.bandpath(npoints=100) # 生成高对称路径，总共100个k点

复制代码

（3）调用ABACUS calculator做NSCF

1. parameters.update({'calculation': 'nscf', 'init_chg': 'file', 'out_chg': 0, 'out_band': 1})  # 更新参数做能带计算

复制代码

（4）画能带

1. bs = BandStructure(path=kpath, energies=energies, reference=fermi_level) # 传入高对称路径，本征值，费米能
   
2. bsp = bs.subtract_reference() # 减去费米能
   
3. bsp.plot(emin=-6, emax=6, filename='band.png') # 画图
   

复制代码

3. 构建分子模型

ABACUS做分子计算需要人为构建真空层，接口提供了两种方法来构建真空层

（1）atoms.center(vacuum=None, axis=(0, 1, 2), about=None)

这个方法可以将atoms放置于Cell的中心位置，用户可以指定任意x, y, z方向上的真空层大小

1. from ase.build import molecule

复制代码

由该方法得到的STRU如下：

1. ATOMIC_SPECIES
   
2. H   1.008         H_ONCV_PBE-1.0.upf
   
3. 
   
4. NUMERICAL_ORBITAL
   
5. dpsi_H.dat
   
6. 
   
7. LATTICE_CONSTANT
   
8. 1.8897261258369282
   
9. 
   
10. LATTICE_VECTORS
    
11. 30.000000000      0.0000000000      0.0000000000      
    
12. 0.0000000000      30.000000000      0.0000000000      
    
13. 0.0000000000      0.0000000000      30.737166000      
    
14. 
    
15. ATOMIC_POSITIONS
    
16. Direct
    
17. 
    
18. H
    
19. 0.0000000000
    
20. 2
    
21. 0.5000000000 0.5000000000 0.5119910000 1 1 1
    
22. 0.5000000000 0.5000000000 0.4880090000 1 1 1

复制代码

（2）atoms.set_cell(cell, scale_atoms=False, apply_constraint=True)

这个方法相当于直接在STRU文件中修改LATTICE_VECTORS。

1. import numpy as np

复制代码

得到的分子可能会处于Cell的边界，此时HR/SR矩阵会分配到不同的cells，而不是集中在(0, 0, 0) cell，但并不影响基于HR/SR矩阵得到的计算结果。

由该方法得到的STRU如下：

1. ATOMIC_SPECIES
   
2. H   1.008         H_ONCV_PBE-1.0.upf
   
3. 
   
4. NUMERICAL_ORBITAL
   
5. dpsi_H.dat
   
6. 
   
7. LATTICE_CONSTANT
   
8. 1.8897261258369282
   
9. 
   
10. LATTICE_VECTORS
    
11. 30.000000000      0.0000000000      0.0000000000      
    
12. 0.0000000000      30.000000000      0.0000000000      
    
13. 0.0000000000      0.0000000000      30.000000000      
    
14. 
    
15. ATOMIC_POSITIONS
    
16. Direct
    
17. 
    
18. H
    
19. 0.0000000000
    
20. 2
    
21. 0.5000000000 0.5000000000 0.5245720000 1 1 1
    
22. 0.5000000000 0.5000000000 0.5000000000 1 1 1

复制代码

4. 基组的Delta测试
假设有C的SZ，DZP，TZDP，QZTP基组，用接口可以计算彼此之间的Delta值

1. import json
   
2. from pathlib import Path
   
3. 
   
4. from ase.collections import dcdft
   
5. from ase.calculators.abacus import Abacus, AbacusProfile
   
6. from ase.io import Trajectory, read
   
7. from ase.eos import EquationOfState as EOS
   
8. from ase.units import kJ
   
9. from ase.utils.deltacodesdft import delta
   
10. 
    
11. # 计算七个能量点，并保存至traj文件
    
12. def cal_EOS(symbol, parameters, labels):
    
13.     for label, orb in labels.items():
    
14.         traj = Trajectory(f'{symbol}-{label}.traj', 'w')
    
15.         if label == 'PW':
    
16.             parameters.update({
    
17.                 'basis_type': 'pw',
    
18.                 'ks_solver': 'cg',
    
19.             })
    
20.         else:
    
21.             basis_dir = f"/home/jiyy/Orb/Gen/C/Orbital_C_{label}"
    
22.             basis = {symbol: orb}
    
23.             parameters.update({
    
24.                 'basis_type': 'lcao',
    
25.                 'ks_solver': 'genelpa',
    
26.                 'basis': basis,
    
27.                 'basis_dir': basis_dir
    
28.             })
    
29.         for s in range(94, 108, 2):
    
30.             atoms = dcdft[symbol]
    
31.             atoms.set_cell(atoms.cell * (s / 100)**(1 / 3), scale_atoms=True)
    
32.             atoms.calc = Abacus(profile=profile, directory=f'{symbol}-{label}',
    
33.                                 scaled=False, **parameters)
    
34.             atoms.get_potential_energy()
    
35.             traj.write(atoms)
    
36. 
    
37. #读取traj文件，用Birch-Murnaghan公式拟合体积模量
    
38. def fit(symbol, label):
    
39.     V = []
    
40.     E = []
    
41.     for atoms in read(f'{symbol}-{label}.traj@:'):
    
42.         V.append(atoms.get_volume() / len(atoms))
    
43.         E.append(atoms.get_potential_energy() / len(atoms))
    
44.     eos = EOS(V, E, 'birchmurnaghan')
    
45.     eos.fit(warn=False)
    
46.     e0, B, Bp, v0 = eos.eos_parameters
    
47.     return e0, v0, B, Bp
    
48. 
    
49. # 将基组结果和Wien2K数据保存到json文件
    
50. def save_json(symbol):
    
51.     data = {f'{symbol}': {}}
    
52.     labels = []
    
53.     for path in Path().glob('*.traj'):
    
54.         filename = path.stem
    
55.         symbol, label = filename.split('-')
    
56.         labels.append(label)
    
57.         e0, v0, B, Bp = fit(symbol, label)
    
58.         sub_data = {f'{label}': {
    
59.             'volume': v0,
    
60.             'B': B * 1e24 / kJ,
    
61.             'Bp': Bp
    
62.         }}
    
63.         data[symbol].update(sub_data)
    
64. 
    
65.     data[symbol]['wien2k'] = {
    
66.         'volume': dcdft.data[symbol]['wien2k_volume'],
    
67.         'B': dcdft.data[symbol]['wien2k_B'],
    
68.         'Bp': dcdft.data[symbol]['wien2k_Bp'],
    
69.     }
    
70. 
    
71.     with open(f'{symbol}.json', 'w') as fp:
    
72.         json.dump(data, fp)
    
73. 
    
74. 
    
75. def get_tuple(data, label):
    
76.     return (data[label]['volume'], data[label]['B'], data[label]['Bp'])
    
77. 
    
78. 
    
79. if __name__ == "__main__":
    
80.     import sys
    
81.     xc = sys.argv[1]
    
82.     mpi = sys.argv[2]
    
83.     omp = sys.argv[3]
    
84. 
    
85.     profile = AbacusProfile(
    
86.         argv=['mpirun', '-env', f'OMP_NUM_THREADS={omp}', '-n', f'{mpi}', '/home/jiyy/software/ABACUS/abacus-develop-2.3.5/bin/abacus'])
    
87. 
    
88.     symbol = 'C'
    
89.     pp = {"C": "C_ONCV_PBE-1.0.upf"}
    
90.     pseudo_dir = "/home/jiyy/PP/SG15"
    
91.     kpts = [8, 8, 8]
    
92.     parameters = {
    
93.         'xc': f'{xc}',
    
94.         'ecutwfc': 100,
    
95.         'smearing_method': 'gaussian',
    
96.         'smearing_sigma': 0.01,
    
97.         'mixing_type': 'pulay',
    
98.         'scf_thr': 1e-8,
    
99.         'out_chg': 1,
    
100.         'kpts': kpts,
     
101.         'pp': pp,
     
102.         'pseudo_dir': pseudo_dir,
     
103.     }
     
104.     labels = {'PW': '',
     
105.               'Z': 'C_gga_8au_100Ry_1s1p.orb',
     
106.               'DZP': 'C_gga_8au_100Ry_2s2p1d.orb',
     
107.               'TZDP': 'C_gga_8au_100Ry_3s3p2d.orb',
     
108.               'QZTP': 'C_gga_8au_100Ry_4s4p3d.orb',
     
109.               }
     
110.     cal_EOS(symbol, parameters, labels)
     
111.     save_json(symbol)
     
112. 
     
113. # 计算彼此delta值
     
114.     with open(f'{symbol}.json', 'r') as fp:
     
115.         data = json.load(fp)[symbol]
     
116.     pw = get_tuple(data, 'PW')
     
117.     wien2k = get_tuple(data, 'wien2k')
     
118.     z = get_tuple(data, 'Z')
     
119.     dzp = get_tuple(data, 'DZP')
     
120.     tzdp = get_tuple(data, 'TZDP')
     
121.     qztp = get_tuple(data, 'QZTP')
     
122.     delta_data = {
     
123.         'pw-wien2k': delta(*pw, *wien2k),
     
124.         'z-pw': delta(*z, *pw),
     
125.         'dzp-pw': delta(*dzp, *pw),
     
126.         'tzdp': delta(*tzdp, *pw),
     
127.         'qztp': delta(*qztp, *pw)
     
128.     }
     
129.     with open(f'delta-{symbol}.json', 'w') as fp:
     
130.         json.dump(delta_data, fp)

复制代码

delta-C.json:

1. {"pw-wien2k": 0.3519292445820627, "z-pw": 6.826961630421406, "dzp-pw": 0.0590395574697385, "tzdp": 0.026290147268529143, "qztp": 0.04396233617693006}

复制代码

量子御坂

在AISI发现很多小伙伴不知道ASE是有GUI的（当然其实最开始我也不知道），所以借着这个板块说一下基于ASE-ABACUS接口进行ABACUS计算结构的可视化
ASE的GUI可以在命令行直接调用。当然也可以在代码写，对应模块为 ase.visualize.view()。
在Window的Anaconda Prompt（如果python环境由conda管理），mac, WSL2, Linux的终端上都可以借助ASE完成对ABACUS输入结构和输出结构的可视化。后两者由于对tkinker可视化库的支持缺失，可能需要安装一些依赖组件，如

1. sudo apt-get install libreadline-dev

复制代码

比如，如下MgO结构是一个STRU文件：

1. ATOMIC_SPECIES
   
2. Mg 24.305 Mg_ONCV_PBE-1.0.upf
   
3. O  15.999 O_ONCV_PBE-1.0.upf
   
4. 
   
5. NUMERICAL_ORBITAL
   
6. Mg_gga_8au_100Ry_4s2p1d.orb
   
7. O_gga_8au_100Ry_2s2p1d.orb
   
8. 
   
9. LATTICE_CONSTANT
   
10. 1.8897259886    # 1.8897259886 Bohr = 1.0 Angstrom
    
11. 
    
12. LATTICE_VECTORS
    
13. 4.27957 0.00000 0.00000
    
14. 0.00000 4.27957 0.00000
    
15. 0.00000 0.00000 4.27957
    
16. 
    
17. ATOMIC_POSITIONS
    
18. Direct
    
19. Mg
    
20. 0.0
    
21. 4
    
22. 0.0 0.0 0.0 1 1 1
    
23. 0.0 0.5 0.5 1 1 1
    
24. 0.5 0.0 0.5 1 1 1
    
25. 0.5 0.5 0.0 1 1 1
    
26. 
    
27. O
    
28. 0.0
    
29. 4
    
30. 0.5 0.0 0.0 1 1 1
    
31. 0.5 0.5 0.5 1 1 1
    
32. 0.0 0.0 0.5 1 1 1
    
33. 0.0 0.5 0.0 1 1 1
    

复制代码

可以通过一行简单的命令可视化

1. ase -T gui STRU

复制代码

我们就能直接看到可视化结果了


可以通过鼠标右键进行面外转动，通过shift+鼠标右键进行平移，鼠标左键选中原子，ctrl+鼠标左键可以选中多个原子，并随着选中的原子数目增多，逐步显示键长/键角/二面角。
可视化窗口还支持直接操作各种选项。

可视化一些更复杂的结构也完全没问题，比如ASE可视化一个三百原子的沸石分子筛体系


输入文件STRU可以可视化，输出文件running*.log当然也可以可视化。比如

1. ase -T gui running_relax.log

复制代码

类似的方法也可以可视化ABACUS的MD轨迹，以及可视化ASE调用ABACUS进行NEB等计算的轨迹。