[Pw_forum] Mean Square Displacement (MSD) of atoms

[Giovanni La Penna]

>
> How to remove the artificial
> translation when processing the coordinates?
>
>
If boxc[] = a,b,c
of an orthorhombic cell and the origin is in 0,0,0
a piece of the C code for storing the variables to be
correlated in a molecular liquid at time t is:

         for(k=0;k<3;k++) {
          boxcm[k]=0.5*boxc[k];
         }

         for(im=0; im<nmol; im++) {
          for(ia=0; ia<natom; ia++) {
           for (k=0;k<3;k++,crdiold++) {
            dri=(*crdi-*crdiold);
            *crdiold=*crdi;
            if(dri >  boxcm[k]) (*jumppi)++;
            if(dri < -boxcm[k]) (*jumpmi)++;
            x = (*crdi++)-(*jumppi-*jumpmi)*boxc[k];
            *stor++ = x;
            jumppi++; jumpmi++;
           }
          }
         }

All arrays are nmol*natom*(x,y,z) arrays. stor is ndt*nmol*natom*(x,y,z),
where ndt is the number of points in the msqd.
The procedure of storing and correlating is that of Allen & Tildesley
book. In the Frenkel's book (Algorithm 8, the code at page 82) stor
corresponds to x0(i,t).
It can be adapted to an arbitrary cell, applying to dri the usual
cartesian-to-cell transformation. If co is the cell written in the xsf
snapshot
(PRIMVEC entry):

/* calculation of cell_to_cart (see wikipedia: Fractional coordinates) */
   a=sqrt(co[0]*co[0]+co[1]*co[1]+co[2]*co[2]);
   b=sqrt(co[3]*co[3]+co[4]*co[4]+co[5]*co[5]);
   c=sqrt(co[6]*co[6]+co[7]*co[7]+co[8]*co[8]);
   cos_alpha = (co[3]*co[6]+co[4]*co[7]+co[5]*co[8])/(b*c);
   cos_beta  = (co[0]*co[6]+co[1]*co[7]+co[2]*co[8])/(a*c);
   cos_gamma = (co[0]*co[3]+co[1]*co[4]+co[2]*co[5])/(a*b);
   sin_gamma = sqrt(1.-cos_gamma*cos_gamma);
   v = sqrt(1.-cos_gamma*cos_gamma
                -cos_alpha*cos_alpha
                -cos_beta*cos_beta
                +2.*cos_alpha*cos_beta*cos_gamma );
   cell_to_cart[1][1] = a;
   cell_to_cart[1][2] = 0.;
   cell_to_cart[1][3] = 0.;
   cell_to_cart[2][1] = b*cos_gamma;
   cell_to_cart[2][2] = b*sin_gamma;
   cell_to_cart[2][3] = 0.;
   cell_to_cart[3][1] = c*cos_beta;
   cell_to_cart[3][2] = c*(cos_alpha - cos_beta*cos_gamma)/sin_gamma;
   cell_to_cart[3][3] = c*v/sin_gamma;
   invert(cell_to_cart,3,3,cart_to_cell);

Once in fractional coordinates, the test on dri is made using 1/2.
I never had non-orthorhombic cells, I work with solution, so I did not
check the second part.
There are special cases where the test above is not safe, but they
should be very rare.

               Giovanni
============================================================
Giovanni La Penna - National research council (Cnr)
Institute for chemistry of organo-metallic compounds (Iccom)
via Madonna del Piano 10,
I-50019 Sesto Fiorentino, Firenze, Italy
tel.: +39 055 522-5264, fax: +39 055 522-5203
e-mail: glapenna at iccom.cnr.it - http://www.iccom.cnr.it/lapenna
skype: giovannilapenna
============================================================
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://www.democritos.it/pipermail/pw_forum/attachments/20110118/2a902753/attachment-0001.htm