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