在ASE中使用xtb-6.3.0(pre.2)做优化或MD的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))

复制代码