espressomd-users
[Top][All Lists]

## Re: [ESPResSo] BOND_ANGLE_HARMONIC

 From: Jon Halverson Subject: Re: [ESPResSo] BOND_ANGLE_HARMONIC Date: Wed, 09 Jul 2008 14:25:18 +0200 User-agent: Thunderbird 2.0.0.12 (X11/20080213)

```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?

```
```
> 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
```
```

```