Modern open-source tools for simulation of NMR spectra

Question

I am interested in predicting the NMR spectrum of small organic compounds. It doesn't matter to me if the prediction is very accurate. I'll eventually be comparing the prediction to experimental results.

It looks as if several free online tools offer the ability to do such predictions:

1. NMR DB is a free website that offers a prediction of 1H, 13C, and supposedly even 2D NMR experiments. Structures are entered via a GUI editor.

2. VeSPA and PyGamma are Python packages (or Python wrappers around C libraries) that supposedly simulate NMR spectra of molecules. However, both packages are very old and have not been substantially updated in some time.

My question is, what is the best available (free) computational tool for predicting NMR spectra? As I mentioned, I understand that predictions are imperfect. What software, if any, are widely used by practicing chemists for NMR simulation? What is the best reference to read about such software, including the heuristics and algorithms it uses to simulate spectra? Is Vespa and/or PyGamma state of the art? The NMR DB site is nice, but I'd very much prefer an interface that is amenable to scripting, and where structures could be supplied computationally rather than via a GUI interface.

(I'm a mass spectrometrist and biochemical engineer with little training in NMR methods.)

Answer 1

In Org Biomol Chem 2016, 14, 3943, Goodman reports the replacement of some "tradition" (expensive) programs with free/open-source ones, as applied to their dp4 method.

Specifically:

• Molecular mechanics with tinker (I believe as a replacement for MacroModel), you may or may not want this step depending on what you're doing. The paper is looking at molecules with stereocentres and as such, they appear to have wanted to ensure the NMR values calculated matched the most likely conformation of the molecule experimentally studied.
• NMR prediction using NWChem (replacing Jaguar as the gold standard) to do the DFT calculations. There is some mention that this was slower, but again, it would depend on what you're interested in studying.

I think that fulfils your requirement for being able to do the prediction using open-source programs.

As for scripting, there is a mention in the paper of wrapping the whole thing up in a python script (and to my knowledge, there isn't actually a GUI for either piece of the aforementioned programs), so it should be amenable to this.

One thing that it doesn't appear to do, is to graphically plot the resulting data generated (I believe they just compare a list of values to a list of experimental values), but you should be able to plot this in something like MatLab, or, if you play enough, it should be importable into TopSpin via something to turn it into a set of suitably written values.

Answer 2

You can submit the structure to the NMR database website using this link.

Concerning the algorithms, there are references on the nmrdb.org website. For $\ce{^1H}$ prediction we are using Spinus (based on the neural network). For $\ce{^13C}$ we are using NMR shift DB that is based on hose code.

To get an overview of the prediction algorithms you may be interested in our last publication.

Answer 3

The problem with computational chemistry software (ORCA, GAMESS etc.), as described in NotEvans' answer and in various comments, is that they only spit out a list of chemical shifts and coupling constants. For example, this is part of a typical ORCA output:

 --------------
 Nucleus  34H :
 --------------

[...]

Total shielding tensor (ppm): 
            28.820         -0.203         1.178
            -2.079         29.083         1.927
             1.392          3.320        24.895


 Diagonalized sT*s matrix:

 sDSO            23.705           26.255           39.984  iso=      29.981
 sPSO            -0.555            2.878           -9.468  iso=      -2.382
        ---------------  ---------------  ---------------
 Total           23.149           29.133           30.516  iso=      27.599

Getting from a list of numbers to something that looks like a spectrum is an entirely different matter altogether. In order to do so, one would have to use these parameters to construct a spin Hamiltonian:

$$\hat{H}_{\text{free}} = \sum_i \Omega_i\hat{I}_{iz} + \sum_{ij} 2\pi J_{ij}\hat{I}_{iz}\hat{I}_{jz}$$

and then perform a quantum mechanical simulation, where a density matrix $\rho$ evolves under this free Hamiltonian, as well as any other Hamiltonians corresponding to pulses, according to the Liouville–von Neumann equation. After the pulse sequence, one can obtain the FID by measuring the expectation value of the operator $\sum_i\hat{I}_{i-}$ at regular intervals corresponding to the NMR dwell time, and then the spectrum by Fourier transforming that.

The best general tool for this is Spinach, which is coded in MATLAB. Spinach itself is open-source, but MATLAB isn't (and the Spinach code does not run on GNU Octave). So, it does not quite meet the criteria in the question, but I thought I should mention it anyway.

---

Specifically for use with ORCA, Stefan Grimme's group have done a lot of work in this area, see Angew. Chem. Int. Ed. 2017, 56 (46), 14763–14769. There are various Python scripts that are capable of doing this spin Hamiltonian simulation as well as plotting the spectrum. Indeed, it should now be possible to interface with all these directly from ORCA, such that with one ORCA input file you can perform all the steps described in the paper (conformer generation, NMR property calculation, spin Hamiltonian simulation, and plotting). I suggest looking in the ORCA manual for more information.

ORCA is not open-source per se, but it is free to use, which is what I think you're broadly aiming for (I'd rather not get into debates about what "free" software means).

---

For what it's worth, if you have the programming knowledge required, it's easy enough to write a (e.g.) Python script that is capable of simulating two- or three-spin (maybe four-spin) systems. At the end of the day, an NMR experiment is "just" a bunch of unitary operators acting on density matrices, i.e. lots and lots of matrix multiplication. For example, you could take a look at some of my MATLAB code here which simulates the first FID of a sensitivity-enhanced HSQC experiment. It is fairly trivial to port this to Python, since numpy provides you with pretty much every function you might need.

[Disclaimer: the code will not necessarily stay there forever.]

The only problem is that this scales exponentially with the number of spins, so unless you do some serious optimisation of the code (like that done in Spinach), it becomes intractable very quickly. (It depends on what kind of spectrum you want, too: a simple 1D proton experiment is easy, but if you want a 2D then you need to simulate a number of FIDs, and the complexity increases rapidly if you are using shaped pulses and especially if there is any kind of spatially-dependent element, like pulsed field gradients.)

But if your molecule can be "decomposed" into several spin systems which are small, then the overall spectrum is just the sum of the spectrum of each spin system, so you might be able to get away with simpler code. I've never thought about how you might do that, although my guess is that you would want to construct a Hamiltonian matrix, then pass it to a function which would call itself recursively if the matrix can be written in block-diagonal form.