在ASE中使用xtb-6.3.0(pre.2)的py脚本 - 计算化学公社

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

复制代码

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

复制代码

#!/usr/bin/python3
def labels_convertion(input_string):
new_string = []
tmp = input_string.split(",")
for i in tmp:
if "-" in i:
start = int(i.split("-")[0]) - 1
end = int(i.split("-")[1])
for j in range(start,end):
new_string.append(j)
else:
new_string.append(int(i) - 1)
return new_string
from ase import Atoms
from ase.units import Hartree
from ase.constraints import FixAtoms
from ase.optimize import BFGS
from ase.io import read,write,iread
from xtb import GFN1
import os,sys,datetime
start_time = datetime.datetime.now()
name = sys.argv[1]
mol = read(name+".xyz", format='xyz')
total_atoms = len(mol.get_atomic_numbers())
# orthorhombic cell parameter
a, b, c = 7.8598, 7.8598, 19.5257
mol.cell = [a, b, c]
# PBC: periodic boundary conditions
mol.pbc = (True, True, True)
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())
# fix atoms (index starts from 1)
fixed_atoms = "1,3-4,6-17,19-20,22-33,35-36,38-49,51-52,54-64"
fix = FixAtoms(indices=labels_convertion(fixed_atoms))
mol.set_constraint(fix)
print('-> The atoms with following labels are fixed:\n{0}'.format(labels_convertion(fixed_atoms)))
# set calculation level, charge, and magnetic moment
total_charges = 0
total_unpaired_electrons = 0
mol.calc = GFN1(accuracy=1, electronic_temperature=300, max_iterations=250)
mol.set_initial_charges(charges=[total_charges/total_atoms,]*total_atoms)
mol.set_initial_magnetic_moments(magmoms=[total_unpaired_electrons/total_atoms,]*total_atoms)
# calculate single point energy (SPE) of the initial geometry
e_mol = mol.get_potential_energy()
print('\nTotal energy of the initial geomertry: \n{0:.6f} eV <==> {1:.8f} Hartree\n'.format(e_mol, e_mol/Hartree))
print('~~~ Geometry optimization started... ~~~')
opt = BFGS(mol, trajectory=name+".traj", restart=name+".pckl")
opt.run(fmax=0.001, steps=800)
print('~~~ Geometry optimization terminated... ~~~')
# calculate single point energy (SPE) of the optimized geometry
e_mol_opted = mol.get_potential_energy()
print('\nTotal energy of the optimized geometry: \n{0:.6f} eV <==> {1:.8f} Hartree\n'.format(e_mol_opted, e_mol_opted/Hartree))
write(name+"_opted.xyz", read(name+".traj"))
write(name+"_traj.xyz", iread(name+".traj"))
os.remove('xtbrestart')
os.remove('wbo')
end_time = datetime.datetime.now()
print('*** ASE program with xTB ended at {0} ***'.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))

复制代码

#!/bin/bash
#PBS -N ASE_xtb
#PBS -l nodes=1:ppn=12
#PBS -l walltime=1440:00:00
#PBS -q GPU
export OMP_NUM_THREADS=12 # CPU cores
export MKL_NUM_THREADS=12 # CPU cores
export OMP_STACKSIZE=3000m # memory
ulimit -s unlimited
source /home/yjy/intel/parallel_studio_xe_2020/psxevars.sh
cd $PBS_O_WORKDIR
touch jobID.$PBS_JOBID
FILENAME=Cell.py # input file name
python3 ${FILENAME} ${FILENAME%.py} | tee ${FILENAME/%py/log} # running python script

复制代码

write(name+"_opted.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')

复制代码

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

复制代码

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

复制代码