|
|
本帖最后由 tjuptz 于 2021-11-10 23:09 编辑
这几个脚本还是去年学习ASE及xtb时，仿照冰释之川老师帖子 在ASE中使用xtb-6.3.0(pre.2)的py脚本 练习写的。由于只简单测试了可以跑，没有实际派上用场一直搁置。与其遗忘在角落，还不如分享出来抛砖引玉。
优化的脚本就是把上贴中的优化算法换换，主体内容就不贴了。MD的分别写了NVT、NPT对应不同控温控压算法的，这里贴出一个Berendsen热浴压浴npt_B.py的内容示例。当时现学现卖写的进度条和计算速度统计。
- import xtb
- from xtb import GFN0
- import ase
- from ase.io import read, iread, write
- from ase.units import fs, kB
- from ase.md.nptberendsen import NPTBerendsen
- from ase.md.velocitydistribution import MaxwellBoltzmannDistribution, Stationary, ZeroRotation
- from ase.io.trajectory import Trajectory
- from ase.md import MDLogger
- import os,sys,datetime
- start_time = datetime.datetime.now()
- # read molecular structure data, here from energy minimization
- basename, suffix = sys.argv[1].split(".")
- name = basename+'_npt'
- mol = read(sys.argv[1])
- total_atoms = mol.get_global_number_of_atoms()
- print('*** ASE program with xTB started at {0} ***\n'.format(start_time.strftime("%Y-%m-%d %H:%M:%S")))
- print('-> Task name: %50s' %name)
- print('-> Orthorhombic cell parameters: %31s' %mol.get_cell())
- print('-> Total number of atoms in a cell: %27s' %total_atoms)
- print('-> PBC dimensions: %45s' %mol.get_pbc())
- # set calculation level, charge, and magnetic moment
- total_charges = 0
- total_unpaired_electrons = 0
- mol.set_initial_charges(charges=[total_charges/total_atoms,]*total_atoms)
- mol.set_initial_magnetic_moments(magmoms=[total_unpaired_electrons/total_atoms,]*total_atoms)
- mol.calc = GFN0(accuracy=1, electronic_temperature=300, max_iterations=250)
- # setup the npt md
- nstep = 100
- dstep = 10
- istep = -10 # -dstep
- dt = 1 # fs
- temp = 400 # kB
- dyn = NPTBerendsen(mol, timestep=dt*fs, temperature=temp, taut=0.1*1000*fs,
- pressure=1.01325, taup=1.0*1000*fs, compressibility=4.57e-5, logfile=name+'.log')
- # Set the momenta corresponding to reference temperature
- MaxwellBoltzmannDistribution(mol, temp*kB, force_temp=True)
- Stationary(mol) # zero linear momentum
- #ZeroRotation(mol) # zero angular momentum
- # save the positions of all atoms after every dstep
- traj = Trajectory(name+'.traj', 'w', mol)
- dyn.attach(traj.write, interval=dstep)
- # save the logfile after every dstep
- dyn.attach(MDLogger(dyn, mol, name+'.log', header=True, stress=True,
- peratom=True, mode="a"), interval=dstep)
- # performance counter
- def perf_counter(dstep=dstep, nstep=nstep, start_time=start_time):
- """Function to counter the performance."""
- global istep
- if istep < 0 :
- istep += dstep
- else :
- istep += dstep
- d_time = (datetime.datetime.now() - start_time).seconds
- scale = 50
- i = int(scale * istep/nstep)
- prdt_end_time = datetime.datetime.now() + datetime.timedelta(seconds=d_time*nstep/istep)
- a = '*' * i
- b = '.' * (scale - i)
- c = (istep/nstep)*100
- print("\r{:^3.0f}%[{}->{}]{:.2f}s\n*** This task will ended at {} ***".format(c,a,b,d_time,prdt_end_time.strftime("%Y-%m-%d %H:%M:%S")), end="")
- dyn.attach(perf_counter, interval=dstep)
- # Now run the dynamics
- print('\n'+' NPT MD started... '.center(60,'~'))
- dyn.run(nstep)
- print('\n'+' NPT MD terminated... '.center(60,'~'))
- # write trajectory and final geometry to file
- write(name+"_ed.xyz", read(name+".traj"), columns=['symbols','positions'],write_results=False)
- write(name+"_traj.xyz", iread(name+".traj"), columns=['symbols','positions'],write_results=False)
- #os.remove('xtbrestart')
- #os.remove('wbo')
- #os.remove('gfn0.out')
- end_time = datetime.datetime.now()
- print('*** ASE program with xTB ended at {0} ***\n'.format(end_time.strftime("%Y-%m-%d %H:%M:%S")))
- tot_seconds = (end_time - start_time).seconds
- days = tot_seconds // 86400
- hours = (tot_seconds % 86400) // 3600
- minutes = (tot_seconds % 86400 % 3600)// 60
- seconds = tot_seconds % 60
- print(">> Elapsed time: {0:2d} days {1:2d} hours {2:2d} minutes {3:4.1f} seconds <<".format(days,hours,minutes,seconds))
- ps_per_day = dt * nstep * 86400 / tot_seconds / 1e3
- hour_per_ps = 24 / ps_per_day
- print(">> Performance: {0:.2f} ps/day {1:.2f} hour/ps <<".format(ps_per_day,hour_per_ps))
注：对于量子化学包括半经验的方法，我不太清楚ASE处理PBC的机制！像第一性原理程序，方法和基组本来就是处理周期性的，稍能理解。所以，不太清楚把这些程序当作力的计算器时候咋处理跨越边界的分子键。待大佬解惑。
另外，知乎有两篇文章是在ASE中使用xtb做优化和过渡态计算的，列出以供学习。
https://zhuanlan.zhihu.com/p/294508681
https://zhuanlan.zhihu.com/p/295273428
以下是脚本附件
nvt_NH.py
(3.62 KB, 下载次数 Times of downloads: 17)
nvt_B.py
(3.62 KB, 下载次数 Times of downloads: 23)
npt_NH_PR.py
(3.69 KB, 下载次数 Times of downloads: 19)
npt_B.py
(3.71 KB, 下载次数 Times of downloads: 19)
opt_Precon_ecf.py
(2.76 KB, 下载次数 Times of downloads: 17)
opt_Precon.py
(2.59 KB, 下载次数 Times of downloads: 16)
opt_BFGSs.py
(2.49 KB, 下载次数 Times of downloads: 16)
|
评分 Rate
-
查看全部评分 View all ratings
|
发表于 Post on 2021-11-10 23:04:31
|
只看该作者 Only view this author
|倒序浏览 Reverse view
|阅读模式 Reading model
收藏 Add to favorites
楼主 Author