# Source code for nmrsim.discrete

"""Non-quantum mechanical solutions for specific second-order patterns.

These are adapted from the routines from WINDNMR _ by Hans Reich,
U. Wisconsin, and equations from Pople, Schneider and Bernstein _. Note that
many of the names for arguments, etc. are non-Pythonic but chosen to match the
WINDNMR interface and source code (for now).

The partials module provides the following functions:

* AB: simulates an AB quartet.

* AB2: simulates an AB2 system.

* ABX: simulates an ABX system.

* ABX3: simulates an ABX3 system.

* AAXX: simulates an AA'XX' system.

* AABB: simulates an AA'BB' system.

References
----------
..  Pople, J.A.; Schneider, W.G.; Bernstein, H.J. *High-Resolution Nuclear Magnetic Resonance.* New York:
McGraw-Hill, 1959.

"""
# TODO: refactor away from legacy argument/variable names to user-friendly,
# pythonic ones.
from math import sqrt

import numpy as np

from nmrsim.firstorder import multiplet
from nmrsim.math import _normalize

[docs]def AB(Jab, Vab, Vcentr, normalize=True):
"""
Calculates the signal frequencies and intensities for two strongly
coupled protons (Ha and Hb).

Parameters
---------
Jab : float
The coupling constant (Hz) between Ha and Hb
Vab : float
The chemical shift difference (Hz) between Ha and Hb in the absence
of coupling.
Vcentr : float
The frequency (Hz) for the center of the AB quartet.
normalize: bool
Whether the signal intensity should be normalized.

Returns
-------
[(float, float)...]
A list of four (frequency, intensity) tuples.
"""
J = Jab
dv = Vab
c = ((dv ** 2 + J ** 2) ** 0.5) / 2
center = Vcentr
v1 = center - c - (J / 2)
v2 = v1 + J
v3 = center + c - (J / 2)
v4 = v3 + J
dI = J / (2 * c)
I1 = 1 - dI
I2 = 1 + dI
I3 = I2
I4 = I1
vList = [v1, v2, v3, v4]
IList = [I1, I2, I3, I4]
if normalize:
_normalize(IList, 2)
return list(zip(vList, IList))

[docs]def AB2(Jab, Vab, Vcentr, normalize=True):
"""
Calculates signal frequencies and intensities for an AB2 spin system.

Parameters
---------
Jab : float
the Ha-Hb coupling constant (Hz).
Vab : float
the difference in the frequencies (Hz). A positive value means vb > va;
negative means va > vb.
Vcentr : float
the frequency (Hz) for the center of the AB2 signal.
normalize: bool
whether the signal intensity should be normalized.

Returns
-------
[(float, float)...]
a list of (frequency, intensity) tuples.
"""
# There is a disconnect between the variable names in the WINDNMR GUI and
# the variable names in this function.
# The following code provides a temporary interface until this is
# refactored.
J, dV, Vab = Jab, Vab, Vcentr

# Also, inconsistencies in WINDNMR GUI, internal WINDNMR code, and Pople
# equations require a conversion.
dV = - dV
va = Vab + (dV / 2)
vb = va - dV

Jmod = J * (3 / 4)
C_plus = sqrt(dV ** 2 + dV * J + (9 / 4) * (J ** 2)) / 2
C_minus = sqrt(dV ** 2 - dV * J + (9 / 4) * (J ** 2)) / 2
cos2theta_plus = (dV / 2 + J / 4) / C_plus
cos2theta_minus = (dV / 2 - J / 4) / C_minus
sintheta_plus = sqrt((1 - cos2theta_plus) / 2)
sintheta_minus = sqrt((1 - cos2theta_minus) / 2)
costheta_plus = sqrt((1 + cos2theta_plus) / 2)
costheta_minus = sqrt((1 + cos2theta_minus) / 2)
sin_dtheta = sintheta_plus * costheta_minus - costheta_plus * sintheta_minus
cos_dtheta = costheta_plus * costheta_minus + sintheta_plus * sintheta_minus

# In Pople, Schneider and Bernstein, Table 6-8:
# V1-V4 are "Origin: A";
# V5-V8 are "Origin: B";
# V9 is "Origin: Comb."
V1 = Vab + Jmod + C_plus
V2 = vb + C_plus + C_minus
V3 = va
V4 = Vab - Jmod + C_minus
V5 = vb + C_plus - C_minus
V6 = Vab + Jmod - C_plus
V7 = vb - C_plus + C_minus
V8 = Vab - Jmod - C_minus
V9 = vb - C_plus - C_minus

I1 = (sqrt(2) * sintheta_plus - costheta_plus) ** 2
I2 = (sqrt(2) * sin_dtheta + costheta_plus * costheta_minus) ** 2
I3 = 1
I4 = (sqrt(2) * sintheta_minus + costheta_minus) ** 2
I5 = (sqrt(2) * cos_dtheta + costheta_plus * sintheta_minus) ** 2
I6 = (sqrt(2) * costheta_plus + sintheta_plus) ** 2
I7 = (sqrt(2) * cos_dtheta - sintheta_plus * costheta_minus) ** 2
I8 = (sqrt(2) * costheta_minus - sintheta_minus) ** 2
I9 = (sqrt(2) * sin_dtheta + sintheta_plus * sintheta_minus) ** 2
vList = [V1, V2, V3, V4, V5, V6, V7, V8, V9]
IList = [I1, I2, I3, I4, I5, I6, I7, I8, I9]

if normalize:
_normalize(IList, 3)
return list(zip(vList, IList))

[docs]def ABX(Jab, Jax, Jbx, Vab, Vcentr, vx, normalize=True):
"""
Non-QM approximation for an ABX spin system. The approximation assumes
that Hx is very far away in chemical shift from Ha/Hb.

Parameters
---------
Jab : float
The Ha-Hb coupling constant (Hz).
Jax : float
The Ha-Hx coupling constant (Hz).
Jbx : float
The Hb-Hx coupling constant (Hz).
Vab : float
The difference in the frequencies (in the absence of
coupling) of Ha and Hb (Hz).
Vcentr : float
The frequency (Hz) for the center of the AB signal.
vx : float
The frequency (Hz) for Hx in the absence of coupling.

normalize: bool (optional)
whether the signal intensity should be normalized. If false, the total
signal intensity happens to be ~12.

Returns
-------
[(float, float)...]
a list of (frequency, intensity) tuples.
"""
# Contradictions in naming between WINDNMR's interface, internal code, and
# Pople/Schneider/Bernstein "fixed" with these reassignments:
dVab = Vab
Vab = Vcentr

dJx = Jax - Jbx
cm = dJx / 2
cp = Jax + Jbx
M = dVab + cm
L = dVab - cm
D_plus = sqrt(M ** 2 + Jab ** 2) / 2
D_minus = sqrt(L ** 2 + Jab ** 2) / 2
sin2phi_plus = Jab / (2 * D_plus)
sin2phi_minus = Jab / (2 * D_minus)
cos2phi_plus = M / (2 * D_plus)
cos2phi_minus = L / (2 * D_minus)
m = (cp + 2 * Jab) / 4
n = (cp - 2 * Jab) / 4
t = cos2phi_plus * cos2phi_minus + sin2phi_plus * sin2phi_minus

# In Pople, Schneider and Bernstein, Table 6-15:
# V1-V4 are "Origin: B" (PSB Table 6-15);
# V5-V8 are "Origin: A";
# V9-V12 are "Origin: X";
# V13-14 are "Origin: Comb. (X)"
V1 = Vab - m - D_minus
V2 = Vab + n - D_plus
V3 = Vab - n - D_minus
V4 = Vab + m - D_plus
V5 = Vab - m + D_minus
V6 = Vab + n + D_plus
V7 = Vab - n + D_minus
V8 = Vab + m + D_plus
V9 = vx - cp / 2
V10 = vx + D_plus - D_minus
V11 = vx - D_plus + D_minus
V12 = vx + cp / 2
V13 = vx - D_plus - D_minus
V14 = vx + D_plus + D_minus
I1 = 1 - sin2phi_minus
I2 = 1 - sin2phi_plus
I3 = 1 + sin2phi_minus
I4 = 1 + sin2phi_plus
I5 = I3
I6 = I4
I7 = I1
I8 = I2
I9 = 1
I10 = (1 + t) / 2
I11 = I10
I12 = 1
I13 = (1 - t) / 2
I14 = I13
VList = [V1, V2, V3, V4, V5, V6, V7, V8, V9, V10, V11, V12, V13, V14]
IList = [I1, I2, I3, I4, I5, I6, I7, I8, I9, I10, I11, I12, I13, I14]
if normalize:
_normalize(IList, 3)
return list(zip(VList, IList))

[docs]def ABX3(Jab, Jax, Jbx, Vab, Vcentr):
"""
Simulation of the AB part of an ABX3 spin system.

Parameters
---------
Jab : float
the Ha-Hb coupling constant (Hz).
Jax : float
the Ha-Hb coupling constant (Hz).
Jbx : float
the Ha-Hb coupling constant (Hz).
Vab : float
the difference in the frequencies (Hz) of Ha and Hb in the absence of
coupling. Positive when vb > va.
Vcentr : float
the frequency (Hz) for the center of the AB signal.

Returns
-------
[(float, float)...]
a list of (frequency, intensity) tuples.
"""
# First: simulate two quartets for va and vb ("Jab turned off")
va = Vcentr - Vab / 2
vb = Vcentr + Vab / 2
a_quartet = multiplet((va, 1), [(Jax, 3)])
b_quartet = multiplet((vb, 1), [(Jbx, 3)])
res = []
# Then: for each pair of a and b singlets in the quartets, calculate an
# AB quartet ("Turn Jab on").
for i in range(4):
dv = b_quartet[i] - a_quartet[i]
abcenter = (b_quartet[i] + a_quartet[i]) / 2
sub_abq = AB(Jab, dv, abcenter, normalize=True)
scale_factor = a_quartet[i]
scaled_sub_abq = [(v, i * scale_factor) for v, i in sub_abq]
res.extend(scaled_sub_abq)
return res

[docs]def AAXX(Jaa, Jxx, Jax, Jax_prime, Vcentr, normalize=True):
"""
Simulates one half ('A' part) of an AA'XX' spin system.

All frequencies are in Hz.

Parameters
---------
Jaa, Jax, Jax, Jax_prime : float
Jaa is the JAA' coupling constant;
Jxx the JXX';
Jax the JAX; and
JAX_prime the JAX'.
Vcentr : float
the frequency for the center of the signal.
normalize: bool
whether the signal intensity should be normalized (to 2).

Returns
-------
[(float, float)...]
a list of (frequency, intensity) tuples.
"""
# Define the constants required to calculate frequencies and intensities

# K, L, M, N are as defined in Pople, Schneider and Bernstein
K = Jaa + Jxx
M = Jaa - Jxx
L = Jax - Jax_prime
N = Jax + Jax_prime
p = sqrt((K ** 2 + L ** 2)) / 2
r = sqrt((M ** 2 + L ** 2)) / 2
sin2theta_s = (1 - K / (2 * p)) / 2
sin2theta_a = (1 - M / (2 * r)) / 2
cos2theta_s = (1 + K / (2 * p)) / 2
cos2theta_a = (1 + M / (2 * r)) / 2

# See PSB Table 6-18. Transitions 1-4 are condensed into V1 and V2.
V1 = Vcentr + N / 2
V2 = Vcentr - N / 2
V3 = Vcentr + K / 2 + p
V4 = Vcentr - K / 2 + p
V5 = Vcentr + K / 2 - p
V6 = Vcentr - K / 2 - p
V7 = Vcentr + M / 2 + r
V8 = Vcentr - M / 2 + r
V9 = Vcentr + M / 2 - r
V10 = Vcentr - M / 2 - r

I1 = 2
I2 = I1
I3 = sin2theta_s
I4 = cos2theta_s
I5 = I4
I6 = I3
I7 = sin2theta_a
I8 = cos2theta_a
I9 = I8
I10 = I7

VList = [V1, V2, V3, V4, V5, V6, V7, V8, V9, V10]
IList = [I1, I2, I3, I4, I5, I6, I7, I8, I9, I10]
if normalize:
_normalize(IList, 2)
return list(zip(VList, IList))

[docs]def AABB(Vab, Jaa, Jbb, Jab, Jab_prime, Vcentr, normalize=True, **kwargs):
"""
A wrapper for a second-order AA'BB' calculation, but using the
same arguments as WINDNMR.

Parameters
---------
Vab : float
the difference in frequency (Hz) between Ha and Hb in the absence of
coupling. A positive number indicates vb > va.
Jaa, Jbb, Jab, Jab_prime : float
Jaa is the JAA' coupling constant;
Jxx the JXX';
Jax the JAX; and
JAX_prime the JAX'.
Vcentr : float
the frequency for the center of the signal.
normalize: bool
whether the signal intensity should be normalized.

Returns
-------
[(float, float)...]
a list of (frequency, intensity) tuples.
"""
from nmrsim.qm import qm_spinsystem
va = Vcentr - Vab / 2
vb = Vcentr + Vab / 2
freqlist = [va, va, vb, vb]
J = np.zeros((4, 4))
J[0, 1] = Jaa
J[0, 2] = Jab
J[0, 3] = Jab_prime
J[1, 2] = Jab_prime
J[1, 3] = Jab
J[2, 3] = Jbb
J = J + J.T

result = qm_spinsystem(freqlist, J, normalize=normalize, sparse=False, **kwargs)
return result