GFN2-xTB and CREST: expanded calculation scope and caffeine ensemble validation

Abstract

The GFN2-xTB implementation now supports additional chemical environments and property outputs. The caffeine study tests whether independently generated CREST archives recover the same structures after refinement. Every final rotamer frame in three retained archives undergoes the same geometry refinement, followed by bidirectional nearest-neighbor comparisons of all-atom RMSD and relative energy. The property calculations use Mayer bond indices, spin populations and orbitals from a common self-consistent state, together with electronic-energy-weighted conformer properties and harmonic IR spectra.

The three caffeine search archives contain 224 final rotamer frames, producing 448 converged native and WASM refinement paths. All nine bidirectional set comparisons meet the prescribed criteria: maximum all-atom RMSD is 0.0835379023 Å and maximum relative-energy difference is 0.0198879192 kJ mol⁻¹, below 0.15 Å and 0.25 kJ mol⁻¹, respectively. A complete set denotes every frame in each retained archive; the comparisons concern these finite search samples.

Refinement of the complete rotamer ensembles recovers geometric coverage omitted by representative sets. Differences between representative and complete ensembles show that cross-implementation comparisons depend on ensemble selection and matching as well as electronic-structure numerics.

Coverage differences in caffeine representative sets

The principal caffeine discrepancy concerns the hydrogen orientation of one methyl group. The heavy-atom scaffolds nearly coincide, while alignment including every hydrogen gives a substantially larger RMSD. Figure 1 uses the old native rank-5 strict endpoint and overlays its nearest neighbors from the old serial representative set and the refined complete rotamer set. The latter restores coverage of this methyl orientation.

The original study comprises nine searches for butane, hexane and caffeine. All 100,717 actual electronic calls from the six WASM searches were replayed natively at the recorded coordinates and in their original order. After agreement of the electronic evaluations was established, the original ensemble comparisons passed for butane and retained coverage differences for hexane and caffeine.

Strict refinement first used 89 exported representatives across the three molecules, producing 178 paths. Six caffeine paths did not reach the force criterion within 300 steps. A separately versioned 1,000-step budget brought all 178 paths to convergence, but the 27 set comparisons still contained 15 passes and 12 failures: six each for hexane and caffeine. The larger budget resolved convergence, leaving ensemble coverage as the remaining problem.

Inspection of the retained archives found geometries close to this old native endpoint in the final rotamers of all three sources. The exported representative sets did not retain their coverage. The H21/H22/H23 methyl torsion also follows a soft vibrational direction in the retained Hessians. These observations motivated refinement of the complete final rotamer archives.

Refinement of all final rotamers

The new caffeine workflow reads the final archives of the native, plain-serial WASM and JSPI-candidate searches, containing 81, 75 and 68 frames, respectively. Each source is processed independently in its original order, preserving group membership and exported representatives. The count of 224 includes repeated and nearly repeated orientations. For the old rank-5 endpoint, each refinement kernel finds 13, 13 and 9 source frames satisfying both geometry and relative-energy gates.

The original searches use the neutral closed-shell gas-phase GFN2 quick configuration of CREST 3.0.2. Each final frame is refined with both electronic implementations using a frozen ASE 3.29.0 BFGS profile: 300 K electronic temperature, accuracy 0.01, at most 500 SCF iterations and 1,000 optimization steps, with maximum atomic force below $10^{-4}\,\mathrm{eV}/\text{Å}$. BFGS uses alpha=70 and maxstep=0.2 Å; each evaluation creates a fresh electronic calculator.

All 448 paths converge and pass review of the electronic arrays, BFGS updates and endpoint structural identities. They comprise 32 exact reuses of old v2 caffeine representative paths, 35 reviewed native paths from the earlier execution and 381 fresh paths starting from original frames. Additional coverage comes from frames previously omitted from the representative output, while the old endpoints remain unchanged. The original search uses ANCOPT stopping criteria; the post-processing uses maximum Cartesian atomic force.

Ensemble comparison and independent review

Each source pair is evaluated through a primary comparison, a comparison with both sets refined natively and a comparison with both sets refined using WASM, giving nine comparisons. Every endpoint is checked in both directions. Geometry selects the nearest neighbor before its relative-energy difference is assessed. Each source uses its own minimum energy as zero. The gates remain 0.15 Å all-atom RMSD, 0.25 kJ mol⁻¹ relative-energy difference and 0.1 kJ mol⁻¹ minimum-energy difference.

The independent review recomputes the complete alignment matrices using 216 full-graph automorphisms and SVD. These inputs specify no stereochemistry, so the established comparison permits proper and improper alignments. All nine complete-set comparisons pass. The maximum bidirectional nearest-neighbor RMSD is 0.0835379023 Å, the maximum relative-energy discrepancy is 0.0198879192 kJ mol⁻¹, and the largest absolute minimum-energy difference is $2.586\times10^{-8}\,\mathrm{kJ/mol}$.

The two comparators agree on every nearest-neighbor index and scientific verdict. Four near-zero, non-nearest matrix cells differ by approximately 1.10–2.16 microångström; the associated precheck and full matrices are retained with the source data. The final pass applies to the combination of the original complete final rotamers and the specified strict refinement. Original representative-set and v2 failures remain unchanged.

Extensions to the GFN2 calculation domain

GFN2 extensions are defined by structural family, electronic state, method and environment. The added organic structures contain one Al/Ga/In/Bi trialkyl center or one Pb tetraalkyl center, with either all methyl ligands or exactly two ethyl ligands and the remaining ligands methyl. Ordinary GFN2 provides single-point calculations and geometry optimization within this input scope for neutral closed-shell gas-phase inputs. The heavy-hydride branch includes PbH₄, BiH₃ and the specified monomers and three PbH₄ pairs, using fixed-geometry single points and reference-distance windows. Branch-specific conditions are listed below.

Open-shell inputs use specified total charge and unpaired-electron count: UHF=1 for neutral H, OH, NO and CH₃; UHF=3 for N; and UHF=2 for O₂. Gas-phase calculations support ordinary GFN2 and GFN2 + SP. The solvated salt and cluster combinations comprise NaCl, KCl, water dimer, Na⁺–water, NH₄Cl and water–ammonia under ordinary GFN2. OH spin–solvent coupling uses GFN2 + SP with ALPB water.

The heavy-hydride reference comparison uses the original HEAVY28 geometries and dissociation reactions. Across the three complete pair reactions, the maximum absolute deviation is 0.31376 kcal mol⁻¹, mean absolute deviation is 0.12679 kcal mol⁻¹ and maximum relative deviation is 25.10 %. The branch was selected after inspection of the complete benchmark, and its reference-accuracy evidence concerns those three reactions. Same-method native comparisons assess numerical migration; both evidence types retain their calculation conditions and reference objects.

GFN2 + SP denotes GFN2-xTB with an explicit spin-polarization term. Open-shell occupations in ordinary GFN2 and spin-dependent energies in SP are evaluated within their respective method definitions; the UHF parameter specifies the number of unpaired electrons. ALPB water represents continuum solvation, and the OH bond lengths and populations refer to the separately optimized endpoints in each environment. [6, 7]

Unified structural admission and same-SCF properties

Structural admission uses the same rules during preflight, execution and result restoration, checking indexed graphs, connected fragments, formal charges, file radical declarations and explicit hydrogens. Composition or coordinate screening is recorded separately from declared graphs, with initial identity retained. Single-point and optimization records read energy, gradients, charges, dipoles, signed atomic spin populations and continuous Mayer indices from the same successful SCF. Every actual optimization evaluation retains its corresponding matrix and geometry.

The GFN2 + SP OH optimizations illustrate the correspondence between these properties. Gas-phase and ALPB-water O–H distances are 0.96296 and 0.97565 Å; total Mayer indices are 0.88708 and 0.79604, and magnetization indices are −0.04339 and −0.05264. The O/H spin populations are +1.06584/−0.06584 in gas and +1.05878/−0.05878 in ALPB water. Each pair sums to the declared one unpaired electron. Each plotted set comes from a single SCF at its own final geometry.

Mayer indices are retained as continuous matrices, with complete exports of total and available magnetization channels. The 13×13 AlMe₃ matrix retains all 169 entries, including 12 negative total entries. Atomic spin uses signed Mulliken alpha-minus-beta populations; the ordinary single-Hamiltonian-channel route records this quantity as unavailable when it is not supplied. Declared input bonds and calculated indices remain separate, and missing quantities in old records retain their unrecorded state.

Orbitals, charges and conformer-weighted properties

HOMO, LUMO and selected-orbital isosurfaces are generated from the original MO coefficients and basis functions, retaining both signs of the wavefunction. Charge and dipole analysis provides per-atom charges, dipole magnitude and vector, and the coordinate-origin convention. Existing charges and dipoles can be reused from a valid GFN2 result. These analyses use gas-phase closed-shell structures containing H, B, C, N, O, F, P, S, Cl, Br or I, at charge −2 to +2 subject to an even-electron check.

Charges, dipoles, HOMO and LUMO energies, and gaps are evaluated for each input conformer and averaged with electronic-energy weights:

Here $E$ is molar electronic energy converted to units consistent with $R$; each exported conformer uses $g_i=1$. This is an electronic-energy-weighted approximation for the selected set. Vibrational, rotational and solvation free-energy contributions enter a fuller thermodynamic population model. At 298.15 K, the three CREST butane conformers have weights of 57.95 %, 21.18 % and 20.87 %, and a weighted mean dipole magnitude of 0.01373 D.

Conformer-weighted IR and isotope comparison

Conformer-weighted IR first optimizes each geometry, then obtains harmonic frequencies and intensities from central-difference Hessians and dipole derivatives. Weights are recalculated from the electronic energies of those optimized endpoints. Figure 5 reconstructs the saved calculations for the three actual CREST butane conformers at 298.15 K, using Gaussian lines with a 20 cm⁻¹ full width at half maximum, and compares hydrogen with full deuterium substitution.

Isotope comparison changes the mass-weighted vibrational problem while reusing the Hessian, dipole derivatives and electronic-energy weights. Temperature, width and line-shape adjustments likewise operate on the retained results. Analysis records preserve the harmonic approximation, electronic-energy weighting, selected conformers and coverage. Imaginary modes, missing records and duplicate optimized structures are handled explicitly, and spectrum graphics and CSV exports use the same conditions and values.

Spectral coverage is the sum of the weights of successfully evaluated members of the input conformer set. Here, 100 % means that all three input conformers are included. The H/D comparison changes atomic masses within the Born–Oppenheimer harmonic approximation, using a common electronic potential-energy surface and the original electronic-energy weights. Isotope-dependent zero-point and thermal free-energy changes can be used for subsequent population corrections.

Structure input, persistence and result reuse

Names, SMILES, drawings and existing three-dimensional geometries share the molecular canvas. Names resolve to definite structures, while textual and drawn inputs prepare coordinates before calculation. Existing geometries reuse their coordinates directly. Optimized geometries can become inputs to subsequent calculations, and trajectory results transfer the actual selected frame. Structure, method and electronic state remain associated with the same object.

Drawings retain a matched KET/MOL/SMILES version. Calculation records preserve atom order, coordinates, method, electronic state, complete results and provenance, with bindings checked again on restoration. Saving, refreshing, reimporting and returning to results have been verified. The shared file controls export complete records, Mayer matrices, structures and spectra. Restored orbital, charge and conformer results remain available for comparison and property analysis.

Data and reproducibility

The public data support reproduction of the nine caffeine set comparisons and the reported property figures. raw-results.json.gzip retains 448 endpoints, complete alignment matrices, bidirectional nearest neighbors and source counts; caffeine-endpoints.csv and summary.csv provide endpoint and comparison summaries. Verification uses the recorded atom mappings, 216 automorphisms, reflection-permitting SVD alignment and each set’s own energy zero to check nearest-neighbor RMSD, relative-energy differences and minimum-energy differences.

The OH, AlMe₃ and butane exports retain the corresponding structures, methods and numerical results, with H/D modes and weighted spectra supplied for IR analysis. Properties and spectra can be reconstructed using the recorded temperature, line shape, full width at half maximum and electronic-energy weights. Repeating searches and refinements uses the versions, original frames, random streams, SCF settings and BFGS conditions specified in Methods; figure reconstruction and electronic-structure recalculation are recorded separately.