Skip to content

Commit 03ad309

Browse files
alongdalongd
committed
Added the Arkane xTB ESS adapter
1 parent a35dd78 commit 03ad309

3 files changed

Lines changed: 321 additions & 1 deletion

File tree

‎arkane/ess/__init__.py‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -39,3 +39,4 @@
3939
from arkane.ess.psi4_parser import Psi4Log
4040
from arkane.ess.qchem import QChemLog
4141
from arkane.ess.terachem import TeraChemLog
42+
from arkane.ess.xtb import XTBLog

‎arkane/ess/factory.py‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -102,9 +102,12 @@ def ess_factory(fullpath: str,
102102
elif 'terachem' in line:
103103
ess_name = 'TeraChemLog'
104104
break
105+
elif 'x t b' in line or 'xtb version' in line:
106+
ess_name = 'XTBLog'
107+
break
105108
line = f.readline().lower()
106109
if ess_name is None:
107110
raise InputError(f'The file at {fullpath} could not be identified as a '
108-
f'Gaussian, Molpro, Orca, Psi4, QChem, or TeraChem log file.')
111+
f'Gaussian, Molpro, Orca, Psi4, QChem, TeraChem, or xTB log file.')
109112

110113
return _registered_ess_adapters[ess_name](path=fullpath, check_for_errors=check_for_errors, scratch_directory=scratch_directory)

‎arkane/ess/xtb.py‎

Lines changed: 316 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,316 @@
1+
#!/usr/bin/env python3
2+
3+
###############################################################################
4+
# #
5+
# RMG - Reaction Mechanism Generator #
6+
# #
7+
# Copyright (c) 2002-2023 Prof. William H. Green (whgreen@mit.edu), #
8+
# Prof. Richard H. West (r.west@neu.edu) and the RMG Team (rmg_dev@mit.edu) #
9+
# #
10+
# Permission is hereby granted, free of charge, to any person obtaining a #
11+
# copy of this software and associated documentation files (the 'Software'), #
12+
# to deal in the Software without restriction, including without limitation #
13+
# the rights to use, copy, modify, merge, publish, distribute, sublicense, #
14+
# and/or sell copies of the Software, and to permit persons to whom the #
15+
# Software is furnished to do so, subject to the following conditions: #
16+
# #
17+
# The above copyright notice and this permission notice shall be included in #
18+
# all copies or substantial portions of the Software. #
19+
# #
20+
# THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EXPRESS OR #
21+
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, #
22+
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE #
23+
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER #
24+
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING #
25+
# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER #
26+
# DEALINGS IN THE SOFTWARE. #
27+
# #
28+
###############################################################################
29+
30+
"""
31+
Arkane ESS adapter for xTB (extended Tight Binding) output files.
32+
Supports geometry, energy, frequencies, and conformer loading.
33+
"""
34+
35+
import logging
36+
import re
37+
38+
import numpy as np
39+
40+
import rmgpy.constants as constants
41+
from rmgpy.statmech import (
42+
Conformer,
43+
HarmonicOscillator,
44+
IdealGasTranslation,
45+
LinearRotor,
46+
NonlinearRotor,
47+
)
48+
49+
from arkane.common import (
50+
get_element_mass,
51+
get_principal_moments_of_inertia,
52+
symbol_by_number,
53+
)
54+
from arkane.exceptions import LogError
55+
from arkane.ess.adapter import ESSAdapter
56+
from arkane.ess.factory import register_ess_adapter
57+
58+
logger = logging.getLogger(__name__)
59+
60+
61+
class XTBLog(ESSAdapter):
62+
"""
63+
Arkane ESS adapter for xTB output files.
64+
65+
Parses output from the xtb program (https://github.com/grimme-lab/xtb),
66+
including geometry, electronic energy, vibrational frequencies, and
67+
thermochemical data.
68+
"""
69+
70+
def check_for_errors(self):
71+
"""Check for common xTB errors in the output file."""
72+
with open(self.path, 'r') as f:
73+
for line in f:
74+
if 'ERROR' in line and 'SETUP' not in line:
75+
raise LogError(f'xTB error found in {self.path}: {line.strip()}')
76+
if 'abnormal termination' in line.lower():
77+
raise LogError(f'xTB job terminated abnormally in {self.path}')
78+
79+
def get_number_of_atoms(self) -> int:
80+
"""Return the number of atoms from the xTB output."""
81+
with open(self.path, 'r') as f:
82+
for line in f:
83+
if 'number of atoms' in line:
84+
return int(line.split()[-1])
85+
raise LogError(f'Could not find the number of atoms in {self.path}')
86+
87+
def load_geometry(self):
88+
"""
89+
Load the molecular geometry from the xTB output file.
90+
91+
For optimization outputs, reads the final V2000 structure block.
92+
For frequency outputs, reads the initial coordinate input section.
93+
94+
Returns:
95+
Tuple of (coord, number, mass):
96+
coord: np.ndarray (n, 3) in Angstroms
97+
number: np.ndarray (n,) atomic numbers
98+
mass: np.ndarray (n,) atomic masses in amu
99+
"""
100+
coord, number, mass = [], [], []
101+
102+
with open(self.path, 'r') as f:
103+
lines = f.readlines()
104+
105+
# Try V2000 format (from geometry optimization output)
106+
for i, line in enumerate(lines):
107+
if 'V2000' in line:
108+
# Parse the counts line: "n_atoms n_bonds 0 ..."
109+
parts = line.split()
110+
n_atoms = int(parts[0])
111+
coord, number, mass = [], [], []
112+
for j in range(i + 1, i + 1 + n_atoms):
113+
tokens = lines[j].split()
114+
x, y, z = float(tokens[0]), float(tokens[1]), float(tokens[2])
115+
symbol = tokens[3]
116+
coord.append([x, y, z])
117+
mass_i, num_i = get_element_mass(symbol)
118+
number.append(num_i)
119+
mass.append(mass_i)
120+
121+
if not coord:
122+
# Fall back: parse from the "Calculation Setup" section
123+
# xTB prints "ID Z sym. atoms" followed by element list,
124+
# and coordinates may come from the input SDF parsed at the start.
125+
# Try parsing the Cartesian coordinate block if present
126+
for i, line in enumerate(lines):
127+
if 'final xyz coordinates' in line.lower() or 'cartesian coordinates' in line.lower():
128+
coord, number, mass = [], [], []
129+
for j in range(i + 1, len(lines)):
130+
tokens = lines[j].split()
131+
if len(tokens) < 4:
132+
break
133+
try:
134+
symbol = tokens[0]
135+
x, y, z = float(tokens[1]), float(tokens[2]), float(tokens[3])
136+
coord.append([x, y, z])
137+
mass_i, num_i = get_element_mass(symbol)
138+
number.append(num_i)
139+
mass.append(mass_i)
140+
except (ValueError, KeyError):
141+
break
142+
143+
if not coord:
144+
raise LogError(f'Could not find geometry in {self.path}')
145+
146+
return np.array(coord, dtype=np.float64), np.array(number, dtype=int), np.array(mass, dtype=np.float64)
147+
148+
def load_energy(self, zpe_scale_factor=1.):
149+
"""
150+
Load the electronic energy in J/mol from the xTB output.
151+
152+
Returns the last "total energy" value found in the SUMMARY block.
153+
The zero-point energy is NOT included.
154+
"""
155+
e_elect = None
156+
with open(self.path, 'r') as f:
157+
for line in f:
158+
if ':: total energy' in line:
159+
# Format: ":: total energy -9.907945789900 Eh ::"
160+
match = re.search(r'(-?\d+\.\d+)\s+Eh', line)
161+
if match:
162+
e_elect = float(match.group(1))
163+
elif 'TOTAL ENERGY' in line and 'Eh' in line:
164+
# Format: "| TOTAL ENERGY -9.907945789911 Eh |"
165+
match = re.search(r'(-?\d+\.\d+)\s+Eh', line)
166+
if match:
167+
e_elect = float(match.group(1))
168+
if e_elect is None:
169+
raise LogError(f'Unable to find energy in xTB output file {self.path}')
170+
return e_elect * constants.E_h * constants.Na
171+
172+
def load_zero_point_energy(self):
173+
"""
174+
Load the zero-point energy in J/mol from the xTB output.
175+
Only available in frequency (--hess) calculations.
176+
"""
177+
zpe = None
178+
with open(self.path, 'r') as f:
179+
for line in f:
180+
if ':: zero point energy' in line:
181+
match = re.search(r'(-?\d+\.\d+)\s+Eh', line)
182+
if match:
183+
zpe = float(match.group(1))
184+
if zpe is None:
185+
raise LogError(f'Unable to find zero-point energy in xTB output file {self.path}')
186+
return zpe * constants.E_h * constants.Na
187+
188+
def _load_frequencies(self):
189+
"""
190+
Load vibrational frequencies from the xTB output.
191+
192+
xTB prints frequencies twice (after Hessian and in Frequency Printout).
193+
We take only the last complete set.
194+
195+
Returns:
196+
List of positive (real) frequencies in cm^-1.
197+
"""
198+
all_blocks = []
199+
current_block = []
200+
with open(self.path, 'r') as f:
201+
for line in f:
202+
if line.strip().startswith('eigval :'):
203+
values = line.split(':')[1].split()
204+
current_block.extend(float(v) for v in values)
205+
elif current_block and not line.strip().startswith('eigval'):
206+
all_blocks.append(current_block)
207+
current_block = []
208+
if current_block:
209+
all_blocks.append(current_block)
210+
211+
if not all_blocks:
212+
return []
213+
214+
# Use the last frequency block
215+
frequencies = all_blocks[-1]
216+
# Filter out translational/rotational modes (near-zero) and imaginary (negative)
217+
return [f for f in frequencies if f > 0.1]
218+
219+
def _load_spin_multiplicity(self):
220+
"""
221+
Determine spin multiplicity from the xTB output.
222+
223+
xTB reports 'spin' as S (not 2S+1), so multiplicity = 2*S + 1.
224+
"""
225+
with open(self.path, 'r') as f:
226+
for line in f:
227+
if 'spin' in line and 'spin' == line.split()[0].strip():
228+
# Format: "spin : 0.5"
229+
s = float(line.split()[-1])
230+
return int(2 * s + 1)
231+
return 1 # default singlet
232+
233+
def load_conformer(self, symmetry=None, spin_multiplicity=0, optical_isomers=None, label=''):
234+
"""
235+
Load the molecular degree of freedom data from an xTB output file.
236+
237+
Returns:
238+
Tuple of (Conformer, unscaled_frequencies).
239+
"""
240+
if optical_isomers is None or symmetry is None:
241+
_optical_isomers, _symmetry, _ = self.get_symmetry_properties()
242+
if optical_isomers is None:
243+
optical_isomers = _optical_isomers
244+
if symmetry is None:
245+
symmetry = _symmetry
246+
247+
if spin_multiplicity == 0:
248+
spin_multiplicity = self._load_spin_multiplicity()
249+
250+
# Load frequencies
251+
unscaled_frequencies = self._load_frequencies()
252+
modes = []
253+
254+
# Translation
255+
coord, number, mass = self.load_geometry()
256+
translation = IdealGasTranslation(mass=(sum(mass), "amu"))
257+
modes.append(translation)
258+
259+
# Rotation
260+
symbols = [symbol_by_number[i] for i in number]
261+
inertia = get_principal_moments_of_inertia(coord, numbers=number, symbols=symbols)
262+
inertia = list(inertia[0])
263+
if inertia and not all(i == 0.0 for i in inertia):
264+
if any(i == 0.0 for i in inertia):
265+
inertia.remove(0.0)
266+
modes.append(LinearRotor(inertia=(inertia[0], "amu*angstrom^2"), symmetry=symmetry))
267+
else:
268+
modes.append(NonlinearRotor(inertia=(inertia, "amu*angstrom^2"), symmetry=symmetry))
269+
270+
# Vibration
271+
if unscaled_frequencies:
272+
modes.append(HarmonicOscillator(frequencies=(unscaled_frequencies, "cm^-1")))
273+
274+
return Conformer(E0=(0.0, "kJ/mol"),
275+
modes=modes,
276+
spin_multiplicity=spin_multiplicity,
277+
optical_isomers=optical_isomers), unscaled_frequencies
278+
279+
def load_negative_frequency(self):
280+
"""
281+
Load the imaginary (negative) frequency in cm^-1 for a transition state.
282+
"""
283+
# Reuse _load_frequencies logic to get the last block, but include negatives
284+
all_blocks = []
285+
current_block = []
286+
with open(self.path, 'r') as f:
287+
for line in f:
288+
if line.strip().startswith('eigval :'):
289+
values = line.split(':')[1].split()
290+
current_block.extend(float(v) for v in values)
291+
elif current_block and not line.strip().startswith('eigval'):
292+
all_blocks.append(current_block)
293+
current_block = []
294+
if current_block:
295+
all_blocks.append(current_block)
296+
frequencies = all_blocks[-1] if all_blocks else []
297+
neg_freqs = [f for f in frequencies if f < -0.1]
298+
if not neg_freqs:
299+
raise LogError(f'No imaginary frequencies found in {self.path}')
300+
return neg_freqs[0]
301+
302+
def load_force_constant_matrix(self):
303+
"""xTB writes a separate hessian file; not parsed from the main output."""
304+
return None
305+
306+
def load_scan_energies(self):
307+
raise NotImplementedError('Rotor scans are not supported by the xTB adapter.')
308+
309+
def load_scan_pivot_atoms(self):
310+
raise NotImplementedError('Rotor scans are not supported by the xTB adapter.')
311+
312+
def load_scan_frozen_atoms(self):
313+
raise NotImplementedError('Rotor scans are not supported by the xTB adapter.')
314+
315+
316+
register_ess_adapter('XTBLog', XTBLog)

0 commit comments

Comments
 (0)