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