2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
44 #include "gromacs/legacyheaders/typedefs.h"
45 #include "gromacs/legacyheaders/macros.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"
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"
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. */
71 /* Peter Tieleman, April 1995 */
72 /* P.J. van Maaren, November 2005 Added tetrahedral stuff */
73 /****************************************************************************/
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[],
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];
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++)
100 snew(r_nn[i], natoms);
103 for (j = 0; (j < natoms); j++)
112 /* Must init pbc every step because of pressure coupling */
113 set_pbc(&pbc, ePBC, box);
115 gmx_rmpbc(gpbc, natoms, box, x);
117 nsgbin = 1 + 1/0.0005;
123 for (i = 0; (i < maxidx); i++) /* loop over index file */
126 for (j = 0; (j < maxidx); j++)
135 pbc_dx(&pbc, x[ix], x[jx], dx);
138 /* set_mat_entry(dmat,i,j,r2); */
140 /* determine the nearest neighbours */
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;
148 else if (r2 < r_nn[1][i])
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;
154 else if (r2 < r_nn[2][i])
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;
159 else if (r2 < r_nn[3][i])
161 r_nn[3][i] = r2; nn[3][i] = j;
166 /* calculate mean distance between nearest neighbours */
168 for (j = 0; (j < 4); j++)
170 r_nn[j][i] = sqrt(r_nn[j][i]);
179 /* Chau1998a eqn 3 */
180 /* angular part tetrahedrality order parameter per atom */
181 for (j = 0; (j < 3); j++)
183 for (k = j+1; (k < 4); k++)
185 pbc_dx(&pbc, x[ix], x[index[nn[k][i]]], rk);
186 pbc_dx(&pbc, x[ix], x[index[nn[j][i]]], rj);
191 cost = iprod(urk, urj) + onethird;
194 /* sgmol[i] += 3*cost2/32; */
197 /* determine distribution */
198 ibin = nsgbin * cost2;
203 /* printf("%d %d %f %d %d\n", j, k, cost * cost, ibin, sgbin[ibin]);*/
209 /* normalize sgmol between 0.0 and 1.0 */
210 sgmol[i] = 3*sgmol[i]/32;
213 /* distance part tetrahedrality order parameter per atom */
214 rmean2 = 4 * 3 * rmean * rmean;
215 for (j = 0; (j < 4); j++)
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]) );
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 */
235 for (i = 0; (i < nslice); i++)
239 sgslice[i] /= sl_count[i];
240 skslice[i] /= sl_count[i];
247 for (i = 0; (i < 4); i++)
255 static void calc_tetra_order_parm(const char *fnNDX, const char *fnTPS,
256 const char *fnTRX, const char *sgfn,
258 int nslice, int slice_dim,
259 const char *sgslfn, const char *skslfn,
260 const output_env_t oenv)
262 FILE *fpsg = NULL, *fpsk = NULL;
265 char title[STRLEN], fn[STRLEN], subtitle[STRLEN];
274 int i, *isize, ng, nframes;
275 real *sg_slice, *sg_slice_tot, *sk_slice, *sk_slice_tot;
276 gmx_rmpbc_t gpbc = NULL;
279 read_tps_conf(fnTPS, title, &top, &ePBC, &xtop, NULL, box, FALSE);
281 snew(sg_slice, nslice);
282 snew(sk_slice, nslice);
283 snew(sg_slice_tot, nslice);
284 snew(sk_slice_tot, nslice);
286 /* get index groups */
287 printf("Select the group that contains the atoms you want to use for the tetrahedrality order parameter calculation:\n");
291 get_index(&top.atoms, fnNDX, ng, isize, index, grpname);
293 /* Analyze trajectory */
294 natoms = read_first_x(oenv, &status, fnTRX, &t, &x, box);
295 if (natoms > top.atoms.nr)
297 gmx_fatal(FARGS, "Topology (%d atoms) does not match trajectory (%d atoms)",
298 top.atoms.nr, natoms);
300 check_index(NULL, ng, index[0], NULL, natoms);
302 fpsg = xvgropen(sgfn, "S\\sg\\N Angle Order Parameter", "Time (ps)", "S\\sg\\N",
304 fpsk = xvgropen(skfn, "S\\sk\\N Distance Order Parameter", "Time (ps)", "S\\sk\\N",
307 /* loop over frames */
308 gpbc = gmx_rmpbc_init(&top.idef, ePBC, natoms);
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++)
316 sg_slice_tot[i] += sg_slice[i];
317 sk_slice_tot[i] += sk_slice[i];
319 fprintf(fpsg, "%f %f\n", t, sg);
320 fprintf(fpsk, "%f %f\n", t, sk);
323 while (read_next_x(oenv, status, &t, x, box));
325 gmx_rmpbc_done(gpbc);
334 fpsg = xvgropen(sgslfn,
335 "S\\sg\\N Angle Order Parameter / Slab", "(nm)", "S\\sg\\N",
337 fpsk = xvgropen(skslfn,
338 "S\\sk\\N Distance Order Parameter / Slab", "(nm)", "S\\sk\\N",
340 for (i = 0; (i < nslice); i++)
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);
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)
358 fprintf(stderr, "Using following groups: \n");
359 for (i = 0; i < ngrps; i++)
361 fprintf(stderr, "Groupname: %s First atomname: %s First atomnr %d\n",
362 groups[i], *(top->atoms.atomname[a[index[i]]]), a[index[i]]);
364 fprintf(stderr, "\n");
367 static void check_length(real length, int a, int b)
371 fprintf(stderr, "WARNING: distance between atoms %d and "
372 "%d > 0.3 nm (%f). Index file might be corrupt.\n",
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,
382 const output_env_t oenv)
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) */
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 */
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;
414 real arcdist, tmpdist;
415 gmx_rmpbc_t gpbc = NULL;
417 /* PBC added for center-of-mass vector*/
418 /* Initiate the pbc structure */
419 memset(&pbc, 0, sizeof(pbc));
421 if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
423 gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
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*/
433 bSliced = FALSE; /*force slices off */
434 fprintf(stderr, "Calculating order parameters for each of %d molecules\n",
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);
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*/
455 if (use_unitvector && bSliced)
457 fprintf(stderr, "Warning: slicing and specified unit vectors are not currently compatible\n");
460 snew(slCount, nslices);
461 snew(*slOrder, nslices);
462 for (i = 0; i < nslices; i++)
464 snew((*slOrder)[i], ngrps);
468 snew(*distvals, nslices);
469 for (i = 0; i < nslices; i++)
471 snew((*distvals)[i], ngrps);
475 snew(slFrameorder, nslices);
480 *slWidth = box[axis][axis]/nslices;
481 fprintf(stderr, "Box divided in %d slices. Initial width of slice: %f\n",
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
495 gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
496 /*********** Start processing trajectory ***********/
501 *slWidth = box[axis][axis]/nslices;
505 set_pbc(&pbc, ePBC, box);
506 gmx_rmpbc_copy(gpbc, natoms, box, x0, x1);
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
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++)
523 rvec_inc(com, x1[comidx[j]]);
525 svmul(1.0/comsize, com, com);
529 dref[XX] = 0.0; dref[YY] = 0.0; dref[ZZ] = 0.0;
530 for (j = 0; j < distsize; j++)
532 rvec_inc(dist, x1[distidx[j]]);
534 svmul(1.0/distsize, dref, dref);
537 pbc_dx(&pbc, dref, com, dvec);
542 for (i = 1; i < ngrps - 1; i++)
544 clear_rvec(frameorder);
546 size = index[i+1] - index[i];
547 if (size != nr_tails)
549 gmx_fatal(FARGS, "grp %d does not have same number of"
550 " elements as grp 1\n", i);
553 for (j = 0; j < size; j++)
556 /*create unit vector*/
558 pbc_dx(&pbc, x1[a[index[i]+j]], com, direction);
559 unitv(direction, direction);
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]);*/
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);
572 check_length(length, a[index[i]+j], a[index[i+1]+j]);
573 svmul(1/length, dist, Sz);
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 */
580 sdbangle += gmx_angle(direction, Sz); /*this can probably be optimized*/
584 sdbangle += acos(Sz[axis]);
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]);
594 svmul(1/length, dist, Sz);
595 /* Sz is now the molecular axis Sz, normalized and all that */
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);
605 /* now we can get Sy from the outer product of Sx and Sz */
607 svmul(1/norm(Sy), Sy, Sy);
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
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));
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]);
626 for (m = 0; m < DIM; m++)
628 frameorder[m] += 0.5 * (3 * cossum[m] - 1);
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
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);
645 z_ave += box[axis][axis];
647 if (z_ave > box[axis][axis])
649 z_ave -= box[axis][axis];
652 slice = (int)(0.5 + (z_ave / (*slWidth))) - 1;
653 slCount[slice]++; /* determine slice, increase count */
655 slFrameorder[slice] += 0.5 * (3 * cossum[axis] - 1);
657 else if (permolecule)
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));
668 /* bin order parameter by arc distance from reference group*/
669 arcdist = gmx_angle(dvec, direction);
670 (*distvals)[j][i] += arcdist;
674 /* Want minimum lateral distance to first group calculated */
675 tmpdist = trace(box); /* should be max value */
676 for (k = 0; k < distsize; k++)
678 pbc_dx(&pbc, x1[distidx[k]], x1[a[index[i]+j]], dvec);
679 /* at the moment, just remove dvec[axis] */
681 tmpdist = min(tmpdist, norm2(dvec));
683 //fprintf(stderr, "Min dist %f; trace %f\n", tmpdist, trace(box));
684 (*distvals)[j][i] += sqrt(tmpdist);
687 } /* end loop j, over all atoms in group */
689 for (m = 0; m < DIM; m++)
691 (*order)[i][m] += (frameorder[m]/size);
695 { /*Skip following if doing per-molecule*/
696 for (k = 0; k < nslices; k++)
698 if (slCount[k]) /* if no elements, nothing has to be added */
700 (*slOrder)[k][i] += slFrameorder[k]/slCount[k];
701 slFrameorder[k] = 0; slCount[k] = 0;
704 } /* end loop i, over all groups in indexfile */
709 while (read_next_x(oenv, status, &t, x0, box));
710 /*********** done with status file **********/
712 fprintf(stderr, "\nRead trajectory. Printing parameters to file\n");
713 gmx_rmpbc_done(gpbc);
715 /* average over frames */
716 for (i = 1; i < ngrps - 1; i++)
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)
723 for (k = 0; k < nslices; k++)
725 (*slOrder)[k][i] /= nr_frames;
730 for (k = 0; k < nslices; k++)
732 (*distvals)[k][i] /= nr_frames;
739 fprintf(stderr, "Average angle between double bond and normal: %f\n",
740 180*sdbangle/(nr_frames * size*M_PI));
743 sfree(x0); /* free memory used by coordinate arrays */
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)
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 */
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++)
777 fprintf(ord, "%12d %12g\n", atom, -1 * (0.6667 * order[atom][XX] +
778 0.333 * order[atom][YY]));
781 for (slice = 0; slice < nslices; slice++)
783 fprintf(slOrd, "%12d\t", slice);
786 fprintf(slOrd, "%12g\t", distvals[slice][1]); /*use distance value at second carbon*/
788 for (atom = 1; atom < ngrps - 1; atom++)
790 fprintf(slOrd, "%12g\t", slOrder[slice][atom]);
792 fprintf(slOrd, "\n");
798 sprintf(buf, "Orderparameters Sz per atom");
799 ord = xvgropen(afile, buf, "Atom", "S", oenv);
800 fprintf(stderr, "ngrps = %d, nslices = %d", ngrps, nslices);
802 sprintf(buf, "Orderparameters per atom per slice");
803 slOrd = xvgropen(bfile, buf, "Slice", "S", oenv);
805 for (atom = 1; atom < ngrps - 1; atom++)
807 fprintf(ord, "%12d %12g\n", atom, order[atom][ZZ]);
810 for (slice = 0; slice < nslices; slice++)
813 for (atom = 1; atom < ngrps - 1; atom++)
815 S += slOrder[slice][atom];
817 fprintf(slOrd, "%12g %12g\n", slice*slWidth, S/atom);
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);
828 for (atom = 1; atom < ngrps - 1; atom++)
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]));
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)
843 /*function to write order parameters as B factors in PDB file using
844 first frame of trajectory*/
847 t_trxframe fr, frout;
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,
864 init_t_atoms(&useatoms, nout, TRUE);
867 /*initialize PDBinfo*/
868 for (i = 0; i < useatoms.nr; ++i)
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;
876 for (j = 0, ctr = 0; j < nslices; j++)
878 for (i = 0; i < ngrps; i++, ctr++)
880 /*iterate along each chain*/
881 useatoms.pdbinfo[ctr].bfac = order[j][i+1];
884 useatoms.pdbinfo[ctr].occup = distvals[j][i+1];
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*/
894 write_sto_conf(opt2fn("-ob", nfile, fnm), "Order parameters", &useatoms, frout.x, NULL, frout.ePBC, frout.box);
897 free_t_atoms(&useatoms, FALSE);
900 int gmx_order(int argc, char *argv[])
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",
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]",
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 */
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" },
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 */
956 axis = 0; /* normal axis */
957 t_topology *top; /* topology */
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 */
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;
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))
989 gmx_fatal(FARGS, "Can not have nslices < 1");
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);
998 if (strcmp(normal_axis[0], "x") == 0)
1002 else if (strcmp(normal_axis[0], "y") == 0)
1006 else if (strcmp(normal_axis[0], "z") == 0)
1012 gmx_fatal(FARGS, "Invalid axis, use x, y or z");
1018 fprintf(stderr, "Taking x axis as normal to the membrane\n");
1021 fprintf(stderr, "Taking y axis as normal to the membrane\n");
1024 fprintf(stderr, "Taking z axis as normal to the membrane\n");
1028 /* tetraheder order parameter */
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),
1037 /* view xvgr files */
1038 do_view(oenv, opt2fn("-Sg", NFILE, fnm), NULL);
1039 do_view(oenv, opt2fn("-Sk", NFILE, fnm), NULL);
1042 do_view(oenv, opt2fn("-Sgsl", NFILE, fnm), NULL);
1043 do_view(oenv, opt2fn("-Sksl", NFILE, fnm), NULL);
1048 /* tail order parameter */
1053 fprintf(stderr, "Dividing box in %d slices.\n\n", nslices);
1058 fprintf(stderr, "Only calculating Sz\n");
1062 fprintf(stderr, "Taking carbons as unsaturated!\n");
1065 top = read_top(ftp2fn(efTPX, NFILE, fnm), &ePBC); /* read topology file */
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 */
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);
1080 fprintf(stderr, "Calculating radial distances\n");
1083 gmx_fatal(FARGS, "Cannot yet output radial distances without permolecule\n");
1087 /* show atomtypes, to check if index file is correct */
1088 print_types(index, a, ngrps, grpname, top);
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);
1096 ngrps--; /*don't print the last group--was used for
1097 center-of-mass determination*/
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);
1103 if (opt2bSet("-ob", NFILE, fnm))
1108 "Won't write B-factors with averaged order parameters; use -permolecule\n");
1112 write_bfactors(fnm, NFILE, index, a, nslices, ngrps, slOrder, top, distvals, oenv);
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 */
1122 if (distvals != NULL)
1124 for (i = 0; i < nslices; ++i)