Re: [ESPResSo] bug in random.tcl in mbtools

[Torsten Stuehn]

Hi all,

to generate Random Walk Polymers (polymer-command), ESPResSo (internally)
uses the following algorithm to generate homogeneously distributed points
on the surface of a sphere with radius=bond_length:

        zz     = (2.0*d_random()-1.0)*bond_length;
        rr     = sqrt(SQR(bond_length)-SQR(zz));
        phi    = 2.0*PI*d_random();
        x      = rr*cos(phi);
        y      = rr*sin(phi);
        z      = zz;

The advantage of this algorithm is, that it needs only two random
numbers and you never have to throw away any numbers.
Experience has shown, that on modern CPUs this algorithm runs faster
than the one Burkhard proposed.
It is based on the (not very intuitive) fact, that the projection of
homogeneously distributed points on a sphere-surface onto any axis gives a
normal distribution of points on that axis (a prove of that can be found
in math books).

Greetings,
Torsten


Tristan Bereau wrote:
> Hello,
>
> I'm not sure whether a gaussian random number generator has already  
> been implemented in the TCL shell of ESPResSo, but if not, there is a  
> couple of methods described at the end of the "Old Testament" of  
> molecular simulations (i.e. Allen and Tildesley), and Burkhard's is  
> one of them.
>
> Best,
>
> Tristan
>
>
> On Mar 19, 2009, at 4:52 AM, Jacob Kirkensgaard wrote:
>
> > Hi again
> >
> > Thanks to Tristan and Burkhard - you are of course right. For my  
> > purposes my suggestion is sufficient but I will see if I have time  
> > to implement a small routine to do it properly as outlined by  
> > Burkhard. Is there a routine for generating Gaussian random numbers?  
> > As Tristan mentioned this would be a quick way to solve the issue.
> >
> > Best,
> > Jacob
> >
> >
> >
> > On Mar 19, 2009, at 9:40 AM, address@hidden wrote:
> >
> > > Hello,
> > >
> > > from Tristan's remarks I understand that there
> > > is a question of how to get random numbers
> > > uniformly in or on a sphere? If so, let me tell
> > > you one possible answer:
> > >
> > > 1. You first draw 3 RNs u_1, u_2, u_3 from [0:1]
> > > 2. Then you transform u_1 = 2 * u_1 - 1 etc.,
> > >  resulting in RNs in [-1,1]
> > > 3. Calculate r2 = u_1**2 + u_2**2 + u_3**2
> > > 4. If r2 > 1 --> throw everything away,
> > >  go to step 1, try again
> > > 5. If r2 < 1: Either you are done (supposing
> > >  that you want RN IN the sphere), or you set
> > >  norm = 1 / sqrt(r2)
> > >  u_1 = u_1 * norm etc.
> > >  By this you project the output onto the
> > >  surface of the unit sphere, so you get
> > >  something ON the sphere.
> > >
> > > If that was trivial: Please accept my apologies
> > > for bothering you.
> > >
> > > Regards
> > > Burkhard.
>
>
> _______________________________________________
> ESPResSo mailing list
> address@hidden
> https://fias.uni-frankfurt.de/mailman/listinfo/espresso


--
Dr. Torsten Stühn
MPI für Polymerforschung
Ackermannweg 10
55128 Mainz / Germany

Tel.  +49-(0)6131-379268
Fax   +49-(0)6131-379100 
EMail address@hidden