src/tools/gmx_order.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <math.h>
  40 #include <ctype.h>
  41
  42 #include "sysstuff.h"
  43 #include "string.h"
  44 #include "typedefs.h"
  45 #include "smalloc.h"
  46 #include "macros.h"
  47 #include "gstat.h"
  48 #include "vec.h"
  49 #include "xvgr.h"
  50 #include "pbc.h"
  51 #include "copyrite.h"
  52 #include "futil.h"
  53 #include "statutil.h"
  54 #include "index.h"
  55 #include "tpxio.h"
  56 #include "confio.h"
  57 #include "cmat.h"
  58 #include "gmx_ana.h"
  59
  60
  61 /****************************************************************************/
  62 /* This program calculates the order parameter per atom for an interface or */
  63 /* bilayer, averaged over time.                                             */
  64 /* S = 1/2 * (3 * cos(i)cos(j) - delta(ij))                                 */
  65 /* It is assumed that the order parameter with respect to a box-axis        */
  66 /* is calculated. In that case i = j = axis, and delta(ij) = 1.             */
  67 /*                                                                          */
  68 /* Peter Tieleman,  April 1995                                              */
  69 /* P.J. van Maaren, November 2005     Added tetrahedral stuff               */
  70 /****************************************************************************/
  71
  72 static void find_nearest_neighbours(t_topology top, int ePBC,
  73                                     int natoms, matrix box,
  74                                     rvec x[],int maxidx,atom_id index[],
  75                                     real time,
  76                                     real *sgmean, real *skmean,
  77                                     int nslice,int slice_dim,
  78                                     real sgslice[],real skslice[],
  79                                     gmx_rmpbc_t gpbc)
  80 {
  81   FILE    *fpoutdist;
  82   char    fnsgdist[32];
  83   int     ix,jx,nsgbin, *sgbin;
  84   int     i1,i2,i,ibin,j,k,l,n,*nn[4];
  85   rvec    dx,dx1,dx2,rj,rk,urk,urj;
  86   real    cost,cost2,*sgmol,*skmol,rmean,rmean2,r2,box2,*r_nn[4];
  87   t_pbc   pbc;
  88   t_mat   *dmat;
  89   t_dist  *d;
  90   int     m1,mm,sl_index;
  91   int    **nnb,*sl_count;
  92   real   onethird=1.0/3.0;
  93   /*  dmat = init_mat(maxidx, FALSE); */
  94   box2 = box[XX][XX] * box[XX][XX];
  95   snew(sl_count,nslice);
  96   for (i=0; (i<4); i++) {
  97     snew(r_nn[i],natoms);
  98     snew(nn[i],natoms);
  99
 100     for (j=0; (j<natoms); j++) {
 101       r_nn[i][j] = box2;
 102     }
 103   }
 104
 105   snew(sgmol,maxidx);
 106   snew(skmol,maxidx);
 107
 108   /* Must init pbc every step because of pressure coupling */
 109   set_pbc(&pbc,ePBC,box);
 110
 111   gmx_rmpbc(gpbc,natoms,box,x);
 112
 113   nsgbin = 1 + 1/0.0005;
 114   snew(sgbin,nsgbin);
 115
 116   *sgmean = 0.0;
 117   *skmean = 0.0;
 118   l=0;
 119   for (i=0; (i<maxidx); i++) { /* loop over index file */
 120     ix = index[i];
 121     for (j=0; (j<maxidx); j++) {
 122       if (i == j) continue;
 123
 124       jx = index[j];
 125
 126       pbc_dx(&pbc,x[ix],x[jx],dx);
 127       r2=iprod(dx,dx);
 128
 129       /* set_mat_entry(dmat,i,j,r2); */
 130
 131       /* determine the nearest neighbours */
 132       if (r2 < r_nn[0][i]) {
 133         r_nn[3][i] = r_nn[2][i]; nn[3][i] = nn[2][i];
 134         r_nn[2][i] = r_nn[1][i]; nn[2][i] = nn[1][i];
 135         r_nn[1][i] = r_nn[0][i]; nn[1][i] = nn[0][i];
 136         r_nn[0][i] = r2;         nn[0][i] = j;
 137       } else if (r2 < r_nn[1][i]) {
 138         r_nn[3][i] = r_nn[2][i]; nn[3][i] = nn[2][i];
 139         r_nn[2][i] = r_nn[1][i]; nn[2][i] = nn[1][i];
 140         r_nn[1][i] = r2;         nn[1][i] = j;
 141       } else if (r2 < r_nn[2][i]) {
 142         r_nn[3][i] = r_nn[2][i]; nn[3][i] = nn[2][i];
 143         r_nn[2][i] = r2;         nn[2][i] = j;
 144       } else if (r2 < r_nn[3][i]) {
 145         r_nn[3][i] = r2;         nn[3][i] = j;
 146       }
 147     }
 148
 149
 150     /* calculate mean distance between nearest neighbours */
 151     rmean = 0;
 152     for (j=0; (j<4); j++) {
 153       r_nn[j][i] = sqrt(r_nn[j][i]);
 154       rmean += r_nn[j][i];
 155     }
 156     rmean /= 4;
 157
 158     n = 0;
 159     sgmol[i] = 0.0;
 160     skmol[i] = 0.0;
 161
 162     /* Chau1998a eqn 3 */
 163     /* angular part tetrahedrality order parameter per atom */
 164     for (j=0; (j<3); j++) {
 165       for (k=j+1; (k<4); k++) {
 166         pbc_dx(&pbc,x[ix],x[index[nn[k][i]]],rk);
 167         pbc_dx(&pbc,x[ix],x[index[nn[j][i]]],rj);
 168
 169         unitv(rk,urk);
 170         unitv(rj,urj);
 171
 172         cost = iprod(urk,urj) + onethird;
 173         cost2 = cost * cost;
 174
 175         /* sgmol[i] += 3*cost2/32;  */
 176            sgmol[i] += cost2;
 177
 178         /* determine distribution */
 179         ibin = nsgbin * cost2;
 180         if (ibin < nsgbin)
 181           sgbin[ibin]++;
 182         /* printf("%d %d %f %d %d\n", j, k, cost * cost, ibin, sgbin[ibin]);*/
 183         l++;
 184         n++;
 185       }
 186     }
 187
 188     /* normalize sgmol between 0.0 and 1.0 */
 189     sgmol[i] = 3*sgmol[i]/32;
 190     *sgmean += sgmol[i];
 191
 192     /* distance part tetrahedrality order parameter per atom */
 193     rmean2 = 4 * 3 * rmean * rmean;
 194     for (j=0; (j<4); j++) {
 195       skmol[i] += (rmean - r_nn[j][i]) * (rmean - r_nn[j][i]) / rmean2;
 196       /*      printf("%d %f (%f %f %f %f) \n",
 197               i, skmol[i], rmean, rmean2, r_nn[j][i], (rmean - r_nn[j][i]) );
 198       */
 199     }
 200
 201     *skmean += skmol[i];
 202
 203     /* Compute sliced stuff */
 204     sl_index = gmx_nint((1+x[i][slice_dim]/box[slice_dim][slice_dim])*nslice) % nslice;
 205     sgslice[sl_index] += sgmol[i];
 206     skslice[sl_index] += skmol[i];
 207     sl_count[sl_index]++;
 208   } /* loop over entries in index file */
 209
 210   *sgmean /= maxidx;
 211   *skmean /= maxidx;
 212
 213   for(i=0; (i<nslice); i++) {
 214     if (sl_count[i] > 0) {
 215       sgslice[i] /= sl_count[i];
 216       skslice[i] /= sl_count[i];
 217     }
 218   }
 219   sfree(sl_count);
 220   sfree(sgbin);
 221   sfree(sgmol);
 222   sfree(skmol);
 223   for (i=0; (i<4); i++) {
 224     sfree(r_nn[i]);
 225     sfree(nn[i]);
 226   }
 227 }
 228
 229
 230 static void calc_tetra_order_parm(const char *fnNDX,const char *fnTPS,
 231                                   const char *fnTRX, const char *sgfn,
 232                                   const char *skfn,
 233                                   int nslice,int slice_dim,
 234                                   const char *sgslfn,const char *skslfn,
 235                                   const output_env_t oenv)
 236 {
 237   FILE       *fpsg=NULL,*fpsk=NULL;
 238   t_topology top;
 239   int        ePBC;
 240   char       title[STRLEN],fn[STRLEN],subtitle[STRLEN];
 241   t_trxstatus *status;
 242   int        natoms;
 243   real       t;
 244   rvec       *xtop,*x;
 245   matrix     box;
 246   real       sg,sk;
 247   atom_id    **index;
 248   char       **grpname;
 249   int        i,*isize,ng,nframes;
 250   real       *sg_slice,*sg_slice_tot,*sk_slice,*sk_slice_tot;
 251   gmx_rmpbc_t  gpbc=NULL;
 252
 253
 254   read_tps_conf(fnTPS,title,&top,&ePBC,&xtop,NULL,box,FALSE);
 255
 256   snew(sg_slice,nslice);
 257   snew(sk_slice,nslice);
 258   snew(sg_slice_tot,nslice);
 259   snew(sk_slice_tot,nslice);
 260   ng = 1;
 261   /* get index groups */
 262   printf("Select the group that contains the atoms you want to use for the tetrahedrality order parameter calculation:\n");
 263   snew(grpname,ng);
 264   snew(index,ng);
 265   snew(isize,ng);
 266   get_index(&top.atoms,fnNDX,ng,isize,index,grpname);
 267
 268   /* Analyze trajectory */
 269   natoms=read_first_x(oenv,&status,fnTRX,&t,&x,box);
 270   if ( natoms > top.atoms.nr )
 271     gmx_fatal(FARGS,"Topology (%d atoms) does not match trajectory (%d atoms)",
 272               top.atoms.nr,natoms);
 273   check_index(NULL,ng,index[0],NULL,natoms);
 274
 275   fpsg=xvgropen(sgfn,"S\\sg\\N Angle Order Parameter","Time (ps)","S\\sg\\N",
 276                 oenv);
 277   fpsk=xvgropen(skfn,"S\\sk\\N Distance Order Parameter","Time (ps)","S\\sk\\N",
 278                 oenv);
 279
 280   /* loop over frames */
 281   gpbc = gmx_rmpbc_init(&top.idef,ePBC,natoms,box);
 282   nframes = 0;
 283   do {
 284     find_nearest_neighbours(top,ePBC,natoms,box,x,isize[0],index[0],t,
 285                             &sg,&sk,nslice,slice_dim,sg_slice,sk_slice,gpbc);
 286     for(i=0; (i<nslice); i++) {
 287       sg_slice_tot[i] += sg_slice[i];
 288       sk_slice_tot[i] += sk_slice[i];
 289     }
 290     fprintf(fpsg,"%f %f\n", t, sg);
 291     fprintf(fpsk,"%f %f\n", t, sk);
 292     nframes++;
 293   } while (read_next_x(oenv,status,&t,natoms,x,box));
 294   close_trj(status);
 295   gmx_rmpbc_done(gpbc);
 296
 297   sfree(grpname);
 298   sfree(index);
 299   sfree(isize);
 300
 301   ffclose(fpsg);
 302   ffclose(fpsk);
 303
 304   fpsg = xvgropen(sgslfn,
 305                   "S\\sg\\N Angle Order Parameter / Slab","(nm)","S\\sg\\N",
 306                    oenv);
 307   fpsk = xvgropen(skslfn,
 308                   "S\\sk\\N Distance Order Parameter / Slab","(nm)","S\\sk\\N",
 309                   oenv);
 310   for(i=0; (i<nslice); i++) {
 311     fprintf(fpsg,"%10g  %10g\n",(i+0.5)*box[slice_dim][slice_dim]/nslice,
 312             sg_slice_tot[i]/nframes);
 313     fprintf(fpsk,"%10g  %10g\n",(i+0.5)*box[slice_dim][slice_dim]/nslice,
 314             sk_slice_tot[i]/nframes);
 315   }
 316   ffclose(fpsg);
 317   ffclose(fpsk);
 318 }
 319
 320
 321 /* Print name of first atom in all groups in index file */
 322 static void print_types(atom_id index[], atom_id a[], int ngrps,
 323                         char *groups[], t_topology *top)
 324 {
 325   int i;
 326
 327   fprintf(stderr,"Using following groups: \n");
 328   for(i = 0; i < ngrps; i++)
 329     fprintf(stderr,"Groupname: %s First atomname: %s First atomnr %u\n",
 330             groups[i], *(top->atoms.atomname[a[index[i]]]), a[index[i]]);
 331   fprintf(stderr,"\n");
 332 }
 333
 334 static void check_length(real length, int a, int b)
 335 {
 336   if (length > 0.3)
 337     fprintf(stderr,"WARNING: distance between atoms %d and "
 338             "%d > 0.3 nm (%f). Index file might be corrupt.\n",
 339             a, b, length);
 340 }
 341
 342 void calc_order(const char *fn, atom_id *index, atom_id *a, rvec **order,
 343                 real ***slOrder, real *slWidth, int nslices, gmx_bool bSliced,
 344                 gmx_bool bUnsat, t_topology *top, int ePBC, int ngrps, int axis,
 345                 gmx_bool permolecule, gmx_bool radial, gmx_bool distcalc, const char *radfn,
 346                 real ***distvals,
 347                 const output_env_t oenv)
 348 {
 349   /* if permolecule = TRUE, order parameters will be calculed per molecule
 350    * and stored in slOrder with #slices = # molecules */
 351   rvec *x0,          /* coordinates with pbc                           */
 352     *x1,             /* coordinates without pbc                        */
 353     dist;            /* vector between two atoms                       */
 354   matrix box;        /* box (3x3)                                      */
 355   t_trxstatus *status;
 356   rvec  cossum,      /* sum of vector angles for three axes            */
 357     Sx, Sy, Sz,      /* the three molecular axes                       */
 358     tmp1, tmp2,      /* temp. rvecs for calculating dot products       */
 359     frameorder;      /* order parameters for one frame                 */
 360   real *slFrameorder; /* order parameter for one frame, per slice      */
 361   real length,       /* total distance between two atoms               */
 362     t,               /* time from trajectory                           */
 363     z_ave,z1,z2;     /* average z, used to det. which slice atom is in */
 364   int natoms,        /* nr. atoms in trj                               */
 365     nr_tails,        /* nr tails, to check if index file is correct    */
 366     size=0,          /* nr. of atoms in group. same as nr_tails        */
 367     i,j,m,k,l,teller = 0,
 368     slice,           /* current slice number                           */
 369     nr_frames = 0,
 370     *slCount;        /* nr. of atoms in one slice                      */
 371    real dbangle = 0, /* angle between double bond and  axis            */
 372         sdbangle = 0;/* sum of these angles                            */
 373    gmx_bool use_unitvector = FALSE; /* use a specified unit vector instead of axis to specify unit normal*/
 374    rvec direction, com,dref,dvec;
 375    int comsize, distsize;
 376    atom_id *comidx=NULL, *distidx=NULL;
 377    char *grpname=NULL;
 378    t_pbc pbc;
 379    real arcdist;
 380    gmx_rmpbc_t  gpbc=NULL;
 381
 382   /* PBC added for center-of-mass vector*/
 383   /* Initiate the pbc structure */
 384   memset(&pbc,0,sizeof(pbc));
 385
 386   if ((natoms = read_first_x(oenv,&status,fn,&t,&x0,box)) == 0)
 387     gmx_fatal(FARGS,"Could not read coordinates from statusfile\n");
 388
 389   nr_tails = index[1] - index[0];
 390   fprintf(stderr,"Number of elements in first group: %d\n",nr_tails);
 391   /* take first group as standard. Not rocksolid, but might catch error in index*/
 392
 393   if (permolecule)
 394   {
 395           nslices=nr_tails;
 396           bSliced=FALSE;  /*force slices off */
 397       fprintf(stderr,"Calculating order parameters for each of %d molecules\n",
 398             nslices);
 399   }
 400
 401   if (radial)
 402   {
 403         use_unitvector=TRUE;
 404         fprintf(stderr,"Select an index group to calculate the radial membrane normal\n");
 405         get_index(&top->atoms,radfn,1,&comsize,&comidx,&grpname);
 406         if (distcalc)
 407         {
 408                 if (grpname!=NULL)
 409                         sfree(grpname);
 410                 fprintf(stderr,"Select an index group to use as distance reference\n");
 411                 get_index(&top->atoms,radfn,1,&distsize,&distidx,&grpname);
 412                 bSliced=FALSE; /*force slices off*/
 413         }
 414   }
 415
 416   if (use_unitvector && bSliced)
 417         fprintf(stderr,"Warning:  slicing and specified unit vectors are not currently compatible\n");
 418
 419   snew(slCount, nslices);
 420   snew(*slOrder, nslices);
 421   for(i = 0; i < nslices; i++)
 422     snew((*slOrder)[i],ngrps);
 423   if (distcalc)
 424   {
 425     snew(*distvals, nslices);
 426     for(i = 0; i < nslices; i++)
 427       snew((*distvals)[i],ngrps);
 428   }
 429   snew(*order,ngrps);
 430   snew(slFrameorder, nslices);
 431   snew(x1, natoms);
 432
 433   if (bSliced) {
 434     *slWidth = box[axis][axis]/nslices;
 435     fprintf(stderr,"Box divided in %d slices. Initial width of slice: %f\n",
 436             nslices, *slWidth);
 437   }
 438
 439 #if 0
 440   nr_tails = index[1] - index[0];
 441   fprintf(stderr,"Number of elements in first group: %d\n",nr_tails);
 442   /* take first group as standard. Not rocksolid, but might catch error
 443      in index*/
 444 #endif
 445
 446   teller = 0;
 447
 448   gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box);
 449   /*********** Start processing trajectory ***********/
 450   do {
 451     if (bSliced)
 452       *slWidth = box[axis][axis]/nslices;
 453     teller++;
 454
 455     set_pbc(&pbc,ePBC,box);
 456     gmx_rmpbc_copy(gpbc,natoms,box,x0,x1);
 457
 458     /* Now loop over all groups. There are ngrps groups, the order parameter can
 459        be calculated for grp 1 to grp ngrps - 1. For each group, loop over all
 460        atoms in group, which is index[i] to (index[i+1] - 1) See block.h. Of
 461        course, in this case index[i+1] -index[i] has to be the same for all
 462        groups, namely the number of tails. i just runs over all atoms in a tail,
 463        so for DPPC ngrps = 16 and i runs from 1 to 14, including 14
 464      */
 465
 466           if (radial)
 467           {
 468                 /*center-of-mass determination*/
 469                 com[XX]=0.0; com[YY]=0.0; com[ZZ]=0.0;
 470                 for (j=0;j<comsize;j++)
 471                         rvec_inc(com,x1[comidx[j]]);
 472                 svmul(1.0/comsize,com,com);
 473           }
 474           if (distcalc)
 475           {
 476                 dref[XX]=0.0; dref[YY]=0.0;dref[ZZ]=0.0;
 477                 for (j=0;j<distsize;j++)
 478                         rvec_inc(dist,x1[distidx[j]]);
 479                 svmul(1.0/distsize,dref,dref);
 480                 pbc_dx(&pbc,dref,com,dvec);
 481                 unitv(dvec,dvec);
 482           }
 483
 484     for (i = 1; i < ngrps - 1; i++) {
 485       clear_rvec(frameorder);
 486
 487       size = index[i+1] - index[i];
 488       if (size != nr_tails)
 489         gmx_fatal(FARGS,"grp %d does not have same number of"
 490                 " elements as grp 1\n",i);
 491
 492       for (j = 0; j < size; j++) {
 493           if (radial)
 494           /*create unit vector*/
 495           {
 496                 pbc_dx(&pbc,x1[a[index[i]+j]],com,direction);
 497                 unitv(direction,direction);
 498                 /*DEBUG*/
 499                 /*if (j==0)
 500                         fprintf(stderr,"X %f %f %f\tcom %f %f %f\tdirection %f %f %f\n",x1[a[index[i]+j]][0],x1[a[index[i]+j]][1],x1[a[index[i]+j]][2],com[0],com[1],com[2],
 501                                 direction[0],direction[1],direction[2]);*/
 502           }
 503
 504         if (bUnsat) {
 505           /* Using convention for unsaturated carbons */
 506           /* first get Sz, the vector from Cn to Cn+1 */
 507           rvec_sub(x1[a[index[i+1]+j]], x1[a[index[i]+j]], dist);
 508           length = norm(dist);
 509           check_length(length, a[index[i]+j], a[index[i+1]+j]);
 510           svmul(1/length, dist, Sz);
 511
 512           /* this is actually the cosine of the angle between the double bond
 513              and axis, because Sz is normalized and the two other components of
 514              the axis on the bilayer are zero */
 515           if (use_unitvector)
 516           {
 517             sdbangle += acos(iprod(direction,Sz));  /*this can probably be optimized*/
 518           }
 519           else
 520             sdbangle += acos(Sz[axis]);
 521         } else {
 522           /* get vector dist(Cn-1,Cn+1) for tail atoms */
 523           rvec_sub(x1[a[index[i+1]+j]], x1[a[index[i-1]+j]], dist);
 524           length = norm(dist);      /* determine distance between two atoms */
 525           check_length(length, a[index[i-1]+j], a[index[i+1]+j]);
 526
 527           svmul(1/length, dist, Sz);
 528           /* Sz is now the molecular axis Sz, normalized and all that */
 529         }
 530
 531         /* now get Sx. Sx is normal to the plane of Cn-1, Cn and Cn+1 so
 532            we can use the outer product of Cn-1->Cn and Cn+1->Cn, I hope */
 533         rvec_sub(x1[a[index[i+1]+j]], x1[a[index[i]+j]], tmp1);
 534         rvec_sub(x1[a[index[i-1]+j]], x1[a[index[i]+j]], tmp2);
 535         cprod(tmp1, tmp2, Sx);
 536         svmul(1/norm(Sx), Sx, Sx);
 537
 538         /* now we can get Sy from the outer product of Sx and Sz   */
 539         cprod(Sz, Sx, Sy);
 540         svmul(1/norm(Sy), Sy, Sy);
 541
 542         /* the square of cosine of the angle between dist and the axis.
 543            Using the innerproduct, but two of the three elements are zero
 544            Determine the sum of the orderparameter of all atoms in group
 545            */
 546         if (use_unitvector)
 547         {
 548         cossum[XX] = sqr(iprod(Sx,direction)); /* this is allowed, since Sa is normalized */
 549         cossum[YY] = sqr(iprod(Sy,direction));
 550         cossum[ZZ] = sqr(iprod(Sz,direction));
 551         }
 552         else
 553         {
 554         cossum[XX] = sqr(Sx[axis]); /* this is allowed, since Sa is normalized */
 555         cossum[YY] = sqr(Sy[axis]);
 556         cossum[ZZ] = sqr(Sz[axis]);
 557     }
 558
 559         for (m = 0; m < DIM; m++)
 560           frameorder[m] += 0.5 * (3 * cossum[m] - 1);
 561
 562         if (bSliced) {
 563           /* get average coordinate in box length for slicing,
 564              determine which slice atom is in, increase count for that
 565              slice. slFrameorder and slOrder are reals, not
 566              rvecs. Only the component [axis] of the order tensor is
 567              kept, until I find it necessary to know the others too
 568            */
 569
 570           z1 = x1[a[index[i-1]+j]][axis];
 571           z2 = x1[a[index[i+1]+j]][axis];
 572           z_ave = 0.5 * (z1 + z2);
 573           if (z_ave < 0)
 574             z_ave += box[axis][axis];
 575           if (z_ave > box[axis][axis])
 576             z_ave -= box[axis][axis];
 577
 578           slice  = (int)(0.5 + (z_ave / (*slWidth))) - 1;
 579           slCount[slice]++;               /* determine slice, increase count */
 580
 581           slFrameorder[slice] += 0.5 * (3 * cossum[axis] - 1);
 582         }
 583         else if (permolecule)
 584         {
 585                 /*  store per-molecule order parameter
 586                  *  To just track single-axis order: (*slOrder)[j][i] += 0.5 * (3 * iprod(cossum,direction) - 1);
 587                  *  following is for Scd order: */
 588                  (*slOrder)[j][i] += -1* (0.3333 * (3 * cossum[XX] - 1) + 0.3333 * 0.5 * (3 * cossum[YY] - 1));
 589         }
 590         if (distcalc)
 591         {
 592                 /* bin order parameter by arc distance from reference group*/
 593                 arcdist=acos(iprod(dvec,direction));
 594                 (*distvals)[j][i]+=arcdist;
 595         }
 596       }   /* end loop j, over all atoms in group */
 597
 598       for (m = 0; m < DIM; m++)
 599         (*order)[i][m] += (frameorder[m]/size);
 600
 601           if (!permolecule)
 602           {  /*Skip following if doing per-molecule*/
 603          for (k = 0; k < nslices; k++) {
 604                if (slCount[k]) {     /* if no elements, nothing has to be added */
 605                   (*slOrder)[k][i] += slFrameorder[k]/slCount[k];
 606                   slFrameorder[k] = 0; slCount[k] = 0;
 607                }
 608           }
 609       }   /* end loop i, over all groups in indexfile */
 610     }
 611     nr_frames++;
 612
 613   } while (read_next_x(oenv,status,&t,natoms,x0,box));
 614   /*********** done with status file **********/
 615
 616   fprintf(stderr,"\nRead trajectory. Printing parameters to file\n");
 617   gmx_rmpbc_done(gpbc);
 618
 619   /* average over frames */
 620   for (i = 1; i < ngrps - 1; i++) {
 621     svmul(1.0/nr_frames, (*order)[i], (*order)[i]);
 622     fprintf(stderr,"Atom %d Tensor: x=%g , y=%g, z=%g\n",i,(*order)[i][XX],
 623             (*order)[i][YY], (*order)[i][ZZ]);
 624     if (bSliced || permolecule) {
 625       for (k = 0; k < nslices; k++)
 626         (*slOrder)[k][i] /= nr_frames;
 627     }
 628         if (distcalc)
 629       for (k = 0; k < nslices; k++)
 630                 (*distvals)[k][i] /= nr_frames;
 631   }
 632
 633   if (bUnsat)
 634     fprintf(stderr,"Average angle between double bond and normal: %f\n",
 635             180*sdbangle/(nr_frames * size*M_PI));
 636
 637   sfree(x0);  /* free memory used by coordinate arrays */
 638   sfree(x1);
 639   if (comidx!=NULL)
 640         sfree(comidx);
 641   if (distidx!=NULL)
 642         sfree(distidx);
 643   if (grpname!=NULL)
 644     sfree(grpname);
 645 }
 646
 647
 648 void order_plot(rvec order[], real *slOrder[], const char *afile, const char *bfile,
 649                 const char *cfile, int ngrps, int nslices, real slWidth, gmx_bool bSzonly,
 650                 gmx_bool permolecule, real **distvals, const output_env_t oenv)
 651 {
 652   FILE       *ord, *slOrd;        /* xvgr files with order parameters  */
 653   int        atom, slice;         /* atom corresponding to order para.*/
 654   char       buf[256];            /* for xvgr title */
 655   real      S;                    /* order parameter averaged over all atoms */
 656
 657   if (permolecule)
 658   {
 659     sprintf(buf,"Scd order parameters");
 660     ord = xvgropen(afile,buf,"Atom","S",oenv);
 661     sprintf(buf, "Orderparameters per atom per slice");
 662     slOrd = xvgropen(bfile, buf, "Molecule", "S",oenv);
 663     for (atom = 1; atom < ngrps - 1; atom++) {
 664       fprintf(ord,"%12d   %12g\n", atom, -1 * (0.6667 * order[atom][XX] +
 665                                                  0.333 * order[atom][YY]));
 666     }
 667
 668     for (slice = 0; slice < nslices; slice++) {
 669           fprintf(slOrd,"%12d\t",slice);
 670           if (distvals)
 671                 fprintf(slOrd,"%12g\t", distvals[slice][1]); /*use distance value at second carbon*/
 672       for (atom = 1; atom < ngrps - 1; atom++)
 673             fprintf(slOrd,"%12g\t", slOrder[slice][atom]);
 674           fprintf(slOrd,"\n");
 675     }
 676
 677   }
 678   else if (bSzonly) {
 679     sprintf(buf,"Orderparameters Sz per atom");
 680     ord = xvgropen(afile,buf,"Atom","S",oenv);
 681     fprintf(stderr,"ngrps = %d, nslices = %d",ngrps, nslices);
 682
 683     sprintf(buf, "Orderparameters per atom per slice");
 684     slOrd = xvgropen(bfile, buf, "Slice", "S",oenv);
 685
 686     for (atom = 1; atom < ngrps - 1; atom++)
 687       fprintf(ord,"%12d       %12g\n", atom, order[atom][ZZ]);
 688
 689     for (slice = 0; slice < nslices; slice++) {
 690       S = 0;
 691       for (atom = 1; atom < ngrps - 1; atom++)
 692         S += slOrder[slice][atom];
 693       fprintf(slOrd,"%12g     %12g\n", slice*slWidth, S/atom);
 694     }
 695
 696   } else {
 697     sprintf(buf,"Order tensor diagonal elements");
 698     ord = xvgropen(afile,buf,"Atom","S",oenv);
 699     sprintf(buf,"Deuterium order parameters");
 700     slOrd = xvgropen(cfile,buf, "Atom", "Scd",oenv);
 701
 702     for (atom = 1; atom < ngrps - 1; atom++) {
 703       fprintf(ord,"%12d   %12g   %12g   %12g\n", atom, order[atom][XX],
 704               order[atom][YY], order[atom][ZZ]);
 705       fprintf(slOrd,"%12d   %12g\n", atom, -1 * (0.6667 * order[atom][XX] +
 706                                                  0.333 * order[atom][YY]));
 707     }
 708
 709     ffclose(ord);
 710     ffclose(slOrd);
 711   }
 712 }
 713
 714 void write_bfactors(t_filenm  *fnm, int nfile, atom_id *index, atom_id *a, int nslices, int ngrps, real **order, t_topology *top, real **distvals,output_env_t oenv)
 715 {
 716         /*function to write order parameters as B factors in PDB file using
 717           first frame of trajectory*/
 718         t_trxstatus *status;
 719         int natoms;
 720         t_trxframe fr, frout;
 721         t_atoms useatoms;
 722         int i,j,ctr,nout;
 723
 724         ngrps-=2;  /*we don't have an order parameter for the first or
 725                      last atom in each chain*/
 726         nout=nslices*ngrps;
 727         natoms=read_first_frame(oenv,&status,ftp2fn(efTRX,nfile,fnm),&fr,
 728                                 TRX_NEED_X);
 729         close_trj(status);
 730         frout = fr;
 731         frout.natoms=nout;
 732         frout.bF=FALSE;
 733         frout.bV=FALSE;
 734         frout.x=0;
 735         snew(frout.x,nout);
 736
 737         init_t_atoms(&useatoms,nout,TRUE);
 738         useatoms.nr=nout;
 739
 740         /*initialize PDBinfo*/
 741         for (i=0;i<useatoms.nr;++i)
 742         {
 743                 useatoms.pdbinfo[i].type=0;
 744                 useatoms.pdbinfo[i].occup=0.0;
 745                 useatoms.pdbinfo[i].bfac=0.0;
 746                 useatoms.pdbinfo[i].bAnisotropic=FALSE;
 747         }
 748
 749         for (j=0,ctr=0;j<nslices;j++)
 750                 for (i=0;i<ngrps;i++,ctr++)
 751                 {
 752                         /*iterate along each chain*/
 753                         useatoms.pdbinfo[ctr].bfac=order[j][i+1];
 754                         if (distvals)
 755                                 useatoms.pdbinfo[ctr].occup=distvals[j][i+1];
 756                         copy_rvec(fr.x[a[index[i+1]+j]],frout.x[ctr]);
 757                         useatoms.atomname[ctr]=top->atoms.atomname[a[index[i+1]+j]];
 758                         useatoms.atom[ctr]=top->atoms.atom[a[index[i+1]+j]];
 759                         useatoms.nres=max(useatoms.nres,useatoms.atom[ctr].resind+1);
 760                         useatoms.resinfo[useatoms.atom[ctr].resind]=top->atoms.resinfo[useatoms.atom[ctr].resind]; /*copy resinfo*/
 761                 }
 762
 763         write_sto_conf(opt2fn("-ob",nfile,fnm),"Order parameters",&useatoms,frout.x,NULL,frout.ePBC,frout.box);
 764
 765         sfree(frout.x);
 766         free_t_atoms(&useatoms,FALSE);
 767 }
 768
 769 int gmx_order(int argc,char *argv[])
 770 {
 771   const char *desc[] = {
 772     "Compute the order parameter per atom for carbon tails. For atom i the",
 773     "vector i-1, i+1 is used together with an axis. ",
 774     "The index file should contain only the groups to be used for calculations,",
 775     "with each group of equivalent carbons along the relevant acyl chain in its own",
 776     "group. There should not be any generic groups (like System, Protein) in the index",
 777     "file to avoid confusing the program (this is not relevant to tetrahedral order",
 778     "parameters however, which only work for water anyway).[PAR]",
 779     "The program can also give all",
 780     "diagonal elements of the order tensor and even calculate the deuterium",
 781     "order parameter Scd (default). If the option [TT]-szonly[tt] is given, only one",
 782     "order tensor component (specified by the [TT]-d[tt] option) is given and the",
 783     "order parameter per slice is calculated as well. If [TT]-szonly[tt] is not",
 784     "selected, all diagonal elements and the deuterium order parameter is",
 785     "given.[PAR]"
 786     "The tetrahedrality order parameters can be determined",
 787     "around an atom. Both angle an distance order parameters are calculated. See",
 788     "P.-L. Chau and A.J. Hardwick, Mol. Phys., 93, (1998), 511-518.",
 789     "for more details.[BR]",
 790     ""
 791   };
 792
 793   static int  nslices = 1;                    /* nr of slices defined       */
 794   static gmx_bool bSzonly = FALSE;                /* True if only Sz is wanted  */
 795   static gmx_bool bUnsat = FALSE;                 /* True if carbons are unsat. */
 796   static const char *normal_axis[] = { NULL, "z", "x", "y", NULL };
 797   static gmx_bool permolecule = FALSE;  /*compute on a per-molecule basis */
 798   static gmx_bool radial = FALSE; /*compute a radial membrane normal */
 799   static gmx_bool distcalc = FALSE; /*calculate distance from a reference group */
 800   t_pargs pa[] = {
 801     { "-d",      FALSE, etENUM, {normal_axis},
 802       "Direction of the normal on the membrane" },
 803     { "-sl",     FALSE, etINT, {&nslices},
 804       "Calculate order parameter as function of boxlength, dividing the box"
 805       " in #nr slices." },
 806     { "-szonly", FALSE, etBOOL,{&bSzonly},
 807       "Only give Sz element of order tensor. (axis can be specified with -d)" },
 808     { "-unsat",  FALSE, etBOOL,{&bUnsat},
 809       "Calculate order parameters for unsaturated carbons. Note that this can"
 810       "not be mixed with normal order parameters." },
 811         { "-permolecule", FALSE, etBOOL,{&permolecule},
 812       "Compute per-molecule Scd order parameters" },
 813         { "-radial", FALSE, etBOOL,{&radial},
 814       "Compute a radial membrane normal" },
 815         { "-calcdist", FALSE, etBOOL,{&distcalc},
 816       "Compute distance from a reference (currently defined only for radial and permolecule)" },
 817   };
 818
 819   rvec      *order;                         /* order par. for each atom   */
 820   real      **slOrder;                      /* same, per slice            */
 821   real      slWidth = 0.0;                  /* width of a slice           */
 822   char      **grpname;                      /* groupnames                 */
 823   int       ngrps,                          /* nr. of groups              */
 824             i,
 825             axis=0;                         /* normal axis                */
 826   t_topology *top;                          /* topology                   */
 827   int       ePBC;
 828   atom_id   *index,                         /* indices for a              */
 829             *a;                             /* atom numbers in each group */
 830   t_blocka  *block;                         /* data from index file       */
 831   t_filenm  fnm[] = {                       /* files for g_order          */
 832     { efTRX, "-f", NULL,  ffREAD },         /* trajectory file            */
 833     { efNDX, "-n", NULL,  ffREAD },         /* index file                 */
 834     { efNDX, "-nr", NULL,  ffREAD },            /* index for radial axis calculation      */
 835     { efTPX, NULL, NULL,  ffREAD },         /* topology file              */
 836     { efXVG,"-o","order", ffWRITE },        /* xvgr output file           */
 837     { efXVG,"-od","deuter", ffWRITE },      /* xvgr output file           */
 838     { efPDB,"-ob",NULL, ffWRITE },          /* write Scd as B factors to PDB if permolecule           */
 839     { efXVG,"-os","sliced", ffWRITE },      /* xvgr output file           */
 840     { efXVG,"-Sg","sg-ang", ffOPTWR },      /* xvgr output file           */
 841     { efXVG,"-Sk","sk-dist", ffOPTWR },     /* xvgr output file           */
 842     { efXVG,"-Sgsl","sg-ang-slice", ffOPTWR },      /* xvgr output file           */
 843     { efXVG,"-Sksl","sk-dist-slice", ffOPTWR },     /* xvgr output file           */
 844   };
 845   gmx_bool      bSliced = FALSE;                /* True if box is sliced      */
 846 #define NFILE asize(fnm)
 847   real **distvals=NULL;
 848   const char *sgfnm,*skfnm,*ndxfnm,*tpsfnm,*trxfnm;
 849   output_env_t oenv;
 850
 851   CopyRight(stderr,argv[0]);
 852
 853   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
 854                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0, NULL,&oenv);
 855   if (nslices < 1)
 856     gmx_fatal(FARGS,"Can not have nslices < 1");
 857   sgfnm = opt2fn_null("-Sg",NFILE,fnm);
 858   skfnm = opt2fn_null("-Sk",NFILE,fnm);
 859   ndxfnm = opt2fn_null("-n",NFILE,fnm);
 860   tpsfnm = ftp2fn(efTPX,NFILE,fnm);
 861   trxfnm = ftp2fn(efTRX,NFILE,fnm);
 862
 863   /* Calculate axis */
 864   if (strcmp(normal_axis[0],"x") == 0) axis = XX;
 865   else if (strcmp(normal_axis[0],"y") == 0) axis = YY;
 866   else if (strcmp(normal_axis[0],"z") == 0) axis = ZZ;
 867   else gmx_fatal(FARGS,"Invalid axis, use x, y or z");
 868
 869   switch (axis) {
 870   case 0:
 871     fprintf(stderr,"Taking x axis as normal to the membrane\n");
 872     break;
 873   case 1:
 874     fprintf(stderr,"Taking y axis as normal to the membrane\n");
 875     break;
 876   case 2:
 877     fprintf(stderr,"Taking z axis as normal to the membrane\n");
 878     break;
 879   }
 880
 881   /* tetraheder order parameter */
 882   if (skfnm || sgfnm) {
 883     /* If either of theoptions is set we compute both */
 884     sgfnm = opt2fn("-Sg",NFILE,fnm);
 885     skfnm = opt2fn("-Sk",NFILE,fnm);
 886     calc_tetra_order_parm(ndxfnm,tpsfnm,trxfnm,sgfnm,skfnm,nslices,axis,
 887                           opt2fn("-Sgsl",NFILE,fnm),opt2fn("-Sksl",NFILE,fnm),
 888                           oenv);
 889     /* view xvgr files */
 890     do_view(oenv,opt2fn("-Sg",NFILE,fnm), NULL);
 891     do_view(oenv,opt2fn("-Sk",NFILE,fnm), NULL);
 892     if (nslices > 1) {
 893       do_view(oenv,opt2fn("-Sgsl",NFILE,fnm), NULL);
 894       do_view(oenv,opt2fn("-Sksl",NFILE,fnm), NULL);
 895     }
 896   }
 897   else {
 898     /* tail order parameter */
 899
 900     if (nslices > 1) {
 901       bSliced = TRUE;
 902       fprintf(stderr,"Dividing box in %d slices.\n\n", nslices);
 903     }
 904
 905     if (bSzonly)
 906       fprintf(stderr,"Only calculating Sz\n");
 907     if (bUnsat)
 908       fprintf(stderr,"Taking carbons as unsaturated!\n");
 909
 910     top = read_top(ftp2fn(efTPX,NFILE,fnm),&ePBC);     /* read topology file */
 911
 912     block = init_index(ftp2fn(efNDX,NFILE,fnm),&grpname);
 913     index = block->index;                       /* get indices from t_block block */
 914     a = block->a;                               /* see block.h                    */
 915     ngrps = block->nr;
 916
 917   if (permolecule)
 918   {
 919     nslices = index[1] - index[0];  /*I think this assumes contiguous lipids in topology*/
 920           fprintf(stderr,"Calculating Scd order parameters for each of %d molecules\n",nslices);
 921   }
 922
 923   if (distcalc)
 924   {
 925           radial=TRUE;
 926           fprintf(stderr,"Calculating radial distances\n");
 927           if (!permolecule)
 928                 gmx_fatal(FARGS,"Cannot yet output radial distances without permolecule\n");
 929   }
 930
 931         /* show atomtypes, to check if index file is correct */
 932     print_types(index, a, ngrps, grpname, top);
 933
 934     calc_order(ftp2fn(efTRX,NFILE,fnm), index, a, &order,
 935                &slOrder, &slWidth, nslices, bSliced, bUnsat,
 936                top, ePBC, ngrps, axis,permolecule,radial,distcalc,opt2fn_null("-nr",NFILE,fnm),&distvals, oenv);
 937
 938         if (radial)
 939                 ngrps--; /*don't print the last group--was used for
 940                            center-of-mass determination*/
 941
 942     order_plot(order, slOrder, opt2fn("-o",NFILE,fnm), opt2fn("-os",NFILE,fnm),
 943                opt2fn("-od",NFILE,fnm), ngrps, nslices, slWidth, bSzonly,permolecule,distvals,oenv);
 944
 945         if (opt2bSet("-ob",NFILE,fnm))
 946         {
 947                 if (!permolecule)
 948                         fprintf(stderr,
 949                                 "Won't write B-factors with averaged order parameters; use -permolecule\n");
 950                 else
 951                         write_bfactors(fnm,NFILE,index,a,nslices,ngrps,slOrder,top,distvals,oenv);
 952         }
 953
 954
 955     do_view(oenv,opt2fn("-o",NFILE,fnm), NULL);      /* view xvgr file */
 956     do_view(oenv,opt2fn("-os",NFILE,fnm), NULL);     /* view xvgr file */
 957     do_view(oenv,opt2fn("-od",NFILE,fnm), NULL);     /* view xvgr file */
 958   }
 959
 960   thanx(stderr);
 961   if (distvals!=NULL)
 962   {
 963         for (i=0;i<nslices;++i)
 964                 sfree(distvals[i]);
 965         sfree(distvals);
 966   }
 967
 968   return 0;
 969 }