在ASE中使用xtb-6.3.0(pre.2)的py脚本

冰释之川

最近想利用xtb优化一些周期性体系，可惜xtb软件的周期性模块比较弱，而且支持的优化引擎也只有“engine=inertial”，因而转向Atomic Simulation Environment(ASE) [https://wiki.fysik.dtu.dk/ase/]，利用其接口调用xtb来进行优化，并且支持的优化算法相对比较多（ BFGS, BFGSLineSearch, LBFGS, LBFGSLineSearch, GPMin, MDMin and FIRE），还能开启 1D，2D和3D的PBC，用起来真心美滋滋。

2020.11.18 新增xtb-6.3.0(pre.2)版对应的xtb-python模块安装教程:

1. 编辑 ~/.bashrc 增加xtb-6.3.0(pre.2)版动态库的环境变量并source ~/.bashrc：

1. export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/mnt/d/xtb_Linux/xtb_6.3.pre2/lib64

复制代码

2. 切换至xtb主体目录下的python子目录并运行:

1. $ python3 setup.py install --prefix=/home/yjy/.local

复制代码

下面提供ASE与xtb联用的py脚本，供有需求的童鞋参考借鉴：
(1) py脚本：

1. #!/usr/bin/python3
   
2. def labels_convertion(input_string):
   
3.   new_string = []
   
4.   tmp = input_string.split(",")
   
5.   for i in tmp:
   
6.     if "-" in i:
   
7.       start = int(i.split("-")[0]) - 1
   
8.       end = int(i.split("-")[1])
   
9.       for j in range(start,end):
   
10.         new_string.append(j)
    
11.     else:
    
12.       new_string.append(int(i) - 1)
    
13.   return new_string
    
14. 
    
15. from ase import Atoms
    
16. from ase.units import Hartree
    
17. from ase.constraints import FixAtoms
    
18. from ase.optimize import BFGS
    
19. from ase.io import read,write,iread
    
20. from xtb import GFN1
    
21. import os,sys,datetime
    
22. 
    
23. start_time = datetime.datetime.now()
    
24. 
    
25. name = sys.argv[1]
    
26. mol = read(name+".xyz", format='xyz')
    
27. total_atoms = len(mol.get_atomic_numbers())
    
28. 
    
29. # orthorhombic cell parameter
    
30. a, b, c = 7.8598, 7.8598, 19.5257
    
31. mol.cell = [a, b, c]
    
32. 
    
33. # PBC: periodic boundary conditions
    
34. mol.pbc = (True, True, True)
    
35. 
    
36. print('***  ASE program with xTB started at {0}  ***\n'.format(start_time.strftime("%Y-%m-%d %H:%M:%S")))
    
37. print('-> Task name:  %50s' %name)
    
38. print('-> Orthorhombic cell parameters:  %31s' %mol.get_cell())
    
39. print('-> Total number of atoms in a cell:  %27s' %total_atoms)
    
40. print('-> PBC dimensions:  %45s' %mol.get_pbc())
    
41. 
    
42. #  fix atoms (index starts from 1)
    
43. fixed_atoms = "1,3-4,6-17,19-20,22-33,35-36,38-49,51-52,54-64"
    
44. fix = FixAtoms(indices=labels_convertion(fixed_atoms))
    
45. mol.set_constraint(fix)
    
46. print('-> The atoms with following labels are fixed:\n{0}'.format(labels_convertion(fixed_atoms)))
    
47. 
    
48. # set calculation level, charge, and magnetic moment
    
49. total_charges = 0
    
50. total_unpaired_electrons = 0
    
51. mol.calc = GFN1(accuracy=1, electronic_temperature=300, max_iterations=250)
    
52. mol.set_initial_charges(charges=[total_charges/total_atoms,]*total_atoms)
    
53. mol.set_initial_magnetic_moments(magmoms=[total_unpaired_electrons/total_atoms,]*total_atoms)
    
54. 
    
55. # calculate single point energy (SPE) of the initial geometry
    
56. e_mol = mol.get_potential_energy()
    
57. print('\nTotal energy of the initial geomertry: \n{0:.6f} eV   <==>   {1:.8f} Hartree\n'.format(e_mol, e_mol/Hartree))
    
58. 
    
59. print('~~~ Geometry optimization started... ~~~')
    
60. opt = BFGS(mol, trajectory=name+".traj", restart=name+".pckl")
    
61. opt.run(fmax=0.001, steps=800)
    
62. print('~~~ Geometry optimization terminated... ~~~')
    
63. 
    
64. # calculate single point energy (SPE) of the optimized geometry
    
65. e_mol_opted = mol.get_potential_energy()
    
66. print('\nTotal energy of the optimized geometry: \n{0:.6f} eV   <==>   {1:.8f} Hartree\n'.format(e_mol_opted, e_mol_opted/Hartree))
    
67. 
    
68. write(name+"_opted.xyz", read(name+".traj"))
    
69. write(name+"_traj.xyz", iread(name+".traj"))
    
70. os.remove('xtbrestart')
    
71. os.remove('wbo')
    
72. 
    
73. end_time = datetime.datetime.now()
    
74. print('***  ASE program with xTB ended at {0}  ***'.format(end_time.strftime("%Y-%m-%d %H:%M:%S")))
    
75. 
    
76. tot_seconds = (end_time - start_time).seconds
    
77. days = tot_seconds // 86400
    
78. hours = (tot_seconds % 86400) // 3600
    
79. minutes = (tot_seconds % 86400 % 3600)// 60
    
80. seconds = tot_seconds % 60
    
81. print(">> Elapsed time: {0:2d} days {1:2d} hours {2:2d} minutes {3:4.1f} seconds <<".format(days,hours,minutes,seconds))
    

复制代码

（2）调用py脚本的PBS作业提交脚本：

1. #!/bin/bash
   
2. #PBS -N ASE_xtb
   
3. #PBS -l nodes=1:ppn=12
   
4. #PBS -l walltime=1440:00:00
   
5. #PBS -q GPU
   
6. 
   
7. export OMP_NUM_THREADS=12     # CPU cores
   
8. export MKL_NUM_THREADS=12     # CPU cores
   
9. export OMP_STACKSIZE=3000m    # memory
   
10. ulimit -s unlimited
    
11. source /home/yjy/intel/parallel_studio_xe_2020/psxevars.sh
    
12. 
    
13. cd $PBS_O_WORKDIR
    
14. touch jobID.$PBS_JOBID
    
15. 
    
16. FILENAME=Cell.py  # input file name
    
17. 
    
18. python3 ${FILENAME} ${FILENAME%.py} | tee ${FILENAME/%py/log} # running python script

复制代码

实例包：
ASE_xtb_test.7z (607.17 KB, 下载次数 Times of downloads: 162)

2020-4-20 22:44 上传 Uploaded

点击下载
Click to download

diamond_test.7z (6.03 KB, 下载次数 Times of downloads: 89)

2020-4-23 09:40 上传 Uploaded

点击下载
Click to download

wei

很好，能给个简单手册和算例吗？

tjuptz

请问您的固体是有机的还是无机（离子or共价or金属）？xtb可靠性如何？

冰释之川

tjuptz 发表于 2020-4-22 22:37
请问您的固体是有机的还是无机（离子or共价or金属）？xtb可靠性如何？

测试的是二氧化硅表面

冰释之川

wei 发表于 2020-4-22 16:43
很好，能给个简单手册和算例吗？

算例见附件压缩包，因为还在摸索中，就暂时不写tutorial了   其实算例和手册翻翻 ASE在线手册就行了，必要时再回头看看xtb在线手册中的python部分示例

aaronzjw

这个并行效率怎么样啊

tjuptz

冰释之川 发表于 2020-4-23 08:19
算例见附件压缩包，因为还在摸索中，就暂时不写tutorial了   其实算例和手册翻翻 ASE在线手册就行了 ...

请教下，为何  diamond的脚本里最后一段

1. write(name+"_opted.xyz", read(name+".traj"), columns=['symbols','positions'], write_results=False)
   
2. write(name+"_traj.xyz", iread(name+".traj"), columns=['symbols','positions'], write_results=False)
   
3. #os.remove('xtbrestart')
   
4. #os.remove('wbo')

复制代码

而SiO2的脚本里是这样？

1. write(name+"_opted.xyz", read(name+".traj"))
   
2. write(name+"_traj.xyz", iread(name+".traj"))
   
3. os.remove('xtbrestart')
   
4. os.remove('wbo')

复制代码

我看了write_xzy的源码，大概理解了二者的差异，不过没搞懂为啥产生的文件不一样。是GFN0和GFN1的输出不同导致的嘛？
PS: 我用SiO2脚本上述语句转换ASE+XTB GFN0的MD轨迹时候遇到了问题，用diamond脚本里的语句就行了。

冰释之川

tjuptz 发表于 2020-6-15 15:20
请教下，为何  diamond的脚本里最后一段而SiO2的脚本里是这样？我看了write_xzy的源码，大概理解了二者的 ...

因为我想让生成的*.xyz信息只包含原子坐标，其他信息不要写入进去

chrinide

赞一个

tjuptz

冰释之川 发表于 2020-6-15 15:48
因为我想让生成的*.xyz信息只包含原子坐标，其他信息不要写入进去

Soga，谢谢您。xyz文件里的信息确实有差异。

Aridea

冰姐，xtb的周期体系的md是真周期嘛

冰释之川

Aridea 发表于 2021-10-10 16:33
冰姐，xtb的周期体系的md是真周期嘛

不管是不是真，xtb里最好不要跑周期性体系，因为xtb里周期性模块比较垃圾。。。
你要用GFN1-xTB 跑周期性，干脆就用CP2K

Aridea

冰释之川 发表于 2021-10-11 08:14
不管是不是真，xtb里最好不要跑周期性体系，因为xtb里周期性模块比较垃圾。。。
你要用GFN1-xTB 跑周期 ...

冰姐，我昨天第一次跑md，用的用就是cp2k xtb结果上来就崩，我，，，￥%##%*@

neocc

请问冰冰姐，运行python ASE_xtb.py diamond 没有任何回显，运行ASE_xtb_new.py和ASE_xtb_MinimaHopping.py
python ASE_xtb_MinimaHopping.py diamond  
报错

1. ~~~ Global geometry optimization started... ~~~
   
2. Traceback (most recent call last):
   
3.   File "ASE_xtb_MinimaHopping.py", line 71, in <module>
   
4.     mh(totalsteps=10)
   
5.   File "/home/xxx/anaconda3/envs/cp2k/lib/python3.8/site-packages/ase/optimize/minimahopping.py", line 58, in __call__
   
6.     self._startup()
   
7.   File "/home/xxx/anaconda3/envs/cp2k/lib/python3.8/site-packages/ase/optimize/minimahopping.py", line 107, in _startup
   
8.     self._optimize()
   
9.   File "/home/xxx/anaconda3/envs/cp2k/lib/python3.8/site-packages/ase/optimize/minimahopping.py", line 255, in _optimize
   
10.     self._atoms.set_momenta(np.zeros(self._atoms.get_momenta().shape))
    
11.   File "/home/xxx/anaconda3/envs/cp2k/lib/python3.8/site-packages/ase/atoms.py", line 586, in set_momenta
    
12.     constraint.adjust_momenta(self, momenta)
    
13.   File "/home/xxx/anaconda3/envs/cp2k/lib/python3.8/site-packages/ase/constraints.py", line 78, in adjust_momenta
    
14.     self.adjust_forces(atoms, momenta)
    
15.   File "/home/xxx/anaconda3/envs/cp2k/lib/python3.8/site-packages/ase/constraints.py", line 162, in adjust_forces
    
16.     forces[self.index] = 0.0
    
17. IndexError: index 8 is out of bounds for axis 0 with size 8
    

复制代码

冰释之川

neocc 发表于 2021-11-17 23:15
请问冰冰姐，运行python ASE_xtb.py diamond 没有任何回显，运行ASE_xtb_new.py和ASE_xtb_MinimaHopping.py ...

好像和ASE版本还有关系，如果是用的最新版，估计要去翻翻ASE手册（https://wiki.fysik.dtu.dk/ase/），调整一下代码了。