Hi Troels,
I think this TP02 model in the 'no Rex' edge cases should return:
R1rho = R1 cos^2(theta) + R1rho' sin^2(theta) (1)
rather than:
R1rho = R1rho' (2)
The problem is in the unit tests in
test_suite/unit_tests/_lib/_dispersion/test_tp02.py, specifically the
calc_r1rho() method. This method should calculate the equation (1)
for R1rho and check against that in all cases, rather than the current
check:
# Check all R1rho values.
if self.kex > 1.e5:
for i in range(self.num_points):
self.assertAlmostEqual(R1rho[i], self.r1, 6)
else:
for i in range(self.num_points):
self.assertAlmostEqual(R1rho[i], self.r1rho_prime)
You can see this in equation 11.68 in the manual
(http://www.nmr-relax.com/manual/TP02_2_site_exchange_R1_model.html),
the TP02 model equation. The values of pA = 1.0 (hence pB = 0.0), kex
= 0.0, and dw = 0.0 will all cause the numerator in the third part of
the equation to be zero, so all that is left is equation (1) above.
Large values of kex will cause the other parts of the third part of
equation 11.68 to be insignificant so you have kex / kex^2 which heads
to zero as kex heads to infinity. So then you have the third part
disappearing, leaving equation (1) again.
So you would need to calculate R1rho as in equation (1) and modify the checks:
# Calculate R1rho in the absence of exchange.
r1rho_no_rex = ...
# Check all R1rho values.
if self.kex > 1.e5:
for i in range(self.num_points):
self.assertAlmostEqual(R1rho[i], r1rho_no_rex, 6)
else:
for i in range(self.num_points):
self.assertAlmostEqual(R1rho[i], r1rho_no_rex)
Cheers,
Edward
On 26 May 2014 23:09, <[email protected]> wrote:
> Author: tlinnet
> Date: Mon May 26 23:09:56 2014
> New Revision: 23445
>
> URL: http://svn.gna.org/viewcvs/relax?rev=23445&view=rev
> Log:
> Critical fix for the math domain catching of model TP02.
>
> This was discovered with the added 8 unit tests demonstrating edge case 'no
> Rex' failures.
>
> This follows from the ideas in the post
> http://article.gmane.org/gmane.science.nmr.relax.devel/5858.
> This is related to: task #7793: (https://gna.org/task/?7793) Speed-up of
> dispersion models.
>
> This is to implement catching of math domain errors, before they occur.
>
> Modified:
> branches/disp_speed/lib/dispersion/tp02.py
>
> Modified: branches/disp_speed/lib/dispersion/tp02.py
> URL:
> http://svn.gna.org/viewcvs/relax/branches/disp_speed/lib/dispersion/tp02.py?rev=23445&r1=23444&r2=23445&view=diff
> ==============================================================================
> --- branches/disp_speed/lib/dispersion/tp02.py (original)
> +++ branches/disp_speed/lib/dispersion/tp02.py Mon May 26 23:09:56 2014
> @@ -123,18 +123,13 @@
>
> # Catch zeros (to avoid pointless mathematical operations).
> # This will result in no exchange, returning flat lines.
> - if min(numer) == 0.0:
> - return R1_cos_theta2 + R1rho_prime_sin_theta2
> + if numer == 0.0:
> + return array([r1rho_prime]*num_points)
>
> # Denominator.
> denom = waeff2 * wbeff2 / weff2 + kex2
> #denom_extended = waeff2*wbeff2/weff2+kex2-2*sin_theta2*pA*pB*dw**2
>
> - # Catch math domain error of dividing with 0.
> - # This is when denom=0.
> - if min(abs(denom)) == 0:
> - return array([1e100]*num_points)
> -
> # R1rho calculation.
> R1rho = R1_cos_theta2 + R1rho_prime_sin_theta2 + sin_theta2 * numer /
> denom
>
>
>
> _______________________________________________
> relax (http://www.nmr-relax.com)
>
> This is the relax-commits mailing list
> [email protected]
>
> To unsubscribe from this list, get a password
> reminder, or change your subscription options,
> visit the list information page at
> https://mail.gna.org/listinfo/relax-commits
_______________________________________________
relax (http://www.nmr-relax.com)
This is the relax-devel mailing list
[email protected]
To unsubscribe from this list, get a password
reminder, or change your subscription options,
visit the list information page at
https://mail.gna.org/listinfo/relax-devel
