Re: [ESPResSo] possible bug in rigid bonds in forces.h ?

[Christoph Junghans]

Hi Ana,

you are right, this is how it looks in the cvs version:

#ifdef BOND_CONSTRAINT
    case BONDED_IA_RIGID_BOND:
      //add_rigid_bond_pair_force(p1,p2, iaparams, force, force2);
      bond_broken = 0;
      force[0]=force[1]=force[2]=0.0;
      break;
#endif

So the same line was already added by me on the April 9th.

I think we (who is that?) just forgot to backport it to 2.0.5l.

Bye,

Christoph

Ana Vila Verde schrieb:
> Hi all,
>
> I noticed that in forces.h (espresso-2.0.5l), array forces is not 
> initialized in the part corresponding to rigid bonds.
> #ifdef BOND_CONSTRAINT
>     case BONDED_IA_RIGID_BOND:
>       //add_rigid_bond_pair_force(p1,p2, iaparams, force, force2);
>       bond_broken = 0;
>       break;
> #endif
>
> Shouldn't it be initialized   to zero to prevent unphysical forces to be 
> added to the particles? (see "case 1" further down in the subroutine)
> #ifdef BOND_CONSTRAINT
>     case BONDED_IA_RIGID_BOND:
>       //add_rigid_bond_pair_force(p1,p2, iaparams, force, force2);
>       bond_broken = 0;
>       f*orce[0]=0.; force[1]=0;force[2]=0; // added by ana*
>       break;
> #endif
>
> I realize that most likely I am missing something that explains why this 
> initialization is not necessary, but I've struggled with this for a 
> while so I figured it was time to ask the experts.  Thank you all for 
> any help you can provide here.
>
> Ana

--
Dipl.-Phys. Christoph Junghans
Max Planck Institute for Polymer Research
Theory Group
POBox 3148
D 55021 Mainz, Germany

Phone: +49 6131 379 335
Web: http://www.mpip-mainz.mpg.de/~junghans