[ESPResSo] RE: [ESPResSo-devel] bugs in polymer.c

[Limbach, Hans Joerg, LAUSANNE, NRC-FS]

Dear Peter,

Thanks for the routines. I'll try to push some people to implement this,
since I also do not have the time to do this.
For the future, please always send the mails to the mailing list.

Hanjo 

> -----Original Message-----
> From: Peter Kosovan [mailto:address@hidden 
> Sent: jeudi, 1. novembre 2007 12:52
> To: Limbach,Hans Joerg,LAUSANNE,NRC-FS
> Subject: RE: [ESPResSo-devel] bugs in polymer.c
> 
> Dear Hanjo
> 
> As for the crude routine, I think it is not so bad, 
> especially if it is not used too often. I am sorry but I do 
> not have enough time to dive into the source to correctly 
> implement the function. I still feel more like a passive 
> reader of the code rather than a developer.
> 
> Yet I can at least provide the C code of functions which are 
> being used in by group in Prague so that someone else can use 
> and implement them. One efficient method for generating the 
> random vector is that of Marsaglia. 
> The code was implemented and tested some time ago by my 
> colleague Filip Uhlik. Filip says that this issue is also 
> discussed in Allen&Tildesley but I do not have the book at 
> hand now to check it.
> 
> void s2_deviates(double *x,double *y,double *z)
> /* uniform deviates on a surface of a unit sphere, Marsaglia 
> method */ {
>    double t,t1,t2;
>    do{
>      t1=1.0-2.0*udev(); t2=1.0-2.0*udev();
>      t=SQR(t1)+SQR(t2);
>    }while(t>1.0);
>    *z=1.0-2.0*t;
>    t=sqrt(1.0-t);
>    *x=2.0*t1*t; *y=2.0*t2*t;
>    return;
> }
> 
> Another option is the one which was pointed out to me by 
> Milan Predota. 
> Probably it can be found in one of his papers. It genereates 
> a random point on a surface of a cylinder and then uses the 
> fact that a ring on the sphere of constant surface area 
> corresponds to a given increment dz, independent of the value of z.
> 
> // assumes result points to a space allocated for double 
> result[3] {x,y,z} void rnd_orient (double *result) {
>    double sintheta;
>    double phi=2.0*PI*d_random();
>    result[2]=1.0-2.0*d_random();
>    sintheta=sqrt(1-result[2]*result[2]);
>    result[1]=sintheta*cos(phi);
>    result[0]=sintheta*sin(phi);
>    return;
> }
> 
> peter
> 
> On Tue, 30 Oct 2007, Limbach,Hans Joerg,LAUSANNE,NRC-FS wrote:
> 
> > Date: Tue, 30 Oct 2007 14:48:06 +0100
> > From: "Limbach,Hans Joerg,LAUSANNE,NRC-FS"
> >     <address@hidden>
> > To: Peter Kosovan <address@hidden>,
> >     address@hidden
> > Subject: RE: [ESPResSo-devel] bugs in polymer.c
> > 
> > Dear Peter, dear all
> >
> > Thank you for pointing out this bug. There is a very crude routine 
> > doing exactly this on the tcl level in scripts/ABHmath.tcl. Maybe 
> > someone knows a better way. I did not spent a whole lot of 
> time when I 
> > wrote that piece of code.
> > It would be great if you could write a corresponding c-function and 
> > incorporate it in random.h as an MDINLINE function (and send us an 
> > updated version of random.h then so that we can place it in 
> the next 
> > version).
> >
> > With best regards,
> > Hanjo
> >
> >> -----Original Message-----
> >> From: address@hidden
> >> [mailto:address@hidden On Behalf Of 
> >> Peter Kosovan
> >> Sent: lundi, 29. octobre 2007 15:11
> >> To: address@hidden
> >> Subject: [ESPResSo-devel] bugs in polymer.c
> >>
> >> Though I am a member of the user mailing list, I think this better 
> >> fits into the developers' list.
> >>
> >> I found out that a call to polymer command with a fixed angle 
> >> produces very different results in the version 1.9.7h and 
> 2.0.2n even 
> >> though I have taken care of the change of the angle definition. 
> >> Examining polymer.c (version 2.0.2n), on the line 440 I found
> >>
> >> phi = 2.0*d_random();
> >>
> >> where should be
> >>
> >> phi = 2.0*PI*d_random();
> >>
> >> Using the polymer command with angle but without angle2 
> thus produces 
> >> wierd helix-like structures with prefered absolute orientation in 
> >> space.
> >> When both angle and angle2 are used, everything is correct.
> >>
> >> Also there is another minor bug repeated several times in the same 
> >> file. Using
> >>
> >> theta  = PI*d_random();
> >> phi    = 2.0*PI*d_random();
> >> pos[0] = poz[0]+bond_length*sin(theta)*cos(phi);
> >> pos[1] = poz[1]+bond_length*sin(theta)*sin(phi);
> >> pos[2] = poz[2]+bond_length*cos(theta);
> >>
> >> does not place a new monomer uniformly on the surface of a unit 
> >> sphere around the previous one, as it is supposed to.
> >> Since the surface of 2*PI*sin(theta)*dtheta corresponds to each 
> >> dtheta, this code actually prefers the new monomer to be placed 
> >> around the poles (theta=0).
> >>
> >> With regards
> >>
> >> Peter
> >>
> >> Peter Kosovan
> >> Department of Physical and Macromolecular Chemistry Faculty of 
> >> Science Charles University in Prague Czech Republic 
> >> address@hidden http://lynette.natur.cuni.cz 
> Tel. +420 
> >> 221 951 290 Currently on a stay at the Wageningen University, the 
> >> Netherlands
> >>
> >>
> >> _______________________________________________
> >> ESPResSo-devel mailing list
> >> address@hidden
> >> https://fias.uni-frankfurt.de/mailman/listinfo/espresso-devel
> >>
> >
>