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/units.h"
55 #include "gromacs/math/utilities.h"
56 #include "gromacs/math/vec.h"
57 #include "gromacs/pbcutil/pbc.h"
58 #include "gromacs/pbcutil/rmpbc.h"
59 #include "gromacs/topology/index.h"
60 #include "gromacs/topology/topology.h"
61 #include "gromacs/trajectory/trajectoryframe.h"
62 #include "gromacs/utility/arraysize.h"
63 #include "gromacs/utility/cstringutil.h"
64 #include "gromacs/utility/fatalerror.h"
65 #include "gromacs/utility/futil.h"
66 #include "gromacs/utility/gmxassert.h"
67 #include "gromacs/utility/smalloc.h"
69 /****************************************************************************/
70 /* This program calculates the order parameter per atom for an interface or */
71 /* bilayer, averaged over time. */
72 /* S = 1/2 * (3 * cos(i)cos(j) - delta(ij)) */
73 /* It is assumed that the order parameter with respect to a box-axis */
74 /* is calculated. In that case i = j = axis, and delta(ij) = 1. */
76 /* Peter Tieleman, April 1995 */
77 /* P.J. van Maaren, November 2005 Added tetrahedral stuff */
78 /****************************************************************************/
80 static void find_nearest_neighbours(PbcType pbcType,
94 int ix, jx, nsgbin, *sgbin;
95 int i, ibin, j, k, l, n, *nn[4];
96 rvec dx, rj, rk, urk, urj;
97 real cost, cost2, *sgmol, *skmol, rmean, rmean2, r2, box2, *r_nn[4];
101 real onethird = 1.0 / 3.0;
102 /* dmat = init_mat(maxidx, FALSE); */
103 box2 = box[XX][XX] * box[XX][XX];
104 snew(sl_count, nslice);
105 for (i = 0; (i < 4); i++)
107 snew(r_nn[i], natoms);
110 for (j = 0; (j < natoms); j++)
119 /* Must init pbc every step because of pressure coupling */
120 set_pbc(&pbc, pbcType, box);
122 gmx_rmpbc(gpbc, natoms, box, x);
124 nsgbin = 2001; // Calculated as (1 + 1/0.0005)
130 for (i = 0; (i < maxidx); i++) /* loop over index file */
133 for (j = 0; (j < maxidx); j++)
142 pbc_dx(&pbc, x[ix], x[jx], dx);
145 /* set_mat_entry(dmat,i,j,r2); */
147 /* determine the nearest neighbours */
150 r_nn[3][i] = r_nn[2][i];
152 r_nn[2][i] = r_nn[1][i];
154 r_nn[1][i] = r_nn[0][i];
159 else if (r2 < r_nn[1][i])
161 r_nn[3][i] = r_nn[2][i];
163 r_nn[2][i] = r_nn[1][i];
168 else if (r2 < r_nn[2][i])
170 r_nn[3][i] = r_nn[2][i];
175 else if (r2 < r_nn[3][i])
183 /* calculate mean distance between nearest neighbours */
185 for (j = 0; (j < 4); j++)
187 r_nn[j][i] = std::sqrt(r_nn[j][i]);
196 /* Chau1998a eqn 3 */
197 /* angular part tetrahedrality order parameter per atom */
198 for (j = 0; (j < 3); j++)
200 for (k = j + 1; (k < 4); k++)
202 pbc_dx(&pbc, x[ix], x[index[nn[k][i]]], rk);
203 pbc_dx(&pbc, x[ix], x[index[nn[j][i]]], rj);
208 cost = iprod(urk, urj) + onethird;
211 /* sgmol[i] += 3*cost2/32; */
214 /* determine distribution */
215 ibin = static_cast<int>(static_cast<real>(nsgbin) * cost2);
220 /* printf("%d %d %f %d %d\n", j, k, cost * cost, ibin, sgbin[ibin]);*/
226 /* normalize sgmol between 0.0 and 1.0 */
227 sgmol[i] = 3 * sgmol[i] / 32;
230 /* distance part tetrahedrality order parameter per atom */
231 rmean2 = 4 * 3 * rmean * rmean;
232 for (j = 0; (j < 4); j++)
234 skmol[i] += (rmean - r_nn[j][i]) * (rmean - r_nn[j][i]) / rmean2;
235 /* printf("%d %f (%f %f %f %f) \n",
236 i, skmol[i], rmean, rmean2, r_nn[j][i], (rmean - r_nn[j][i]) );
242 /* Compute sliced stuff */
243 sl_index = static_cast<int>(std::round((1 + x[i][slice_dim] / box[slice_dim][slice_dim])
244 * static_cast<real>(nslice)))
246 sgslice[sl_index] += sgmol[i];
247 skslice[sl_index] += skmol[i];
248 sl_count[sl_index]++;
249 } /* loop over entries in index file */
251 *sgmean /= static_cast<real>(maxidx);
252 *skmean /= static_cast<real>(maxidx);
254 for (i = 0; (i < nslice); i++)
258 sgslice[i] /= sl_count[i];
259 skslice[i] /= sl_count[i];
266 for (i = 0; (i < 4); i++)
274 static void calc_tetra_order_parm(const char* fnNDX,
283 const gmx_output_env_t* oenv)
285 FILE * fpsg = nullptr, *fpsk = nullptr;
296 int i, *isize, ng, nframes;
297 real * sg_slice, *sg_slice_tot, *sk_slice, *sk_slice_tot;
298 gmx_rmpbc_t gpbc = nullptr;
301 read_tps_conf(fnTPS, &top, &pbcType, &xtop, nullptr, box, FALSE);
303 snew(sg_slice, nslice);
304 snew(sk_slice, nslice);
305 snew(sg_slice_tot, nslice);
306 snew(sk_slice_tot, nslice);
308 /* get index groups */
309 printf("Select the group that contains the atoms you want to use for the tetrahedrality order "
310 "parameter calculation:\n");
314 get_index(&top.atoms, fnNDX, ng, isize, index, grpname);
316 /* Analyze trajectory */
317 natoms = read_first_x(oenv, &status, fnTRX, &t, &x, box);
318 if (natoms > top.atoms.nr)
320 gmx_fatal(FARGS, "Topology (%d atoms) does not match trajectory (%d atoms)", top.atoms.nr, natoms);
322 check_index(nullptr, ng, index[0], nullptr, natoms);
324 fpsg = xvgropen(sgfn, "S\\sg\\N Angle Order Parameter", "Time (ps)", "S\\sg\\N", oenv);
325 fpsk = xvgropen(skfn, "S\\sk\\N Distance Order Parameter", "Time (ps)", "S\\sk\\N", oenv);
327 /* loop over frames */
328 gpbc = gmx_rmpbc_init(&top.idef, pbcType, natoms);
332 find_nearest_neighbours(
333 pbcType, natoms, box, x, isize[0], index[0], &sg, &sk, nslice, slice_dim, sg_slice, sk_slice, gpbc);
334 for (i = 0; (i < nslice); i++)
336 sg_slice_tot[i] += sg_slice[i];
337 sk_slice_tot[i] += sk_slice[i];
339 fprintf(fpsg, "%f %f\n", t, sg);
340 fprintf(fpsk, "%f %f\n", t, sk);
342 } while (read_next_x(oenv, status, &t, x, box));
344 gmx_rmpbc_done(gpbc);
353 fpsg = xvgropen(sgslfn, "S\\sg\\N Angle Order Parameter / Slab", "(nm)", "S\\sg\\N", oenv);
354 fpsk = xvgropen(skslfn, "S\\sk\\N Distance Order Parameter / Slab", "(nm)", "S\\sk\\N", oenv);
355 for (i = 0; (i < nslice); i++)
359 (i + 0.5) * box[slice_dim][slice_dim] / nslice,
360 sg_slice_tot[i] / static_cast<real>(nframes));
363 (i + 0.5) * box[slice_dim][slice_dim] / nslice,
364 sk_slice_tot[i] / static_cast<real>(nframes));
371 /* Print name of first atom in all groups in index file */
372 static void print_types(const int index[], int a[], int ngrps, char* groups[], const t_topology* top)
376 fprintf(stderr, "Using following groups: \n");
377 for (i = 0; i < ngrps; i++)
380 "Groupname: %s First atomname: %s First atomnr %d\n",
382 *(top->atoms.atomname[a[index[i]]]),
385 fprintf(stderr, "\n");
388 static void check_length(real length, int a, int b)
393 "WARNING: distance between atoms %d and "
394 "%d > 0.3 nm (%f). Index file might be corrupt.\n",
401 static void calc_order(const char* fn,
410 const t_topology* top,
414 gmx_bool permolecule,
419 const gmx_output_env_t* oenv)
421 /* if permolecule = TRUE, order parameters will be calculed per molecule
422 * and stored in slOrder with #slices = # molecules */
423 rvec *x0, /* coordinates with pbc */
424 *x1; /* coordinates without pbc */
425 matrix box; /* box (3x3) */
427 rvec cossum, /* sum of vector angles for three axes */
428 Sx, Sy, Sz, /* the three molecular axes */
429 tmp1, tmp2, /* temp. rvecs for calculating dot products */
430 frameorder; /* order parameters for one frame */
431 real* slFrameorder; /* order parameter for one frame, per slice */
432 real length, /* total distance between two atoms */
433 t, /* time from trajectory */
434 z_ave, z1, z2; /* average z, used to det. which slice atom is in */
435 int natoms, /* nr. atoms in trj */
436 nr_tails, /* nr tails, to check if index file is correct */
437 size = 0, /* nr. of atoms in group. same as nr_tails */
438 i, j, m, k, teller = 0, slice; /* current slice number */
440 int* slCount; /* nr. of atoms in one slice */
441 real sdbangle = 0; /* sum of these angles */
442 gmx_bool use_unitvector = FALSE; /* use a specified unit vector instead of axis to specify unit normal*/
444 int comsize, distsize;
445 int * comidx = nullptr, *distidx = nullptr;
446 char* grpname = nullptr;
448 real arcdist, tmpdist;
449 gmx_rmpbc_t gpbc = nullptr;
451 /* PBC added for center-of-mass vector*/
452 /* Initiate the pbc structure */
453 std::memset(&pbc, 0, sizeof(pbc));
455 if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
457 gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
460 nr_tails = index[1] - index[0];
461 fprintf(stderr, "Number of elements in first group: %d\n", nr_tails);
462 /* take first group as standard. Not rocksolid, but might catch error in index*/
467 bSliced = FALSE; /*force slices off */
468 fprintf(stderr, "Calculating order parameters for each of %d molecules\n", nslices);
473 use_unitvector = TRUE;
474 fprintf(stderr, "Select an index group to calculate the radial membrane normal\n");
475 get_index(&top->atoms, radfn, 1, &comsize, &comidx, &grpname);
479 if (grpname != nullptr)
483 fprintf(stderr, "Select an index group to use as distance reference\n");
484 get_index(&top->atoms, radfn, 1, &distsize, &distidx, &grpname);
485 bSliced = FALSE; /*force slices off*/
488 if (use_unitvector && bSliced)
491 "Warning: slicing and specified unit vectors are not currently compatible\n");
494 snew(slCount, nslices);
495 snew(*slOrder, nslices);
496 for (i = 0; i < nslices; i++)
498 snew((*slOrder)[i], ngrps);
502 snew(*distvals, nslices);
503 for (i = 0; i < nslices; i++)
505 snew((*distvals)[i], ngrps);
509 snew(slFrameorder, nslices);
514 *slWidth = box[axis][axis] / static_cast<real>(nslices);
515 fprintf(stderr, "Box divided in %d slices. Initial width of slice: %f\n", nslices, *slWidth);
520 nr_tails = index[1] - index[0];
521 fprintf(stderr, "Number of elements in first group: %d\n", nr_tails);
522 /* take first group as standard. Not rocksolid, but might catch error
528 gpbc = gmx_rmpbc_init(&top->idef, pbcType, natoms);
529 /*********** Start processing trajectory ***********/
534 *slWidth = box[axis][axis] / static_cast<real>(nslices);
538 set_pbc(&pbc, pbcType, box);
539 gmx_rmpbc_copy(gpbc, natoms, box, x0, x1);
541 /* Now loop over all groups. There are ngrps groups, the order parameter can
542 be calculated for grp 1 to grp ngrps - 1. For each group, loop over all
543 atoms in group, which is index[i] to (index[i+1] - 1) See block.h. Of
544 course, in this case index[i+1] -index[i] has to be the same for all
545 groups, namely the number of tails. i just runs over all atoms in a tail,
546 so for DPPC ngrps = 16 and i runs from 1 to 14, including 14
552 /*center-of-mass determination*/
556 for (j = 0; j < comsize; j++)
558 rvec_inc(com, x1[comidx[j]]);
560 svmul(1.0 / comsize, com, com);
562 rvec displacementFromReference;
565 rvec dref = { 0.0, 0.0, 0.0 };
566 for (j = 0; j < distsize; j++)
568 rvec_inc(dref, x1[distidx[j]]);
570 svmul(1.0 / distsize, dref, dref);
573 pbc_dx(&pbc, dref, com, displacementFromReference);
574 unitv(displacementFromReference, displacementFromReference);
578 for (i = 1; i < ngrps - 1; i++)
580 clear_rvec(frameorder);
582 size = index[i + 1] - index[i];
583 if (size != nr_tails)
586 "grp %d does not have same number of"
587 " elements as grp 1\n",
591 for (j = 0; j < size; j++)
594 /*create unit vector*/
596 pbc_dx(&pbc, x1[a[index[i] + j]], com, direction);
597 unitv(direction, direction);
600 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],
601 direction[0],direction[1],direction[2]);*/
607 /* Using convention for unsaturated carbons */
608 /* first get Sz, the vector from Cn to Cn+1 */
609 rvec_sub(x1[a[index[i + 1] + j]], x1[a[index[i] + j]], dist);
611 check_length(length, a[index[i] + j], a[index[i + 1] + j]);
612 svmul(1.0 / length, dist, Sz);
614 /* this is actually the cosine of the angle between the double bond
615 and axis, because Sz is normalized and the two other components of
616 the axis on the bilayer are zero */
619 sdbangle += gmx_angle(direction, Sz); /*this can probably be optimized*/
623 sdbangle += std::acos(Sz[axis]);
629 /* get vector dist(Cn-1,Cn+1) for tail atoms */
630 rvec_sub(x1[a[index[i + 1] + j]], x1[a[index[i - 1] + j]], dist);
631 length = norm(dist); /* determine distance between two atoms */
632 check_length(length, a[index[i - 1] + j], a[index[i + 1] + j]);
634 svmul(1.0 / length, dist, Sz);
635 /* Sz is now the molecular axis Sz, normalized and all that */
638 /* now get Sx. Sx is normal to the plane of Cn-1, Cn and Cn+1 so
639 we can use the outer product of Cn-1->Cn and Cn+1->Cn, I hope */
640 rvec_sub(x1[a[index[i + 1] + j]], x1[a[index[i] + j]], tmp1);
641 rvec_sub(x1[a[index[i - 1] + j]], x1[a[index[i] + j]], tmp2);
642 cprod(tmp1, tmp2, Sx);
643 svmul(1.0 / norm(Sx), Sx, Sx);
645 /* now we can get Sy from the outer product of Sx and Sz */
647 svmul(1.0 / norm(Sy), Sy, Sy);
649 /* the square of cosine of the angle between dist and the axis.
650 Using the innerproduct, but two of the three elements are zero
651 Determine the sum of the orderparameter of all atoms in group
655 cossum[XX] = gmx::square(iprod(Sx, direction)); /* this is allowed, since Sa is normalized */
656 cossum[YY] = gmx::square(iprod(Sy, direction));
657 cossum[ZZ] = gmx::square(iprod(Sz, direction));
661 cossum[XX] = gmx::square(Sx[axis]); /* this is allowed, since Sa is normalized */
662 cossum[YY] = gmx::square(Sy[axis]);
663 cossum[ZZ] = gmx::square(Sz[axis]);
666 for (m = 0; m < DIM; m++)
668 frameorder[m] += 0.5 * (3.0 * cossum[m] - 1.0);
673 /* get average coordinate in box length for slicing,
674 determine which slice atom is in, increase count for that
675 slice. slFrameorder and slOrder are reals, not
676 rvecs. Only the component [axis] of the order tensor is
677 kept, until I find it necessary to know the others too
680 z1 = x1[a[index[i - 1] + j]][axis];
681 z2 = x1[a[index[i + 1] + j]][axis];
682 z_ave = 0.5 * (z1 + z2);
683 slice = static_cast<int>((static_cast<real>(nslices) * z_ave) / box[axis][axis]);
686 slice += static_cast<real>(nslices);
688 slice = slice % nslices;
690 slCount[slice]++; /* determine slice, increase count */
692 slFrameorder[slice] += 0.5 * (3 * cossum[axis] - 1);
694 else if (permolecule)
696 /* store per-molecule order parameter
697 * To just track single-axis order: (*slOrder)[j][i] += 0.5 * (3 *
698 * iprod(cossum,direction) - 1); following is for Scd order: */
699 (*slOrder)[j][i] += -1
700 * (1.0 / 3.0 * (3 * cossum[XX] - 1)
701 + 1.0 / 3.0 * 0.5 * (3.0 * cossum[YY] - 1));
707 /* bin order parameter by arc distance from reference group*/
708 arcdist = gmx_angle(displacementFromReference, direction);
709 (*distvals)[j][i] += arcdist;
713 /* Want minimum lateral distance to first group calculated */
714 tmpdist = trace(box); /* should be max value */
715 for (k = 0; k < distsize; k++)
718 pbc_dx(&pbc, x1[distidx[k]], x1[a[index[i] + j]], displacement);
719 /* at the moment, just remove displacement[axis] */
720 displacement[axis] = 0;
721 tmpdist = std::min(tmpdist, norm2(displacement));
723 // fprintf(stderr, "Min dist %f; trace %f\n", tmpdist, trace(box));
724 (*distvals)[j][i] += std::sqrt(tmpdist);
727 } /* end loop j, over all atoms in group */
729 for (m = 0; m < DIM; m++)
731 (*order)[i][m] += (frameorder[m] / static_cast<real>(size));
735 { /*Skip following if doing per-molecule*/
736 for (k = 0; k < nslices; k++)
738 if (slCount[k]) /* if no elements, nothing has to be added */
740 (*slOrder)[k][i] += slFrameorder[k] / static_cast<real>(slCount[k]);
745 } /* end loop i, over all groups in indexfile */
749 } while (read_next_x(oenv, status, &t, x0, box));
750 /*********** done with status file **********/
752 fprintf(stderr, "\nRead trajectory. Printing parameters to file\n");
753 gmx_rmpbc_done(gpbc);
755 /* average over frames */
756 for (i = 1; i < ngrps - 1; i++)
758 svmul(1.0 / nr_frames, (*order)[i], (*order)[i]);
760 "Atom %d Tensor: x=%g , y=%g, z=%g\n",
765 if (bSliced || permolecule)
767 for (k = 0; k < nslices; k++)
769 (*slOrder)[k][i] /= nr_frames;
774 for (k = 0; k < nslices; k++)
776 (*distvals)[k][i] /= nr_frames;
784 "Average angle between double bond and normal: %f\n",
785 180 * sdbangle / (nr_frames * static_cast<real>(size) * M_PI));
788 sfree(x0); /* free memory used by coordinate arrays */
790 if (comidx != nullptr)
794 if (distidx != nullptr)
798 if (grpname != nullptr)
805 static void order_plot(rvec order[],
814 gmx_bool permolecule,
816 const gmx_output_env_t* oenv)
818 FILE *ord, *slOrd; /* xvgr files with order parameters */
819 int atom, slice; /* atom corresponding to order para.*/
820 char buf[256]; /* for xvgr title */
821 real S; /* order parameter averaged over all atoms */
825 sprintf(buf, "Scd order parameters");
826 ord = xvgropen(afile, buf, "Atom", "S", oenv);
827 sprintf(buf, "Orderparameters per atom per slice");
828 slOrd = xvgropen(bfile, buf, "Molecule", "S", oenv);
829 for (atom = 1; atom < ngrps - 1; atom++)
834 -1.0 * (2.0 / 3.0 * order[atom][XX] + 1.0 / 3.0 * order[atom][YY]));
837 for (slice = 0; slice < nslices; slice++)
839 fprintf(slOrd, "%12d\t", slice);
842 fprintf(slOrd, "%12g\t", distvals[slice][1]); /*use distance value at second carbon*/
844 for (atom = 1; atom < ngrps - 1; atom++)
846 fprintf(slOrd, "%12g\t", slOrder[slice][atom]);
848 fprintf(slOrd, "\n");
853 sprintf(buf, "Orderparameters Sz per atom");
854 ord = xvgropen(afile, buf, "Atom", "S", oenv);
855 fprintf(stderr, "ngrps = %d, nslices = %d", ngrps, nslices);
857 sprintf(buf, "Orderparameters per atom per slice");
858 slOrd = xvgropen(bfile, buf, "Slice", "S", oenv);
860 for (atom = 1; atom < ngrps - 1; atom++)
862 fprintf(ord, "%12d %12g\n", atom, order[atom][ZZ]);
865 for (slice = 0; slice < nslices; slice++)
868 for (atom = 1; atom < ngrps - 1; atom++)
870 S += slOrder[slice][atom];
872 fprintf(slOrd, "%12g %12g\n", static_cast<real>(slice) * slWidth, S / static_cast<real>(atom));
877 sprintf(buf, "Order tensor diagonal elements");
878 ord = xvgropen(afile, buf, "Atom", "S", oenv);
879 sprintf(buf, "Deuterium order parameters");
880 slOrd = xvgropen(cfile, buf, "Atom", "Scd", oenv);
882 for (atom = 1; atom < ngrps - 1; atom++)
884 fprintf(ord, "%12d %12g %12g %12g\n", atom, order[atom][XX], order[atom][YY], order[atom][ZZ]);
888 -1.0 * (2.0 / 3.0 * order[atom][XX] + 1.0 / 3.0 * order[atom][YY]));
895 static void write_bfactors(t_filenm* fnm,
902 const t_topology* top,
904 gmx_output_env_t* oenv)
906 /*function to write order parameters as B factors in PDB file using
907 first frame of trajectory*/
909 t_trxframe fr, frout;
913 ngrps -= 2; /*we don't have an order parameter for the first or
914 last atom in each chain*/
915 nout = nslices * ngrps;
916 read_first_frame(oenv, &status, ftp2fn(efTRX, nfile, fnm), &fr, TRX_NEED_X);
926 init_t_atoms(&useatoms, nout, TRUE);
929 /*initialize PDBinfo*/
930 for (i = 0; i < useatoms.nr; ++i)
932 useatoms.pdbinfo[i].type = PdbRecordType::Atom;
933 useatoms.pdbinfo[i].occup = 0.0;
934 useatoms.pdbinfo[i].bfac = 0.0;
935 useatoms.pdbinfo[i].bAnisotropic = FALSE;
938 for (j = 0, ctr = 0; j < nslices; j++)
940 for (i = 0; i < ngrps; i++, ctr++)
942 /*iterate along each chain*/
943 useatoms.pdbinfo[ctr].bfac = order[j][i + 1];
946 useatoms.pdbinfo[ctr].occup = distvals[j][i + 1];
948 copy_rvec(fr.x[a[index[i + 1] + j]], frout.x[ctr]);
949 useatoms.atomname[ctr] = top->atoms.atomname[a[index[i + 1] + j]];
950 useatoms.atom[ctr] = top->atoms.atom[a[index[i + 1] + j]];
951 useatoms.nres = std::max(useatoms.nres, useatoms.atom[ctr].resind + 1);
952 useatoms.resinfo[useatoms.atom[ctr].resind] =
953 top->atoms.resinfo[useatoms.atom[ctr].resind]; /*copy resinfo*/
958 opt2fn("-ob", nfile, fnm), "Order parameters", &useatoms, frout.x, nullptr, frout.pbcType, frout.box);
961 done_atom(&useatoms);
964 int gmx_order(int argc, char* argv[])
966 const char* desc[] = {
967 "[THISMODULE] computes the order parameter per atom for carbon tails. For atom i the",
968 "vector i-1, i+1 is used together with an axis. ",
969 "The index file should contain only the groups to be used for calculations,",
970 "with each group of equivalent carbons along the relevant acyl chain in its own",
971 "group. There should not be any generic groups (like System, Protein) in the index",
972 "file to avoid confusing the program (this is not relevant to tetrahedral order",
973 "parameters however, which only work for water anyway).[PAR]",
974 "[THISMODULE] can also give all",
975 "diagonal elements of the order tensor and even calculate the deuterium",
976 "order parameter Scd (default). If the option [TT]-szonly[tt] is given, only one",
977 "order tensor component (specified by the [TT]-d[tt] option) is given and the",
978 "order parameter per slice is calculated as well. If [TT]-szonly[tt] is not",
979 "selected, all diagonal elements and the deuterium order parameter is",
981 "The tetrahedrality order parameters can be determined",
982 "around an atom. Both angle an distance order parameters are calculated. See",
983 "P.-L. Chau and A.J. Hardwick, Mol. Phys., 93, (1998), 511-518.",
987 static int nslices = 1; /* nr of slices defined */
988 static gmx_bool bSzonly = FALSE; /* True if only Sz is wanted */
989 static gmx_bool bUnsat = FALSE; /* True if carbons are unsat. */
990 static const char* normal_axis[] = { nullptr, "z", "x", "y", nullptr };
991 static gmx_bool permolecule = FALSE; /*compute on a per-molecule basis */
992 static gmx_bool radial = FALSE; /*compute a radial membrane normal */
993 static gmx_bool distcalc = FALSE; /*calculate distance from a reference group */
995 { "-d", FALSE, etENUM, { normal_axis }, "Direction of the normal on the membrane" },
1000 "Calculate order parameter as function of box length, dividing the box"
1001 " into this number of slices." },
1006 "Only give Sz element of order tensor. (axis can be specified with [TT]-d[tt])" },
1011 "Calculate order parameters for unsaturated carbons. Note that this can"
1012 "not be mixed with normal order parameters." },
1017 "Compute per-molecule Scd order parameters" },
1018 { "-radial", FALSE, etBOOL, { &radial }, "Compute a radial membrane normal" },
1019 { "-calcdist", FALSE, etBOOL, { &distcalc }, "Compute distance from a reference" },
1022 rvec* order; /* order par. for each atom */
1023 real** slOrder; /* same, per slice */
1024 real slWidth = 0.0; /* width of a slice */
1025 char** grpname; /* groupnames */
1026 int ngrps, /* nr. of groups */
1027 i, axis = 0; /* normal axis */
1028 t_topology* top; /* topology */
1029 PbcType pbcType; /* type of periodic boundary conditions */
1030 int * index, /* indices for a */
1031 *a; /* atom numbers in each group */
1032 t_blocka* block; /* data from index file */
1034 /* files for g_order */
1035 { efTRX, "-f", nullptr, ffREAD }, /* trajectory file */
1036 { efNDX, "-n", nullptr, ffREAD }, /* index file */
1037 { efNDX, "-nr", nullptr, ffOPTRD }, /* index for radial axis calculation */
1038 { efTPR, nullptr, nullptr, ffREAD }, /* topology file */
1039 { efXVG, "-o", "order", ffWRITE }, /* xvgr output file */
1040 { efXVG, "-od", "deuter", ffWRITE }, /* xvgr output file */
1041 { efPDB, "-ob", nullptr, ffOPTWR }, /* write Scd as B factors to PDB if permolecule */
1042 { efXVG, "-os", "sliced", ffWRITE }, /* xvgr output file */
1043 { efXVG, "-Sg", "sg-ang", ffOPTWR }, /* xvgr output file */
1044 { efXVG, "-Sk", "sk-dist", ffOPTWR }, /* xvgr output file */
1045 { efXVG, "-Sgsl", "sg-ang-slice", ffOPTWR }, /* xvgr output file */
1046 { efXVG, "-Sksl", "sk-dist-slice", ffOPTWR }, /* xvgr output file */
1048 gmx_bool bSliced = FALSE; /* True if box is sliced */
1049 #define NFILE asize(fnm)
1050 real** distvals = nullptr;
1051 const char * sgfnm, *skfnm, *ndxfnm, *tpsfnm, *trxfnm;
1052 gmx_output_env_t* oenv;
1054 if (!parse_common_args(
1055 &argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME, NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
1061 gmx_fatal(FARGS, "Can not have nslices < 1");
1063 sgfnm = opt2fn_null("-Sg", NFILE, fnm);
1064 skfnm = opt2fn_null("-Sk", NFILE, fnm);
1065 ndxfnm = opt2fn_null("-n", NFILE, fnm);
1066 tpsfnm = ftp2fn(efTPR, NFILE, fnm);
1067 trxfnm = ftp2fn(efTRX, NFILE, fnm);
1069 /* Calculate axis */
1070 GMX_RELEASE_ASSERT(normal_axis[0] != nullptr, "Options inconsistency; normal_axis[0] is NULL");
1071 if (std::strcmp(normal_axis[0], "x") == 0)
1075 else if (std::strcmp(normal_axis[0], "y") == 0)
1079 else if (std::strcmp(normal_axis[0], "z") == 0)
1085 gmx_fatal(FARGS, "Invalid axis, use x, y or z");
1090 case 0: fprintf(stderr, "Taking x axis as normal to the membrane\n"); break;
1091 case 1: fprintf(stderr, "Taking y axis as normal to the membrane\n"); break;
1092 case 2: fprintf(stderr, "Taking z axis as normal to the membrane\n"); break;
1095 /* tetraheder order parameter */
1098 /* If either of theoptions is set we compute both */
1099 sgfnm = opt2fn("-Sg", NFILE, fnm);
1100 skfnm = opt2fn("-Sk", NFILE, fnm);
1101 calc_tetra_order_parm(ndxfnm,
1108 opt2fn("-Sgsl", NFILE, fnm),
1109 opt2fn("-Sksl", NFILE, fnm),
1111 /* view xvgr files */
1112 do_view(oenv, opt2fn("-Sg", NFILE, fnm), nullptr);
1113 do_view(oenv, opt2fn("-Sk", NFILE, fnm), nullptr);
1116 do_view(oenv, opt2fn("-Sgsl", NFILE, fnm), nullptr);
1117 do_view(oenv, opt2fn("-Sksl", NFILE, fnm), nullptr);
1122 /* tail order parameter */
1127 fprintf(stderr, "Dividing box in %d slices.\n\n", nslices);
1132 fprintf(stderr, "Only calculating Sz\n");
1136 fprintf(stderr, "Taking carbons as unsaturated!\n");
1139 top = read_top(ftp2fn(efTPR, NFILE, fnm), &pbcType); /* read topology file */
1141 block = init_index(ftp2fn(efNDX, NFILE, fnm), &grpname);
1142 index = block->index; /* get indices from t_block block */
1143 a = block->a; /* see block.h */
1148 nslices = index[1] - index[0]; /*I think this assumes contiguous lipids in topology*/
1149 fprintf(stderr, "Calculating Scd order parameters for each of %d molecules\n", nslices);
1154 fprintf(stderr, "Calculating radial distances\n");
1157 gmx_fatal(FARGS, "Cannot yet output radial distances without permolecule\n");
1161 /* show atomtypes, to check if index file is correct */
1162 print_types(index, a, ngrps, grpname, top);
1164 calc_order(ftp2fn(efTRX, NFILE, fnm),
1180 opt2fn_null("-nr", NFILE, fnm),
1186 ngrps--; /*don't print the last group--was used for
1187 center-of-mass determination*/
1191 opt2fn("-o", NFILE, fnm),
1192 opt2fn("-os", NFILE, fnm),
1193 opt2fn("-od", NFILE, fnm),
1202 if (opt2bSet("-ob", NFILE, fnm))
1207 "Won't write B-factors with averaged order parameters; use -permolecule\n");
1211 write_bfactors(fnm, NFILE, index, a, nslices, ngrps, slOrder, top, distvals, oenv);
1216 do_view(oenv, opt2fn("-o", NFILE, fnm), nullptr); /* view xvgr file */
1217 do_view(oenv, opt2fn("-os", NFILE, fnm), nullptr); /* view xvgr file */
1218 do_view(oenv, opt2fn("-od", NFILE, fnm), nullptr); /* view xvgr file */
1221 if (distvals != nullptr)
1223 for (i = 0; i < nslices; ++i)