Cluster solver

Introduction

Similar to the direct solver, the cluster solver calculates transitions in a 1D NMR spectrum through the resolvent approach. However, it can calculate the spectral function for much larger molecules, as it automatically divides the molecule into small clusters that are only weakly coupled and solves them as being independent. The solver can either be specified as "cluster_nmr_solver" or "cluster_nmr_solver (symmetry)". If the symmetry version is chosen, the solver is allowed to exploit local SU(2) symmetry. While this does allow to go to larger system sizes, it is not necessarily faster for small molecules, due to additional computational overhead in the evaluation.

Although this is the fastest solver for large molecules, it can give inaccurate results for spin contributions of spins at the boundary between two clusters. To improve accuracy, check out the spin dependent cluster solver.

The setup is the same as for the direct solver. Meaning based on a struqture Hamiltonian input, the spectral function is computed from the correlation function ${\hat{M}}^{-}{\hat{M}}^{+}$, where $\hat{M}$ corresponds to the spin operator including the gyromagnetic factors gamma, ${\hat{M}}^{\alpha }={\sum }_{\ell }{\gamma }_{\ell }{\hat{I}}_{\ell }^{\alpha }$, ${\hat{I}}^{\alpha }={\hat{S}}^{\alpha }/\hslash$.

Math

As for the direct solver the resolvent representation (see chapter math) is used. The correlation function is given as

$$\begin{array}{rl}{\mathcal{G}}_{\eta }(\omega )& =\sum _n{\mathrm{e}}^{-\beta E_n}⟨n|{\hat{M}}^{-}\frac{1}{\omega -\left(\hat{H}-E_n\right)+i\eta }{\hat{M}}^{+}|n⟩/Z\\ & \\ {\mathcal{G}}_{\eta }(\ell ,\omega )& =\sum _n{\mathrm{e}}^{-\beta E_n}{\gamma }_{\ell }⟨n|{\hat{I}}_{\ell }^{-}\frac{1}{\omega -\left(\hat{H}-E_n\right)+i\eta }{\hat{M}}^{+}|n⟩/Z\end{array}$$

However, it is only evaluated on each individual cluster. The clusters are specified in multiple steps:

Partitioning

First, a partitioning of the system into strictly independent parts (from now on called partitions) is performed. This is achieved by diagonalizing the Laplacian matrix based on the J-coupling matrix of the Hamiltonian. By counting the number of zero eigenvalues, the number of partitions in the system can be easily determined. To identify them explicitly, one can use an analysis of the Fiedler vector, which is normally defined as the eigenvector corresponding to the second-lowest eigenvalue. In an ideal case a partition is then given by the zero entries in the vector leading to a bipartitioning of the system. By recursively applying this method on each of the two identified partitions, one can split the system into all of its individual parts. However, special care has to be taken as there can be no zeros in the Fiedler vector although the number of zero eigenvalues is larger than one. This can happen as the Fiedler vector is actually not well-defined if the zero eigenvalues, which are also the lowest eigenvalues, are degenerate. In these cases a linear combination of two eigenvectors with zero eigenvalues can be created such that at least one entry of the resulting vector is zero. This vector can then be used as the Fiedler vector instead. After performing the partitioning, each partition is treated independently, however since these are truly non interacting parts of the system, no approximation was made here.

Symmetry Considerations

As a next stage, if the symmetry option was specified, symmetric groups are identified in each partition. A symmetric group is defined as a group of spins, where each individual spin has the same chemical shift and couples in the exact same way to the rest of the system. Identifying these groups is advantageous, as one can combine them into higher order spin representations. Exploiting the SU(2) symmetry then leads to a significant reduction in computational space. As an example consider propane, which has eight hydrogens. By identifying all symmetrically coupled groups in this molecule, the number of spins can be reduced to two. One being the CH2 group which has a combined spin representation of one and the second being the two methyl groups each representing a spin 3/2 and adding up to a total spin representation of spin 3.

Clustering

Once all of the groups in a partition have been identified, a final clustering of the partition into its weakly coupling parts can be performed. Note that up to here, everything was still exact and approximations are only introduced through clustering. To perform the clustering lower triangular matrices are created, representing each partition as a weighted graph with the spins being the nodes and the edge weights being ${\Delta }_{x,y}$ based on assumptions from perturbation theory:

$${\Delta }_{x,y}=\frac{J_{xy}^2}{({\gamma }_x{\delta }_x-{\gamma }_y{\delta }_y)B}$$

Where $J_{x,y}$ are the entries in the J coupling matrix connecting sites x and y, $\gamma$ the gyromagnetic ratio, $\delta$ the chemical shift and $B$ the magnetic field strength. After defining this graph, a Stör-Wagner splitting of the graph is used to obtain individual clusters. Note that the groups which were previously identified in the partition are always preserved by this splitting as they are treated as a single spin.

Extra Arguments

eta: float Artificial broadening.

The broadening $\eta$ is formally necessary to ensure the convergence of the Fourier integral. In practice it corresponds to the unknown or neglected noise, being intrinsic to the resolution limitation of the spectrometer. Note that the complete spectral function is positive, $A(\omega )\ge 0$. However, the spin-resolved spectral functions can be negative.

beta: float Inverse temperature

$\beta =1/(k_{\mathrm{B}}T)$.

Due to the low energy of nuclear spin excitations the temperature in real measurements is practically infinite, i.e. $\beta =0$. For theoretical studies or comparison to very low temperature NMR measurements one can apply finite temperatures.

Note: Currently, only the case $\beta =0$ is implemented.

max_cluster_size: int Maximum cluster size.

This has to be set according to restraints imposed by the computational resources available. Note that symmetric groups are viewed as one site, with an effective number of sites corresponding to the spin representation of the group.

tolerance_couplings: float Tolerance for the J-coupling when identifying symmetry groups

When the symmetry groups are identified, the coupling to the outside does not have to be perfectly symmetric, as there might be small deviations in the input parameters. Hence this parameter allows to define a tolerance for the J-coupling until which spins are assumed to couple symmetric.

tolerance_shifts: float Tolerance for the chemical shifts when identifying symmetry groups

Similarly to the tolerance for the J-couplings the chemical shifts don't have to be perfectly identical. This parameter allows in the same way to define a tolerance until which they are assumed to be the same.

delta: float

Small parameter to avoid division by zero, when creating the weight matrix for the clustering.

threshold_matrix_elements: float

In the final evaluation a lot of matrix elements $⟨n|I_l^{+}|m⟩$ evaluate to zero and hence do not need to be included, when performing the summation for each frequency. We therefore remove all matrix elements below this threshold beforehand.

verbose: int Verbosity of output

Specifically when set to one, the partitioning, as well as the identified groups and the clustering is printed. For further debug information it may be set higher.