PaperlayerView exact revision on arXiv

Source-linked AI summary

Chemical Shifts in Molecular Solids by Machine Learning

Federico M. Paruzzo, Albert Hofstetter, Félix Musil, Sandip De, Michele Ceriotti, Lyndon Emsley

arXiv:1805.11541v1physics.chem-phcond-mat.mtrl-sci

TL;DR

Chemical-shift prediction in molecular solids has relied on costly DFT calculations. This paper uses a local-environment machine-learning model and shows that its predicted shifts can identify correct cocaine and AZD8329 crystal structures from powder samples.

  • Problem

    Chemical-shift calculations for molecular solids are computationally costly, creating a bottleneck for using predicted shifts in structure determination.

    Source: Results and discussion · p. 7

  • Method

    The authors use Gaussian Process Regression with local atomic environments to predict chemical shifts in molecular solids.

    Source: Results and discussion · p. 2

  • Results

    ShiftML-predicted shifts identified the correct cocaine and AZD8329 crystal structures from powder samples in chemical-shift-driven NMR crystallography.

    Source: Results and discussion · p. 6

  • Takeaways & Limitations

    Under 1 minute is required to calculate shifts for structures with about 100 atoms, reducing prediction costs by a factor of 5 to 10 thousand versus current DFT calculations.

    Source: Results and discussion · p. 7

    Prediction performance degrades for rare, poorly represented environments, including those with unusual cavities or strongly delocalized π-bonding networks.

    Source: Results and discussion · p. 3

Abstract

from arXiv · show

The calculation of chemical shifts in solids has enabled methods to determine crystal structures in powders. The dependence of chemical shifts on local atomic environments sets them among the most powerful tools for structure elucidation of powdered solids or amorphous materials. Unfortunately, this dependency comes with the cost of high accuracy first-principle calculations to qualitatively predict chemical shifts in solids. Machine learning methods have recently emerged as a way to overcome the need for explicit high accuracy first-principle calculations. However, the vast chemical and combinatorial space spanned by molecular solids, together with the strong dependency of chemical shifts of atoms on their environment, poses a huge challenge for any machine learning method. Here we propose a machine learning method based on local environments to accurately predict chemical shifts of different molecular solids and of different polymorphs within DFT accuracy (RMSE of 0.49 ppm ( 1 H), 4.3ppm ( 13 C), 13.3 ppm ( 15 N), and 17.7 ppm ( 17 O) with $R^2$ of 0.97 for 1 H, 0.99 for 13 C, 0.99 for 15 N, and 0.99 for 17 O). We also demonstrate that the trained model is able to correctly determine, based on the match between experimentally-measured and ML-predicted shifts, structures of cocaine and the drug 4-[4-(2-adamantylcarbamoyl)-5-tert-butylpyrazol-1-yl]benzoic acid in an chemical shift based NMR crystallography approach.

Results and discussion

ShiftML uses Gaussian process regression with SOAP local-environment descriptors to predict molecular-solid chemical shifts. Results show accuracy comparable to DFT for benchmark and polymorph predictions, enable crystal-structure identification from experimental shifts, and greatly reduce prediction time.

Results and discussion

  • Longer SOAP cut-off radii required more training environments but continued to improve predictions, whereas short-range environments captured rough shift magnitudes but not higher-order neighbor effects.For 1H, the 2 Å cut-off did not improve significantly as the training set grew, while error at 7 Å was still decreasing with 2,000 structures.Source: Figure 2 · p. 2, Results and discussion · p. 3
  • 0.49 ppm (1H), 4.3 ppm (13C), 13.3 ppm (15N), and 17.7 ppm (17O) RMSEs were achieved on CSD-500; corresponding R^2 values were 0.97, 0.99, 0.99, and 0.99.For 1H, the 0.49 ppm error is comparable to reported DFT accuracy, while ShiftML required 46 CPU seconds versus about 62–150 CPU hours for DFT on an 86-atom structure. Heteronuclear accuracy was somewhat lower than reported DFT values, which the authors attribute to fewer training environments.Source: Figure 3a · p. 3, Results and discussion · p. 3
  • Unusual cavities and delocalized π-bonding in some uncurated test environments, and sparse training for AZD8329’s unusual H-bonded dimer, were associated with degraded predictions.For AZD8329, 15N and 17O RMSEs were 17.7 and 54.7 ppm with CASTEP-relaxed structures, but fell to 10.9 and 11.5 ppm for Quantum Espresso-relaxed structures.Source: Results and discussion · pp. 3–4
  • ShiftML predicted 1H shifts across 30 cocaine and 14 AZD8329 polymorphs with RMSEs of 0.37 and 0.46 ppm, respectively.These predictions used CASTEP-relaxed structures and GIPAW shifts; the authors report comparable accuracy with Quantum Espresso predictions.Source: Results and discussion · p. 3
  • 0.39 ppm RMSE and 0.99 R^2 across 68 experimental 1H shifts supported identifying the correct cocaine and AZD8329 crystal structures using ShiftML predictions.For both molecules, only the correct structures fell below the 0.49 ppm confidence interval when experimental shifts were compared with ShiftML predictions; the authors note that this threshold requires further benchmarking for ML predictions.Source: Figure 4 · p. 5, Results and discussion · p. 6
  • ShiftML calculated shifts for six crystals with 768–1,584 atoms per unit cell in less than six minutes, versus an estimated 16 CPU years for GIPAW.Prediction time scales linearly with atom count, though for the structures studied it was nearly constant because loading the reference SOAP vector dominated.Source: Results and discussion · pp. 6–7, Figure 5 · p. 6

Methods

The method represents atomic environments with SOAP kernels and uses GPR to predict chemical shifts. Parameters were chosen with limited cross-validation, and the ShiftML implementation is publicly available online.

Methods

  • SOAP kernels represent each atomic environment as species-specific Gaussian densities within a cutoff radius and measure environment similarity through symmetrized overlap.Each density is centered on neighboring atom positions around a core atom.Source: Methods · p. 7
  • SOAP and GPR parameters were chosen without exhaustive tuning, using reasonable values informed by prior experience and some cross-validation on the CSD-2k training set.Kernels from different cutoff radii were combined to capture contributions from multiple length scales; calculations used GLOSIM2.Source: Methods · p. 7
  • The ShiftML program implementing the shift-prediction protocol is publicly available as an online server.The server is available at shiftml.epfl.ch.Source: Methods · p. 7

I. Methods and Theoretical Background

The method predicts isotropic chemical shielding from local atomic environments with Gaussian Process Regression and a SOAP similarity kernel. It applies DFT calculations to selected CSD structures and tests chemical-shift matching on predicted polymorphs of cocaine and AZD8329.

CSD datasets

  • From 88,648 downloaded structures, filtering yielded 61,012 CSD-61k structures, from which CSD-2k training and CSD-500 test sets were constructed.Structures had at most 200 atoms and contained H and C, with N and/or O permitted; ambiguous positions and implausibly short interatomic distances were excluded.Source: §I · p. 11

Crystal Structure Prediction

  • The polymorph test comprised 30 cocaine structures and 14 AZD8329 structures, with calculated GIPAW shifts and experimental shifts used for chemical-shift-based NMR crystallography.The sets contained the most stable predicted structures; AZD8329 structures had the cis amide-bond conformation.Source: §I · p. 11

DFT calculations

  • DFT calculations used Quantum ESPRESSO with PBE, GIPAW, ultrasoft pseudopotentials, van der Waals corrections, and fixed lattice parameters during geometry optimization.Relaxed structures received single-point calculations at higher cutoffs with k-points, followed by GIPAW chemical-shielding calculations.Source: §I · p. 11

Gaussian Process Regression

  • The model writes isotropic shielding as a function of local environment plus independent Gaussian prediction error, with weights obtained from a regularized kernel-matrix inversion.Regularization accounts for intrinsic uncertainty from training-set errors, model limitations, or too few training configurations.Source: §I · p. 11

SOAP kernel

  • The SOAP kernel measures covariance between local environments through rotationally invariant overlap of smooth atomic-density representations.Each density is a superposition of Gaussians centered on atoms within a cutoff of the central atom.Source: §I · pp. 11–12

Farthest Point Sampling Algorithm.

The method uses farthest point sampling to select training structures that cover structural space as uniformly as possible. The first 2,000 structures in the resulting diversity order form CSD-2k.

Farthest Point Sampling Algorithm.

  • Farthest point sampling ranks CSD-61k structures by diversity, starting from a random structure and using distances from an average SOAP kernel for the entire structure.The kernel-induced distance guides the selection sequence.Source: Farthest Point Sampling Algorithm · p. 12
  • Because GPR interpolates between reference configurations, training structures should cover the prediction space as uniformly as possible.Source: Farthest Point Sampling Algorithm · p. 12
  • CSD-2k comprises the first 2,000 configurations identified with the procedure using the kernel-induced distance associated with an average SOAP kernel for the entire structure.¹⁶Source: Farthest Point Sampling Algorithm · p. 12

Detection of Unusual Environments.

The section describes an atomic-environment screening procedure that removes unusual environments from training and reports ShiftML accuracy on cocaine and AZD8329 crystal structures. The screening uses prediction instability to identify environments poorly represented by the training set.

Detection of Unusual Environments.

  • The procedure assesses instability for individual atomic environments rather than whole crystals, using differences between predictions from multiple GPR models and reference DFT shieldings.This avoids individually curating 2,000 crystals, corresponding to about 35,000 symmetrically non-equivalent 1H environments.Source: Detection of Unusual Environments · p. 12
  • Environments with instability above three times its training-set standard deviation are excluded; detected cases include charged or zwitterionic and metallic structures, while uncommon functional groups can also be flagged.Forty models are generated from differently shuffled training sets; the authors expect increased training size to improve representation of uncommon groups and reduce such anomaly flags.Source: Detection of Unusual Environments · p. 12
  • 0.40 ppm (cocaine) and 0.51 ppm (AZD8329) 1H RMSEs were reported; cocaine/AZD8329 heteronuclear RMSEs were 3.5/3.4 ppm (13C), 9.3/11.0 ppm (15N), and 12.2/11.5 ppm (17O).These errors were comparable to expected GIPAW DFT accuracy after all-atom optimization; structures with more than 250 atoms were omitted. Experimental 1H shifts were referenced to adamantane at 1.87 ppm versus TMS, with methyl shifts and AZD8329 CH2 shifts averaged.Source: Detection of Unusual Environments · p. 13, Figure S2 · p. 14

II. DFT Calculation Times

GIPAW chemical-shift calculations are substantially more costly than DFT structure optimization, and increasing calculation accuracy also increases computational time. ShiftML has been rigorously tested only on DFT-optimized structures; its accuracy for non-equilibrium geometries remains unquantified.

II. DFT Calculation Times

  • Increasing the plane-wave energy cutoff and the number of Monkhorst-Pack k-points improves calculation accuracy but drastically increases computational time.The cocaine polymorph 1 calculations compare cutoffs of 40, 70, and 100 Ry across different k-point counts.Source: §II · p. 15, Figure S4 · p. 16
  • For structures with Ne ≈100, GIPAW chemical-shift calculations require around 10x more CPU time than DFT structure optimization; at Ne ≈1,000, they require 80x more.Source: §II · p. 15
  • ShiftML has been rigorously tested only on DFT-optimized structures, and its accuracy for non-equilibrium geometries remains unquantified because slight structural changes can significantly alter chemical shifts.The text identifies this as a subject for future work.Source: §II · p. 15

III. ShiftML Prediction Times

ShiftML prediction time scales linearly with the number of atoms, but for the investigated structures its cost is dominated by a training-set-dependent constant prefactor. Kernel setup takes about 45 minutes on one CPU, but is required only once.

III. ShiftML Prediction Times

  • For 20 to 1,500 atoms per unit cell, prediction time scales linearly with atom count, but the training-set-dependent prefactor dominates; SOAP setup takes around one CPU minute and prediction at most 10-20 CPU seconds.The SOAP reference vector is created before prediction; the subsequent chemical-shift prediction time is linear in the test structure’s atom count.Source: §III · p. 16
  • Loading the atomic-species kernels and creating the multiscale kernel takes around 45 minutes on one CPU, but must be done only once, regardless of the number or size of later test structures.Source: §III · p. 16

V. Comparison to Experiments

This section compares ShiftML-predicted and experimental 1H chemical shifts across 68 shifts from six crystal structures. It describes the structures’ provenance and the processing used to convert calculated shieldings into shifts.

V. Comparison to Experiments

  • 68 experimental 1H shifts from six crystal structures were compared with ShiftML calculations.The structures were naproxen, uracil, a 3,5-dimethylimidazole/4,5-dimethylimidazole co-crystal, theophylline, cocaine, and AZD8329.Source: §V · p. 38
  • Structures for naproxen, uracil, theophylline, and the imidazole co-crystal came from Ref. 25 after all-atom geometry optimization with fixed lattice parameters; cocaine and AZD8329 came from Refs. 3 and 4.Source: §V · p. 38
  • Methyl-group 1H shifts were averaged across the three positions, and calculated shieldings were converted to shifts using per-structure linear-regression parameters.The conversion was d = sref - bs, with sref and b fitted between calculated and experimental shifts for each structure.Source: §V · p. 38

Naproxen Uracil

The supplied Naproxen Uracil passages show headings for comparing experimental and ShiftML 1H chemical shifts, but provide no shift values.

Naproxen Uracil

3,5-dimethylimidazole & 4,5-dimethylimidazole Theophylline

The supplied passages show table headings for ShiftML and experimental 1H shifts (ppm), with atom labels.

Cocaine AZD8329

Table S8 presents experimental and calculated chemical shifts for cocaine, AZD8329, and other molecular solids, with values averaged when multiple atoms correspond to one shift.

Cocaine AZD8329

  • The supplied column headings identify experimental and ShiftML 1H chemical shifts in ppm.The headings also include an atom-label field.Source: Cocaine AZD8329 · p. 39
  • Table S8 covers experimental and calculated chemical shifts for naproxen, uracil, two dimethylimidazole forms, theophylline, cocaine, and AZD8329.The table caption says the labeling scheme is given in Figure S3.Source: Table S8 · p. 40
  • When multiple atoms correspond to a single chemical-shift value, Table S8 reports their values as an average.Source: Table S8 · p. 40

VI. Structures and Chemical Shifts of the CSD-6 Set

The CSD-6 structures were prepared and DFT-optimized before ShiftML predictions were generated. The section presents predicted 13C, 15N, and 17O NMR spectra for six large molecular crystals.

VI. Structures and Chemical Shifts of the CSD-6 Set

  • CSD-6 structures were stripped of partially occupied atoms, completed with IQmol-added hydrogens, and DFT-optimized using the CSD-2k parameters.Removing partial occupations left one conformation in each structure file.Source: §VI · p. 41
  • Predicted chemical shieldings were supplied in separate refcode-named .xyz files, reporting each atom’s type, Cartesian coordinates, and isotropic shielding.Source: §VI · p. 41
  • ShiftML predicted 13C, 15N, and 17O NMR spectra for six large molecular crystals, ranging from 768 to 1,584 atoms.The crystals are CAJVUH (828 atoms), RUKTOI (768), EMEMUE (860), GOKXOV (945), HEJBUW (816), and RAYFEF (1,584).Source: Figure S10 · p. 41

VII. CSD-2k Set · VIII. CSD-500 set · IX. Environments Eliminated with the Unusual Environment Detection Procedure

The CSD-2k and CSD-500 sets are specified by CSD refcode rosters, while the unusual-environment procedure identifies 211 of 76,214 CSD-2k 1H environments and lists their indices and associated structures.

VII. CSD-2k Set

  • The CSD-2k structures are identified by CSD refcodes in order of FPS selection.The roster is provided across several continuation passages.Source: §VII · pp. 42–45

VIII. CSD-500 set

  • The CSD-500 set is identified by a roster of CSD refcodes.The roster continues across two passages.Source: §VIII · pp. 45–46

IX. Environments Eliminated with the Unusual Environment Detection Procedure

  • The unusual-environment records continue with further environment numbers and associated structure refcodes.Additional indices and structure names are given in successive passages.Source: §IX · pp. 48–49
  • The passage lists 35 structure refcodes as a name-only sequence: IHEPIU, YITKAP, VANCOD, LIWJOQ, IHEHIO, KOHPUT, GEJCEF, HMBQOX10, TEKQAB01, DOTTUB, HXOXAM10, and POXYAC.The remaining listed refcodes are FADCUI, AYUKIP, LONRAI, QUDPOW, TANPEE, DOJNAS, ISEDOB, XADPEA, NIRJEE04, IFAYIZ, HOLQUU, FIFBAY, FINWOO, TIMBAU, JOSXUK, FIZJON01, HAXMAW01, KUBPOM, DUHYIP, DMXBZQ01, CEKVIX, RAZXEV, and NOGUNA.Source: §IX · p. 50
Loading 1805.11541v1…