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.
46 #include "gromacs/utility/smalloc.h"
52 #include "gromacs/fileio/futil.h"
53 #include "gromacs/commandline/pargs.h"
55 #include "gromacs/fileio/tpxio.h"
56 #include "gromacs/fileio/trxio.h"
57 #include "gromacs/fileio/confio.h"
62 /****************************************************************************/
63 /* This program calculates the order parameter per atom for an interface or */
64 /* bilayer, averaged over time. */
65 /* S = 1/2 * (3 * cos(i)cos(j) - delta(ij)) */
66 /* It is assumed that the order parameter with respect to a box-axis */
67 /* is calculated. In that case i = j = axis, and delta(ij) = 1. */
69 /* Peter Tieleman, April 1995 */
70 /* P.J. van Maaren, November 2005 Added tetrahedral stuff */
71 /****************************************************************************/
73 static void find_nearest_neighbours(int ePBC,
74 int natoms, matrix box,
75 rvec x[], int maxidx, atom_id index[],
76 real *sgmean, real *skmean,
77 int nslice, int slice_dim,
78 real sgslice[], real skslice[],
83 int ix, jx, nsgbin, *sgbin;
84 int i1, i2, i, ibin, j, k, l, n, *nn[4];
85 rvec dx, dx1, dx2, rj, rk, urk, urj;
86 real cost, cost2, *sgmol, *skmol, rmean, rmean2, r2, box2, *r_nn[4];
92 real onethird = 1.0/3.0;
93 /* dmat = init_mat(maxidx, FALSE); */
94 box2 = box[XX][XX] * box[XX][XX];
95 snew(sl_count, nslice);
96 for (i = 0; (i < 4); i++)
98 snew(r_nn[i], natoms);
101 for (j = 0; (j < natoms); j++)
110 /* Must init pbc every step because of pressure coupling */
111 set_pbc(&pbc, ePBC, box);
113 gmx_rmpbc(gpbc, natoms, box, x);
115 nsgbin = 1 + 1/0.0005;
121 for (i = 0; (i < maxidx); i++) /* loop over index file */
124 for (j = 0; (j < maxidx); j++)
133 pbc_dx(&pbc, x[ix], x[jx], dx);
136 /* set_mat_entry(dmat,i,j,r2); */
138 /* determine the nearest neighbours */
141 r_nn[3][i] = r_nn[2][i]; nn[3][i] = nn[2][i];
142 r_nn[2][i] = r_nn[1][i]; nn[2][i] = nn[1][i];
143 r_nn[1][i] = r_nn[0][i]; nn[1][i] = nn[0][i];
144 r_nn[0][i] = r2; nn[0][i] = j;
146 else if (r2 < r_nn[1][i])
148 r_nn[3][i] = r_nn[2][i]; nn[3][i] = nn[2][i];
149 r_nn[2][i] = r_nn[1][i]; nn[2][i] = nn[1][i];
150 r_nn[1][i] = r2; nn[1][i] = j;
152 else if (r2 < r_nn[2][i])
154 r_nn[3][i] = r_nn[2][i]; nn[3][i] = nn[2][i];
155 r_nn[2][i] = r2; nn[2][i] = j;
157 else if (r2 < r_nn[3][i])
159 r_nn[3][i] = r2; nn[3][i] = j;
164 /* calculate mean distance between nearest neighbours */
166 for (j = 0; (j < 4); j++)
168 r_nn[j][i] = sqrt(r_nn[j][i]);
177 /* Chau1998a eqn 3 */
178 /* angular part tetrahedrality order parameter per atom */
179 for (j = 0; (j < 3); j++)
181 for (k = j+1; (k < 4); k++)
183 pbc_dx(&pbc, x[ix], x[index[nn[k][i]]], rk);
184 pbc_dx(&pbc, x[ix], x[index[nn[j][i]]], rj);
189 cost = iprod(urk, urj) + onethird;
192 /* sgmol[i] += 3*cost2/32; */
195 /* determine distribution */
196 ibin = nsgbin * cost2;
201 /* printf("%d %d %f %d %d\n", j, k, cost * cost, ibin, sgbin[ibin]);*/
207 /* normalize sgmol between 0.0 and 1.0 */
208 sgmol[i] = 3*sgmol[i]/32;
211 /* distance part tetrahedrality order parameter per atom */
212 rmean2 = 4 * 3 * rmean * rmean;
213 for (j = 0; (j < 4); j++)
215 skmol[i] += (rmean - r_nn[j][i]) * (rmean - r_nn[j][i]) / rmean2;
216 /* printf("%d %f (%f %f %f %f) \n",
217 i, skmol[i], rmean, rmean2, r_nn[j][i], (rmean - r_nn[j][i]) );
223 /* Compute sliced stuff */
224 sl_index = gmx_nint((1+x[i][slice_dim]/box[slice_dim][slice_dim])*nslice) % nslice;
225 sgslice[sl_index] += sgmol[i];
226 skslice[sl_index] += skmol[i];
227 sl_count[sl_index]++;
228 } /* loop over entries in index file */
233 for (i = 0; (i < nslice); i++)
237 sgslice[i] /= sl_count[i];
238 skslice[i] /= sl_count[i];
245 for (i = 0; (i < 4); i++)
253 static void calc_tetra_order_parm(const char *fnNDX, const char *fnTPS,
254 const char *fnTRX, const char *sgfn,
256 int nslice, int slice_dim,
257 const char *sgslfn, const char *skslfn,
258 const output_env_t oenv)
260 FILE *fpsg = NULL, *fpsk = NULL;
263 char title[STRLEN], fn[STRLEN], subtitle[STRLEN];
272 int i, *isize, ng, nframes;
273 real *sg_slice, *sg_slice_tot, *sk_slice, *sk_slice_tot;
274 gmx_rmpbc_t gpbc = NULL;
277 read_tps_conf(fnTPS, title, &top, &ePBC, &xtop, NULL, box, FALSE);
279 snew(sg_slice, nslice);
280 snew(sk_slice, nslice);
281 snew(sg_slice_tot, nslice);
282 snew(sk_slice_tot, nslice);
284 /* get index groups */
285 printf("Select the group that contains the atoms you want to use for the tetrahedrality order parameter calculation:\n");
289 get_index(&top.atoms, fnNDX, ng, isize, index, grpname);
291 /* Analyze trajectory */
292 natoms = read_first_x(oenv, &status, fnTRX, &t, &x, box);
293 if (natoms > top.atoms.nr)
295 gmx_fatal(FARGS, "Topology (%d atoms) does not match trajectory (%d atoms)",
296 top.atoms.nr, natoms);
298 check_index(NULL, ng, index[0], NULL, natoms);
300 fpsg = xvgropen(sgfn, "S\\sg\\N Angle Order Parameter", "Time (ps)", "S\\sg\\N",
302 fpsk = xvgropen(skfn, "S\\sk\\N Distance Order Parameter", "Time (ps)", "S\\sk\\N",
305 /* loop over frames */
306 gpbc = gmx_rmpbc_init(&top.idef, ePBC, natoms);
310 find_nearest_neighbours(ePBC, natoms, box, x, isize[0], index[0],
311 &sg, &sk, nslice, slice_dim, sg_slice, sk_slice, gpbc);
312 for (i = 0; (i < nslice); i++)
314 sg_slice_tot[i] += sg_slice[i];
315 sk_slice_tot[i] += sk_slice[i];
317 fprintf(fpsg, "%f %f\n", t, sg);
318 fprintf(fpsk, "%f %f\n", t, sk);
321 while (read_next_x(oenv, status, &t, x, box));
323 gmx_rmpbc_done(gpbc);
332 fpsg = xvgropen(sgslfn,
333 "S\\sg\\N Angle Order Parameter / Slab", "(nm)", "S\\sg\\N",
335 fpsk = xvgropen(skslfn,
336 "S\\sk\\N Distance Order Parameter / Slab", "(nm)", "S\\sk\\N",
338 for (i = 0; (i < nslice); i++)
340 fprintf(fpsg, "%10g %10g\n", (i+0.5)*box[slice_dim][slice_dim]/nslice,
341 sg_slice_tot[i]/nframes);
342 fprintf(fpsk, "%10g %10g\n", (i+0.5)*box[slice_dim][slice_dim]/nslice,
343 sk_slice_tot[i]/nframes);
350 /* Print name of first atom in all groups in index file */
351 static void print_types(atom_id index[], atom_id a[], int ngrps,
352 char *groups[], t_topology *top)
356 fprintf(stderr, "Using following groups: \n");
357 for (i = 0; i < ngrps; i++)
359 fprintf(stderr, "Groupname: %s First atomname: %s First atomnr %u\n",
360 groups[i], *(top->atoms.atomname[a[index[i]]]), a[index[i]]);
362 fprintf(stderr, "\n");
365 static void check_length(real length, int a, int b)
369 fprintf(stderr, "WARNING: distance between atoms %d and "
370 "%d > 0.3 nm (%f). Index file might be corrupt.\n",
375 void calc_order(const char *fn, atom_id *index, atom_id *a, rvec **order,
376 real ***slOrder, real *slWidth, int nslices, gmx_bool bSliced,
377 gmx_bool bUnsat, t_topology *top, int ePBC, int ngrps, int axis,
378 gmx_bool permolecule, gmx_bool radial, gmx_bool distcalc, const char *radfn,
380 const output_env_t oenv)
382 /* if permolecule = TRUE, order parameters will be calculed per molecule
383 * and stored in slOrder with #slices = # molecules */
384 rvec *x0, /* coordinates with pbc */
385 *x1, /* coordinates without pbc */
386 dist; /* vector between two atoms */
387 matrix box; /* box (3x3) */
389 rvec cossum, /* sum of vector angles for three axes */
390 Sx, Sy, Sz, /* the three molecular axes */
391 tmp1, tmp2, /* temp. rvecs for calculating dot products */
392 frameorder; /* order parameters for one frame */
393 real *slFrameorder; /* order parameter for one frame, per slice */
394 real length, /* total distance between two atoms */
395 t, /* time from trajectory */
396 z_ave, z1, z2; /* average z, used to det. which slice atom is in */
397 int natoms, /* nr. atoms in trj */
398 nr_tails, /* nr tails, to check if index file is correct */
399 size = 0, /* nr. of atoms in group. same as nr_tails */
400 i, j, m, k, l, teller = 0,
401 slice, /* current slice number */
403 int *slCount; /* nr. of atoms in one slice */
404 real dbangle = 0, /* angle between double bond and axis */
405 sdbangle = 0; /* sum of these angles */
406 gmx_bool use_unitvector = FALSE; /* use a specified unit vector instead of axis to specify unit normal*/
407 rvec direction, com, dref, dvec;
408 int comsize, distsize;
409 atom_id *comidx = NULL, *distidx = NULL;
410 char *grpname = NULL;
412 real arcdist, tmpdist;
413 gmx_rmpbc_t gpbc = NULL;
415 /* PBC added for center-of-mass vector*/
416 /* Initiate the pbc structure */
417 memset(&pbc, 0, sizeof(pbc));
419 if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
421 gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
424 nr_tails = index[1] - index[0];
425 fprintf(stderr, "Number of elements in first group: %d\n", nr_tails);
426 /* take first group as standard. Not rocksolid, but might catch error in index*/
431 bSliced = FALSE; /*force slices off */
432 fprintf(stderr, "Calculating order parameters for each of %d molecules\n",
438 use_unitvector = TRUE;
439 fprintf(stderr, "Select an index group to calculate the radial membrane normal\n");
440 get_index(&top->atoms, radfn, 1, &comsize, &comidx, &grpname);
448 fprintf(stderr, "Select an index group to use as distance reference\n");
449 get_index(&top->atoms, radfn, 1, &distsize, &distidx, &grpname);
450 bSliced = FALSE; /*force slices off*/
453 if (use_unitvector && bSliced)
455 fprintf(stderr, "Warning: slicing and specified unit vectors are not currently compatible\n");
458 snew(slCount, nslices);
459 snew(*slOrder, nslices);
460 for (i = 0; i < nslices; i++)
462 snew((*slOrder)[i], ngrps);
466 snew(*distvals, nslices);
467 for (i = 0; i < nslices; i++)
469 snew((*distvals)[i], ngrps);
473 snew(slFrameorder, nslices);
478 *slWidth = box[axis][axis]/nslices;
479 fprintf(stderr, "Box divided in %d slices. Initial width of slice: %f\n",
485 nr_tails = index[1] - index[0];
486 fprintf(stderr, "Number of elements in first group: %d\n", nr_tails);
487 /* take first group as standard. Not rocksolid, but might catch error
493 gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
494 /*********** Start processing trajectory ***********/
499 *slWidth = box[axis][axis]/nslices;
503 set_pbc(&pbc, ePBC, box);
504 gmx_rmpbc_copy(gpbc, natoms, box, x0, x1);
506 /* Now loop over all groups. There are ngrps groups, the order parameter can
507 be calculated for grp 1 to grp ngrps - 1. For each group, loop over all
508 atoms in group, which is index[i] to (index[i+1] - 1) See block.h. Of
509 course, in this case index[i+1] -index[i] has to be the same for all
510 groups, namely the number of tails. i just runs over all atoms in a tail,
511 so for DPPC ngrps = 16 and i runs from 1 to 14, including 14
517 /*center-of-mass determination*/
518 com[XX] = 0.0; com[YY] = 0.0; com[ZZ] = 0.0;
519 for (j = 0; j < comsize; j++)
521 rvec_inc(com, x1[comidx[j]]);
523 svmul(1.0/comsize, com, com);
527 dref[XX] = 0.0; dref[YY] = 0.0; dref[ZZ] = 0.0;
528 for (j = 0; j < distsize; j++)
530 rvec_inc(dist, x1[distidx[j]]);
532 svmul(1.0/distsize, dref, dref);
535 pbc_dx(&pbc, dref, com, dvec);
540 for (i = 1; i < ngrps - 1; i++)
542 clear_rvec(frameorder);
544 size = index[i+1] - index[i];
545 if (size != nr_tails)
547 gmx_fatal(FARGS, "grp %d does not have same number of"
548 " elements as grp 1\n", i);
551 for (j = 0; j < size; j++)
554 /*create unit vector*/
556 pbc_dx(&pbc, x1[a[index[i]+j]], com, direction);
557 unitv(direction, direction);
560 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],
561 direction[0],direction[1],direction[2]);*/
566 /* Using convention for unsaturated carbons */
567 /* first get Sz, the vector from Cn to Cn+1 */
568 rvec_sub(x1[a[index[i+1]+j]], x1[a[index[i]+j]], dist);
570 check_length(length, a[index[i]+j], a[index[i+1]+j]);
571 svmul(1/length, dist, Sz);
573 /* this is actually the cosine of the angle between the double bond
574 and axis, because Sz is normalized and the two other components of
575 the axis on the bilayer are zero */
578 sdbangle += gmx_angle(direction, Sz); /*this can probably be optimized*/
582 sdbangle += acos(Sz[axis]);
587 /* get vector dist(Cn-1,Cn+1) for tail atoms */
588 rvec_sub(x1[a[index[i+1]+j]], x1[a[index[i-1]+j]], dist);
589 length = norm(dist); /* determine distance between two atoms */
590 check_length(length, a[index[i-1]+j], a[index[i+1]+j]);
592 svmul(1/length, dist, Sz);
593 /* Sz is now the molecular axis Sz, normalized and all that */
596 /* now get Sx. Sx is normal to the plane of Cn-1, Cn and Cn+1 so
597 we can use the outer product of Cn-1->Cn and Cn+1->Cn, I hope */
598 rvec_sub(x1[a[index[i+1]+j]], x1[a[index[i]+j]], tmp1);
599 rvec_sub(x1[a[index[i-1]+j]], x1[a[index[i]+j]], tmp2);
600 cprod(tmp1, tmp2, Sx);
601 svmul(1/norm(Sx), Sx, Sx);
603 /* now we can get Sy from the outer product of Sx and Sz */
605 svmul(1/norm(Sy), Sy, Sy);
607 /* the square of cosine of the angle between dist and the axis.
608 Using the innerproduct, but two of the three elements are zero
609 Determine the sum of the orderparameter of all atoms in group
613 cossum[XX] = sqr(iprod(Sx, direction)); /* this is allowed, since Sa is normalized */
614 cossum[YY] = sqr(iprod(Sy, direction));
615 cossum[ZZ] = sqr(iprod(Sz, direction));
619 cossum[XX] = sqr(Sx[axis]); /* this is allowed, since Sa is normalized */
620 cossum[YY] = sqr(Sy[axis]);
621 cossum[ZZ] = sqr(Sz[axis]);
624 for (m = 0; m < DIM; m++)
626 frameorder[m] += 0.5 * (3 * cossum[m] - 1);
631 /* get average coordinate in box length for slicing,
632 determine which slice atom is in, increase count for that
633 slice. slFrameorder and slOrder are reals, not
634 rvecs. Only the component [axis] of the order tensor is
635 kept, until I find it necessary to know the others too
638 z1 = x1[a[index[i-1]+j]][axis];
639 z2 = x1[a[index[i+1]+j]][axis];
640 z_ave = 0.5 * (z1 + z2);
643 z_ave += box[axis][axis];
645 if (z_ave > box[axis][axis])
647 z_ave -= box[axis][axis];
650 slice = (int)(0.5 + (z_ave / (*slWidth))) - 1;
651 slCount[slice]++; /* determine slice, increase count */
653 slFrameorder[slice] += 0.5 * (3 * cossum[axis] - 1);
655 else if (permolecule)
657 /* store per-molecule order parameter
658 * To just track single-axis order: (*slOrder)[j][i] += 0.5 * (3 * iprod(cossum,direction) - 1);
659 * following is for Scd order: */
660 (*slOrder)[j][i] += -1* (0.3333 * (3 * cossum[XX] - 1) + 0.3333 * 0.5 * (3 * cossum[YY] - 1));
666 /* bin order parameter by arc distance from reference group*/
667 arcdist = gmx_angle(dvec, direction);
668 (*distvals)[j][i] += arcdist;
672 /* Want minimum lateral distance to first group calculated */
673 tmpdist = trace(box); /* should be max value */
674 for (k = 0; k < distsize; k++)
676 pbc_dx(&pbc, x1[distidx[k]], x1[a[index[i]+j]], dvec);
677 /* at the moment, just remove dvec[axis] */
679 tmpdist = min(tmpdist, norm2(dvec));
681 //fprintf(stderr, "Min dist %f; trace %f\n", tmpdist, trace(box));
682 (*distvals)[j][i] += sqrt(tmpdist);
685 } /* end loop j, over all atoms in group */
687 for (m = 0; m < DIM; m++)
689 (*order)[i][m] += (frameorder[m]/size);
693 { /*Skip following if doing per-molecule*/
694 for (k = 0; k < nslices; k++)
696 if (slCount[k]) /* if no elements, nothing has to be added */
698 (*slOrder)[k][i] += slFrameorder[k]/slCount[k];
699 slFrameorder[k] = 0; slCount[k] = 0;
702 } /* end loop i, over all groups in indexfile */
707 while (read_next_x(oenv, status, &t, x0, box));
708 /*********** done with status file **********/
710 fprintf(stderr, "\nRead trajectory. Printing parameters to file\n");
711 gmx_rmpbc_done(gpbc);
713 /* average over frames */
714 for (i = 1; i < ngrps - 1; i++)
716 svmul(1.0/nr_frames, (*order)[i], (*order)[i]);
717 fprintf(stderr, "Atom %d Tensor: x=%g , y=%g, z=%g\n", i, (*order)[i][XX],
718 (*order)[i][YY], (*order)[i][ZZ]);
719 if (bSliced || permolecule)
721 for (k = 0; k < nslices; k++)
723 (*slOrder)[k][i] /= nr_frames;
728 for (k = 0; k < nslices; k++)
730 (*distvals)[k][i] /= nr_frames;
737 fprintf(stderr, "Average angle between double bond and normal: %f\n",
738 180*sdbangle/(nr_frames * size*M_PI));
741 sfree(x0); /* free memory used by coordinate arrays */
758 void order_plot(rvec order[], real *slOrder[], const char *afile, const char *bfile,
759 const char *cfile, int ngrps, int nslices, real slWidth, gmx_bool bSzonly,
760 gmx_bool permolecule, real **distvals, const output_env_t oenv)
762 FILE *ord, *slOrd; /* xvgr files with order parameters */
763 int atom, slice; /* atom corresponding to order para.*/
764 char buf[256]; /* for xvgr title */
765 real S; /* order parameter averaged over all atoms */
769 sprintf(buf, "Scd order parameters");
770 ord = xvgropen(afile, buf, "Atom", "S", oenv);
771 sprintf(buf, "Orderparameters per atom per slice");
772 slOrd = xvgropen(bfile, buf, "Molecule", "S", oenv);
773 for (atom = 1; atom < ngrps - 1; atom++)
775 fprintf(ord, "%12d %12g\n", atom, -1 * (0.6667 * order[atom][XX] +
776 0.333 * order[atom][YY]));
779 for (slice = 0; slice < nslices; slice++)
781 fprintf(slOrd, "%12d\t", slice);
784 fprintf(slOrd, "%12g\t", distvals[slice][1]); /*use distance value at second carbon*/
786 for (atom = 1; atom < ngrps - 1; atom++)
788 fprintf(slOrd, "%12g\t", slOrder[slice][atom]);
790 fprintf(slOrd, "\n");
796 sprintf(buf, "Orderparameters Sz per atom");
797 ord = xvgropen(afile, buf, "Atom", "S", oenv);
798 fprintf(stderr, "ngrps = %d, nslices = %d", ngrps, nslices);
800 sprintf(buf, "Orderparameters per atom per slice");
801 slOrd = xvgropen(bfile, buf, "Slice", "S", oenv);
803 for (atom = 1; atom < ngrps - 1; atom++)
805 fprintf(ord, "%12d %12g\n", atom, order[atom][ZZ]);
808 for (slice = 0; slice < nslices; slice++)
811 for (atom = 1; atom < ngrps - 1; atom++)
813 S += slOrder[slice][atom];
815 fprintf(slOrd, "%12g %12g\n", slice*slWidth, S/atom);
821 sprintf(buf, "Order tensor diagonal elements");
822 ord = xvgropen(afile, buf, "Atom", "S", oenv);
823 sprintf(buf, "Deuterium order parameters");
824 slOrd = xvgropen(cfile, buf, "Atom", "Scd", oenv);
826 for (atom = 1; atom < ngrps - 1; atom++)
828 fprintf(ord, "%12d %12g %12g %12g\n", atom, order[atom][XX],
829 order[atom][YY], order[atom][ZZ]);
830 fprintf(slOrd, "%12d %12g\n", atom, -1 * (0.6667 * order[atom][XX] +
831 0.333 * order[atom][YY]));
839 void write_bfactors(t_filenm *fnm, int nfile, atom_id *index, atom_id *a, int nslices, int ngrps, real **order, t_topology *top, real **distvals, output_env_t oenv)
841 /*function to write order parameters as B factors in PDB file using
842 first frame of trajectory*/
845 t_trxframe fr, frout;
849 ngrps -= 2; /*we don't have an order parameter for the first or
850 last atom in each chain*/
851 nout = nslices*ngrps;
852 natoms = read_first_frame(oenv, &status, ftp2fn(efTRX, nfile, fnm), &fr,
862 init_t_atoms(&useatoms, nout, TRUE);
865 /*initialize PDBinfo*/
866 for (i = 0; i < useatoms.nr; ++i)
868 useatoms.pdbinfo[i].type = 0;
869 useatoms.pdbinfo[i].occup = 0.0;
870 useatoms.pdbinfo[i].bfac = 0.0;
871 useatoms.pdbinfo[i].bAnisotropic = FALSE;
874 for (j = 0, ctr = 0; j < nslices; j++)
876 for (i = 0; i < ngrps; i++, ctr++)
878 /*iterate along each chain*/
879 useatoms.pdbinfo[ctr].bfac = order[j][i+1];
882 useatoms.pdbinfo[ctr].occup = distvals[j][i+1];
884 copy_rvec(fr.x[a[index[i+1]+j]], frout.x[ctr]);
885 useatoms.atomname[ctr] = top->atoms.atomname[a[index[i+1]+j]];
886 useatoms.atom[ctr] = top->atoms.atom[a[index[i+1]+j]];
887 useatoms.nres = max(useatoms.nres, useatoms.atom[ctr].resind+1);
888 useatoms.resinfo[useatoms.atom[ctr].resind] = top->atoms.resinfo[useatoms.atom[ctr].resind]; /*copy resinfo*/
892 write_sto_conf(opt2fn("-ob", nfile, fnm), "Order parameters", &useatoms, frout.x, NULL, frout.ePBC, frout.box);
895 free_t_atoms(&useatoms, FALSE);
898 int gmx_order(int argc, char *argv[])
900 const char *desc[] = {
901 "[THISMODULE] computes the order parameter per atom for carbon tails. For atom i the",
902 "vector i-1, i+1 is used together with an axis. ",
903 "The index file should contain only the groups to be used for calculations,",
904 "with each group of equivalent carbons along the relevant acyl chain in its own",
905 "group. There should not be any generic groups (like System, Protein) in the index",
906 "file to avoid confusing the program (this is not relevant to tetrahedral order",
907 "parameters however, which only work for water anyway).[PAR]",
908 "[THISMODULE] can also give all",
909 "diagonal elements of the order tensor and even calculate the deuterium",
910 "order parameter Scd (default). If the option [TT]-szonly[tt] is given, only one",
911 "order tensor component (specified by the [TT]-d[tt] option) is given and the",
912 "order parameter per slice is calculated as well. If [TT]-szonly[tt] is not",
913 "selected, all diagonal elements and the deuterium order parameter is",
915 "The tetrahedrality order parameters can be determined",
916 "around an atom. Both angle an distance order parameters are calculated. See",
917 "P.-L. Chau and A.J. Hardwick, Mol. Phys., 93, (1998), 511-518.",
918 "for more details.[BR]",
922 static int nslices = 1; /* nr of slices defined */
923 static gmx_bool bSzonly = FALSE; /* True if only Sz is wanted */
924 static gmx_bool bUnsat = FALSE; /* True if carbons are unsat. */
925 static const char *normal_axis[] = { NULL, "z", "x", "y", NULL };
926 static gmx_bool permolecule = FALSE; /*compute on a per-molecule basis */
927 static gmx_bool radial = FALSE; /*compute a radial membrane normal */
928 static gmx_bool distcalc = FALSE; /*calculate distance from a reference group */
930 { "-d", FALSE, etENUM, {normal_axis},
931 "Direction of the normal on the membrane" },
932 { "-sl", FALSE, etINT, {&nslices},
933 "Calculate order parameter as function of box length, dividing the box"
934 " into this number of slices." },
935 { "-szonly", FALSE, etBOOL, {&bSzonly},
936 "Only give Sz element of order tensor. (axis can be specified with [TT]-d[tt])" },
937 { "-unsat", FALSE, etBOOL, {&bUnsat},
938 "Calculate order parameters for unsaturated carbons. Note that this can"
939 "not be mixed with normal order parameters." },
940 { "-permolecule", FALSE, etBOOL, {&permolecule},
941 "Compute per-molecule Scd order parameters" },
942 { "-radial", FALSE, etBOOL, {&radial},
943 "Compute a radial membrane normal" },
944 { "-calcdist", FALSE, etBOOL, {&distcalc},
945 "Compute distance from a reference" },
948 rvec *order; /* order par. for each atom */
949 real **slOrder; /* same, per slice */
950 real slWidth = 0.0; /* width of a slice */
951 char **grpname; /* groupnames */
952 int ngrps, /* nr. of groups */
954 axis = 0; /* normal axis */
955 t_topology *top; /* topology */
957 atom_id *index, /* indices for a */
958 *a; /* atom numbers in each group */
959 t_blocka *block; /* data from index file */
960 t_filenm fnm[] = { /* files for g_order */
961 { efTRX, "-f", NULL, ffREAD }, /* trajectory file */
962 { efNDX, "-n", NULL, ffREAD }, /* index file */
963 { efNDX, "-nr", NULL, ffREAD }, /* index for radial axis calculation */
964 { efTPX, NULL, NULL, ffREAD }, /* topology file */
965 { efXVG, "-o", "order", ffWRITE }, /* xvgr output file */
966 { efXVG, "-od", "deuter", ffWRITE }, /* xvgr output file */
967 { efPDB, "-ob", NULL, ffWRITE }, /* write Scd as B factors to PDB if permolecule */
968 { efXVG, "-os", "sliced", ffWRITE }, /* xvgr output file */
969 { efXVG, "-Sg", "sg-ang", ffOPTWR }, /* xvgr output file */
970 { efXVG, "-Sk", "sk-dist", ffOPTWR }, /* xvgr output file */
971 { efXVG, "-Sgsl", "sg-ang-slice", ffOPTWR }, /* xvgr output file */
972 { efXVG, "-Sksl", "sk-dist-slice", ffOPTWR }, /* xvgr output file */
974 gmx_bool bSliced = FALSE; /* True if box is sliced */
975 #define NFILE asize(fnm)
976 real **distvals = NULL;
977 const char *sgfnm, *skfnm, *ndxfnm, *tpsfnm, *trxfnm;
980 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
981 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
987 gmx_fatal(FARGS, "Can not have nslices < 1");
989 sgfnm = opt2fn_null("-Sg", NFILE, fnm);
990 skfnm = opt2fn_null("-Sk", NFILE, fnm);
991 ndxfnm = opt2fn_null("-n", NFILE, fnm);
992 tpsfnm = ftp2fn(efTPX, NFILE, fnm);
993 trxfnm = ftp2fn(efTRX, NFILE, fnm);
996 if (strcmp(normal_axis[0], "x") == 0)
1000 else if (strcmp(normal_axis[0], "y") == 0)
1004 else if (strcmp(normal_axis[0], "z") == 0)
1010 gmx_fatal(FARGS, "Invalid axis, use x, y or z");
1016 fprintf(stderr, "Taking x axis as normal to the membrane\n");
1019 fprintf(stderr, "Taking y axis as normal to the membrane\n");
1022 fprintf(stderr, "Taking z axis as normal to the membrane\n");
1026 /* tetraheder order parameter */
1029 /* If either of theoptions is set we compute both */
1030 sgfnm = opt2fn("-Sg", NFILE, fnm);
1031 skfnm = opt2fn("-Sk", NFILE, fnm);
1032 calc_tetra_order_parm(ndxfnm, tpsfnm, trxfnm, sgfnm, skfnm, nslices, axis,
1033 opt2fn("-Sgsl", NFILE, fnm), opt2fn("-Sksl", NFILE, fnm),
1035 /* view xvgr files */
1036 do_view(oenv, opt2fn("-Sg", NFILE, fnm), NULL);
1037 do_view(oenv, opt2fn("-Sk", NFILE, fnm), NULL);
1040 do_view(oenv, opt2fn("-Sgsl", NFILE, fnm), NULL);
1041 do_view(oenv, opt2fn("-Sksl", NFILE, fnm), NULL);
1046 /* tail order parameter */
1051 fprintf(stderr, "Dividing box in %d slices.\n\n", nslices);
1056 fprintf(stderr, "Only calculating Sz\n");
1060 fprintf(stderr, "Taking carbons as unsaturated!\n");
1063 top = read_top(ftp2fn(efTPX, NFILE, fnm), &ePBC); /* read topology file */
1065 block = init_index(ftp2fn(efNDX, NFILE, fnm), &grpname);
1066 index = block->index; /* get indices from t_block block */
1067 a = block->a; /* see block.h */
1072 nslices = index[1] - index[0]; /*I think this assumes contiguous lipids in topology*/
1073 fprintf(stderr, "Calculating Scd order parameters for each of %d molecules\n", nslices);
1078 fprintf(stderr, "Calculating radial distances\n");
1081 gmx_fatal(FARGS, "Cannot yet output radial distances without permolecule\n");
1085 /* show atomtypes, to check if index file is correct */
1086 print_types(index, a, ngrps, grpname, top);
1088 calc_order(ftp2fn(efTRX, NFILE, fnm), index, a, &order,
1089 &slOrder, &slWidth, nslices, bSliced, bUnsat,
1090 top, ePBC, ngrps, axis, permolecule, radial, distcalc, opt2fn_null("-nr", NFILE, fnm), &distvals, oenv);
1094 ngrps--; /*don't print the last group--was used for
1095 center-of-mass determination*/
1098 order_plot(order, slOrder, opt2fn("-o", NFILE, fnm), opt2fn("-os", NFILE, fnm),
1099 opt2fn("-od", NFILE, fnm), ngrps, nslices, slWidth, bSzonly, permolecule, distvals, oenv);
1101 if (opt2bSet("-ob", NFILE, fnm))
1106 "Won't write B-factors with averaged order parameters; use -permolecule\n");
1110 write_bfactors(fnm, NFILE, index, a, nslices, ngrps, slOrder, top, distvals, oenv);
1115 do_view(oenv, opt2fn("-o", NFILE, fnm), NULL); /* view xvgr file */
1116 do_view(oenv, opt2fn("-os", NFILE, fnm), NULL); /* view xvgr file */
1117 do_view(oenv, opt2fn("-od", NFILE, fnm), NULL); /* view xvgr file */
1120 if (distvals != NULL)
1122 for (i = 0; i < nslices; ++i)