src/gromacs/gmxana/gmx_order.c

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