QCxMS2 Publication Reproducibility Issue

[hq9449]

We are validating QCxMS2 against the test set in Gorges, Engeser, Grimme (2025),
Evaluation of the QCxMS2 method for the calculation of collision-induced
dissociation spectra. Most cases reproduce well, but leucine ([M+H]+, CID)
comes out inverted relative to the paper: the -NH3 loss (m/z 116) is our base
peak and the diagnostic immonium ion (m/z 86) is essentially absent — the paper
reports ESS 0.934 with the immonium as a dominant peak.

We have ruled out every setting we can access and now suspect a CREST /
MSREACT version difference (paper: CREST 1.1.3; ours: CREST 3.0.2). We would
be grateful for your guidance.

Environment

Component	Ours	Paper
QCxMS2	1.2.0 (4db5131), compiled 2025-07-10	(please confirm)
CREST	3.0.2 (2025-08-05)	1.1.3
xtb	6.7.1 (edcfbbe), 2026-06-10	(please confirm)
ORCA	6.1.0	6.0.0
OS	WSL2 (Ubuntu) on Windows; qcxms2 run in native Linux dir	—

Command (composite level = paper's leucine method)

qcxms2 input.xyz -cid -chrg 1 -esiatom 0.250 -pthr 5 \
       -geolevel gfn2 -tslevel wb97x3c \
       -iplevel gfn2 -ip2level wb97x3c -mthr 10 -T 20

Observed vs expected (leucine [M+H]+, m/z 132.10)

m/z	assignment	our intensity (base=100)	paper
116.07	-NH3	100 (base)	minor
132.10	[M+H]+ precursor	~68	precursor
86.10	immonium (diagnostic)	~2 (near 0)	dominant

So -NH3 loss (116) dominates and the immonium (86) is essentially absent — the
inverse of the published spectrum.

What we have already ruled out

• Protomer choice: forced protomer index 1 (O-protonated); also tried the
  automatic Boltzmann ranking (r2SCAN-3c Opt Freq + mRRHO 548 K + SMD water).
• -pthr: used 5 (your paper value). At the qcxms2 default of 1% the immonium
  is further fragmented (to ~m/z 43), so 5% is clearly correct.
• esiatom: 0.25 (also tried higher).
• -msnoiso: enabling it made things WORSE — it removed the immonium (m/z 86
  -> 0) while leaving -NH3 (116) as base peak. (The immonium channel is itself
  rearrangement-mediated, so -msnoiso suppresses the peak we want.)
• Full-DFT geometries (extra check, beyond the paper's leucine level):
  -geolevel wb97x3c -tslevel wb97x3c did NOT complete within a 72-hour wall
  limit for this 22-atom molecule (still in the NEB transition-state phase at
  72 h). So we cannot resolve this by raising the geometry level of theory; it
  points upstream to MSREACT fragment generation.

We could NOT test the paper's CREST version — and here is why

• CREST 1.1.3 (the paper's version) is not on conda-forge, which packages
  only 2.11.1, 2.11.2, 2.12, and 3.0.2.
• We installed the closest available (CREST 2.12 + matched xtb 6.4.0) in
  an isolated env and invoked QCxMS2 1.2.0 with it on the same input. QCxMS2
  aborted at startup, before any fragmentation:

  forrtl: severe (124): Invalid command supplied to EXECUTE_COMMAND_LINE
  Image  PC        Routine            Line   Source
  qcxms2 ...       utility_mp_check_  1260   utility.f90
  qcxms2 ...       argparser_mp_chec  630    argparser.f90
  qcxms2 ...       MAIN__             91     main.f90
  

• On investigation, CREST 2.12 has no MSREACT capability at all — it does not
  recognize --msreact / --msfulliso / --msnoiso, and its --help lists no
  reaction/fragmentation options. MSREACT post-dates the conda-forge 2.x line.
• Therefore the only conda-forge CREST that can drive QCxMS2 1.2.0 is 3.x, and
  we cannot reproduce the paper's CREST 1.1.3 + MSREACT toolchain locally.

Questions

1. Which exact CREST and QCxMS2 versions produced the leucine result in the
   paper (ESS 0.934)? Is our discrepancy consistent with an MSREACT change
   between CREST 1.1.3 and 3.0.2?
2. Could you point us to the CREST 1.1.3 build used for the paper (it is not
   on conda-forge), and confirm which QCxMS2 version it pairs with? We would like
   to reproduce the published toolchain.
3. Are there MSREACT settings (reaction-type restrictions, number of
   fragmentation attempts, isomerization handling) used for the paper that are
   not current defaults?
4. Is -NH3 dominance over the immonium a known artifact for amino acids at
   this version, and is there a mitigation that does NOT also suppress the
   immonium (i.e. not -msnoiso)?
5. Given full-DFT geometries are impractical here, do you consider the issue to
   be upstream in fragment generation (MSREACT) rather than in the barriers?

Files we can provide

• input.xyz (protonated protomer-1 geometry)
• qcxms2.log (composite run)
• allpeaks.dat / peaks.csv
• qcxms2_crest212.log (the severe-124 startup abort with CREST 2.12)

Happy to attach any of these or run additional diagnostics. Thank you very much
for QCxMS2 and for any guidance.

Verbatim external-program commands used (from our run logs)

# Protomer generation:
crest seed.xyz -protonate -ewin 60 -gfn2 -T 20

# Protomer dedup:
molbar protomer_<k>.xyz -c 1 -m mb

# Protomer Gibbs ranking (ORCA):
! r2SCAN-3c Opt Freq
%cpcm smd true smdsolvent "water" end
%freq Temp 548 end

# Conformer search (verbatim):
crest best_protomer.xyz -gfn2 -T 20 -chrg 1

# Fragment NEB-TS (ORCA, verbatim from a fragment dir):
! NEB
! UKS wB97X-3c
! LOOSESCF
%maxcore 8000
%pal nprocs 1 end
%NEB NEB_END_XYZFILE "end.xyz"
     NImages 7
     MaxIter 500
end
*xyzfile 1 1 start.xyz

# Main fragmentation (composite):
qcxms2 input.xyz -cid -chrg 1 -esiatom 0.250 -pthr 5 \
       -geolevel gfn2 -tslevel wb97x3c -iplevel gfn2 -ip2level wb97x3c \
       -mthr 10 -T 20

Any insight is much appreciated.

Thanks,
Michael Mawn

[gorges97]

Hi,

Thank you for the detailed issue report. Here are a few clarifications regarding the setup and toolchain used in our publication (Gorges, Engeser, Grimme, 2025):

1. Software Versions & xtb Prerequisites:

   • CREST: We used CREST version 3.0.2 (the version 1.1.3 you referenced belongs to Molbar, not CREST).
   • QCxMS2 & xtb: We can confirm that QCxMS2 1.2.0 was used. For xtb, we used a development version, as two specific PRs are required for the workflow:
     • xtb PR #1004: Implements the -bhess option in the xtb thermo submodule.
     • tblite PR #1327: Enables the required tblite integration for GFN2-xTB.
       However, a new release version of xtb, that includes these changes is coming soon.
   • MSREACT: Default settings were used throughout. The discrepancy is therefore not rooted in MSREACT settings. You could check this, by searching in the msreact.out output files for the fragment with mass of 86.

2. GFN2-xTB Implementation (tblite vs. xtb):
   For the CID spectra calculations in the paper, we used the tblite implementation of GFN2-xTB via the -geolevel gfn2_tblite and -iplevel gfn2_tblite flags. Unlike the standard xtb implementation, tblite provides a full spin-unrestricted SCF treatment with two independent spin channels, which yields a better physical description of two-electron bond cleavage during fragmentation.

   The exact command to reproduce the paper's leucine results is:

   qcxms2 input.xyz -cid -chrg 1 -esiatom 0.250 -pthr 5 \
          -geolevel gfn2_tblite -tslevel wb97x3c \
          -iplevel gfn2_tblite -ip2level wb97x3c -mthr 10 -T 20

3. ORCA Version Compatibility & Merged Fix:
   The calculation works as intended with ORCA versions like ORCA 6.0.0 (which was used for the paper). However, newer ORCA 6.1.x versions suppress the orca.g98.out file that QCxMS2 historically required for Hessian parsing. We have just merged Pull Request Add the Option to read orca.hess files #24 into main, which resolves this issue for ORCA 6.1.x by adding a fallback to read orca.hess directly. Additionally, I observed on some system architectures that ORCA 6.1.1 returns unexpected exit codes if an xtb calculation does not converge during a NEB run, which can lead to execution crashes in QCxMS2.

---

Next Steps & Diagnostics

I re-ran the calculations locally with QCxMS2 1.2.0 (including the merged PR #24), xtb (with PR #1004 and tblite PR #1327), and ORCA 6.0.1. I still obtain the diagnostic peak at m/z 86 with both -geolevel gfn2 and -geolevel gfn2_tblite, suggesting a technical issue or silent failure might be playing a role in your setup.

To help troubleshoot:

• Could you please attach your qcxms2.log and output files so I can inspect if everything completed without silent technical failures?
• Please update your QCxMS2 build to the latest main branch (which includes the fix for ORCA 6.1.x Hessian parsing) and retry the calculation.

Best regards,
Johannes Gorges