Hi Troels, Note that you should do this for all R1rho model types. Exceptions are the 'M61' and 'M61 skew' models which are on-resonance, hence they should return arrays of R1rho' as theta is always 0.0.
Regards, Edward On 27 May 2014 09:14, Edward d'Auvergne <[email protected]> wrote: > 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
