chemistrykit.spectro #

Returns:

ndarray, shape matching x

extra: dict#

intensities: ndarray#

labels: list#

normalized()[source]#

Return a copy with intensities rescaled so the maximum is 1.

Return type:: Spectrum
Returns:: Spectrum

Examples

>>> spectrum = Spectrum(positions=np.array([1.0, 2.0]), intensities=np.array([4.0, 2.0]))
>>> spectrum.normalized().intensities.tolist()
[1.0, 0.5]

positions: ndarray#

class chemistrykit.spectro.TriatomicNormalModes(masses, equilibrium_coordinates, k_r1, k_r2, k_theta)[source]#

Bases: object

Normal-mode vibrational frequencies of an A-B-A triatomic, via the Wilson GF-matrix method.

A genuine (not tabulated) normal-mode calculation: an internal- coordinate force field is combined with the molecule’s actual equilibrium geometry and atomic masses through Wilson’s G (kinetic energy) and F (force constant) matrices, and the resulting generalized eigenvalue problem is solved for the true normal-mode frequencies (Wilson, Decius & Cross, Molecular Vibrations, 1955, Ch. 4). Three internal coordinates are used: the two bond stretches \(r_1\), \(r_2\), and the bond angle (or, for a linear equilibrium geometry, a linearized bend coordinate – see below).

The `G` matrix (mass-weighted, purely geometric) is built here numerically: the Wilson B-matrix, \(B_{ki}=\partial S_k/\partial x_i\) (internal coordinate k with respect to Cartesian displacement i), is obtained by central finite differences of the internal-coordinate functions, then \(G=BM^{-1}B^T\) with M the diagonal atomic-mass matrix – an exact numerical evaluation of the standard analytic B-matrix formulas, not an independent approximation.

The `F` matrix is a simple (diagonal) valence force field: \(F=\mathrm{diag}(k_{r_1},k_{r_2},k_\theta)\), with no stretch-stretch or stretch-bend interaction constants. Flagged as an approximation: a full experimental force field for a real molecule (e.g. CO2, H2O) generally does include small off-diagonal interaction terms, fit to reproduce observed frequencies to high precision; the diagonal approximation used here still gives the correct qualitative mode pattern (symmetric stretch, antisymmetric stretch, bend) and the right order of magnitude for reasonable literature-typical k_r/k_theta, but not spectroscopic-accuracy absolute wavenumbers.

Linear equilibrium geometries (is_linear=True) need a bend coordinate that stays well-defined exactly at 180 degrees, where the ordinary \(\theta=\arccos(\dots)\) bond-angle formula’s derivative diverges. The standard fix (Wilson, Decius & Cross, Ch. 4) is used: the bend coordinate is instead the linearized transverse (in-plane) displacement combination \(S_{\text{bend}}=(x_0-2x_1+x_2)/r_{eq}\), well-defined at exactly linear. A linear triatomic’s bend is genuinely doubly degenerate (in-plane and out-of-plane bending are equivalent by the molecule’s cylindrical symmetry); this module solves only the in-plane 3-internal-coordinate problem, so the reported bend wavenumber should be understood as representing both degenerate components (Herzberg, Molecular Spectra and Molecular Structure II, 1945, Ch. I.3).

The resulting symmetric generalized eigenvalue problem, GFL=L\Lambda, is solved here as the equivalent symmetric form \(Fc=\lambda G^{-1}c\) via chemistrykit.quantum.utils.secular_equation.solve_secular_equation() (with F playing the role of the Hamiltonian and \(G^{-1}\) the role of the overlap matrix) – reusing the same shared-eigenvalue- solve machinery chemistrykit.quantum’s variational models use, since both are literally the same kind of generalized symmetric eigenproblem, \(\lambda=\omega^2\) in SI units, converted to a wavenumber via \(\tilde\nu=\omega/(2\pi c)\).

Parameters:

masses (array-like of float, shape (3,)) – Atomic masses (atom0, central atom1, atom2), in kg.
equilibrium_coordinates (array-like of float, shape (3, 3)) – Cartesian equilibrium coordinates, in meters.
k_r1 (float) – Stretch force constants for bonds (atom0-atom1) and (atom2-atom1), in N/m.
k_r2 (float) – Stretch force constants for bonds (atom0-atom1) and (atom2-atom1), in N/m.
k_theta (float) – Bend force constant, in J/rad^2 (equivalently N m/rad^2).

Examples

See linear() and bent() for the usual way to construct one; solve() returns three positive wavenumbers for a stable equilibrium geometry.

classmethod bent(mass_terminal, mass_central, bond_length, bond_angle_degrees, k_r, k_theta)[source]#

Build a symmetric bent A-B-A triatomic (e.g. H2O-like) in the xz plane.

Parameters:

mass_terminal (float) – Mass of each terminal atom (A), in kg.
mass_central (float) – Mass of the central atom (B), in kg.
bond_length (float) – Equilibrium A-B bond length, in meters.
bond_angle_degrees (float) – Equilibrium A-B-A bond angle, in degrees (< 180).
k_r (float) – Shared stretch force constant for both bonds, in N/m.
k_theta (float) – Bend force constant, in J/rad^2.

Return type:

TriatomicNormalModes

Returns:

TriatomicNormalModes

classmethod linear(mass_terminal, mass_central, bond_length, k_r, k_theta)[source]#

Build a symmetric linear A-B-A triatomic (e.g. CO2-like) along the z axis.

Parameters:

mass_terminal (float) – Mass of each terminal atom (A), in kg.
mass_central (float) – Mass of the central atom (B), in kg.
bond_length (float) – Equilibrium A-B bond length, in meters.
k_r (float) – Shared stretch force constant for both bonds, in N/m.
k_theta (float) – Bend force constant, in J/rad^2.

Return type:

TriatomicNormalModes

Returns:

TriatomicNormalModes

solve()[source]#

Diagonalize the GF-matrix eigenproblem and return the normal-mode wavenumbers.

Return type:: NormalModeResult
Returns:: NormalModeResult

Examples

A bent triatomic (H2O-like) has 3 real, positive vibrational wavenumbers, with the bend lower in energy than either stretch – the correct qualitative ordering for a real bent AX2 molecule:

>>> import scipy.constants as sc
>>> water = TriatomicNormalModes.bent(
...     mass_terminal=1.008 * sc.atomic_mass,
...     mass_central=15.999 * sc.atomic_mass,
...     bond_length=95.8e-12,
...     bond_angle_degrees=104.5,
...     k_r=770.0,
...     k_theta=0.7e-18,
... )
>>> result = water.solve()
>>> bool(np.all(result.wavenumbers > 0))
True
>>> bool(result.wavenumbers[0] < result.wavenumbers[1] < result.wavenumbers[2])
True

chemistrykit.spectro.absorbance(molar_absorptivity, concentration, path_length)[source]#

The Beer-Lambert law: \(A=\varepsilon c l\).

Linear in concentration by construction – the exact relationship assumed to hold for a dilute, non-interacting, monochromatic-light measurement (Skoog, Holler & Crouch, Principles of Instrumental Analysis, 7th ed., Ch. 13.2). See apparent_absorbance_with_stray_light() for one well-known real instrumental effect that breaks this linearity at high absorbance.

Parameters:

molar_absorptivity (float) – Molar absorptivity \(\varepsilon\), in L mol^-1 cm^-1.
concentration (float or array-like of float) – Concentration c, in mol/L.
path_length (float) – Path length l, in cm.

Returns:

float or ndarray – Absorbance (dimensionless).

Examples

>>> round(absorbance(molar_absorptivity=5000.0, concentration=2.0e-5, path_length=1.0), 4)
0.1
>>> absorbance(5000.0, 4.0e-5, 1.0) == 2.0 * absorbance(5000.0, 2.0e-5, 1.0)
True

chemistrykit.spectro.anharmonicity_from_overtones(wavenumber_01, wavenumber_02)[source]#

Recover the Morse parameters \(\omega_e\), \(\omega_ex_e\) from two measured band positions.

The Morse-oscillator term-value formula gives the fundamental and first-overtone wavenumbers as

\[\tilde\nu_{0\to1} = \omega_e - 2\omega_ex_e, \qquad \tilde\nu_{0\to2} = 2\omega_e - 6\omega_ex_e\]

(Atkins & de Paula, Physical Chemistry, 11th ed., Ch. 12.4; Herzberg, Molecular Spectra and Molecular Structure I, 2nd ed., Ch. III), a linear system solved here exactly (no fitting) for the two spectroscopic constants \(\omega_e\) (the harmonic wavenumber) and \(\omega_ex_e\) (the anharmonicity constant) from two observed transitions – the standard way real vibrational anharmonicity constants are extracted from an observed overtone spectrum.

Parameters:

wavenumber_01 (float) – Observed \(0\to1\) (fundamental) wavenumber, in cm^-1.
wavenumber_02 (float) – Observed \(0\to2\) (first overtone) wavenumber, in cm^-1.

Return type:

tuple

Returns:

omega_e (float) – Harmonic wavenumber, in cm^-1.
omega_e_xe (float) – Anharmonicity constant, in cm^-1.

Examples

Round-trip against a MorseOscillator built from known omega_e/omega_e_xe-equivalent parameters recovers them from its own predicted band positions:

>>> morse = MorseOscillator(mass=1.6e-27, force_constant=500.0, dissociation_energy=7.0e-19)
>>> nu01, nu02 = morse_transition_wavenumbers(morse, v_max=2)
>>> omega_e, omega_e_xe = anharmonicity_from_overtones(nu01, nu02)
>>> bool(omega_e > nu01 > 0)  # the harmonic wavenumber exceeds the (red-shifted) observed fundamental
True
>>> bool(omega_e_xe > 0)
True

chemistrykit.spectro.apparent_absorbance_with_stray_light(molar_absorptivity, concentration, path_length, stray_light_fraction)[source]#

The apparent absorbance measured in the presence of instrumental stray light.

Approximation flagged explicitly: real Beer-Lambert-law deviations at high concentration have several distinct physical origins (a concentration-dependent chemical equilibrium shifting the true absorbing species’ concentration, use of polychromatic rather than strictly monochromatic radiation, and instrumental stray light); this function models only the last of those, one specific, well-characterized, purely instrumental effect (Skoog, Holler & Crouch, Principles of Instrumental Analysis, 7th ed., Ch. 13.4).

A small fraction \(p_s\) of the light reaching the detector is “stray” radiation that never actually passed through the fully absorbing sample path (e.g. light scattered around the sample cell), so it is not attenuated by the sample the way the main beam is. The detector reads the combined transmitted-plus-stray intensity, so the apparent absorbance is

\[A_{\text{apparent}} = -\log_{10}\!\left(10^{-A_{\text{true}}} + p_s\right)\]

(normalized so that at stray_light_fraction=0, or in the dilute limit \(A_{\text{true}}\to0\), this reduces exactly to the true Beer-Lambert absorbance). Because the stray-light term does not shrink as the sample absorbs more strongly, it dominates the detector reading at high true absorbance, causing the classic negative (sublinear, “rolling over”) deviation from Beer-Lambert linearity seen in real instruments at high concentration.

Parameters:

molar_absorptivity (float) – In L mol^-1 cm^-1.
concentration (float) – In mol/L.
path_length (float) – In cm.
stray_light_fraction (float) – Fraction of incident intensity that reaches the detector as stray light, \(p_s\in[0,1)\).

Return type:

Returns:

float – Apparent absorbance, always \(\le A_{\text{true}}\).

Raises:

ValueError – If stray_light_fraction is not in \([0, 1)\).

Examples

With no stray light, the apparent absorbance is exactly the true one:

>>> round(apparent_absorbance_with_stray_light(5000.0, 2.0e-5, 1.0, stray_light_fraction=0.0), 6)
0.1

At low concentration (low true absorbance) the deviation is negligible, but it grows sharply at high concentration – Beer- Lambert linearity breaks down exactly where the textbook warns it does:

>>> low = apparent_absorbance_with_stray_light(5000.0, 1.0e-6, 1.0, stray_light_fraction=0.001)
>>> low_true = absorbance(5000.0, 1.0e-6, 1.0)
>>> high = apparent_absorbance_with_stray_light(5000.0, 1.0e-3, 1.0, stray_light_fraction=0.001)
>>> high_true = absorbance(5000.0, 1.0e-3, 1.0)
>>> bool(abs(low - low_true) < abs(high - high_true))
True

chemistrykit.spectro.broaden_stick_spectrum(centers, intensities, x, shape='gaussian', fwhm=1.0, fwhm_lorentzian=None)[source]#

Sum a lineshape over every (center, intensity) stick to build a continuous spectrum.

Parameters:

centers (array-like of float) – Peak positions.
intensities (array-like of float) – Relative intensity (area) of each peak, same length as centers.
x (array-like of float) – Grid to evaluate the broadened spectrum on.
shape (str) – Lineshape to use for every peak.
fwhm (float) – FWHM (or the Gaussian-component FWHM, for shape=”voigt”).
fwhm_lorentzian (float, optional) – Lorentzian-component FWHM, required (and only used) for shape=”voigt”.

Returns:

ndarray, shape matching x

Raises:

ValueError – For an unrecognized shape, or a missing fwhm_lorentzian when shape=”voigt”.

Examples

A single stick recovers the plain lineshape function, scaled by its intensity:

>>> import numpy as np
>>> x = np.linspace(-10, 10, 501)
>>> spectrum = broaden_stick_spectrum([0.0], [2.0], x, shape="gaussian", fwhm=1.0)
>>> bool(np.allclose(spectrum, 2.0 * gaussian(x, 0.0, 1.0)))
True

chemistrykit.spectro.chemical_shift_ppm(frequency_hz, reference_frequency_hz)[source]#

Chemical shift \(\delta=10^6(\nu-\nu_\text{ref})/\nu_\text{ref}\), in ppm.

Dividing by the reference frequency removes the (field-proportional) Larmor frequency, so \(\delta\) is the same on every spectrometer – the reason NMR spectra are reported in ppm.

Parameters:

frequency_hz (float or array-like of float) – Resonance frequency of the nucleus of interest, in Hz.
reference_frequency_hz (float) – Resonance frequency of the reference compound (e.g. TMS), in Hz.

Returns:

float or ndarray

Examples

A proton resonating 2800 Hz above TMS on a 400 MHz instrument sits at 7.0 ppm:

>>> round(chemical_shift_ppm(400.0e6 + 2800.0, 400.0e6), 6)
7.0

chemistrykit.spectro.concentration_from_absorbance(absorbance_value, molar_absorptivity, path_length)[source]#

Invert the Beer-Lambert law for concentration: \(c=A/(\varepsilon l)\).

Parameters:

absorbance_value (float or array-like of float)
molar_absorptivity (float) – In L mol^-1 cm^-1.
path_length (float) – In cm.

Returns:

float or ndarray – Concentration, in mol/L.

Examples

>>> c = concentration_from_absorbance(absorbance_value=0.1, molar_absorptivity=5000.0, path_length=1.0)
>>> round(c, 8)
2e-05

chemistrykit.spectro.energy_to_wavenumber(energy_j)[source]#

Convert an energy (in J) to a spectroscopic wavenumber \(\tilde\nu=E/(hc)\), in cm^-1.

Parameters:: energy_j (float or array-like of float)
Return type:: float
Returns:: float or ndarray

Examples

>>> round(float(energy_to_wavenumber(H * C * 100.0)), 6)
1.0

chemistrykit.spectro.fid_to_spectrum(fid, dwell_time)[source]#

Fourier-transform an FID into a frequency-domain (absorption-mode) spectrum.

Uses the discrete Fourier transform with the first point halved (the standard correction that removes a constant baseline offset) and returns the real (absorption) part, frequency axis centered on zero. A line of transverse relaxation time \(T_2\) becomes a Lorentzian of FWHM \(1/(\pi T_2)\) Hz.

Parameters:

fid (array-like of complex) – Uniformly sampled FID, starting at \(t=0\).
dwell_time (float) – Sampling interval, in s.

Return type:

tuple

Returns:

frequencies_hz (ndarray) – Ascending frequency axis, in Hz.
spectrum (ndarray) – Real (absorption-mode) spectrum, in units of FID amplitude times seconds.

Examples

One resonance at +50 Hz Fourier-transforms into a single peak at +50 Hz:

>>> t = np.arange(4096) * 1e-3
>>> freqs, spec = fid_to_spectrum(free_induction_decay(t, [50.0], [1.0], t2=0.2), 1e-3)
>>> round(float(freqs[np.argmax(spec)]), 1)
50.0

chemistrykit.spectro.first_order_multiplet(chemical_shift_ppm, j_coupling_hz, n_neighbors, spectrometer_frequency_mhz)[source]#

Simulate a simple multiplet from coupling to one set of n equivalent neighbors.

Line positions are evenly spaced by \(J/\nu_0\) in ppm (J in Hz converted to ppm by dividing by the spectrometer frequency in MHz, since 1 ppm on the chemical-shift scale corresponds to \(\nu_0\,[\text{Hz}]/10^6\) at operating frequency \(\nu_0\)), symmetric about chemical_shift_ppm, with the Pascal’s-triangle relative intensities from pascals_triangle_intensities().

Parameters:

chemical_shift_ppm (float) – Unperturbed chemical shift, in ppm.
j_coupling_hz (float) – Coupling constant J, in Hz.
n_neighbors (int) – Number of equivalent coupled neighbors.
spectrometer_frequency_mhz (float) – Spectrometer (proton, or whichever nucleus) operating frequency, in MHz – needed to convert J from Hz to the ppm scale.

Return type:

Returns:

Spectrum – positions in ppm, intensities the relative (unnormalized) binomial line intensities.

Examples

A triplet, symmetric about the chemical shift, spaced by J:

>>> spectrum = first_order_multiplet(chemical_shift_ppm=1.2, j_coupling_hz=7.0, n_neighbors=2, spectrometer_frequency_mhz=400.0)
>>> len(spectrum.positions)
3
>>> round(float(np.mean(spectrum.positions)), 6) == round(1.2, 6)
True
>>> spectrum.intensities.tolist()
[1.0, 2.0, 1.0]

chemistrykit.spectro.franck_condon_factor(v, S)[source]#

The \(0\to v\) Franck-Condon factor for two displaced, equal-frequency harmonic oscillators.

\[\text{FC}(0\to v) = \frac{e^{-S}S^v}{v!}\]

(a Poisson distribution in v with mean S; see the module docstring for the underlying approximation).

Parameters:

v (int) – Final-state vibrational quantum number, >= 0.
S (float) – Huang-Rhys parameter, >= 0.

Return type:

Returns:

float – In \([0, 1]\).

Examples

Zero displacement (S=0) puts all intensity in the 0->0 transition (no vibronic structure at all):

>>> round(franck_condon_factor(0, S=0.0), 6)
1.0
>>> round(franck_condon_factor(1, S=0.0), 6)
0.0

chemistrykit.spectro.franck_condon_progression(v_max, S)[source]#

Tabulate the Franck-Condon factors for \(v=0,1,\dots,v_{max}\).

Parameters:

v_max (int)
S (float) – Huang-Rhys parameter.

Return type:

Returns:

ndarray, shape (v_max + 1,)

Examples

The Franck-Condon factors are (approximately) a complete probability distribution over the final vibrational level, so they sum to (approximately) 1 once v_max comfortably exceeds the Poisson distribution’s mean S – the required normalization check for any genuine Franck-Condon calculation:

>>> round(float(np.sum(franck_condon_progression(v_max=40, S=3.0))), 9)
1.0

The most probable (highest-intensity) vibronic transition lands near \(v'\approx S\) – the qualitative rule of thumb for reading off a vibronic progression’s displacement from its band shape:

>>> progression = franck_condon_progression(v_max=20, S=5.0)
>>> bool(abs(int(np.argmax(progression)) - 5) <= 1)
True

chemistrykit.spectro.franck_condon_spectrum(origin_wavenumber, vibrational_wavenumber, S, v_max)[source]#

Build a vibronic-progression Spectrum.

Stick positions are the electronic origin (the \(0\to0\) transition) plus v’ quanta of the excited-state vibrational wavenumber; stick intensities are the corresponding Franck-Condon factors (Atkins & de Paula, Physical Chemistry, 11th ed., Ch. 12.7(b)).

Parameters:

origin_wavenumber (float) – The \(0\to0\) (electronic origin) transition wavenumber, in cm^-1.
vibrational_wavenumber (float) – Excited-state vibrational spacing, in cm^-1.
S (float) – Huang-Rhys parameter.
v_max (int) – Highest excited-state vibrational level to include.

Return type:

Returns:

Spectrum – positions in cm^-1, intensities the Franck-Condon factors.

Examples

>>> spectrum = franck_condon_spectrum(origin_wavenumber=25000.0, vibrational_wavenumber=1500.0, S=1.5, v_max=10)
>>> round(float(spectrum.positions[0]), 1)
25000.0
>>> round(float(spectrum.positions[2] - spectrum.positions[1]), 1)
1500.0

chemistrykit.spectro.free_induction_decay(t, offsets_hz, amplitudes, t2)[source]#

Complex free-induction decay (FID) after a single 90-degree pulse.

Each resonance, at offset \(\Delta\nu_k\) from the rotating-frame reference, contributes a decaying complex exponential,

\[s(t)=\sum_k a_k\,e^{2\pi i\Delta\nu_k t}\,e^{-t/T_2},\]

the time-domain signal whose Fourier transform is the spectrum (R. R. Ernst & W. A. Anderson, Rev. Sci. Instrum. 37, 93 (1966)).

Parameters:

t (array-like of float) – Sample times, in s.
offsets_hz (array-like of float) – Resonance offsets, in Hz.
amplitudes (array-like of float) – Relative amplitude of each resonance.
t2 (float) – Transverse relaxation time \(T_2\), in s (shared by all lines).

Return type:

Returns:

ndarray of complex, shape matching t

Examples

At \(t=0\) every component is in phase, so the signal is the sum of the amplitudes:

>>> fid = free_induction_decay([0.0, 0.1], offsets_hz=[10.0, -25.0], amplitudes=[1.0, 2.0], t2=0.5)
>>> complex(fid[0])
(3+0j)

chemistrykit.spectro.gaussian(x, x0, fwhm)[source]#

A normalized Gaussian lineshape.

\[g(x) = \frac{2}{\text{FWHM}}\sqrt{\frac{\ln 2}{\pi}} \exp\!\left[-4\ln 2\left(\frac{x-x_0}{\text{FWHM}}\right)^2\right]\]

Normalized so \(\int_{-\infty}^{\infty}g(x)\,dx=1\) (Hollas, Modern Spectroscopy, 4th ed., Ch. 2.3).

Parameters:

x (float or array-like of float) – Evaluation point(s).
x0 (float) – Peak center.
fwhm (float) – Full width at half maximum, same units as x.

Returns:

float or ndarray

Examples

The peak height at x0 and the normalization both follow the closed-form prefactor:

>>> import numpy as np
>>> expected_peak = 2.0 / 1.0 * np.sqrt(np.log(2) / np.pi)
>>> round(float(gaussian(0.0, 0.0, fwhm=1.0)), 6) == round(float(expected_peak), 6)
True
>>> x = np.linspace(-20, 20, 200001)
>>> round(float(np.trapezoid(gaussian(x, 0.0, fwhm=2.0), x)), 3)
1.0

chemistrykit.spectro.harmonic_fundamental_wavenumber(force_constant, reduced_mass)[source]#

The harmonic-oscillator fundamental (\(v=0\to1\)) IR band position.

\[\tilde\nu = \frac{1}{2\pi c}\sqrt{\frac{k}{\mu}}\]

Every harmonic-oscillator transition sits at exactly this same wavenumber (evenly spaced levels), unlike the Morse oscillator’s overtones (Atkins & de Paula, Physical Chemistry, 11th ed., Ch. 12.4).

Parameters:

force_constant (float) – Bond force constant k, in N/m.
reduced_mass (float) – Vibrational reduced mass, in kg.

Return type:

Returns:

float – Wavenumber, in cm^-1.

Examples

>>> import scipy.constants as sc
>>> mu = (1.008 * 34.97) / (1.008 + 34.97) * sc.atomic_mass  # HCl-like reduced mass
>>> nu = harmonic_fundamental_wavenumber(force_constant=480.0, reduced_mass=mu)
>>> bool(2500.0 < nu < 3200.0)  # HCl's real fundamental is 2886 cm^-1
True

chemistrykit.spectro.huang_rhys_factor(displacement, mass, angular_frequency)[source]#

The dimensionless Huang-Rhys displacement parameter \(S=\frac12\frac{m\omega}{\hbar}\Delta Q^2\).

Parameters:

displacement (float) – Excited-minus-ground-state equilibrium displacement along the normal coordinate, \(\Delta Q\), in m (for a mass-weighted coordinate) or the natural length unit of the coordinate used.
mass (float) – Oscillator (reduced) mass, in kg.
angular_frequency (float) – Shared vibrational angular frequency \(\omega\), in rad/s.

Return type:

Returns:

float

Examples

>>> round(huang_rhys_factor(displacement=0.0, mass=1.6e-27, angular_frequency=5.0e13), 6)
0.0

chemistrykit.spectro.isotope_shift_ratio(mass1a, mass2, mass1b)[source]#

Ratio of rotational constants (and transition wavenumbers) between two isotopologues.

Substituting atom 1 (e.g. H) for a heavier isotope (e.g. D) changes only the reduced mass \(\mu=m_1m_2/(m_1+m_2)\), and the rotational constant \(B\propto1/I=1/(\mu r^2)\) scales inversely with it (bond length r is, to a very good approximation, unchanged by isotopic substitution, since it depends on the electronic structure, not the nuclear mass – the Born-Oppenheimer approximation). So

\[\frac{B(\text{isotopologue }b)}{B(\text{isotopologue }a)} = \frac{\mu_a}{\mu_b}\]

(Atkins & de Paula, Physical Chemistry, 11th ed., Ch. 12.3). Every rotational transition wavenumber scales by this same ratio, since they are all proportional to B.

Parameters:

mass1a (float) – Mass of the lighter isotope of the substituted atom (e.g. H), in any consistent unit.
mass2 (float) – Mass of the other (unsubstituted) atom, same units.
mass1b (float) – Mass of the heavier isotope of the substituted atom (e.g. D), same units.

Return type:

Returns:

float – \(B_b/B_a\) (and \(\tilde\nu_b/\tilde\nu_a\)), < 1 since substituting a heavier isotope increases the reduced mass and moment of inertia, lowering B.

Examples

Substituting D for H in HCl roughly halves the rotational constant (deuterium is about twice as heavy as hydrogen, and Cl is much heavier than either, so the reduced mass roughly doubles):

>>> import scipy.constants as sc
>>> ratio = isotope_shift_ratio(mass1a=1.008 * sc.atomic_mass, mass2=34.97 * sc.atomic_mass, mass1b=2.014 * sc.atomic_mass)
>>> bool(0.45 < ratio < 0.55)
True

Directly matches a from-scratch RigidRotor comparison:

>>> r = 127.5e-12
>>> rotor_hcl = RigidRotor.from_diatomic(1.008 * sc.atomic_mass, 34.97 * sc.atomic_mass, r)
>>> rotor_dcl = RigidRotor.from_diatomic(2.014 * sc.atomic_mass, 34.97 * sc.atomic_mass, r)
>>> ratio_direct = rotor_dcl.rotational_constant / rotor_hcl.rotational_constant
>>> bool(abs(ratio_direct - ratio) < 1e-9)
True

chemistrykit.spectro.karplus_coupling(dihedral_degrees, coefficients=None)[source]#

Vicinal H-C-C-H coupling constant \(^3J\) from the dihedral angle: the Karplus relation.

With coefficients=None this is Karplus’s original 1959 valence-bond result (M. Karplus, J. Chem. Phys. 30, 11 (1959)),

\[\begin{split}^3J(\phi) = \begin{cases} 8.5\cos^2\phi - 0.28, & 0^\circ\le\phi\le90^\circ \\ 9.5\cos^2\phi - 0.28, & 90^\circ<\phi\le180^\circ \end{cases}\ \text{Hz};\end{split}\]

otherwise the general three-term form \(^3J=A\cos^2\phi+B\cos\phi+C\) (M. Karplus, J. Am. Chem. Soc. 85, 2870 (1963)) with coefficients=(A, B, C) in Hz.

Parameters:

dihedral_degrees (float or array-like of float) – H-C-C-H dihedral angle \(\phi\), in degrees.
coefficients (tuple of float, optional) – (A, B, C) in Hz for the three-term form.

Returns:

float or ndarray – \(^3J\), in Hz.

Examples

Anti (180 degrees) protons couple strongly, gauche (60 degrees) weakly, and perpendicular (90 degrees) protons hardly at all:

>>> round(karplus_coupling(180.0), 2), round(karplus_coupling(60.0), 3), round(karplus_coupling(90.0), 2)
(9.22, 1.845, -0.28)

chemistrykit.spectro.larmor_frequency(gyromagnetic_ratio, field_tesla)[source]#

Nuclear Larmor (resonance) frequency \(\nu_0=\gamma B_0/(2\pi)\), in Hz.

The frequency at which Bloch’s and Purcell’s 1946 experiments found a sharp resonant response from nuclear spins in a static field \(B_0\) (F. Bloch, Phys. Rev. 70, 460 (1946); Atkins & de Paula, Physical Chemistry, 11th ed., Ch. 12A).

Parameters:

gyromagnetic_ratio (float) – Nuclear gyromagnetic ratio \(\gamma\), in rad s^-1 T^-1 (e.g. scipy.constants.physical_constants["proton gyromag. ratio"][0]).
field_tesla (float or array-like of float) – Static field \(B_0\), in tesla.

Returns:

float or ndarray – Resonance frequency, in Hz.

Examples

A “400 MHz” spectrometer’s magnet is about 9.4 T:

>>> import scipy.constants as sc
>>> gamma_h = sc.physical_constants["proton gyromag. ratio"][0]
>>> round(larmor_frequency(gamma_h, 9.3947) / 1e6, 1)
400.0

chemistrykit.spectro.lorentzian(x, x0, fwhm)[source]#

A normalized Lorentzian lineshape.

\[L(x) = \frac{1}{\pi}\frac{\text{FWHM}/2}{(x-x_0)^2+(\text{FWHM}/2)^2}\]

Normalized so \(\int_{-\infty}^{\infty}L(x)\,dx=1\) (Hollas, Modern Spectroscopy, 4th ed., Ch. 2.3).

Parameters:

x (float or array-like of float)
x0 (float) – Peak center.
fwhm (float) – Full width at half maximum, same units as x.

Returns:

float or ndarray

Examples

>>> import numpy as np
>>> x = np.linspace(-2000, 2000, 4_000_001)
>>> round(float(np.trapezoid(lorentzian(x, 0.0, fwhm=2.0), x)), 3)
1.0

By definition, the value has dropped to half its peak at \(x_0\pm\text{FWHM}/2\):

>>> peak = lorentzian(0.0, 0.0, fwhm=2.0)
>>> half_max_point = lorentzian(1.0, 0.0, fwhm=2.0)
>>> round(float(half_max_point), 9) == round(float(peak) / 2.0, 9)
True

chemistrykit.spectro.morse_transition_wavenumbers(morse, v_max)[source]#

Wavenumbers of the \(0\to v\) Morse-oscillator “overtone” transitions, \(v=1,\dots,v_{max}\).

Unlike the harmonic oscillator, successive Morse transitions (fundamental \(0\to1\), first overtone \(0\to2\), …) are not simple multiples of each other – the level spacing shrinks with v (Atkins & de Paula, Physical Chemistry, 11th ed., Ch. 12.4).

Parameters:

morse (MorseOscillator)
v_max (int) – Highest final vibrational level to include.

Return type:

Returns:

ndarray, shape (v_max,) – Wavenumbers of the \(0\to1,0\to2,\dots,0\to v_{max}\) transitions, in cm^-1.

Examples

The first overtone (\(0\to2\)) falls short of exactly twice the fundamental (\(0\to1\)) – the textbook anharmonicity signature:

>>> morse = MorseOscillator(mass=1.6e-27, force_constant=500.0, dissociation_energy=7.0e-19)
>>> wavenumbers = morse_transition_wavenumbers(morse, v_max=2)
>>> bool(wavenumbers[1] < 2.0 * wavenumbers[0])
True

chemistrykit.spectro.multi_coupling_multiplet(chemical_shift_ppm, couplings, spectrometer_frequency_mhz)[source]#

Simulate a multiplet from independent first-order coupling to several inequivalent neighbor sets.

Each coupling constant splits every line already present (the couplings act independently and multiplicatively), giving \(\prod_k(n_k+1)\) lines in the general case (fewer if some coincide) – e.g. a doublet of triplets (dt) from one neighbor with \(J_1\) and two equivalent neighbors with \(J_2\ne J_1\).

Parameters:

chemical_shift_ppm (float) – Unperturbed chemical shift, in ppm.
couplings (sequence of tuple(float, int)) – (J_hz, n_equivalent_neighbors) for each independent coupling partner set.
spectrometer_frequency_mhz (float)

Return type:

Returns:

Spectrum – positions in ppm (sorted ascending), intensities the combined relative intensities (coincident positions within numerical tolerance are merged, summing their intensities).

Examples

A doublet of doublets (dd) with equal J’s degenerates into a triplet with 1:2:1 intensities (two different coupling partners with the same J are indistinguishable from one set of two equivalent partners):

>>> spectrum = multi_coupling_multiplet(chemical_shift_ppm=5.0, couplings=[(7.0, 1), (7.0, 1)], spectrometer_frequency_mhz=400.0)
>>> len(spectrum.positions)
3
>>> spectrum.intensities.tolist()
[1.0, 2.0, 1.0]

A genuine doublet of triplets (\(J_1=12\) Hz to one neighbor, \(J_2=5\) Hz to two equivalent neighbors) has 6 distinct lines:

>>> dt = multi_coupling_multiplet(chemical_shift_ppm=3.0, couplings=[(12.0, 1), (5.0, 2)], spectrometer_frequency_mhz=400.0)
>>> len(dt.positions)
6
>>> round(float(np.sum(dt.intensities)), 6)
8.0

chemistrykit.spectro.multiplicity(n_equivalent_neighbors)[source]#

The n+1 rule: n equivalent spin-1/2 neighbors split a resonance into n+1 lines.

Parameters:: n_equivalent_neighbors (int) – Number of magnetically equivalent, spin-1/2, first-order-coupled neighboring nuclei.
Return type:: int
Returns:: int

Examples

Ethanol’s -CH2- protons, coupled to the 3 equivalent -CH3- protons, form a quartet:

>>> multiplicity(3)
4

chemistrykit.spectro.pascals_triangle_intensities(n_equivalent_neighbors)[source]#

Relative line intensities for n equivalent spin-1/2 neighbors: the binomial coefficients \(\binom{n}{k}\).

Parameters:: n_equivalent_neighbors (int)
Return type:: ndarray
Returns:: ndarray, shape (n_equivalent_neighbors + 1,)

Examples

A quartet’s classic 1:3:3:1 intensity pattern (e.g. ethanol’s CH2, coupled to 3 equivalent CH3 protons):

>>> pascals_triangle_intensities(3)
array([1., 3., 3., 1.])

chemistrykit.spectro.rotational_line_wavenumbers(rotor, J_max)[source]#

Wavenumbers of the \(J\to J+1\) rotational absorption lines, \(J=0,\dots,J_{max}\).

Each line sits at \(\tilde\nu_J=2B(J+1)\) (in wavenumber units), evenly spaced by \(2B\) – the defining signature of a rigid- rotor rotational spectrum (Atkins & de Paula, Physical Chemistry, 11th ed., Ch. 12.2).

Parameters:

rotor (RigidRotor)
J_max (int) – Highest initial J to include.

Return type:

Returns:

ndarray, shape (J_max + 1,) – Wavenumbers, in cm^-1.

Examples

>>> import scipy.constants as sc
>>> rotor = RigidRotor.from_diatomic(mass1=1.008 * sc.atomic_mass, mass2=34.97 * sc.atomic_mass, bond_length=127.5e-12)
>>> lines = rotational_line_wavenumbers(rotor, J_max=4)
>>> spacings = np.diff(lines)
>>> bool(np.allclose(spacings, spacings[0]))
True

chemistrykit.spectro.rotational_spectrum(rotor, J_max, temperature)[source]#

Build a rotational-spectrum Spectrum, with relative intensities.

Approximation flagged explicitly: the relative intensity of the \(J\to J+1\) line is modeled here as proportional to the thermal (Boltzmann) population of the initial level J, weighted by its \((2J+1)\)-fold degeneracy, times a \((J+1)\) line-strength factor coming from the \(\Delta J=\pm1\) transition-dipole matrix element (Atkins & de Paula, Physical Chemistry, 11th ed., Ch. 12.2(b)):

\[I_J \propto (J+1)(2J+1)\exp\!\left[-\frac{E_J}{k_BT}\right]\]

This reproduces the qualitative textbook feature of a rotational spectrum – intensity rising from J=0, peaking at some intermediate J, then falling off – but is not a first-principles transition- dipole calculation for any particular real molecule.

Parameters:

rotor (RigidRotor)
J_max (int) – Highest initial J to include.
temperature (float) – Temperature, in K, for the Boltzmann population weighting.

Return type: