fixed kww calculation with array of betas

This commit is contained in:
dominik 2022-03-23 15:32:56 +01:00
parent 6cd630245c
commit 22217f0544
4 changed files with 81 additions and 6 deletions

View File

@ -514,7 +514,7 @@ class NMRMainWindow(QtWidgets.QMainWindow, Ui_BaseWindow):
def _select_t1tauwidget(self, onoff: bool, pick_required: bool, block_window: bool): def _select_t1tauwidget(self, onoff: bool, pick_required: bool, block_window: bool):
if onoff: # tau from t1 if onoff: # tau from t1
if self.current_graph_widget is None: if self.current_graph_widget is None:
return return pick_required, block_window
idx = self.tabWidget.indexOf(self.t1tauwidget) idx = self.tabWidget.indexOf(self.t1tauwidget)
if len(self.current_graph_widget) != 1: if len(self.current_graph_widget) != 1:
@ -522,7 +522,7 @@ class NMRMainWindow(QtWidgets.QMainWindow, Ui_BaseWindow):
'Only one T1 curve can be evaluated at once') 'Only one T1 curve can be evaluated at once')
self.tabWidget.removeTab(idx) self.tabWidget.removeTab(idx)
self.tabWidget.setCurrentIndex(0) self.tabWidget.setCurrentIndex(0)
return return pick_required, block_window
self.tabWidget.setTabText(idx, 'T1 min (%s)' % {self.current_graph_widget.title}) self.tabWidget.setTabText(idx, 'T1 min (%s)' % {self.current_graph_widget.title})

View File

@ -1,7 +1,7 @@
from typing import List, Tuple from typing import List, Tuple
from ...nmr.coupling import * from ...nmr.coupling import *
from ...distributions import ColeCole, ColeDavidson, HavriliakNegami from ...distributions import ColeCole, ColeDavidson, HavriliakNegami, KWW
from ...utils import pi from ...utils import pi
from ...utils.text import convert from ...utils.text import convert
@ -19,7 +19,7 @@ class QRelaxCalc(QtWidgets.QDialog, Ui_Dialog):
self.graphs = {} self.graphs = {}
self.specdens = [ColeCole, ColeDavidson, HavriliakNegami] self.specdens = [ColeCole, ColeDavidson, HavriliakNegami, KWW]
self.coupling = [Quadrupolar, Czjzek, HomoDipolar] self.coupling = [Quadrupolar, Czjzek, HomoDipolar]
self.tau_parameter = [] self.tau_parameter = []
@ -162,9 +162,9 @@ class QRelaxCalc(QtWidgets.QDialog, Ui_Dialog):
if isinstance(p, Widget): if isinstance(p, Widget):
parameter.append(p.value) parameter.append(p.value)
elif isinstance(p, SelectionWidget): elif isinstance(p, SelectionWidget):
kwargs[p.argname] = p.value kwargs.update(p.value)
else: else:
raise TypeError('WTF?: Unknown widget', p) raise TypeError('WTF: Unknown widget', p)
dic['cp_param'] = (parameter, kwargs) dic['cp_param'] = (parameter, kwargs)

View File

@ -204,6 +204,24 @@ def _kwws_lim_hig(b):
def kww_cos(w, tau, beta): def kww_cos(w, tau, beta):
beta = np.asanyarray(beta)
if beta.ndim == 0:
return _kww_cos_single(w, tau, beta)
else:
tau = np.asanyarray(tau)
if tau.size != beta.size:
raise ValueError('KWW mismatched number of beta and tau')
else:
beta.reshape(tau.shape)
ret_val = np.zeros((w.size, tau.size))
for i, (tau_i, beta_i) in enumerate(zip(tau.ravel(), beta.ravel())):
ret_val[:, i] = _kww_cos_single(w, tau_i, beta_i)
return ret_val.reshape(-1, *tau.shape)
def _kww_cos_single(w, tau, beta):
""" """
Args: Args:
@ -246,6 +264,24 @@ def kww_cos(w, tau, beta):
# \int_0^\infty dt sin(w*t) exp(-(t/tau)^beta) # \int_0^\infty dt sin(w*t) exp(-(t/tau)^beta)
def kww_sin(w, tau, beta): def kww_sin(w, tau, beta):
beta = np.asanyarray(beta)
if beta.ndim == 0:
return _kww_sin_single(w, tau, beta)
else:
tau = np.asanyarray(tau)
if tau.size != beta.size:
raise ValueError('KWW mismatched number of beta and tau')
else:
beta.reshape(tau.shape)
ret_val = np.zeros((w.size, tau.size))
for i, (tau_i, beta_i) in enumerate(zip(tau.ravel(), beta.ravel())):
ret_val[:, i] = _kww_sin_single(w, tau_i, beta_i)
return ret_val.reshape(-1, *tau.shape)
def _kww_sin_single(w, tau, beta):
# check input data # check input data
if not (0.1 <= beta <= 2.): if not (0.1 <= beta <= 2.):
raise ValueError('KWW beta is not inside range 0.1-2') raise ValueError('KWW beta is not inside range 0.1-2')

View File

@ -104,6 +104,8 @@ class Relaxation:
return self.t1_bpp(omega, tau, *specdens_args, **kwargs) return self.t1_bpp(omega, tau, *specdens_args, **kwargs)
elif mode == 'dipolar': elif mode == 'dipolar':
return self.t1_dipolar(omega, tau, *specdens_args, **kwargs) return self.t1_dipolar(omega, tau, *specdens_args, **kwargs)
else:
return self.t1_csa(omega, tau, *specdens_args, **kwargs)
def t1_dipolar(self, omega: ArrayLike, tau: ArrayLike, *specdens_args: Any, inverse: bool = True, def t1_dipolar(self, omega: ArrayLike, tau: ArrayLike, *specdens_args: Any, inverse: bool = True,
prefactor: float = None, omega_coup: ArrayLike = None, prefactor: float = None, omega_coup: ArrayLike = None,
@ -302,6 +304,43 @@ class Relaxation:
else: else:
return rate return rate
def t2(self, omega: ArrayLike, tau: ArrayLike, *specdens_args: Any,
mode: str = 'bpp', **kwargs) -> Union[np.ndarray, float]:
r"""
Convenience function
Args:
omega (array-like): Frequency in 1/s (not Hz)
tau (array-like): Correlation times in s
*specdens_args: Additional parameter for spectral density.
If given this will replace previously set arguments during calculation
mode (str, {`bpp`, `dipolar`, `csa`}): Type of relaxation
- `bpp` : Homonuclear dipolar/quadrupolar interaction :func:`(see here) <nmreval.nmr.Relaxation.t1_bpp>`
- `dipolar` : Heteronuclear dipolar interaction :func:`(see here) <nmreval.nmr.Relaxation.t1_dipolar>`
- `csa` : Chemical shift interaction :func:`(see here) <nmreval.nmr.Relaxation.t1_csa>`
Default is `bpp`
**kwargs:
Returns:
float or ndarray:
A k by n array where k is length of ``omega`` and n is length of ``tau``.
If both are of length 1 a number is returned instead of array.
"""
if mode not in ['bpp', 'dipolar', 'csa']:
raise ValueError(f'Unknown mode {mode} not `bpp`, `dipolar`, `csa`.')
# defaults to BPP
if mode == 'bpp':
return self.t2_bpp(omega, tau, *specdens_args, **kwargs)
elif mode == 'dipolar':
return self.t2_dipolar(omega, tau, *specdens_args, **kwargs)
else:
return self.t2_csa(omega, tau, *specdens_args, **kwargs)
def t2_bpp(self, omega: ArrayLike, tau: ArrayLike, *specdens_args: Any, def t2_bpp(self, omega: ArrayLike, tau: ArrayLike, *specdens_args: Any,
inverse: bool = True, prefactor: float = None) -> Union[np.ndarray, float]: inverse: bool = True, prefactor: float = None) -> Union[np.ndarray, float]:
r"""Calculate T2 under homonuclear dipolar coupling or quadrupolar coupling. r"""Calculate T2 under homonuclear dipolar coupling or quadrupolar coupling.