AMBER Archive (2008)

Subject: Re: AMBER: note on the molsurf and "Assertion `sarg1 >= 0.0' failed" error

From: Kateryna Miroshnychenko (kateryna_mirosh_at_ire.kharkov.ua)
Date: Thu May 29 2008 - 02:22:17 CDT

Funny, but it does work, at least in my particular case. Don't know why ;)
But you are right of course, the sarg should be equaled to zero before sqrt.
See below the last variant of molsurf modification. Please find attached the
pqr file, that caused the problem. (sarg1 = -3.66e-16)

With best regards,
Katya

*** /usr/local/amber9/src/mm_pbsa/molsurf.c_bak 2008-05-29 09:58:45.000000000
+0300
--- /usr/local/amber9/src/mm_pbsa/molsurf.c 2008-05-29 10:20:24.000000000
+0300
*************** static int is_buried (ai, aj, ak, probe_
*** 1030,1036 ****
    POINT t_ij, t_ik, t_diff;
    REAL_T b_factor, h_ijk, xtest;
    REAL_T w_ijk, w_sin, r_pk2;
! REAL_T sarg1, sarg2;
    POINT b_ijk, r_ib;
    void cross ();
    int i;
--- 1030,1036 ----
    POINT t_ij, t_ik, t_diff;
    REAL_T b_factor, h_ijk, xtest;
    REAL_T w_ijk, w_sin, r_pk2;
! REAL_T sarg1, sarg2, eps=0.001;
    POINT b_ijk, r_ib;
    void cross ();
    int i;
*************** static int is_buried (ai, aj, ak, probe_
*** 1052,1057 ****
--- 1052,1062 ----
    sarg1 = (ai->rad + aj->rad + probe_diam) * (ai->rad + aj->rad +
                  probe_diam) - d_ij * d_ij;
    sarg2 = d_ij * d_ij - (ai->rad - aj->rad) * (ai->rad - aj->rad);
+
+ if( sarg1 < 0.0 && fabs(sarg1)<=eps)
+ sarg1=0.0;
+ if( sarg2 < 0.0 && fabs(sarg2)<=eps)
+ sarg2=0.0;
    assert( sarg1 >= 0.0 );
    assert( sarg2 >= 0.0 );
    t_rad = 0.5 * sqrt (sarg1) * sqrt (sarg2) / d_ij;

On Wednesday 28 May 2008 20:27:31 David A. Case wrote:
> Thanks for your detailed analysis, but I really don't understand why this
> would work. In your case, you are allowing sarg1 to be (slightly)
> negative, and then immediately computing sqrt(sarg1). The square root will
> fail unless sarg >=0 (as in the original assert).
>
> It is possible that we should allow for very slightly negative values, but
> then the "solution" would have to involve setting these to zero before
> proceeding to the sqrt() code (at least as far as I can see).
>
> Can you post the structure that is giving the problem?
>
> ...thx...dac
>
> -----------------------------------------------------------------------
> The AMBER Mail Reflector
> To post, send mail to amber_at_scripps.edu
> To unsubscribe, send "unsubscribe amber" (in the *body* of the email)
> to majordomo_at_scripps.edu


-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber_at_scripps.edu
To unsubscribe, send "unsubscribe amber" (in the *body* of the email)
      to majordomo_at_scripps.edu