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, <tlinnet@xxxxxxxxxxxxx> 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 relax-commits@xxxxxxx 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