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