Re: [ESPResSo] BOND_ANGLE_HARMONIC

[Jon Halverson]

Hi all,

I agree with Mehmet. The same holds for the dihedral angle. I ran into 
this problem with another code when studying hydrocarbons which were 
initialized in the all-trans conformation where the angle is 180 degrees 
and cos(angle) may give something outside of /pm 1 due to round-off.

According to Berk Hess, GROMACS does what Mehmet is proposing. 
Furthermore, bonded interactions are inexpensive.

Jon


Mehmet Sayar wrote:
> Hi all,
>
> In the bond_angle_harmonic routine shouldn't there be a check for
> cosine > 1.0. In other words, I think the code should be corrected as:
>
> #ifdef BOND_ANGLE_HARMONIC
>   {
>     double phi,sinphi;
>     if ( cosine >  TINY_COS_VALUE ) cosine = TINY_COS_VALUE;
>     if ( cosine < -TINY_COS_VALUE)  cosine = -TINY_COS_VALUE;
>     phi =  acos(-cosine);
>     sinphi = sin(phi);
>     if ( sinphi < TINY_SIN_VALUE ) sinphi = TINY_SIN_VALUE;
>     fac *= (phi - iaparams->p.angle.phi0)/sinphi;
>   }
> #endif
>
> If everybody agrees, I will add this to the cvs.
>
> P.S. I originally send this message to address@hidden But since 
> that email did not show up in the email-thread, I am re-sending the message.
> Is address@hidden still in use?
>
> And I already got an answer from Hanjo to the email sent to address@hidden, 
> which was:
>
>  > So far we have the checks only for functions with a singularity when
>  > cosine is aproaching 1 (and sinphi is aproaching 0). Can the scalar
>  > product yield values >1 ( <-1)? (which would yield an error in the acos
>  > function). If yes then add
>  >      if ( cosine >  1 ) cosine = 1;
>  >      if ( cosine < -1)  cosine = -1;
>  > To the corresponding force and energy routines.
>  >
>  > Greetings Hanjo
>
> Does everyone agree? Any objections?
>
> Greetings from Istanbul,
>
> Mehmet