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 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
45 #include "gromacs/commandline/pargs.h"
46 #include "gromacs/commandline/viewit.h"
47 #include "gromacs/fileio/confio.h"
48 #include "gromacs/fileio/trxio.h"
49 #include "gromacs/fileio/xvgr.h"
50 #include "gromacs/gmxana/cmat.h"
51 #include "gromacs/gmxana/gmx_ana.h"
52 #include "gromacs/gmxana/gstat.h"
53 #include "gromacs/math/functions.h"
54 #include "gromacs/math/utilities.h"
55 #include "gromacs/math/vec.h"
56 #include "gromacs/pbcutil/pbc.h"
57 #include "gromacs/pbcutil/rmpbc.h"
58 #include "gromacs/topology/index.h"
59 #include "gromacs/topology/topology.h"
60 #include "gromacs/trajectory/trajectoryframe.h"
61 #include "gromacs/utility/arraysize.h"
62 #include "gromacs/utility/cstringutil.h"
63 #include "gromacs/utility/fatalerror.h"
64 #include "gromacs/utility/futil.h"
65 #include "gromacs/utility/gmxassert.h"
66 #include "gromacs/utility/smalloc.h"
68 /****************************************************************************/
69 /* This program calculates the order parameter per atom for an interface or */
70 /* bilayer, averaged over time. */
71 /* S = 1/2 * (3 * cos(i)cos(j) - delta(ij)) */
72 /* It is assumed that the order parameter with respect to a box-axis */
73 /* is calculated. In that case i = j = axis, and delta(ij) = 1. */
75 /* Peter Tieleman, April 1995 */
76 /* P.J. van Maaren, November 2005 Added tetrahedral stuff */
77 /****************************************************************************/
79 static void find_nearest_neighbours(PbcType pbcType,
93 int ix, jx, nsgbin, *sgbin;
94 int i, ibin, j, k, l, n, *nn[4];
95 rvec dx, rj, rk, urk, urj;
96 real cost, cost2, *sgmol, *skmol, rmean, rmean2, r2, box2, *r_nn[4];
100 real onethird = 1.0 / 3.0;
101 /* dmat = init_mat(maxidx, FALSE); */
102 box2 = box[XX][XX] * box[XX][XX];
103 snew(sl_count, nslice);
104 for (i = 0; (i < 4); i++)
106 snew(r_nn[i], natoms);
109 for (j = 0; (j < natoms); j++)
118 /* Must init pbc every step because of pressure coupling */
119 set_pbc(&pbc, pbcType, box);
121 gmx_rmpbc(gpbc, natoms, box, x);
123 nsgbin = 2001; // Calculated as (1 + 1/0.0005)
129 for (i = 0; (i < maxidx); i++) /* loop over index file */
132 for (j = 0; (j < maxidx); j++)
141 pbc_dx(&pbc, x[ix], x[jx], dx);
144 /* set_mat_entry(dmat,i,j,r2); */
146 /* determine the nearest neighbours */
149 r_nn[3][i] = r_nn[2][i];
151 r_nn[2][i] = r_nn[1][i];
153 r_nn[1][i] = r_nn[0][i];
158 else if (r2 < r_nn[1][i])
160 r_nn[3][i] = r_nn[2][i];
162 r_nn[2][i] = r_nn[1][i];
167 else if (r2 < r_nn[2][i])
169 r_nn[3][i] = r_nn[2][i];
174 else if (r2 < r_nn[3][i])
182 /* calculate mean distance between nearest neighbours */
184 for (j = 0; (j < 4); j++)
186 r_nn[j][i] = std::sqrt(r_nn[j][i]);
195 /* Chau1998a eqn 3 */
196 /* angular part tetrahedrality order parameter per atom */
197 for (j = 0; (j < 3); j++)
199 for (k = j + 1; (k < 4); k++)
201 pbc_dx(&pbc, x[ix], x[index[nn[k][i]]], rk);
202 pbc_dx(&pbc, x[ix], x[index[nn[j][i]]], rj);
207 cost = iprod(urk, urj) + onethird;
210 /* sgmol[i] += 3*cost2/32; */
213 /* determine distribution */
214 ibin = static_cast<int>(static_cast<real>(nsgbin) * cost2);
219 /* printf("%d %d %f %d %d\n", j, k, cost * cost, ibin, sgbin[ibin]);*/
225 /* normalize sgmol between 0.0 and 1.0 */
226 sgmol[i] = 3 * sgmol[i] / 32;
229 /* distance part tetrahedrality order parameter per atom */
230 rmean2 = 4 * 3 * rmean * rmean;
231 for (j = 0; (j < 4); j++)
233 skmol[i] += (rmean - r_nn[j][i]) * (rmean - r_nn[j][i]) / rmean2;
234 /* printf("%d %f (%f %f %f %f) \n",
235 i, skmol[i], rmean, rmean2, r_nn[j][i], (rmean - r_nn[j][i]) );
241 /* Compute sliced stuff */
242 sl_index = static_cast<int>(std::round((1 + x[i][slice_dim] / box[slice_dim][slice_dim])
243 * static_cast<real>(nslice)))
245 sgslice[sl_index] += sgmol[i];
246 skslice[sl_index] += skmol[i];
247 sl_count[sl_index]++;
248 } /* loop over entries in index file */
250 *sgmean /= static_cast<real>(maxidx);
251 *skmean /= static_cast<real>(maxidx);
253 for (i = 0; (i < nslice); i++)
257 sgslice[i] /= sl_count[i];
258 skslice[i] /= sl_count[i];
265 for (i = 0; (i < 4); i++)
273 static void calc_tetra_order_parm(const char* fnNDX,
282 const gmx_output_env_t* oenv)
284 FILE * fpsg = nullptr, *fpsk = nullptr;
295 int i, *isize, ng, nframes;
296 real * sg_slice, *sg_slice_tot, *sk_slice, *sk_slice_tot;
297 gmx_rmpbc_t gpbc = nullptr;
300 read_tps_conf(fnTPS, &top, &pbcType, &xtop, nullptr, box, FALSE);
302 snew(sg_slice, nslice);
303 snew(sk_slice, nslice);
304 snew(sg_slice_tot, nslice);
305 snew(sk_slice_tot, nslice);
307 /* get index groups */
308 printf("Select the group that contains the atoms you want to use for the tetrahedrality order "
309 "parameter calculation:\n");
313 get_index(&top.atoms, fnNDX, ng, isize, index, grpname);
315 /* Analyze trajectory */
316 natoms = read_first_x(oenv, &status, fnTRX, &t, &x, box);
317 if (natoms > top.atoms.nr)
319 gmx_fatal(FARGS, "Topology (%d atoms) does not match trajectory (%d atoms)", top.atoms.nr, natoms);
321 check_index(nullptr, ng, index[0], nullptr, natoms);
323 fpsg = xvgropen(sgfn, "S\\sg\\N Angle Order Parameter", "Time (ps)", "S\\sg\\N", oenv);
324 fpsk = xvgropen(skfn, "S\\sk\\N Distance Order Parameter", "Time (ps)", "S\\sk\\N", oenv);
326 /* loop over frames */
327 gpbc = gmx_rmpbc_init(&top.idef, pbcType, natoms);
331 find_nearest_neighbours(
332 pbcType, natoms, box, x, isize[0], index[0], &sg, &sk, nslice, slice_dim, sg_slice, sk_slice, gpbc);
333 for (i = 0; (i < nslice); i++)
335 sg_slice_tot[i] += sg_slice[i];
336 sk_slice_tot[i] += sk_slice[i];
338 fprintf(fpsg, "%f %f\n", t, sg);
339 fprintf(fpsk, "%f %f\n", t, sk);
341 } while (read_next_x(oenv, status, &t, x, box));
343 gmx_rmpbc_done(gpbc);
352 fpsg = xvgropen(sgslfn, "S\\sg\\N Angle Order Parameter / Slab", "(nm)", "S\\sg\\N", oenv);
353 fpsk = xvgropen(skslfn, "S\\sk\\N Distance Order Parameter / Slab", "(nm)", "S\\sk\\N", oenv);
354 for (i = 0; (i < nslice); i++)
358 (i + 0.5) * box[slice_dim][slice_dim] / nslice,
359 sg_slice_tot[i] / static_cast<real>(nframes));
362 (i + 0.5) * box[slice_dim][slice_dim] / nslice,
363 sk_slice_tot[i] / static_cast<real>(nframes));
370 /* Print name of first atom in all groups in index file */
371 static void print_types(const int index[], int a[], int ngrps, char* groups[], const t_topology* top)
375 fprintf(stderr, "Using following groups: \n");
376 for (i = 0; i < ngrps; i++)
379 "Groupname: %s First atomname: %s First atomnr %d\n",
381 *(top->atoms.atomname[a[index[i]]]),
384 fprintf(stderr, "\n");
387 static void check_length(real length, int a, int b)
392 "WARNING: distance between atoms %d and "
393 "%d > 0.3 nm (%f). Index file might be corrupt.\n",
400 static void calc_order(const char* fn,
409 const t_topology* top,
413 gmx_bool permolecule,
418 const gmx_output_env_t* oenv)
420 /* if permolecule = TRUE, order parameters will be calculed per molecule
421 * and stored in slOrder with #slices = # molecules */
422 rvec *x0, /* coordinates with pbc */
423 *x1; /* coordinates without pbc */
424 matrix box; /* box (3x3) */
426 rvec cossum, /* sum of vector angles for three axes */
427 Sx, Sy, Sz, /* the three molecular axes */
428 tmp1, tmp2, /* temp. rvecs for calculating dot products */
429 frameorder; /* order parameters for one frame */
430 real* slFrameorder; /* order parameter for one frame, per slice */
431 real length, /* total distance between two atoms */
432 t, /* time from trajectory */
433 z_ave, z1, z2; /* average z, used to det. which slice atom is in */
434 int natoms, /* nr. atoms in trj */
435 nr_tails, /* nr tails, to check if index file is correct */
436 size = 0, /* nr. of atoms in group. same as nr_tails */
437 i, j, m, k, teller = 0, slice; /* current slice number */
439 int* slCount; /* nr. of atoms in one slice */
440 real sdbangle = 0; /* sum of these angles */
441 gmx_bool use_unitvector = FALSE; /* use a specified unit vector instead of axis to specify unit normal*/
443 int comsize, distsize;
444 int * comidx = nullptr, *distidx = nullptr;
445 char* grpname = nullptr;
447 real arcdist, tmpdist;
448 gmx_rmpbc_t gpbc = nullptr;
450 /* PBC added for center-of-mass vector*/
451 /* Initiate the pbc structure */
452 std::memset(&pbc, 0, sizeof(pbc));
454 if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
456 gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
459 nr_tails = index[1] - index[0];
460 fprintf(stderr, "Number of elements in first group: %d\n", nr_tails);
461 /* take first group as standard. Not rocksolid, but might catch error in index*/
466 bSliced = FALSE; /*force slices off */
467 fprintf(stderr, "Calculating order parameters for each of %d molecules\n", nslices);
472 use_unitvector = TRUE;
473 fprintf(stderr, "Select an index group to calculate the radial membrane normal\n");
474 get_index(&top->atoms, radfn, 1, &comsize, &comidx, &grpname);
478 if (grpname != nullptr)
482 fprintf(stderr, "Select an index group to use as distance reference\n");
483 get_index(&top->atoms, radfn, 1, &distsize, &distidx, &grpname);
484 bSliced = FALSE; /*force slices off*/
487 if (use_unitvector && bSliced)
490 "Warning: slicing and specified unit vectors are not currently compatible\n");
493 snew(slCount, nslices);
494 snew(*slOrder, nslices);
495 for (i = 0; i < nslices; i++)
497 snew((*slOrder)[i], ngrps);
501 snew(*distvals, nslices);
502 for (i = 0; i < nslices; i++)
504 snew((*distvals)[i], ngrps);
508 snew(slFrameorder, nslices);
513 *slWidth = box[axis][axis] / static_cast<real>(nslices);
514 fprintf(stderr, "Box divided in %d slices. Initial width of slice: %f\n", nslices, *slWidth);
519 nr_tails = index[1] - index[0];
520 fprintf(stderr, "Number of elements in first group: %d\n", nr_tails);
521 /* take first group as standard. Not rocksolid, but might catch error
527 gpbc = gmx_rmpbc_init(&top->idef, pbcType, natoms);
528 /*********** Start processing trajectory ***********/
533 *slWidth = box[axis][axis] / static_cast<real>(nslices);
537 set_pbc(&pbc, pbcType, box);
538 gmx_rmpbc_copy(gpbc, natoms, box, x0, x1);
540 /* Now loop over all groups. There are ngrps groups, the order parameter can
541 be calculated for grp 1 to grp ngrps - 1. For each group, loop over all
542 atoms in group, which is index[i] to (index[i+1] - 1) See block.h. Of
543 course, in this case index[i+1] -index[i] has to be the same for all
544 groups, namely the number of tails. i just runs over all atoms in a tail,
545 so for DPPC ngrps = 16 and i runs from 1 to 14, including 14
551 /*center-of-mass determination*/
555 for (j = 0; j < comsize; j++)
557 rvec_inc(com, x1[comidx[j]]);
559 svmul(1.0 / comsize, com, com);
561 rvec displacementFromReference;
564 rvec dref = { 0.0, 0.0, 0.0 };
565 for (j = 0; j < distsize; j++)
567 rvec_inc(dref, x1[distidx[j]]);
569 svmul(1.0 / distsize, dref, dref);
572 pbc_dx(&pbc, dref, com, displacementFromReference);
573 unitv(displacementFromReference, displacementFromReference);
577 for (i = 1; i < ngrps - 1; i++)
579 clear_rvec(frameorder);
581 size = index[i + 1] - index[i];
582 if (size != nr_tails)
585 "grp %d does not have same number of"
586 " elements as grp 1\n",
590 for (j = 0; j < size; j++)
593 /*create unit vector*/
595 pbc_dx(&pbc, x1[a[index[i] + j]], com, direction);
596 unitv(direction, direction);
599 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],
600 direction[0],direction[1],direction[2]);*/
606 /* Using convention for unsaturated carbons */
607 /* first get Sz, the vector from Cn to Cn+1 */
608 rvec_sub(x1[a[index[i + 1] + j]], x1[a[index[i] + j]], dist);
610 check_length(length, a[index[i] + j], a[index[i + 1] + j]);
611 svmul(1.0 / length, dist, Sz);
613 /* this is actually the cosine of the angle between the double bond
614 and axis, because Sz is normalized and the two other components of
615 the axis on the bilayer are zero */
618 sdbangle += gmx_angle(direction, Sz); /*this can probably be optimized*/
622 sdbangle += std::acos(Sz[axis]);
628 /* get vector dist(Cn-1,Cn+1) for tail atoms */
629 rvec_sub(x1[a[index[i + 1] + j]], x1[a[index[i - 1] + j]], dist);
630 length = norm(dist); /* determine distance between two atoms */
631 check_length(length, a[index[i - 1] + j], a[index[i + 1] + j]);
633 svmul(1.0 / length, dist, Sz);
634 /* Sz is now the molecular axis Sz, normalized and all that */
637 /* now get Sx. Sx is normal to the plane of Cn-1, Cn and Cn+1 so
638 we can use the outer product of Cn-1->Cn and Cn+1->Cn, I hope */
639 rvec_sub(x1[a[index[i + 1] + j]], x1[a[index[i] + j]], tmp1);
640 rvec_sub(x1[a[index[i - 1] + j]], x1[a[index[i] + j]], tmp2);
641 cprod(tmp1, tmp2, Sx);
642 svmul(1.0 / norm(Sx), Sx, Sx);
644 /* now we can get Sy from the outer product of Sx and Sz */
646 svmul(1.0 / norm(Sy), Sy, Sy);
648 /* the square of cosine of the angle between dist and the axis.
649 Using the innerproduct, but two of the three elements are zero
650 Determine the sum of the orderparameter of all atoms in group
654 cossum[XX] = gmx::square(iprod(Sx, direction)); /* this is allowed, since Sa is normalized */
655 cossum[YY] = gmx::square(iprod(Sy, direction));
656 cossum[ZZ] = gmx::square(iprod(Sz, direction));
660 cossum[XX] = gmx::square(Sx[axis]); /* this is allowed, since Sa is normalized */
661 cossum[YY] = gmx::square(Sy[axis]);
662 cossum[ZZ] = gmx::square(Sz[axis]);
665 for (m = 0; m < DIM; m++)
667 frameorder[m] += 0.5 * (3.0 * cossum[m] - 1.0);
672 /* get average coordinate in box length for slicing,
673 determine which slice atom is in, increase count for that
674 slice. slFrameorder and slOrder are reals, not
675 rvecs. Only the component [axis] of the order tensor is
676 kept, until I find it necessary to know the others too
679 z1 = x1[a[index[i - 1] + j]][axis];
680 z2 = x1[a[index[i + 1] + j]][axis];
681 z_ave = 0.5 * (z1 + z2);
682 slice = static_cast<int>((static_cast<real>(nslices) * z_ave) / box[axis][axis]);
685 slice += static_cast<real>(nslices);
687 slice = slice % nslices;
689 slCount[slice]++; /* determine slice, increase count */
691 slFrameorder[slice] += 0.5 * (3 * cossum[axis] - 1);
693 else if (permolecule)
695 /* store per-molecule order parameter
696 * To just track single-axis order: (*slOrder)[j][i] += 0.5 * (3 *
697 * iprod(cossum,direction) - 1); following is for Scd order: */
698 (*slOrder)[j][i] += -1
699 * (1.0 / 3.0 * (3 * cossum[XX] - 1)
700 + 1.0 / 3.0 * 0.5 * (3.0 * cossum[YY] - 1));
706 /* bin order parameter by arc distance from reference group*/
707 arcdist = gmx_angle(displacementFromReference, direction);
708 (*distvals)[j][i] += arcdist;
712 /* Want minimum lateral distance to first group calculated */
713 tmpdist = trace(box); /* should be max value */
714 for (k = 0; k < distsize; k++)
717 pbc_dx(&pbc, x1[distidx[k]], x1[a[index[i] + j]], displacement);
718 /* at the moment, just remove displacement[axis] */
719 displacement[axis] = 0;
720 tmpdist = std::min(tmpdist, norm2(displacement));
722 // fprintf(stderr, "Min dist %f; trace %f\n", tmpdist, trace(box));
723 (*distvals)[j][i] += std::sqrt(tmpdist);
726 } /* end loop j, over all atoms in group */
728 for (m = 0; m < DIM; m++)
730 (*order)[i][m] += (frameorder[m] / static_cast<real>(size));
734 { /*Skip following if doing per-molecule*/
735 for (k = 0; k < nslices; k++)
737 if (slCount[k]) /* if no elements, nothing has to be added */
739 (*slOrder)[k][i] += slFrameorder[k] / static_cast<real>(slCount[k]);
744 } /* end loop i, over all groups in indexfile */
748 } while (read_next_x(oenv, status, &t, x0, box));
749 /*********** done with status file **********/
751 fprintf(stderr, "\nRead trajectory. Printing parameters to file\n");
752 gmx_rmpbc_done(gpbc);
754 /* average over frames */
755 for (i = 1; i < ngrps - 1; i++)
757 svmul(1.0 / nr_frames, (*order)[i], (*order)[i]);
759 "Atom %d Tensor: x=%g , y=%g, z=%g\n",
764 if (bSliced || permolecule)
766 for (k = 0; k < nslices; k++)
768 (*slOrder)[k][i] /= nr_frames;
773 for (k = 0; k < nslices; k++)
775 (*distvals)[k][i] /= nr_frames;
783 "Average angle between double bond and normal: %f\n",
784 180 * sdbangle / (nr_frames * static_cast<real>(size) * M_PI));
787 sfree(x0); /* free memory used by coordinate arrays */
789 if (comidx != nullptr)
793 if (distidx != nullptr)
797 if (grpname != nullptr)
804 static void order_plot(rvec order[],
813 gmx_bool permolecule,
815 const gmx_output_env_t* oenv)
817 FILE *ord, *slOrd; /* xvgr files with order parameters */
818 int atom, slice; /* atom corresponding to order para.*/
819 char buf[256]; /* for xvgr title */
820 real S; /* order parameter averaged over all atoms */
824 sprintf(buf, "Scd order parameters");
825 ord = xvgropen(afile, buf, "Atom", "S", oenv);
826 sprintf(buf, "Orderparameters per atom per slice");
827 slOrd = xvgropen(bfile, buf, "Molecule", "S", oenv);
828 for (atom = 1; atom < ngrps - 1; atom++)
833 -1.0 * (2.0 / 3.0 * order[atom][XX] + 1.0 / 3.0 * order[atom][YY]));
836 for (slice = 0; slice < nslices; slice++)
838 fprintf(slOrd, "%12d\t", slice);
841 fprintf(slOrd, "%12g\t", distvals[slice][1]); /*use distance value at second carbon*/
843 for (atom = 1; atom < ngrps - 1; atom++)
845 fprintf(slOrd, "%12g\t", slOrder[slice][atom]);
847 fprintf(slOrd, "\n");
852 sprintf(buf, "Orderparameters Sz per atom");
853 ord = xvgropen(afile, buf, "Atom", "S", oenv);
854 fprintf(stderr, "ngrps = %d, nslices = %d", ngrps, nslices);
856 sprintf(buf, "Orderparameters per atom per slice");
857 slOrd = xvgropen(bfile, buf, "Slice", "S", oenv);
859 for (atom = 1; atom < ngrps - 1; atom++)
861 fprintf(ord, "%12d %12g\n", atom, order[atom][ZZ]);
864 for (slice = 0; slice < nslices; slice++)
867 for (atom = 1; atom < ngrps - 1; atom++)
869 S += slOrder[slice][atom];
871 fprintf(slOrd, "%12g %12g\n", static_cast<real>(slice) * slWidth, S / static_cast<real>(atom));
876 sprintf(buf, "Order tensor diagonal elements");
877 ord = xvgropen(afile, buf, "Atom", "S", oenv);
878 sprintf(buf, "Deuterium order parameters");
879 slOrd = xvgropen(cfile, buf, "Atom", "Scd", oenv);
881 for (atom = 1; atom < ngrps - 1; atom++)
883 fprintf(ord, "%12d %12g %12g %12g\n", atom, order[atom][XX], order[atom][YY], order[atom][ZZ]);
887 -1.0 * (2.0 / 3.0 * order[atom][XX] + 1.0 / 3.0 * order[atom][YY]));
894 static void write_bfactors(t_filenm* fnm,
901 const t_topology* top,
903 gmx_output_env_t* oenv)
905 /*function to write order parameters as B factors in PDB file using
906 first frame of trajectory*/
908 t_trxframe fr, frout;
912 ngrps -= 2; /*we don't have an order parameter for the first or
913 last atom in each chain*/
914 nout = nslices * ngrps;
915 read_first_frame(oenv, &status, ftp2fn(efTRX, nfile, fnm), &fr, TRX_NEED_X);
925 init_t_atoms(&useatoms, nout, TRUE);
928 /*initialize PDBinfo*/
929 for (i = 0; i < useatoms.nr; ++i)
931 useatoms.pdbinfo[i].type = PdbRecordType::Atom;
932 useatoms.pdbinfo[i].occup = 0.0;
933 useatoms.pdbinfo[i].bfac = 0.0;
934 useatoms.pdbinfo[i].bAnisotropic = FALSE;
937 for (j = 0, ctr = 0; j < nslices; j++)
939 for (i = 0; i < ngrps; i++, ctr++)
941 /*iterate along each chain*/
942 useatoms.pdbinfo[ctr].bfac = order[j][i + 1];
945 useatoms.pdbinfo[ctr].occup = distvals[j][i + 1];
947 copy_rvec(fr.x[a[index[i + 1] + j]], frout.x[ctr]);
948 useatoms.atomname[ctr] = top->atoms.atomname[a[index[i + 1] + j]];
949 useatoms.atom[ctr] = top->atoms.atom[a[index[i + 1] + j]];
950 useatoms.nres = std::max(useatoms.nres, useatoms.atom[ctr].resind + 1);
951 useatoms.resinfo[useatoms.atom[ctr].resind] =
952 top->atoms.resinfo[useatoms.atom[ctr].resind]; /*copy resinfo*/
957 opt2fn("-ob", nfile, fnm), "Order parameters", &useatoms, frout.x, nullptr, frout.pbcType, frout.box);
960 done_atom(&useatoms);
963 int gmx_order(int argc, char* argv[])
965 const char* desc[] = {
966 "[THISMODULE] computes the order parameter per atom for carbon tails. For atom i the",
967 "vector i-1, i+1 is used together with an axis. ",
968 "The index file should contain only the groups to be used for calculations,",
969 "with each group of equivalent carbons along the relevant acyl chain in its own",
970 "group. There should not be any generic groups (like System, Protein) in the index",
971 "file to avoid confusing the program (this is not relevant to tetrahedral order",
972 "parameters however, which only work for water anyway).[PAR]",
973 "[THISMODULE] can also give all",
974 "diagonal elements of the order tensor and even calculate the deuterium",
975 "order parameter Scd (default). If the option [TT]-szonly[tt] is given, only one",
976 "order tensor component (specified by the [TT]-d[tt] option) is given and the",
977 "order parameter per slice is calculated as well. If [TT]-szonly[tt] is not",
978 "selected, all diagonal elements and the deuterium order parameter is",
980 "The tetrahedrality order parameters can be determined",
981 "around an atom. Both angle an distance order parameters are calculated. See",
982 "P.-L. Chau and A.J. Hardwick, Mol. Phys., 93, (1998), 511-518.",
986 static int nslices = 1; /* nr of slices defined */
987 static gmx_bool bSzonly = FALSE; /* True if only Sz is wanted */
988 static gmx_bool bUnsat = FALSE; /* True if carbons are unsat. */
989 static const char* normal_axis[] = { nullptr, "z", "x", "y", nullptr };
990 static gmx_bool permolecule = FALSE; /*compute on a per-molecule basis */
991 static gmx_bool radial = FALSE; /*compute a radial membrane normal */
992 static gmx_bool distcalc = FALSE; /*calculate distance from a reference group */
994 { "-d", FALSE, etENUM, { normal_axis }, "Direction of the normal on the membrane" },
999 "Calculate order parameter as function of box length, dividing the box"
1000 " into this number of slices." },
1005 "Only give Sz element of order tensor. (axis can be specified with [TT]-d[tt])" },
1010 "Calculate order parameters for unsaturated carbons. Note that this can"
1011 "not be mixed with normal order parameters." },
1016 "Compute per-molecule Scd order parameters" },
1017 { "-radial", FALSE, etBOOL, { &radial }, "Compute a radial membrane normal" },
1018 { "-calcdist", FALSE, etBOOL, { &distcalc }, "Compute distance from a reference" },
1021 rvec* order; /* order par. for each atom */
1022 real** slOrder; /* same, per slice */
1023 real slWidth = 0.0; /* width of a slice */
1024 char** grpname; /* groupnames */
1025 int ngrps, /* nr. of groups */
1026 i, axis = 0; /* normal axis */
1027 t_topology* top; /* topology */
1028 PbcType pbcType; /* type of periodic boundary conditions */
1029 int * index, /* indices for a */
1030 *a; /* atom numbers in each group */
1031 t_blocka* block; /* data from index file */
1033 /* files for g_order */
1034 { efTRX, "-f", nullptr, ffREAD }, /* trajectory file */
1035 { efNDX, "-n", nullptr, ffREAD }, /* index file */
1036 { efNDX, "-nr", nullptr, ffOPTRD }, /* index for radial axis calculation */
1037 { efTPR, nullptr, nullptr, ffREAD }, /* topology file */
1038 { efXVG, "-o", "order", ffWRITE }, /* xvgr output file */
1039 { efXVG, "-od", "deuter", ffWRITE }, /* xvgr output file */
1040 { efPDB, "-ob", nullptr, ffOPTWR }, /* write Scd as B factors to PDB if permolecule */
1041 { efXVG, "-os", "sliced", ffWRITE }, /* xvgr output file */
1042 { efXVG, "-Sg", "sg-ang", ffOPTWR }, /* xvgr output file */
1043 { efXVG, "-Sk", "sk-dist", ffOPTWR }, /* xvgr output file */
1044 { efXVG, "-Sgsl", "sg-ang-slice", ffOPTWR }, /* xvgr output file */
1045 { efXVG, "-Sksl", "sk-dist-slice", ffOPTWR }, /* xvgr output file */
1047 gmx_bool bSliced = FALSE; /* True if box is sliced */
1048 #define NFILE asize(fnm)
1049 real** distvals = nullptr;
1050 const char * sgfnm, *skfnm, *ndxfnm, *tpsfnm, *trxfnm;
1051 gmx_output_env_t* oenv;
1053 if (!parse_common_args(
1054 &argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME, NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
1060 gmx_fatal(FARGS, "Can not have nslices < 1");
1062 sgfnm = opt2fn_null("-Sg", NFILE, fnm);
1063 skfnm = opt2fn_null("-Sk", NFILE, fnm);
1064 ndxfnm = opt2fn_null("-n", NFILE, fnm);
1065 tpsfnm = ftp2fn(efTPR, NFILE, fnm);
1066 trxfnm = ftp2fn(efTRX, NFILE, fnm);
1068 /* Calculate axis */
1069 GMX_RELEASE_ASSERT(normal_axis[0] != nullptr, "Options inconsistency; normal_axis[0] is NULL");
1070 if (std::strcmp(normal_axis[0], "x") == 0)
1074 else if (std::strcmp(normal_axis[0], "y") == 0)
1078 else if (std::strcmp(normal_axis[0], "z") == 0)
1084 gmx_fatal(FARGS, "Invalid axis, use x, y or z");
1089 case 0: fprintf(stderr, "Taking x axis as normal to the membrane\n"); break;
1090 case 1: fprintf(stderr, "Taking y axis as normal to the membrane\n"); break;
1091 case 2: fprintf(stderr, "Taking z axis as normal to the membrane\n"); break;
1094 /* tetraheder order parameter */
1097 /* If either of theoptions is set we compute both */
1098 sgfnm = opt2fn("-Sg", NFILE, fnm);
1099 skfnm = opt2fn("-Sk", NFILE, fnm);
1100 calc_tetra_order_parm(ndxfnm,
1107 opt2fn("-Sgsl", NFILE, fnm),
1108 opt2fn("-Sksl", NFILE, fnm),
1110 /* view xvgr files */
1111 do_view(oenv, opt2fn("-Sg", NFILE, fnm), nullptr);
1112 do_view(oenv, opt2fn("-Sk", NFILE, fnm), nullptr);
1115 do_view(oenv, opt2fn("-Sgsl", NFILE, fnm), nullptr);
1116 do_view(oenv, opt2fn("-Sksl", NFILE, fnm), nullptr);
1121 /* tail order parameter */
1126 fprintf(stderr, "Dividing box in %d slices.\n\n", nslices);
1131 fprintf(stderr, "Only calculating Sz\n");
1135 fprintf(stderr, "Taking carbons as unsaturated!\n");
1138 top = read_top(ftp2fn(efTPR, NFILE, fnm), &pbcType); /* read topology file */
1140 block = init_index(ftp2fn(efNDX, NFILE, fnm), &grpname);
1141 index = block->index; /* get indices from t_block block */
1142 a = block->a; /* see block.h */
1147 nslices = index[1] - index[0]; /*I think this assumes contiguous lipids in topology*/
1148 fprintf(stderr, "Calculating Scd order parameters for each of %d molecules\n", nslices);
1153 fprintf(stderr, "Calculating radial distances\n");
1156 gmx_fatal(FARGS, "Cannot yet output radial distances without permolecule\n");
1160 /* show atomtypes, to check if index file is correct */
1161 print_types(index, a, ngrps, grpname, top);
1163 calc_order(ftp2fn(efTRX, NFILE, fnm),
1179 opt2fn_null("-nr", NFILE, fnm),
1185 ngrps--; /*don't print the last group--was used for
1186 center-of-mass determination*/
1190 opt2fn("-o", NFILE, fnm),
1191 opt2fn("-os", NFILE, fnm),
1192 opt2fn("-od", NFILE, fnm),
1201 if (opt2bSet("-ob", NFILE, fnm))
1206 "Won't write B-factors with averaged order parameters; use -permolecule\n");
1210 write_bfactors(fnm, NFILE, index, a, nslices, ngrps, slOrder, top, distvals, oenv);
1215 do_view(oenv, opt2fn("-o", NFILE, fnm), nullptr); /* view xvgr file */
1216 do_view(oenv, opt2fn("-os", NFILE, fnm), nullptr); /* view xvgr file */
1217 do_view(oenv, opt2fn("-od", NFILE, fnm), nullptr); /* view xvgr file */
1220 if (distvals != nullptr)
1222 for (i = 0; i < nslices; ++i)