计算化学公社

 找回密码 Forget password
 注册 Register
Views: 9874|回复 Reply: 17

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

[复制链接 Copy URL]

1061

帖子

16

威望

5784

eV
积分
7165

Level 6 (一方通行)

計算化学の社畜

发表于 Post on 2020-4-20 22:44:22 | 显示全部楼层 Show all |阅读模式 Reading model
本帖最后由 冰释之川 于 2020-12-29 22:11 编辑

最近想利用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. from ase import Atoms
  15. from ase.units import Hartree
  16. from ase.constraints import FixAtoms
  17. from ase.optimize import BFGS
  18. from ase.io import read,write,iread
  19. from xtb import GFN1
  20. import os,sys,datetime

  21. start_time = datetime.datetime.now()

  22. name = sys.argv[1]
  23. mol = read(name+".xyz", format='xyz')
  24. total_atoms = len(mol.get_atomic_numbers())

  25. # orthorhombic cell parameter
  26. a, b, c = 7.8598, 7.8598, 19.5257
  27. mol.cell = [a, b, c]

  28. # PBC: periodic boundary conditions
  29. mol.pbc = (True, True, True)

  30. print('***  ASE program with xTB started at {0}  ***\n'.format(start_time.strftime("%Y-%m-%d %H:%M:%S")))
  31. print('-> Task name:  %50s' %name)
  32. print('-> Orthorhombic cell parameters:  %31s' %mol.get_cell())
  33. print('-> Total number of atoms in a cell:  %27s' %total_atoms)
  34. print('-> PBC dimensions:  %45s' %mol.get_pbc())

  35. #  fix atoms (index starts from 1)
  36. fixed_atoms = "1,3-4,6-17,19-20,22-33,35-36,38-49,51-52,54-64"
  37. fix = FixAtoms(indices=labels_convertion(fixed_atoms))
  38. mol.set_constraint(fix)
  39. print('-> The atoms with following labels are fixed:\n{0}'.format(labels_convertion(fixed_atoms)))

  40. # set calculation level, charge, and magnetic moment
  41. total_charges = 0
  42. total_unpaired_electrons = 0
  43. mol.calc = GFN1(accuracy=1, electronic_temperature=300, max_iterations=250)
  44. mol.set_initial_charges(charges=[total_charges/total_atoms,]*total_atoms)
  45. mol.set_initial_magnetic_moments(magmoms=[total_unpaired_electrons/total_atoms,]*total_atoms)

  46. # calculate single point energy (SPE) of the initial geometry
  47. e_mol = mol.get_potential_energy()
  48. print('\nTotal energy of the initial geomertry: \n{0:.6f} eV   <==>   {1:.8f} Hartree\n'.format(e_mol, e_mol/Hartree))

  49. print('~~~ Geometry optimization started... ~~~')
  50. opt = BFGS(mol, trajectory=name+".traj", restart=name+".pckl")
  51. opt.run(fmax=0.001, steps=800)
  52. print('~~~ Geometry optimization terminated... ~~~')

  53. # calculate single point energy (SPE) of the optimized geometry
  54. e_mol_opted = mol.get_potential_energy()
  55. print('\nTotal energy of the optimized geometry: \n{0:.6f} eV   <==>   {1:.8f} Hartree\n'.format(e_mol_opted, e_mol_opted/Hartree))

  56. write(name+"_opted.xyz", read(name+".traj"))
  57. write(name+"_traj.xyz", iread(name+".traj"))
  58. os.remove('xtbrestart')
  59. os.remove('wbo')

  60. end_time = datetime.datetime.now()
  61. print('***  ASE program with xTB ended at {0}  ***'.format(end_time.strftime("%Y-%m-%d %H:%M:%S")))

  62. tot_seconds = (end_time - start_time).seconds
  63. days = tot_seconds // 86400
  64. hours = (tot_seconds % 86400) // 3600
  65. minutes = (tot_seconds % 86400 % 3600)// 60
  66. seconds = tot_seconds % 60
  67. 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. export OMP_NUM_THREADS=12     # CPU cores
  7. export MKL_NUM_THREADS=12     # CPU cores
  8. export OMP_STACKSIZE=3000m    # memory
  9. ulimit -s unlimited
  10. source /home/yjy/intel/parallel_studio_xe_2020/psxevars.sh

  11. cd $PBS_O_WORKDIR
  12. touch jobID.$PBS_JOBID

  13. FILENAME=Cell.py  # input file name

  14. python3 ${FILENAME} ${FILENAME%.py} | tee ${FILENAME/%py/log} # running python script
复制代码
实例包:
ASE_xtb_test.7z (607.17 KB, 下载次数 Times of downloads: 148)

评分 Rate

参与人数
Participants 11
eV +54 收起 理由
Reason
Aridea + 5 GJ!
zsu007 + 5 牛!
chrinide + 5 GJ!
biogon + 5 とてもいい!
bobosiji + 4 不明觉厉
tjuptz + 5 好物!
hebrewsnabla + 5
Jaydu1996 + 2 好物!
978142355 + 5 GJ!
itpfeng + 5 赞!
sobereva + 8

查看全部评分 View all ratings

Stand on the shoulders of giants

90

帖子

0

威望

1653

eV
积分
1743

Level 5 (御坂)

发表于 Post on 2020-4-22 16:43:56 | 显示全部楼层 Show all
很好,能给个简单手册和算例吗?

489

帖子

1

威望

3222

eV
积分
3731

Level 5 (御坂)

发表于 Post on 2020-4-22 22:37:41 | 显示全部楼层 Show all
请问您的固体是有机的还是无机(离子or共价or金属)?xtb可靠性如何?

1061

帖子

16

威望

5784

eV
积分
7165

Level 6 (一方通行)

計算化学の社畜

 楼主 Author| 发表于 Post on 2020-4-23 08:16:50 | 显示全部楼层 Show all
tjuptz 发表于 2020-4-22 22:37
请问您的固体是有机的还是无机(离子or共价or金属)?xtb可靠性如何?

测试的是二氧化硅表面
Stand on the shoulders of giants

1061

帖子

16

威望

5784

eV
积分
7165

Level 6 (一方通行)

計算化学の社畜

 楼主 Author| 发表于 Post on 2020-4-23 08:19:17 | 显示全部楼层 Show all
wei 发表于 2020-4-22 16:43
很好,能给个简单手册和算例吗?

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

18

帖子

0

威望

313

eV
积分
331

Level 3 能力者

发表于 Post on 2020-5-20 21:57:59 | 显示全部楼层 Show all
这个并行效率怎么样啊

489

帖子

1

威望

3222

eV
积分
3731

Level 5 (御坂)

发表于 Post on 2020-6-15 15:20:43 | 显示全部楼层 Show all
本帖最后由 tjuptz 于 2020-6-15 15:25 编辑
冰释之川 发表于 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脚本里的语句就行了。

1061

帖子

16

威望

5784

eV
积分
7165

Level 6 (一方通行)

計算化学の社畜

 楼主 Author| 发表于 Post on 2020-6-15 15:48:41 | 显示全部楼层 Show all
tjuptz 发表于 2020-6-15 15:20
请教下,为何  diamond的脚本里最后一段而SiO2的脚本里是这样?我看了write_xzy的源码,大概理解了二者的 ...

因为我想让生成的*.xyz信息只包含原子坐标,其他信息不要写入进去
Stand on the shoulders of giants

272

帖子

0

威望

3941

eV
积分
4213

Level 6 (一方通行)

发表于 Post on 2020-6-15 16:50:44 | 显示全部楼层 Show all
赞一个

489

帖子

1

威望

3222

eV
积分
3731

Level 5 (御坂)

发表于 Post on 2020-6-15 19:10:22 | 显示全部楼层 Show all
冰释之川 发表于 2020-6-15 15:48
因为我想让生成的*.xyz信息只包含原子坐标,其他信息不要写入进去

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

132

帖子

0

威望

701

eV
积分
833

Level 4 (黑子)

发表于 Post on 2021-10-10 16:33:53 | 显示全部楼层 Show all
冰姐,xtb的周期体系的md是真周期嘛

1061

帖子

16

威望

5784

eV
积分
7165

Level 6 (一方通行)

計算化学の社畜

 楼主 Author| 发表于 Post on 2021-10-11 08:14:21 | 显示全部楼层 Show all
Aridea 发表于 2021-10-10 16:33
冰姐,xtb的周期体系的md是真周期嘛

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

132

帖子

0

威望

701

eV
积分
833

Level 4 (黑子)

发表于 Post on 2021-10-11 13:38:50 | 显示全部楼层 Show all
冰释之川 发表于 2021-10-11 08:14
不管是不是真,xtb里最好不要跑周期性体系,因为xtb里周期性模块比较垃圾。。。
你要用GFN1-xTB 跑周期 ...

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

114

帖子

0

威望

2435

eV
积分
2549

Level 5 (御坂)

发表于 Post on 2021-11-17 23:15:46 | 显示全部楼层 Show all
请问冰冰姐,运行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
复制代码

1061

帖子

16

威望

5784

eV
积分
7165

Level 6 (一方通行)

計算化学の社畜

 楼主 Author| 发表于 Post on 2021-11-18 08:08:15 | 显示全部楼层 Show all
本帖最后由 冰释之川 于 2021-11-18 08:09 编辑
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/),调整一下代码了。
Stand on the shoulders of giants

本版积分规则 Credits rule

手机版 Mobile version|北京科音自然科学研究中心 Beijing Kein Research Center for Natural Sciences|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )|网站地图

GMT+8, 2023-2-5 08:19 , Processed in 0.238820 second(s), 25 queries .

快速回复 返回顶部 返回列表 Return to list