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,2015,2016,2017,2018,2019, 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/commandline/pargs.h"
45 #include "gromacs/commandline/viewit.h"
46 #include "gromacs/fileio/confio.h"
47 #include "gromacs/fileio/trxio.h"
48 #include "gromacs/fileio/xvgr.h"
49 #include "gromacs/gmxana/cmat.h"
50 #include "gromacs/gmxana/gmx_ana.h"
51 #include "gromacs/gmxana/gstat.h"
52 #include "gromacs/math/functions.h"
53 #include "gromacs/math/utilities.h"
54 #include "gromacs/math/vec.h"
55 #include "gromacs/pbcutil/pbc.h"
56 #include "gromacs/pbcutil/rmpbc.h"
57 #include "gromacs/topology/index.h"
58 #include "gromacs/topology/topology.h"
59 #include "gromacs/trajectory/trajectoryframe.h"
60 #include "gromacs/utility/arraysize.h"
61 #include "gromacs/utility/cstringutil.h"
62 #include "gromacs/utility/fatalerror.h"
63 #include "gromacs/utility/futil.h"
64 #include "gromacs/utility/gmxassert.h"
65 #include "gromacs/utility/smalloc.h"
67 /****************************************************************************/
68 /* This program calculates the order parameter per atom for an interface or */
69 /* bilayer, averaged over time. */
70 /* S = 1/2 * (3 * cos(i)cos(j) - delta(ij)) */
71 /* It is assumed that the order parameter with respect to a box-axis */
72 /* is calculated. In that case i = j = axis, and delta(ij) = 1. */
74 /* Peter Tieleman, April 1995 */
75 /* P.J. van Maaren, November 2005 Added tetrahedral stuff */
76 /****************************************************************************/
78 static void find_nearest_neighbours(int ePBC,
92 int ix, jx, nsgbin, *sgbin;
93 int i, ibin, j, k, l, n, *nn[4];
94 rvec dx, rj, rk, urk, urj;
95 real cost, cost2, *sgmol, *skmol, rmean, rmean2, r2, box2, *r_nn[4];
99 real onethird = 1.0 / 3.0;
100 /* dmat = init_mat(maxidx, FALSE); */
101 box2 = box[XX][XX] * box[XX][XX];
102 snew(sl_count, nslice);
103 for (i = 0; (i < 4); i++)
105 snew(r_nn[i], natoms);
108 for (j = 0; (j < natoms); j++)
117 /* Must init pbc every step because of pressure coupling */
118 set_pbc(&pbc, ePBC, box);
120 gmx_rmpbc(gpbc, natoms, box, x);
122 nsgbin = 2001; // Calculated as (1 + 1/0.0005)
128 for (i = 0; (i < maxidx); i++) /* loop over index file */
131 for (j = 0; (j < maxidx); j++)
140 pbc_dx(&pbc, x[ix], x[jx], dx);
143 /* set_mat_entry(dmat,i,j,r2); */
145 /* determine the nearest neighbours */
148 r_nn[3][i] = r_nn[2][i];
150 r_nn[2][i] = r_nn[1][i];
152 r_nn[1][i] = r_nn[0][i];
157 else if (r2 < r_nn[1][i])
159 r_nn[3][i] = r_nn[2][i];
161 r_nn[2][i] = r_nn[1][i];
166 else if (r2 < r_nn[2][i])
168 r_nn[3][i] = r_nn[2][i];
173 else if (r2 < r_nn[3][i])
181 /* calculate mean distance between nearest neighbours */
183 for (j = 0; (j < 4); j++)
185 r_nn[j][i] = std::sqrt(r_nn[j][i]);
194 /* Chau1998a eqn 3 */
195 /* angular part tetrahedrality order parameter per atom */
196 for (j = 0; (j < 3); j++)
198 for (k = j + 1; (k < 4); k++)
200 pbc_dx(&pbc, x[ix], x[index[nn[k][i]]], rk);
201 pbc_dx(&pbc, x[ix], x[index[nn[j][i]]], rj);
206 cost = iprod(urk, urj) + onethird;
209 /* sgmol[i] += 3*cost2/32; */
212 /* determine distribution */
213 ibin = static_cast<int>(static_cast<real>(nsgbin) * cost2);
218 /* printf("%d %d %f %d %d\n", j, k, cost * cost, ibin, sgbin[ibin]);*/
224 /* normalize sgmol between 0.0 and 1.0 */
225 sgmol[i] = 3 * sgmol[i] / 32;
228 /* distance part tetrahedrality order parameter per atom */
229 rmean2 = 4 * 3 * rmean * rmean;
230 for (j = 0; (j < 4); j++)
232 skmol[i] += (rmean - r_nn[j][i]) * (rmean - r_nn[j][i]) / rmean2;
233 /* printf("%d %f (%f %f %f %f) \n",
234 i, skmol[i], rmean, rmean2, r_nn[j][i], (rmean - r_nn[j][i]) );
240 /* Compute sliced stuff */
241 sl_index = static_cast<int>(std::round((1 + x[i][slice_dim] / box[slice_dim][slice_dim])
242 * static_cast<real>(nslice)))
244 sgslice[sl_index] += sgmol[i];
245 skslice[sl_index] += skmol[i];
246 sl_count[sl_index]++;
247 } /* loop over entries in index file */
249 *sgmean /= static_cast<real>(maxidx);
250 *skmean /= static_cast<real>(maxidx);
252 for (i = 0; (i < nslice); i++)
256 sgslice[i] /= sl_count[i];
257 skslice[i] /= sl_count[i];
264 for (i = 0; (i < 4); i++)
272 static void calc_tetra_order_parm(const char* fnNDX,
281 const gmx_output_env_t* oenv)
283 FILE * fpsg = nullptr, *fpsk = nullptr;
294 int i, *isize, ng, nframes;
295 real * sg_slice, *sg_slice_tot, *sk_slice, *sk_slice_tot;
296 gmx_rmpbc_t gpbc = nullptr;
299 read_tps_conf(fnTPS, &top, &ePBC, &xtop, nullptr, box, FALSE);
301 snew(sg_slice, nslice);
302 snew(sk_slice, nslice);
303 snew(sg_slice_tot, nslice);
304 snew(sk_slice_tot, nslice);
306 /* get index groups */
307 printf("Select the group that contains the atoms you want to use for the tetrahedrality order "
308 "parameter calculation:\n");
312 get_index(&top.atoms, fnNDX, ng, isize, index, grpname);
314 /* Analyze trajectory */
315 natoms = read_first_x(oenv, &status, fnTRX, &t, &x, box);
316 if (natoms > top.atoms.nr)
318 gmx_fatal(FARGS, "Topology (%d atoms) does not match trajectory (%d atoms)", top.atoms.nr, natoms);
320 check_index(nullptr, ng, index[0], nullptr, natoms);
322 fpsg = xvgropen(sgfn, "S\\sg\\N Angle Order Parameter", "Time (ps)", "S\\sg\\N", oenv);
323 fpsk = xvgropen(skfn, "S\\sk\\N Distance Order Parameter", "Time (ps)", "S\\sk\\N", oenv);
325 /* loop over frames */
326 gpbc = gmx_rmpbc_init(&top.idef, ePBC, natoms);
330 find_nearest_neighbours(ePBC, natoms, box, x, isize[0], index[0], &sg, &sk, nslice,
331 slice_dim, sg_slice, sk_slice, gpbc);
332 for (i = 0; (i < nslice); i++)
334 sg_slice_tot[i] += sg_slice[i];
335 sk_slice_tot[i] += sk_slice[i];
337 fprintf(fpsg, "%f %f\n", t, sg);
338 fprintf(fpsk, "%f %f\n", t, sk);
340 } while (read_next_x(oenv, status, &t, x, box));
342 gmx_rmpbc_done(gpbc);
351 fpsg = xvgropen(sgslfn, "S\\sg\\N Angle Order Parameter / Slab", "(nm)", "S\\sg\\N", oenv);
352 fpsk = xvgropen(skslfn, "S\\sk\\N Distance Order Parameter / Slab", "(nm)", "S\\sk\\N", oenv);
353 for (i = 0; (i < nslice); i++)
355 fprintf(fpsg, "%10g %10g\n", (i + 0.5) * box[slice_dim][slice_dim] / nslice,
356 sg_slice_tot[i] / static_cast<real>(nframes));
357 fprintf(fpsk, "%10g %10g\n", (i + 0.5) * box[slice_dim][slice_dim] / nslice,
358 sk_slice_tot[i] / static_cast<real>(nframes));
365 /* Print name of first atom in all groups in index file */
366 static void print_types(const int index[], int a[], int ngrps, char* groups[], const t_topology* top)
370 fprintf(stderr, "Using following groups: \n");
371 for (i = 0; i < ngrps; i++)
373 fprintf(stderr, "Groupname: %s First atomname: %s First atomnr %d\n", groups[i],
374 *(top->atoms.atomname[a[index[i]]]), a[index[i]]);
376 fprintf(stderr, "\n");
379 static void check_length(real length, int a, int b)
384 "WARNING: distance between atoms %d and "
385 "%d > 0.3 nm (%f). Index file might be corrupt.\n",
390 static void calc_order(const char* fn,
399 const t_topology* top,
403 gmx_bool permolecule,
408 const gmx_output_env_t* oenv)
410 /* if permolecule = TRUE, order parameters will be calculed per molecule
411 * and stored in slOrder with #slices = # molecules */
412 rvec *x0, /* coordinates with pbc */
413 *x1; /* coordinates without pbc */
414 matrix box; /* box (3x3) */
416 rvec cossum, /* sum of vector angles for three axes */
417 Sx, Sy, Sz, /* the three molecular axes */
418 tmp1, tmp2, /* temp. rvecs for calculating dot products */
419 frameorder; /* order parameters for one frame */
420 real* slFrameorder; /* order parameter for one frame, per slice */
421 real length, /* total distance between two atoms */
422 t, /* time from trajectory */
423 z_ave, z1, z2; /* average z, used to det. which slice atom is in */
424 int natoms, /* nr. atoms in trj */
425 nr_tails, /* nr tails, to check if index file is correct */
426 size = 0, /* nr. of atoms in group. same as nr_tails */
427 i, j, m, k, teller = 0, slice; /* current slice number */
429 int* slCount; /* nr. of atoms in one slice */
430 real sdbangle = 0; /* sum of these angles */
431 gmx_bool use_unitvector = FALSE; /* use a specified unit vector instead of axis to specify unit normal*/
433 int comsize, distsize;
434 int * comidx = nullptr, *distidx = nullptr;
435 char* grpname = nullptr;
437 real arcdist, tmpdist;
438 gmx_rmpbc_t gpbc = nullptr;
440 /* PBC added for center-of-mass vector*/
441 /* Initiate the pbc structure */
442 std::memset(&pbc, 0, sizeof(pbc));
444 if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
446 gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
449 nr_tails = index[1] - index[0];
450 fprintf(stderr, "Number of elements in first group: %d\n", nr_tails);
451 /* take first group as standard. Not rocksolid, but might catch error in index*/
456 bSliced = FALSE; /*force slices off */
457 fprintf(stderr, "Calculating order parameters for each of %d molecules\n", nslices);
462 use_unitvector = TRUE;
463 fprintf(stderr, "Select an index group to calculate the radial membrane normal\n");
464 get_index(&top->atoms, radfn, 1, &comsize, &comidx, &grpname);
468 if (grpname != nullptr)
472 fprintf(stderr, "Select an index group to use as distance reference\n");
473 get_index(&top->atoms, radfn, 1, &distsize, &distidx, &grpname);
474 bSliced = FALSE; /*force slices off*/
477 if (use_unitvector && bSliced)
480 "Warning: slicing and specified unit vectors are not currently compatible\n");
483 snew(slCount, nslices);
484 snew(*slOrder, nslices);
485 for (i = 0; i < nslices; i++)
487 snew((*slOrder)[i], ngrps);
491 snew(*distvals, nslices);
492 for (i = 0; i < nslices; i++)
494 snew((*distvals)[i], ngrps);
498 snew(slFrameorder, nslices);
503 *slWidth = box[axis][axis] / static_cast<real>(nslices);
504 fprintf(stderr, "Box divided in %d slices. Initial width of slice: %f\n", nslices, *slWidth);
509 nr_tails = index[1] - index[0];
510 fprintf(stderr, "Number of elements in first group: %d\n", nr_tails);
511 /* take first group as standard. Not rocksolid, but might catch error
517 gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
518 /*********** Start processing trajectory ***********/
523 *slWidth = box[axis][axis] / static_cast<real>(nslices);
527 set_pbc(&pbc, ePBC, box);
528 gmx_rmpbc_copy(gpbc, natoms, box, x0, x1);
530 /* Now loop over all groups. There are ngrps groups, the order parameter can
531 be calculated for grp 1 to grp ngrps - 1. For each group, loop over all
532 atoms in group, which is index[i] to (index[i+1] - 1) See block.h. Of
533 course, in this case index[i+1] -index[i] has to be the same for all
534 groups, namely the number of tails. i just runs over all atoms in a tail,
535 so for DPPC ngrps = 16 and i runs from 1 to 14, including 14
541 /*center-of-mass determination*/
545 for (j = 0; j < comsize; j++)
547 rvec_inc(com, x1[comidx[j]]);
549 svmul(1.0 / comsize, com, com);
551 rvec displacementFromReference;
554 rvec dref = { 0.0, 0.0, 0.0 };
555 for (j = 0; j < distsize; j++)
557 rvec_inc(dref, x1[distidx[j]]);
559 svmul(1.0 / distsize, dref, dref);
562 pbc_dx(&pbc, dref, com, displacementFromReference);
563 unitv(displacementFromReference, displacementFromReference);
567 for (i = 1; i < ngrps - 1; i++)
569 clear_rvec(frameorder);
571 size = index[i + 1] - index[i];
572 if (size != nr_tails)
575 "grp %d does not have same number of"
576 " elements as grp 1\n",
580 for (j = 0; j < size; j++)
583 /*create unit vector*/
585 pbc_dx(&pbc, x1[a[index[i] + j]], com, direction);
586 unitv(direction, direction);
589 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],
590 direction[0],direction[1],direction[2]);*/
596 /* Using convention for unsaturated carbons */
597 /* first get Sz, the vector from Cn to Cn+1 */
598 rvec_sub(x1[a[index[i + 1] + j]], x1[a[index[i] + j]], dist);
600 check_length(length, a[index[i] + j], a[index[i + 1] + j]);
601 svmul(1.0 / length, dist, Sz);
603 /* this is actually the cosine of the angle between the double bond
604 and axis, because Sz is normalized and the two other components of
605 the axis on the bilayer are zero */
608 sdbangle += gmx_angle(direction, Sz); /*this can probably be optimized*/
612 sdbangle += std::acos(Sz[axis]);
618 /* get vector dist(Cn-1,Cn+1) for tail atoms */
619 rvec_sub(x1[a[index[i + 1] + j]], x1[a[index[i - 1] + j]], dist);
620 length = norm(dist); /* determine distance between two atoms */
621 check_length(length, a[index[i - 1] + j], a[index[i + 1] + j]);
623 svmul(1.0 / length, dist, Sz);
624 /* Sz is now the molecular axis Sz, normalized and all that */
627 /* now get Sx. Sx is normal to the plane of Cn-1, Cn and Cn+1 so
628 we can use the outer product of Cn-1->Cn and Cn+1->Cn, I hope */
629 rvec_sub(x1[a[index[i + 1] + j]], x1[a[index[i] + j]], tmp1);
630 rvec_sub(x1[a[index[i - 1] + j]], x1[a[index[i] + j]], tmp2);
631 cprod(tmp1, tmp2, Sx);
632 svmul(1.0 / norm(Sx), Sx, Sx);
634 /* now we can get Sy from the outer product of Sx and Sz */
636 svmul(1.0 / norm(Sy), Sy, Sy);
638 /* the square of cosine of the angle between dist and the axis.
639 Using the innerproduct, but two of the three elements are zero
640 Determine the sum of the orderparameter of all atoms in group
644 cossum[XX] = gmx::square(iprod(Sx, direction)); /* this is allowed, since Sa is normalized */
645 cossum[YY] = gmx::square(iprod(Sy, direction));
646 cossum[ZZ] = gmx::square(iprod(Sz, direction));
650 cossum[XX] = gmx::square(Sx[axis]); /* this is allowed, since Sa is normalized */
651 cossum[YY] = gmx::square(Sy[axis]);
652 cossum[ZZ] = gmx::square(Sz[axis]);
655 for (m = 0; m < DIM; m++)
657 frameorder[m] += 0.5 * (3.0 * cossum[m] - 1.0);
662 /* get average coordinate in box length for slicing,
663 determine which slice atom is in, increase count for that
664 slice. slFrameorder and slOrder are reals, not
665 rvecs. Only the component [axis] of the order tensor is
666 kept, until I find it necessary to know the others too
669 z1 = x1[a[index[i - 1] + j]][axis];
670 z2 = x1[a[index[i + 1] + j]][axis];
671 z_ave = 0.5 * (z1 + z2);
672 slice = static_cast<int>((static_cast<real>(nslices) * z_ave) / box[axis][axis]);
675 slice += static_cast<real>(nslices);
677 slice = slice % nslices;
679 slCount[slice]++; /* determine slice, increase count */
681 slFrameorder[slice] += 0.5 * (3 * cossum[axis] - 1);
683 else if (permolecule)
685 /* store per-molecule order parameter
686 * To just track single-axis order: (*slOrder)[j][i] += 0.5 * (3 *
687 * iprod(cossum,direction) - 1); following is for Scd order: */
688 (*slOrder)[j][i] += -1
689 * (1.0 / 3.0 * (3 * cossum[XX] - 1)
690 + 1.0 / 3.0 * 0.5 * (3.0 * cossum[YY] - 1));
696 /* bin order parameter by arc distance from reference group*/
697 arcdist = gmx_angle(displacementFromReference, direction);
698 (*distvals)[j][i] += arcdist;
702 /* Want minimum lateral distance to first group calculated */
703 tmpdist = trace(box); /* should be max value */
704 for (k = 0; k < distsize; k++)
707 pbc_dx(&pbc, x1[distidx[k]], x1[a[index[i] + j]], displacement);
708 /* at the moment, just remove displacement[axis] */
709 displacement[axis] = 0;
710 tmpdist = std::min(tmpdist, norm2(displacement));
712 // fprintf(stderr, "Min dist %f; trace %f\n", tmpdist, trace(box));
713 (*distvals)[j][i] += std::sqrt(tmpdist);
716 } /* end loop j, over all atoms in group */
718 for (m = 0; m < DIM; m++)
720 (*order)[i][m] += (frameorder[m] / static_cast<real>(size));
724 { /*Skip following if doing per-molecule*/
725 for (k = 0; k < nslices; k++)
727 if (slCount[k]) /* if no elements, nothing has to be added */
729 (*slOrder)[k][i] += slFrameorder[k] / static_cast<real>(slCount[k]);
734 } /* end loop i, over all groups in indexfile */
738 } while (read_next_x(oenv, status, &t, x0, box));
739 /*********** done with status file **********/
741 fprintf(stderr, "\nRead trajectory. Printing parameters to file\n");
742 gmx_rmpbc_done(gpbc);
744 /* average over frames */
745 for (i = 1; i < ngrps - 1; i++)
747 svmul(1.0 / nr_frames, (*order)[i], (*order)[i]);
748 fprintf(stderr, "Atom %d Tensor: x=%g , y=%g, z=%g\n", i, (*order)[i][XX], (*order)[i][YY],
750 if (bSliced || permolecule)
752 for (k = 0; k < nslices; k++)
754 (*slOrder)[k][i] /= nr_frames;
759 for (k = 0; k < nslices; k++)
761 (*distvals)[k][i] /= nr_frames;
768 fprintf(stderr, "Average angle between double bond and normal: %f\n",
769 180 * sdbangle / (nr_frames * static_cast<real>(size) * M_PI));
772 sfree(x0); /* free memory used by coordinate arrays */
774 if (comidx != nullptr)
778 if (distidx != nullptr)
782 if (grpname != nullptr)
789 static void order_plot(rvec order[],
798 gmx_bool permolecule,
800 const gmx_output_env_t* oenv)
802 FILE *ord, *slOrd; /* xvgr files with order parameters */
803 int atom, slice; /* atom corresponding to order para.*/
804 char buf[256]; /* for xvgr title */
805 real S; /* order parameter averaged over all atoms */
809 sprintf(buf, "Scd order parameters");
810 ord = xvgropen(afile, buf, "Atom", "S", oenv);
811 sprintf(buf, "Orderparameters per atom per slice");
812 slOrd = xvgropen(bfile, buf, "Molecule", "S", oenv);
813 for (atom = 1; atom < ngrps - 1; atom++)
815 fprintf(ord, "%12d %12g\n", atom,
816 -1.0 * (2.0 / 3.0 * order[atom][XX] + 1.0 / 3.0 * order[atom][YY]));
819 for (slice = 0; slice < nslices; slice++)
821 fprintf(slOrd, "%12d\t", slice);
824 fprintf(slOrd, "%12g\t", distvals[slice][1]); /*use distance value at second carbon*/
826 for (atom = 1; atom < ngrps - 1; atom++)
828 fprintf(slOrd, "%12g\t", slOrder[slice][atom]);
830 fprintf(slOrd, "\n");
835 sprintf(buf, "Orderparameters Sz per atom");
836 ord = xvgropen(afile, buf, "Atom", "S", oenv);
837 fprintf(stderr, "ngrps = %d, nslices = %d", ngrps, nslices);
839 sprintf(buf, "Orderparameters per atom per slice");
840 slOrd = xvgropen(bfile, buf, "Slice", "S", oenv);
842 for (atom = 1; atom < ngrps - 1; atom++)
844 fprintf(ord, "%12d %12g\n", atom, order[atom][ZZ]);
847 for (slice = 0; slice < nslices; slice++)
850 for (atom = 1; atom < ngrps - 1; atom++)
852 S += slOrder[slice][atom];
854 fprintf(slOrd, "%12g %12g\n", static_cast<real>(slice) * slWidth,
855 S / static_cast<real>(atom));
860 sprintf(buf, "Order tensor diagonal elements");
861 ord = xvgropen(afile, buf, "Atom", "S", oenv);
862 sprintf(buf, "Deuterium order parameters");
863 slOrd = xvgropen(cfile, buf, "Atom", "Scd", oenv);
865 for (atom = 1; atom < ngrps - 1; atom++)
867 fprintf(ord, "%12d %12g %12g %12g\n", atom, order[atom][XX], order[atom][YY],
869 fprintf(slOrd, "%12d %12g\n", atom,
870 -1.0 * (2.0 / 3.0 * order[atom][XX] + 1.0 / 3.0 * order[atom][YY]));
877 static void write_bfactors(t_filenm* fnm,
884 const t_topology* top,
886 gmx_output_env_t* oenv)
888 /*function to write order parameters as B factors in PDB file using
889 first frame of trajectory*/
891 t_trxframe fr, frout;
895 ngrps -= 2; /*we don't have an order parameter for the first or
896 last atom in each chain*/
897 nout = nslices * ngrps;
898 read_first_frame(oenv, &status, ftp2fn(efTRX, nfile, fnm), &fr, TRX_NEED_X);
908 init_t_atoms(&useatoms, nout, TRUE);
911 /*initialize PDBinfo*/
912 for (i = 0; i < useatoms.nr; ++i)
914 useatoms.pdbinfo[i].type = 0;
915 useatoms.pdbinfo[i].occup = 0.0;
916 useatoms.pdbinfo[i].bfac = 0.0;
917 useatoms.pdbinfo[i].bAnisotropic = FALSE;
920 for (j = 0, ctr = 0; j < nslices; j++)
922 for (i = 0; i < ngrps; i++, ctr++)
924 /*iterate along each chain*/
925 useatoms.pdbinfo[ctr].bfac = order[j][i + 1];
928 useatoms.pdbinfo[ctr].occup = distvals[j][i + 1];
930 copy_rvec(fr.x[a[index[i + 1] + j]], frout.x[ctr]);
931 useatoms.atomname[ctr] = top->atoms.atomname[a[index[i + 1] + j]];
932 useatoms.atom[ctr] = top->atoms.atom[a[index[i + 1] + j]];
933 useatoms.nres = std::max(useatoms.nres, useatoms.atom[ctr].resind + 1);
934 useatoms.resinfo[useatoms.atom[ctr].resind] =
935 top->atoms.resinfo[useatoms.atom[ctr].resind]; /*copy resinfo*/
939 write_sto_conf(opt2fn("-ob", nfile, fnm), "Order parameters", &useatoms, frout.x, nullptr,
940 frout.ePBC, frout.box);
943 done_atom(&useatoms);
946 int gmx_order(int argc, char* argv[])
948 const char* desc[] = {
949 "[THISMODULE] computes the order parameter per atom for carbon tails. For atom i the",
950 "vector i-1, i+1 is used together with an axis. ",
951 "The index file should contain only the groups to be used for calculations,",
952 "with each group of equivalent carbons along the relevant acyl chain in its own",
953 "group. There should not be any generic groups (like System, Protein) in the index",
954 "file to avoid confusing the program (this is not relevant to tetrahedral order",
955 "parameters however, which only work for water anyway).[PAR]",
956 "[THISMODULE] can also give all",
957 "diagonal elements of the order tensor and even calculate the deuterium",
958 "order parameter Scd (default). If the option [TT]-szonly[tt] is given, only one",
959 "order tensor component (specified by the [TT]-d[tt] option) is given and the",
960 "order parameter per slice is calculated as well. If [TT]-szonly[tt] is not",
961 "selected, all diagonal elements and the deuterium order parameter is",
963 "The tetrahedrality order parameters can be determined",
964 "around an atom. Both angle an distance order parameters are calculated. See",
965 "P.-L. Chau and A.J. Hardwick, Mol. Phys., 93, (1998), 511-518.",
969 static int nslices = 1; /* nr of slices defined */
970 static gmx_bool bSzonly = FALSE; /* True if only Sz is wanted */
971 static gmx_bool bUnsat = FALSE; /* True if carbons are unsat. */
972 static const char* normal_axis[] = { nullptr, "z", "x", "y", nullptr };
973 static gmx_bool permolecule = FALSE; /*compute on a per-molecule basis */
974 static gmx_bool radial = FALSE; /*compute a radial membrane normal */
975 static gmx_bool distcalc = FALSE; /*calculate distance from a reference group */
977 { "-d", FALSE, etENUM, { normal_axis }, "Direction of the normal on the membrane" },
982 "Calculate order parameter as function of box length, dividing the box"
983 " into this number of slices." },
988 "Only give Sz element of order tensor. (axis can be specified with [TT]-d[tt])" },
993 "Calculate order parameters for unsaturated carbons. Note that this can"
994 "not be mixed with normal order parameters." },
999 "Compute per-molecule Scd order parameters" },
1000 { "-radial", FALSE, etBOOL, { &radial }, "Compute a radial membrane normal" },
1001 { "-calcdist", FALSE, etBOOL, { &distcalc }, "Compute distance from a reference" },
1004 rvec* order; /* order par. for each atom */
1005 real** slOrder; /* same, per slice */
1006 real slWidth = 0.0; /* width of a slice */
1007 char** grpname; /* groupnames */
1008 int ngrps, /* nr. of groups */
1009 i, axis = 0; /* normal axis */
1010 t_topology* top; /* topology */
1012 int * index, /* indices for a */
1013 *a; /* atom numbers in each group */
1014 t_blocka* block; /* data from index file */
1016 /* files for g_order */
1017 { efTRX, "-f", nullptr, ffREAD }, /* trajectory file */
1018 { efNDX, "-n", nullptr, ffREAD }, /* index file */
1019 { efNDX, "-nr", nullptr, ffOPTRD }, /* index for radial axis calculation */
1020 { efTPR, nullptr, nullptr, ffREAD }, /* topology file */
1021 { efXVG, "-o", "order", ffWRITE }, /* xvgr output file */
1022 { efXVG, "-od", "deuter", ffWRITE }, /* xvgr output file */
1023 { efPDB, "-ob", nullptr, ffOPTWR }, /* write Scd as B factors to PDB if permolecule */
1024 { efXVG, "-os", "sliced", ffWRITE }, /* xvgr output file */
1025 { efXVG, "-Sg", "sg-ang", ffOPTWR }, /* xvgr output file */
1026 { efXVG, "-Sk", "sk-dist", ffOPTWR }, /* xvgr output file */
1027 { efXVG, "-Sgsl", "sg-ang-slice", ffOPTWR }, /* xvgr output file */
1028 { efXVG, "-Sksl", "sk-dist-slice", ffOPTWR }, /* xvgr output file */
1030 gmx_bool bSliced = FALSE; /* True if box is sliced */
1031 #define NFILE asize(fnm)
1032 real** distvals = nullptr;
1033 const char * sgfnm, *skfnm, *ndxfnm, *tpsfnm, *trxfnm;
1034 gmx_output_env_t* oenv;
1036 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME, NFILE, fnm, asize(pa), pa,
1037 asize(desc), desc, 0, nullptr, &oenv))
1043 gmx_fatal(FARGS, "Can not have nslices < 1");
1045 sgfnm = opt2fn_null("-Sg", NFILE, fnm);
1046 skfnm = opt2fn_null("-Sk", NFILE, fnm);
1047 ndxfnm = opt2fn_null("-n", NFILE, fnm);
1048 tpsfnm = ftp2fn(efTPR, NFILE, fnm);
1049 trxfnm = ftp2fn(efTRX, NFILE, fnm);
1051 /* Calculate axis */
1052 GMX_RELEASE_ASSERT(normal_axis[0] != nullptr, "Options inconsistency; normal_axis[0] is NULL");
1053 if (std::strcmp(normal_axis[0], "x") == 0)
1057 else if (std::strcmp(normal_axis[0], "y") == 0)
1061 else if (std::strcmp(normal_axis[0], "z") == 0)
1067 gmx_fatal(FARGS, "Invalid axis, use x, y or z");
1072 case 0: fprintf(stderr, "Taking x axis as normal to the membrane\n"); break;
1073 case 1: fprintf(stderr, "Taking y axis as normal to the membrane\n"); break;
1074 case 2: fprintf(stderr, "Taking z axis as normal to the membrane\n"); break;
1077 /* tetraheder order parameter */
1080 /* If either of theoptions is set we compute both */
1081 sgfnm = opt2fn("-Sg", NFILE, fnm);
1082 skfnm = opt2fn("-Sk", NFILE, fnm);
1083 calc_tetra_order_parm(ndxfnm, tpsfnm, trxfnm, sgfnm, skfnm, nslices, axis,
1084 opt2fn("-Sgsl", NFILE, fnm), opt2fn("-Sksl", NFILE, fnm), oenv);
1085 /* view xvgr files */
1086 do_view(oenv, opt2fn("-Sg", NFILE, fnm), nullptr);
1087 do_view(oenv, opt2fn("-Sk", NFILE, fnm), nullptr);
1090 do_view(oenv, opt2fn("-Sgsl", NFILE, fnm), nullptr);
1091 do_view(oenv, opt2fn("-Sksl", NFILE, fnm), nullptr);
1096 /* tail order parameter */
1101 fprintf(stderr, "Dividing box in %d slices.\n\n", nslices);
1106 fprintf(stderr, "Only calculating Sz\n");
1110 fprintf(stderr, "Taking carbons as unsaturated!\n");
1113 top = read_top(ftp2fn(efTPR, NFILE, fnm), &ePBC); /* read topology file */
1115 block = init_index(ftp2fn(efNDX, NFILE, fnm), &grpname);
1116 index = block->index; /* get indices from t_block block */
1117 a = block->a; /* see block.h */
1122 nslices = index[1] - index[0]; /*I think this assumes contiguous lipids in topology*/
1123 fprintf(stderr, "Calculating Scd order parameters for each of %d molecules\n", nslices);
1128 fprintf(stderr, "Calculating radial distances\n");
1131 gmx_fatal(FARGS, "Cannot yet output radial distances without permolecule\n");
1135 /* show atomtypes, to check if index file is correct */
1136 print_types(index, a, ngrps, grpname, top);
1138 calc_order(ftp2fn(efTRX, NFILE, fnm), index, a, &order, &slOrder, &slWidth, nslices,
1139 bSliced, bUnsat, top, ePBC, ngrps, axis, permolecule, radial, distcalc,
1140 opt2fn_null("-nr", NFILE, fnm), &distvals, oenv);
1144 ngrps--; /*don't print the last group--was used for
1145 center-of-mass determination*/
1147 order_plot(order, slOrder, opt2fn("-o", NFILE, fnm), opt2fn("-os", NFILE, fnm),
1148 opt2fn("-od", NFILE, fnm), ngrps, nslices, slWidth, bSzonly, permolecule,
1151 if (opt2bSet("-ob", NFILE, fnm))
1156 "Won't write B-factors with averaged order parameters; use -permolecule\n");
1160 write_bfactors(fnm, NFILE, index, a, nslices, ngrps, slOrder, top, distvals, oenv);
1165 do_view(oenv, opt2fn("-o", NFILE, fnm), nullptr); /* view xvgr file */
1166 do_view(oenv, opt2fn("-os", NFILE, fnm), nullptr); /* view xvgr file */
1167 do_view(oenv, opt2fn("-od", NFILE, fnm), nullptr); /* view xvgr file */
1170 if (distvals != nullptr)
1172 for (i = 0; i < nslices; ++i)