Diffraction Signals

For CAS-type wavefunctions (CASCI, CASSCF, and FOMO-CASCI), TeraChem can evaluate elastic and inelastic X-ray and electron (UED) scattering intensities, state by state. The implementation follows the formulation worked out by Yang et al.¹ (see the supplementary information for the full derivation); below we summarize what TeraChem actually computes and prints.

Theory

Scattering operator

For an incoming probe with momentum \(\mathbf{k}_0\) and an outgoing probe with momentum \(\mathbf{k}_s\), the scattering vector (momentum transfer) is \(\mathbf{s} = \mathbf{k}_s - \mathbf{k}_0\). Within the first Born approximation, the differential cross section into final electronic state \(|\Psi_i\rangle\) is proportional to

\[ I(\mathbf{s},\epsilon_s,i) \;\propto\; \frac{\epsilon_s}{\epsilon_0}\,\big|\langle\Psi_0|\hat{L}|\Psi_i\rangle\big|^2\, \delta(E_0+\epsilon_0-E_i-\epsilon_s), \]

with the scattering operator depending on the probe:

\[ \hat{L}_{\text{electron}}(\mathbf{s}) \;=\; \frac{1}{s^2}\sum_\alpha N_\alpha\,e^{i\mathbf{s}\!\cdot\!\mathbf{R}_\alpha} \;-\;\frac{1}{s^2}\sum_j e^{i\mathbf{s}\!\cdot\!\mathbf{r}_j}, \qquad \hat{L}_{\text{X-ray}}(\mathbf{s}) \;=\; -\sum_j e^{i\mathbf{s}\!\cdot\!\mathbf{r}_j}. \]

Here \(N_\alpha\) and \(\mathbf{R}_\alpha\) are the charge and position of nucleus \(\alpha\). The X-ray operator is obtained from the electron operator by setting all \(N_\alpha\to 0\) and dropping the \(s^{-2}\) Mott prefactor.

Elastic and inelastic intensities

Because experiments typically do not resolve the final state's energy, the intensity at scattering vector \(\mathbf{s}\) is summed over final states \(i\) and (for inelastic) integrated over \(\epsilon_s\). Using the closure relation \(\sum_i|\Psi_i\rangle\langle\Psi_i|=\hat 1\), the total intensity splits as \(I = I_{\text{elastic}} + I_{\text{inelastic}}\) with

\[ I_{\text{elastic}}(\mathbf{s}) \;=\; \big|\langle\Psi_0|\hat{L}|\Psi_0\rangle\big|^2, \qquad I_{\text{inelastic}}(\mathbf{s}) \;=\; \langle\Psi_0|\hat{L}^\dagger\hat{L}|\Psi_0\rangle \;-\; \big|\langle\Psi_0|\hat{L}|\Psi_0\rangle\big|^2. \]

Introducing the Fourier transforms of the one- and two-electron densities

\[ f(\mathbf{s}) \;\equiv\; \int e^{i\mathbf{s}\!\cdot\!\mathbf{r}}\,\rho(\mathbf{r})\,d\mathbf{r}, \qquad P(\mathbf{s}) \;\equiv\; \int e^{i\mathbf{s}\!\cdot\!(\mathbf{r}-\mathbf{r}')}\,\rho^{(2)}(\mathbf{r},\mathbf{r}')\,d\mathbf{r}\,d\mathbf{r}', \]

and the nuclear form factor \(f_N(\mathbf{s}) \;=\; \sum_\alpha N_\alpha\,e^{i\mathbf{s}\!\cdot\!\mathbf{R}_\alpha}\), the printed quantities take the closed forms

\[ \begin{aligned} I^{\text{X-ray}}_{\text{elastic}}(\mathbf{s}) &= \big|f(\mathbf{s})\big|^2, \\[2pt] I^{\text{UED}}_{\text{elastic}}(\mathbf{s}) &= \big|f_N(\mathbf{s})-f(\mathbf{s})\big|^2, \\[2pt] I_{\text{inelastic}}(\mathbf{s}) &= n \;+\; P(\mathbf{s}) \;-\; \big|f(\mathbf{s})\big|^2, \end{aligned} \]

where \(n\) is the total number of electrons. The full electron-scattering cross section is obtained by multiplying the UED columns by \(s^{-4}\) (the Mott prefactor); TeraChem leaves that factor out of the printout so that the same column can be inspected without an \(s\to 0\) divergence.

Evaluation in a CAS wavefunction

The integrals \(f(\mathbf{s})\) and \(P(\mathbf{s})\) are evaluated directly in the MO basis using

\[ f(\mathbf{s}) = \sum_{pq}\gamma_{pq}\,M_{pq}(\mathbf{s}), \qquad P(\mathbf{s}) = \sum_{pqrs}\Gamma_{pqrs}\,M_{pq}(\mathbf{s})\,M_{rs}(-\mathbf{s}), \]

where \(\gamma_{pq}\) and \(\Gamma_{pqrs}\) are the state's one- and two-particle reduced density matrices and \(M_{pq}(\mathbf{s}) = \int e^{i\mathbf{s}\!\cdot\!\mathbf{r}}\,\phi_p^*(\mathbf{r})\,\phi_q(\mathbf{r})\,d\mathbf{r}\) is a Fourier-of-overlap MO integral. The latter has exactly the structure of an overlap integral once the substitutions of Yang et al. Eqs. (S19)-(S20) are applied, which lets TeraChem reuse its overlap-integral code. Closed-shell, active, and active-active blocks are pre-computed and contracted with the per-state RDMs.

Orientational averaging

A gas-phase experiment averages over molecular orientations. The scattering vector is decomposed as \(\mathbf{s}=s\,\hat{\mathbf{s}}\), and the angular average over \(\hat{\mathbf{s}}\) is done by Lebedev-Laikov quadrature of order diff_lvalue:

\[ \langle I(s)\rangle \;=\; \sum_g w_g\, W(\hat{\mathbf{s}}_g\!\cdot\!\hat{\mathbf{a}})\, I(s\,\hat{\mathbf{s}}_g), \]

where \(w_g\) are the Lebedev weights, \(\hat{\mathbf{a}}\) is the alignment axis set by diff_axis_{x,y,z}, and the additional weight \(W(\cos\theta)\) encodes any pre-alignment of the molecules with respect to \(\hat{\mathbf{a}}\):

`diff_rotation_averaging`	\(W(\cos\theta)\)	When to use
0	1 (isotropic)	Truly isotropic gas; alignment axis ignored.
1	\(\cos^2\theta\)	Molecules pre-aligned parallel to \(\hat{\mathbf{a}}\) (e.g. by a polarized pump pulse).
2	\(\sin^2\theta\)	Molecules pre-aligned perpendicular to \(\hat{\mathbf{a}}\).

Setting diff_lvalue to -1 (the default) disables averaging entirely and evaluates the signal at the single orientation \(\hat{\mathbf{s}}=\hat{\mathbf{a}}\).

Output

The diffraction block is written to standard output (the .out file) as one table per spin manifold (Singlet, Doublet, ...):

Singlet state scattering amplitudes:  axis=  0.0000   0.0000   1.0000
  s range(A-1):      0.0000 to    10.0000
  State       s(A-1)     I(X-Ray)        I(UED)      I(inelastic)
-----------------------------------------------------------------
 Singlet  0   0.0000   100.0000000     0.0000000    -0.0000000
 Singlet  0   0.1000    97.3014448     0.0005328     0.0145411
 ...

The State column indexes the eigenstates within the spin manifold (zero-based); the s(A-1) column is the magnitude of the scattering vector in \(\mathring{A}^{-1}\); the three remaining columns hold the quantities defined above. The s range(A-1): header line just above the table reports the same bounds in \(\mathring{A}^{-1}\), so it matches both the per-row s(A-1) column and the diff_s_start/diff_s_end values you supplied in the input file.

Inputs and example

The NH\(_3\) regression test runs CASSCF(8,9) on the ground and first excited singlet and integrates the signal isotropically over a 20th-order Lebedev grid for \(s\in[0,10]\,\mathring{A}^{-1}\):

ued_nh3.in

basis            ./def2-svp-nrydberg
coordinates      geom.xyz
charge           0
spinmult         1
method           hf
run              energy
precision        double
convthre         1.0e-6
threall          1.0e-14

closed           1
active           9
cassinglets      2
casscf           yes
casguess         c0.casscf

diffraction              yes
diff_lvalue              20
diff_s_start             0.0
diff_s_end              10.0
diff_delta_s             0.1
diff_axis_x              0.0
diff_axis_y              0.0
diff_axis_z              1.0
diff_rotation_averaging  0

Keyword reference

All scattering-vector ranges are specified in \(\mathring{A}^{-1}\). (TeraChem converts them to \(\mathrm{bohr}^{-1}\) internally, but both the s range(A-1): header line and the per-row data column are reported back in \(\mathring{A}^{-1}\).)

Keyword	Type	Default	Description
`diffraction`	bool	`no`	Enable the diffraction-signal evaluation.
`diff_lvalue`	int	`-1`	Order of the Lebedev-Laikov quadrature used for orientational averaging. `-1` skips the average and evaluates at the single direction \(\hat{\mathbf{a}}\) given by `diff_axis_*`.
`diff_s_start`	float	`0.0`	Lower bound of the scattering-vector magnitude \(s\) (\(\mathring{A}^{-1}\)).
`diff_s_end`	float	`10.0`	Upper bound of \(s\) (\(\mathring{A}^{-1}\)).
`diff_delta_s`	float	`0.1`	Step in \(s\) (\(\mathring{A}^{-1}\)).
`diff_axis_x`, `diff_axis_y`, `diff_axis_z`	float	`0, 0, 1`	Cartesian components of the alignment axis \(\hat{\mathbf{a}}\). Renormalized internally.
`diff_rotation_averaging`	int	`0`	\(\cos\theta\) weighting between \(\hat{\mathbf{s}}\) and \(\hat{\mathbf{a}}\): `0` isotropic, `1` \(\cos^2\theta\), `2` \(\sin^2\theta\). Has no effect when `diff_lvalue < 0`.

Limitations and notes

CAS-only. In the present source, the print_diffraction routine is called only from casci.cpp, casscf.cpp, and fomodft_casci.cpp. Although the underlying formalism (Yang et al.¹) applies equally to single-determinant SCF (HF, RKS) and CIS wavefunctions, those entry points are not wired up here. Use casci yes or casscf yes to access the diffraction calculation; for an HF or DFT reference, run a one-state CI on top of it.
Mott \(s^{-4}\) prefactor is omitted from UED. The I(UED) column is \(\lvert f_N - f\rvert^2\); multiply by \(s^{-4}\) to obtain the full differential cross section (and likewise for the inelastic column when comparing to electron-scattering data).
Output is printed, not saved separately. Capture the relevant block from the .out file (everything between the per-manifold header line and the next blank line).
Spin manifolds. The implementation iterates over Singlet through Septet manifolds; the counts come from cassinglets, casdoublets, ..., casseptets. Empty manifolds are skipped.
\(s=0\) sanity check. At \(s=0\) the X-ray column should equal \(n^2\) (square of total electron count) and the UED column should be zero (the nuclear and electronic Fourier amplitudes cancel for a neutral molecule). Both checks are visible in the NH\(_3\) test output (\(n=10\), \(I_{\text{X-ray}}(0)=100\), \(I_{\text{UED}}(0)=0\)).

J. Yang, X. Zhu, J. P. F. Nunes, J. K. Yu, R. M. Parrish, T. J. A. Wolf, M. Centurion, M. Gühr, R. Li, Y. Liu, B. Moore, M. Niebuhr, S. Park, X. Shen, S. Weathersby, T. Weinacht, T. J. Martínez, X. Wang, Science 368, 885 (2020). doi:10.1126/science.abb2235. See Supplementary Materials Eqs. (S5)-(S27) for the full derivation summarized here. ↩↩