3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
60 /****************************************************************************/
61 /* This program calculates the order parameter per atom for an interface or */
62 /* bilayer, averaged over time. */
63 /* S = 1/2 * (3 * cos(i)cos(j) - delta(ij)) */
64 /* It is assumed that the order parameter with respect to a box-axis */
65 /* is calculated. In that case i = j = axis, and delta(ij) = 1. */
67 /* Peter Tieleman, April 1995 */
68 /* P.J. van Maaren, November 2005 Added tetrahedral stuff */
69 /****************************************************************************/
71 static void find_nearest_neighbours(int ePBC,
72 int natoms, matrix box,
73 rvec x[], int maxidx, atom_id index[],
74 real *sgmean, real *skmean,
75 int nslice, int slice_dim,
76 real sgslice[], real skslice[],
81 int ix, jx, nsgbin, *sgbin;
82 int i1, i2, i, ibin, j, k, l, n, *nn[4];
83 rvec dx, dx1, dx2, rj, rk, urk, urj;
84 real cost, cost2, *sgmol, *skmol, rmean, rmean2, r2, box2, *r_nn[4];
90 real onethird = 1.0/3.0;
91 /* dmat = init_mat(maxidx, FALSE); */
92 box2 = box[XX][XX] * box[XX][XX];
93 snew(sl_count, nslice);
94 for (i = 0; (i < 4); i++)
96 snew(r_nn[i], natoms);
99 for (j = 0; (j < natoms); j++)
108 /* Must init pbc every step because of pressure coupling */
109 set_pbc(&pbc, ePBC, box);
111 gmx_rmpbc(gpbc, natoms, box, x);
113 nsgbin = 1 + 1/0.0005;
119 for (i = 0; (i < maxidx); i++) /* loop over index file */
122 for (j = 0; (j < maxidx); j++)
131 pbc_dx(&pbc, x[ix], x[jx], dx);
134 /* set_mat_entry(dmat,i,j,r2); */
136 /* determine the nearest neighbours */
139 r_nn[3][i] = r_nn[2][i]; nn[3][i] = nn[2][i];
140 r_nn[2][i] = r_nn[1][i]; nn[2][i] = nn[1][i];
141 r_nn[1][i] = r_nn[0][i]; nn[1][i] = nn[0][i];
142 r_nn[0][i] = r2; nn[0][i] = j;
144 else if (r2 < r_nn[1][i])
146 r_nn[3][i] = r_nn[2][i]; nn[3][i] = nn[2][i];
147 r_nn[2][i] = r_nn[1][i]; nn[2][i] = nn[1][i];
148 r_nn[1][i] = r2; nn[1][i] = j;
150 else if (r2 < r_nn[2][i])
152 r_nn[3][i] = r_nn[2][i]; nn[3][i] = nn[2][i];
153 r_nn[2][i] = r2; nn[2][i] = j;
155 else if (r2 < r_nn[3][i])
157 r_nn[3][i] = r2; nn[3][i] = j;
162 /* calculate mean distance between nearest neighbours */
164 for (j = 0; (j < 4); j++)
166 r_nn[j][i] = sqrt(r_nn[j][i]);
175 /* Chau1998a eqn 3 */
176 /* angular part tetrahedrality order parameter per atom */
177 for (j = 0; (j < 3); j++)
179 for (k = j+1; (k < 4); k++)
181 pbc_dx(&pbc, x[ix], x[index[nn[k][i]]], rk);
182 pbc_dx(&pbc, x[ix], x[index[nn[j][i]]], rj);
187 cost = iprod(urk, urj) + onethird;
190 /* sgmol[i] += 3*cost2/32; */
193 /* determine distribution */
194 ibin = nsgbin * cost2;
199 /* printf("%d %d %f %d %d\n", j, k, cost * cost, ibin, sgbin[ibin]);*/
205 /* normalize sgmol between 0.0 and 1.0 */
206 sgmol[i] = 3*sgmol[i]/32;
209 /* distance part tetrahedrality order parameter per atom */
210 rmean2 = 4 * 3 * rmean * rmean;
211 for (j = 0; (j < 4); j++)
213 skmol[i] += (rmean - r_nn[j][i]) * (rmean - r_nn[j][i]) / rmean2;
214 /* printf("%d %f (%f %f %f %f) \n",
215 i, skmol[i], rmean, rmean2, r_nn[j][i], (rmean - r_nn[j][i]) );
221 /* Compute sliced stuff */
222 sl_index = gmx_nint((1+x[i][slice_dim]/box[slice_dim][slice_dim])*nslice) % nslice;
223 sgslice[sl_index] += sgmol[i];
224 skslice[sl_index] += skmol[i];
225 sl_count[sl_index]++;
226 } /* loop over entries in index file */
231 for (i = 0; (i < nslice); i++)
235 sgslice[i] /= sl_count[i];
236 skslice[i] /= sl_count[i];
243 for (i = 0; (i < 4); i++)
251 static void calc_tetra_order_parm(const char *fnNDX, const char *fnTPS,
252 const char *fnTRX, const char *sgfn,
254 int nslice, int slice_dim,
255 const char *sgslfn, const char *skslfn,
256 const output_env_t oenv)
258 FILE *fpsg = NULL, *fpsk = NULL;
261 char title[STRLEN], fn[STRLEN], subtitle[STRLEN];
270 int i, *isize, ng, nframes;
271 real *sg_slice, *sg_slice_tot, *sk_slice, *sk_slice_tot;
272 gmx_rmpbc_t gpbc = NULL;
275 read_tps_conf(fnTPS, title, &top, &ePBC, &xtop, NULL, box, FALSE);
277 snew(sg_slice, nslice);
278 snew(sk_slice, nslice);
279 snew(sg_slice_tot, nslice);
280 snew(sk_slice_tot, nslice);
282 /* get index groups */
283 printf("Select the group that contains the atoms you want to use for the tetrahedrality order parameter calculation:\n");
287 get_index(&top.atoms, fnNDX, ng, isize, index, grpname);
289 /* Analyze trajectory */
290 natoms = read_first_x(oenv, &status, fnTRX, &t, &x, box);
291 if (natoms > top.atoms.nr)
293 gmx_fatal(FARGS, "Topology (%d atoms) does not match trajectory (%d atoms)",
294 top.atoms.nr, natoms);
296 check_index(NULL, ng, index[0], NULL, natoms);
298 fpsg = xvgropen(sgfn, "S\\sg\\N Angle Order Parameter", "Time (ps)", "S\\sg\\N",
300 fpsk = xvgropen(skfn, "S\\sk\\N Distance Order Parameter", "Time (ps)", "S\\sk\\N",
303 /* loop over frames */
304 gpbc = gmx_rmpbc_init(&top.idef, ePBC, natoms, box);
308 find_nearest_neighbours(ePBC, natoms, box, x, isize[0], index[0],
309 &sg, &sk, nslice, slice_dim, sg_slice, sk_slice, gpbc);
310 for (i = 0; (i < nslice); i++)
312 sg_slice_tot[i] += sg_slice[i];
313 sk_slice_tot[i] += sk_slice[i];
315 fprintf(fpsg, "%f %f\n", t, sg);
316 fprintf(fpsk, "%f %f\n", t, sk);
319 while (read_next_x(oenv, status, &t, natoms, x, box));
321 gmx_rmpbc_done(gpbc);
330 fpsg = xvgropen(sgslfn,
331 "S\\sg\\N Angle Order Parameter / Slab", "(nm)", "S\\sg\\N",
333 fpsk = xvgropen(skslfn,
334 "S\\sk\\N Distance Order Parameter / Slab", "(nm)", "S\\sk\\N",
336 for (i = 0; (i < nslice); i++)
338 fprintf(fpsg, "%10g %10g\n", (i+0.5)*box[slice_dim][slice_dim]/nslice,
339 sg_slice_tot[i]/nframes);
340 fprintf(fpsk, "%10g %10g\n", (i+0.5)*box[slice_dim][slice_dim]/nslice,
341 sk_slice_tot[i]/nframes);
348 /* Print name of first atom in all groups in index file */
349 static void print_types(atom_id index[], atom_id a[], int ngrps,
350 char *groups[], t_topology *top)
354 fprintf(stderr, "Using following groups: \n");
355 for (i = 0; i < ngrps; i++)
357 fprintf(stderr, "Groupname: %s First atomname: %s First atomnr %u\n",
358 groups[i], *(top->atoms.atomname[a[index[i]]]), a[index[i]]);
360 fprintf(stderr, "\n");
363 static void check_length(real length, int a, int b)
367 fprintf(stderr, "WARNING: distance between atoms %d and "
368 "%d > 0.3 nm (%f). Index file might be corrupt.\n",
373 void calc_order(const char *fn, atom_id *index, atom_id *a, rvec **order,
374 real ***slOrder, real *slWidth, int nslices, gmx_bool bSliced,
375 gmx_bool bUnsat, t_topology *top, int ePBC, int ngrps, int axis,
376 gmx_bool permolecule, gmx_bool radial, gmx_bool distcalc, const char *radfn,
378 const output_env_t oenv)
380 /* if permolecule = TRUE, order parameters will be calculed per molecule
381 * and stored in slOrder with #slices = # molecules */
382 rvec *x0, /* coordinates with pbc */
383 *x1, /* coordinates without pbc */
384 dist; /* vector between two atoms */
385 matrix box; /* box (3x3) */
387 rvec cossum, /* sum of vector angles for three axes */
388 Sx, Sy, Sz, /* the three molecular axes */
389 tmp1, tmp2, /* temp. rvecs for calculating dot products */
390 frameorder; /* order parameters for one frame */
391 real *slFrameorder; /* order parameter for one frame, per slice */
392 real length, /* total distance between two atoms */
393 t, /* time from trajectory */
394 z_ave, z1, z2; /* average z, used to det. which slice atom is in */
395 int natoms, /* nr. atoms in trj */
396 nr_tails, /* nr tails, to check if index file is correct */
397 size = 0, /* nr. of atoms in group. same as nr_tails */
398 i, j, m, k, l, teller = 0,
399 slice, /* current slice number */
401 int *slCount; /* nr. of atoms in one slice */
402 real dbangle = 0, /* angle between double bond and axis */
403 sdbangle = 0; /* sum of these angles */
404 gmx_bool use_unitvector = FALSE; /* use a specified unit vector instead of axis to specify unit normal*/
405 rvec direction, com, dref, dvec;
406 int comsize, distsize;
407 atom_id *comidx = NULL, *distidx = NULL;
408 char *grpname = NULL;
410 real arcdist, tmpdist;
411 gmx_rmpbc_t gpbc = NULL;
413 /* PBC added for center-of-mass vector*/
414 /* Initiate the pbc structure */
415 memset(&pbc, 0, sizeof(pbc));
417 if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
419 gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
422 nr_tails = index[1] - index[0];
423 fprintf(stderr, "Number of elements in first group: %d\n", nr_tails);
424 /* take first group as standard. Not rocksolid, but might catch error in index*/
429 bSliced = FALSE; /*force slices off */
430 fprintf(stderr, "Calculating order parameters for each of %d molecules\n",
436 use_unitvector = TRUE;
437 fprintf(stderr, "Select an index group to calculate the radial membrane normal\n");
438 get_index(&top->atoms, radfn, 1, &comsize, &comidx, &grpname);
446 fprintf(stderr, "Select an index group to use as distance reference\n");
447 get_index(&top->atoms, radfn, 1, &distsize, &distidx, &grpname);
448 bSliced = FALSE; /*force slices off*/
451 if (use_unitvector && bSliced)
453 fprintf(stderr, "Warning: slicing and specified unit vectors are not currently compatible\n");
456 snew(slCount, nslices);
457 snew(*slOrder, nslices);
458 for (i = 0; i < nslices; i++)
460 snew((*slOrder)[i], ngrps);
464 snew(*distvals, nslices);
465 for (i = 0; i < nslices; i++)
467 snew((*distvals)[i], ngrps);
471 snew(slFrameorder, nslices);
476 *slWidth = box[axis][axis]/nslices;
477 fprintf(stderr, "Box divided in %d slices. Initial width of slice: %f\n",
483 nr_tails = index[1] - index[0];
484 fprintf(stderr, "Number of elements in first group: %d\n", nr_tails);
485 /* take first group as standard. Not rocksolid, but might catch error
491 gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms, box);
492 /*********** Start processing trajectory ***********/
497 *slWidth = box[axis][axis]/nslices;
501 set_pbc(&pbc, ePBC, box);
502 gmx_rmpbc_copy(gpbc, natoms, box, x0, x1);
504 /* Now loop over all groups. There are ngrps groups, the order parameter can
505 be calculated for grp 1 to grp ngrps - 1. For each group, loop over all
506 atoms in group, which is index[i] to (index[i+1] - 1) See block.h. Of
507 course, in this case index[i+1] -index[i] has to be the same for all
508 groups, namely the number of tails. i just runs over all atoms in a tail,
509 so for DPPC ngrps = 16 and i runs from 1 to 14, including 14
515 /*center-of-mass determination*/
516 com[XX] = 0.0; com[YY] = 0.0; com[ZZ] = 0.0;
517 for (j = 0; j < comsize; j++)
519 rvec_inc(com, x1[comidx[j]]);
521 svmul(1.0/comsize, com, com);
525 dref[XX] = 0.0; dref[YY] = 0.0; dref[ZZ] = 0.0;
526 for (j = 0; j < distsize; j++)
528 rvec_inc(dist, x1[distidx[j]]);
530 svmul(1.0/distsize, dref, dref);
533 pbc_dx(&pbc, dref, com, dvec);
538 for (i = 1; i < ngrps - 1; i++)
540 clear_rvec(frameorder);
542 size = index[i+1] - index[i];
543 if (size != nr_tails)
545 gmx_fatal(FARGS, "grp %d does not have same number of"
546 " elements as grp 1\n", i);
549 for (j = 0; j < size; j++)
552 /*create unit vector*/
554 pbc_dx(&pbc, x1[a[index[i]+j]], com, direction);
555 unitv(direction, direction);
558 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],
559 direction[0],direction[1],direction[2]);*/
564 /* Using convention for unsaturated carbons */
565 /* first get Sz, the vector from Cn to Cn+1 */
566 rvec_sub(x1[a[index[i+1]+j]], x1[a[index[i]+j]], dist);
568 check_length(length, a[index[i]+j], a[index[i+1]+j]);
569 svmul(1/length, dist, Sz);
571 /* this is actually the cosine of the angle between the double bond
572 and axis, because Sz is normalized and the two other components of
573 the axis on the bilayer are zero */
576 sdbangle += gmx_angle(direction, Sz); /*this can probably be optimized*/
580 sdbangle += acos(Sz[axis]);
585 /* get vector dist(Cn-1,Cn+1) for tail atoms */
586 rvec_sub(x1[a[index[i+1]+j]], x1[a[index[i-1]+j]], dist);
587 length = norm(dist); /* determine distance between two atoms */
588 check_length(length, a[index[i-1]+j], a[index[i+1]+j]);
590 svmul(1/length, dist, Sz);
591 /* Sz is now the molecular axis Sz, normalized and all that */
594 /* now get Sx. Sx is normal to the plane of Cn-1, Cn and Cn+1 so
595 we can use the outer product of Cn-1->Cn and Cn+1->Cn, I hope */
596 rvec_sub(x1[a[index[i+1]+j]], x1[a[index[i]+j]], tmp1);
597 rvec_sub(x1[a[index[i-1]+j]], x1[a[index[i]+j]], tmp2);
598 cprod(tmp1, tmp2, Sx);
599 svmul(1/norm(Sx), Sx, Sx);
601 /* now we can get Sy from the outer product of Sx and Sz */
603 svmul(1/norm(Sy), Sy, Sy);
605 /* the square of cosine of the angle between dist and the axis.
606 Using the innerproduct, but two of the three elements are zero
607 Determine the sum of the orderparameter of all atoms in group
611 cossum[XX] = sqr(iprod(Sx, direction)); /* this is allowed, since Sa is normalized */
612 cossum[YY] = sqr(iprod(Sy, direction));
613 cossum[ZZ] = sqr(iprod(Sz, direction));
617 cossum[XX] = sqr(Sx[axis]); /* this is allowed, since Sa is normalized */
618 cossum[YY] = sqr(Sy[axis]);
619 cossum[ZZ] = sqr(Sz[axis]);
622 for (m = 0; m < DIM; m++)
624 frameorder[m] += 0.5 * (3 * cossum[m] - 1);
629 /* get average coordinate in box length for slicing,
630 determine which slice atom is in, increase count for that
631 slice. slFrameorder and slOrder are reals, not
632 rvecs. Only the component [axis] of the order tensor is
633 kept, until I find it necessary to know the others too
636 z1 = x1[a[index[i-1]+j]][axis];
637 z2 = x1[a[index[i+1]+j]][axis];
638 z_ave = 0.5 * (z1 + z2);
641 z_ave += box[axis][axis];
643 if (z_ave > box[axis][axis])
645 z_ave -= box[axis][axis];
648 slice = (int)(0.5 + (z_ave / (*slWidth))) - 1;
649 slCount[slice]++; /* determine slice, increase count */
651 slFrameorder[slice] += 0.5 * (3 * cossum[axis] - 1);
653 else if (permolecule)
655 /* store per-molecule order parameter
656 * To just track single-axis order: (*slOrder)[j][i] += 0.5 * (3 * iprod(cossum,direction) - 1);
657 * following is for Scd order: */
658 (*slOrder)[j][i] += -1* (0.3333 * (3 * cossum[XX] - 1) + 0.3333 * 0.5 * (3 * cossum[YY] - 1));
664 /* bin order parameter by arc distance from reference group*/
665 arcdist = gmx_angle(dvec, direction);
666 (*distvals)[j][i] += arcdist;
670 /* Want minimum lateral distance to first group calculated */
671 tmpdist = trace(box); /* should be max value */
672 for (k = 0; k < distsize; k++)
674 pbc_dx(&pbc, x1[distidx[k]], x1[a[index[i]+j]], dvec);
675 /* at the moment, just remove dvec[axis] */
677 tmpdist = min(tmpdist, norm2(dvec));
679 //fprintf(stderr, "Min dist %f; trace %f\n", tmpdist, trace(box));
680 (*distvals)[j][i] += sqrt(tmpdist);
683 } /* end loop j, over all atoms in group */
685 for (m = 0; m < DIM; m++)
687 (*order)[i][m] += (frameorder[m]/size);
691 { /*Skip following if doing per-molecule*/
692 for (k = 0; k < nslices; k++)
694 if (slCount[k]) /* if no elements, nothing has to be added */
696 (*slOrder)[k][i] += slFrameorder[k]/slCount[k];
697 slFrameorder[k] = 0; slCount[k] = 0;
700 } /* end loop i, over all groups in indexfile */
705 while (read_next_x(oenv, status, &t, natoms, x0, box));
706 /*********** done with status file **********/
708 fprintf(stderr, "\nRead trajectory. Printing parameters to file\n");
709 gmx_rmpbc_done(gpbc);
711 /* average over frames */
712 for (i = 1; i < ngrps - 1; i++)
714 svmul(1.0/nr_frames, (*order)[i], (*order)[i]);
715 fprintf(stderr, "Atom %d Tensor: x=%g , y=%g, z=%g\n", i, (*order)[i][XX],
716 (*order)[i][YY], (*order)[i][ZZ]);
717 if (bSliced || permolecule)
719 for (k = 0; k < nslices; k++)
721 (*slOrder)[k][i] /= nr_frames;
726 for (k = 0; k < nslices; k++)
728 (*distvals)[k][i] /= nr_frames;
735 fprintf(stderr, "Average angle between double bond and normal: %f\n",
736 180*sdbangle/(nr_frames * size*M_PI));
739 sfree(x0); /* free memory used by coordinate arrays */
756 void order_plot(rvec order[], real *slOrder[], const char *afile, const char *bfile,
757 const char *cfile, int ngrps, int nslices, real slWidth, gmx_bool bSzonly,
758 gmx_bool permolecule, real **distvals, const output_env_t oenv)
760 FILE *ord, *slOrd; /* xvgr files with order parameters */
761 int atom, slice; /* atom corresponding to order para.*/
762 char buf[256]; /* for xvgr title */
763 real S; /* order parameter averaged over all atoms */
767 sprintf(buf, "Scd order parameters");
768 ord = xvgropen(afile, buf, "Atom", "S", oenv);
769 sprintf(buf, "Orderparameters per atom per slice");
770 slOrd = xvgropen(bfile, buf, "Molecule", "S", oenv);
771 for (atom = 1; atom < ngrps - 1; atom++)
773 fprintf(ord, "%12d %12g\n", atom, -1 * (0.6667 * order[atom][XX] +
774 0.333 * order[atom][YY]));
777 for (slice = 0; slice < nslices; slice++)
779 fprintf(slOrd, "%12d\t", slice);
782 fprintf(slOrd, "%12g\t", distvals[slice][1]); /*use distance value at second carbon*/
784 for (atom = 1; atom < ngrps - 1; atom++)
786 fprintf(slOrd, "%12g\t", slOrder[slice][atom]);
788 fprintf(slOrd, "\n");
794 sprintf(buf, "Orderparameters Sz per atom");
795 ord = xvgropen(afile, buf, "Atom", "S", oenv);
796 fprintf(stderr, "ngrps = %d, nslices = %d", ngrps, nslices);
798 sprintf(buf, "Orderparameters per atom per slice");
799 slOrd = xvgropen(bfile, buf, "Slice", "S", oenv);
801 for (atom = 1; atom < ngrps - 1; atom++)
803 fprintf(ord, "%12d %12g\n", atom, order[atom][ZZ]);
806 for (slice = 0; slice < nslices; slice++)
809 for (atom = 1; atom < ngrps - 1; atom++)
811 S += slOrder[slice][atom];
813 fprintf(slOrd, "%12g %12g\n", slice*slWidth, S/atom);
819 sprintf(buf, "Order tensor diagonal elements");
820 ord = xvgropen(afile, buf, "Atom", "S", oenv);
821 sprintf(buf, "Deuterium order parameters");
822 slOrd = xvgropen(cfile, buf, "Atom", "Scd", oenv);
824 for (atom = 1; atom < ngrps - 1; atom++)
826 fprintf(ord, "%12d %12g %12g %12g\n", atom, order[atom][XX],
827 order[atom][YY], order[atom][ZZ]);
828 fprintf(slOrd, "%12d %12g\n", atom, -1 * (0.6667 * order[atom][XX] +
829 0.333 * order[atom][YY]));
837 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)
839 /*function to write order parameters as B factors in PDB file using
840 first frame of trajectory*/
843 t_trxframe fr, frout;
847 ngrps -= 2; /*we don't have an order parameter for the first or
848 last atom in each chain*/
849 nout = nslices*ngrps;
850 natoms = read_first_frame(oenv, &status, ftp2fn(efTRX, nfile, fnm), &fr,
860 init_t_atoms(&useatoms, nout, TRUE);
863 /*initialize PDBinfo*/
864 for (i = 0; i < useatoms.nr; ++i)
866 useatoms.pdbinfo[i].type = 0;
867 useatoms.pdbinfo[i].occup = 0.0;
868 useatoms.pdbinfo[i].bfac = 0.0;
869 useatoms.pdbinfo[i].bAnisotropic = FALSE;
872 for (j = 0, ctr = 0; j < nslices; j++)
874 for (i = 0; i < ngrps; i++, ctr++)
876 /*iterate along each chain*/
877 useatoms.pdbinfo[ctr].bfac = order[j][i+1];
880 useatoms.pdbinfo[ctr].occup = distvals[j][i+1];
882 copy_rvec(fr.x[a[index[i+1]+j]], frout.x[ctr]);
883 useatoms.atomname[ctr] = top->atoms.atomname[a[index[i+1]+j]];
884 useatoms.atom[ctr] = top->atoms.atom[a[index[i+1]+j]];
885 useatoms.nres = max(useatoms.nres, useatoms.atom[ctr].resind+1);
886 useatoms.resinfo[useatoms.atom[ctr].resind] = top->atoms.resinfo[useatoms.atom[ctr].resind]; /*copy resinfo*/
890 write_sto_conf(opt2fn("-ob", nfile, fnm), "Order parameters", &useatoms, frout.x, NULL, frout.ePBC, frout.box);
893 free_t_atoms(&useatoms, FALSE);
896 int gmx_order(int argc, char *argv[])
898 const char *desc[] = {
899 "Compute the order parameter per atom for carbon tails. For atom i the",
900 "vector i-1, i+1 is used together with an axis. ",
901 "The index file should contain only the groups to be used for calculations,",
902 "with each group of equivalent carbons along the relevant acyl chain in its own",
903 "group. There should not be any generic groups (like System, Protein) in the index",
904 "file to avoid confusing the program (this is not relevant to tetrahedral order",
905 "parameters however, which only work for water anyway).[PAR]",
906 "The program can also give all",
907 "diagonal elements of the order tensor and even calculate the deuterium",
908 "order parameter Scd (default). If the option [TT]-szonly[tt] is given, only one",
909 "order tensor component (specified by the [TT]-d[tt] option) is given and the",
910 "order parameter per slice is calculated as well. If [TT]-szonly[tt] is not",
911 "selected, all diagonal elements and the deuterium order parameter is",
913 "The tetrahedrality order parameters can be determined",
914 "around an atom. Both angle an distance order parameters are calculated. See",
915 "P.-L. Chau and A.J. Hardwick, Mol. Phys., 93, (1998), 511-518.",
916 "for more details.[BR]",
920 static int nslices = 1; /* nr of slices defined */
921 static gmx_bool bSzonly = FALSE; /* True if only Sz is wanted */
922 static gmx_bool bUnsat = FALSE; /* True if carbons are unsat. */
923 static const char *normal_axis[] = { NULL, "z", "x", "y", NULL };
924 static gmx_bool permolecule = FALSE; /*compute on a per-molecule basis */
925 static gmx_bool radial = FALSE; /*compute a radial membrane normal */
926 static gmx_bool distcalc = FALSE; /*calculate distance from a reference group */
928 { "-d", FALSE, etENUM, {normal_axis},
929 "Direction of the normal on the membrane" },
930 { "-sl", FALSE, etINT, {&nslices},
931 "Calculate order parameter as function of box length, dividing the box"
932 " into this number of slices." },
933 { "-szonly", FALSE, etBOOL, {&bSzonly},
934 "Only give Sz element of order tensor. (axis can be specified with [TT]-d[tt])" },
935 { "-unsat", FALSE, etBOOL, {&bUnsat},
936 "Calculate order parameters for unsaturated carbons. Note that this can"
937 "not be mixed with normal order parameters." },
938 { "-permolecule", FALSE, etBOOL, {&permolecule},
939 "Compute per-molecule Scd order parameters" },
940 { "-radial", FALSE, etBOOL, {&radial},
941 "Compute a radial membrane normal" },
942 { "-calcdist", FALSE, etBOOL, {&distcalc},
943 "Compute distance from a reference" },
946 rvec *order; /* order par. for each atom */
947 real **slOrder; /* same, per slice */
948 real slWidth = 0.0; /* width of a slice */
949 char **grpname; /* groupnames */
950 int ngrps, /* nr. of groups */
952 axis = 0; /* normal axis */
953 t_topology *top; /* topology */
955 atom_id *index, /* indices for a */
956 *a; /* atom numbers in each group */
957 t_blocka *block; /* data from index file */
958 t_filenm fnm[] = { /* files for g_order */
959 { efTRX, "-f", NULL, ffREAD }, /* trajectory file */
960 { efNDX, "-n", NULL, ffREAD }, /* index file */
961 { efNDX, "-nr", NULL, ffREAD }, /* index for radial axis calculation */
962 { efTPX, NULL, NULL, ffREAD }, /* topology file */
963 { efXVG, "-o", "order", ffWRITE }, /* xvgr output file */
964 { efXVG, "-od", "deuter", ffWRITE }, /* xvgr output file */
965 { efPDB, "-ob", NULL, ffWRITE }, /* write Scd as B factors to PDB if permolecule */
966 { efXVG, "-os", "sliced", ffWRITE }, /* xvgr output file */
967 { efXVG, "-Sg", "sg-ang", ffOPTWR }, /* xvgr output file */
968 { efXVG, "-Sk", "sk-dist", ffOPTWR }, /* xvgr output file */
969 { efXVG, "-Sgsl", "sg-ang-slice", ffOPTWR }, /* xvgr output file */
970 { efXVG, "-Sksl", "sk-dist-slice", ffOPTWR }, /* xvgr output file */
972 gmx_bool bSliced = FALSE; /* True if box is sliced */
973 #define NFILE asize(fnm)
974 real **distvals = NULL;
975 const char *sgfnm, *skfnm, *ndxfnm, *tpsfnm, *trxfnm;
978 parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
979 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv);
982 gmx_fatal(FARGS, "Can not have nslices < 1");
984 sgfnm = opt2fn_null("-Sg", NFILE, fnm);
985 skfnm = opt2fn_null("-Sk", NFILE, fnm);
986 ndxfnm = opt2fn_null("-n", NFILE, fnm);
987 tpsfnm = ftp2fn(efTPX, NFILE, fnm);
988 trxfnm = ftp2fn(efTRX, NFILE, fnm);
991 if (strcmp(normal_axis[0], "x") == 0)
995 else if (strcmp(normal_axis[0], "y") == 0)
999 else if (strcmp(normal_axis[0], "z") == 0)
1005 gmx_fatal(FARGS, "Invalid axis, use x, y or z");
1011 fprintf(stderr, "Taking x axis as normal to the membrane\n");
1014 fprintf(stderr, "Taking y axis as normal to the membrane\n");
1017 fprintf(stderr, "Taking z axis as normal to the membrane\n");
1021 /* tetraheder order parameter */
1024 /* If either of theoptions is set we compute both */
1025 sgfnm = opt2fn("-Sg", NFILE, fnm);
1026 skfnm = opt2fn("-Sk", NFILE, fnm);
1027 calc_tetra_order_parm(ndxfnm, tpsfnm, trxfnm, sgfnm, skfnm, nslices, axis,
1028 opt2fn("-Sgsl", NFILE, fnm), opt2fn("-Sksl", NFILE, fnm),
1030 /* view xvgr files */
1031 do_view(oenv, opt2fn("-Sg", NFILE, fnm), NULL);
1032 do_view(oenv, opt2fn("-Sk", NFILE, fnm), NULL);
1035 do_view(oenv, opt2fn("-Sgsl", NFILE, fnm), NULL);
1036 do_view(oenv, opt2fn("-Sksl", NFILE, fnm), NULL);
1041 /* tail order parameter */
1046 fprintf(stderr, "Dividing box in %d slices.\n\n", nslices);
1051 fprintf(stderr, "Only calculating Sz\n");
1055 fprintf(stderr, "Taking carbons as unsaturated!\n");
1058 top = read_top(ftp2fn(efTPX, NFILE, fnm), &ePBC); /* read topology file */
1060 block = init_index(ftp2fn(efNDX, NFILE, fnm), &grpname);
1061 index = block->index; /* get indices from t_block block */
1062 a = block->a; /* see block.h */
1067 nslices = index[1] - index[0]; /*I think this assumes contiguous lipids in topology*/
1068 fprintf(stderr, "Calculating Scd order parameters for each of %d molecules\n", nslices);
1073 fprintf(stderr, "Calculating radial distances\n");
1076 gmx_fatal(FARGS, "Cannot yet output radial distances without permolecule\n");
1080 /* show atomtypes, to check if index file is correct */
1081 print_types(index, a, ngrps, grpname, top);
1083 calc_order(ftp2fn(efTRX, NFILE, fnm), index, a, &order,
1084 &slOrder, &slWidth, nslices, bSliced, bUnsat,
1085 top, ePBC, ngrps, axis, permolecule, radial, distcalc, opt2fn_null("-nr", NFILE, fnm), &distvals, oenv);
1089 ngrps--; /*don't print the last group--was used for
1090 center-of-mass determination*/
1093 order_plot(order, slOrder, opt2fn("-o", NFILE, fnm), opt2fn("-os", NFILE, fnm),
1094 opt2fn("-od", NFILE, fnm), ngrps, nslices, slWidth, bSzonly, permolecule, distvals, oenv);
1096 if (opt2bSet("-ob", NFILE, fnm))
1101 "Won't write B-factors with averaged order parameters; use -permolecule\n");
1105 write_bfactors(fnm, NFILE, index, a, nslices, ngrps, slOrder, top, distvals, oenv);
1110 do_view(oenv, opt2fn("-o", NFILE, fnm), NULL); /* view xvgr file */
1111 do_view(oenv, opt2fn("-os", NFILE, fnm), NULL); /* view xvgr file */
1112 do_view(oenv, opt2fn("-od", NFILE, fnm), NULL); /* view xvgr file */
1115 if (distvals != NULL)
1117 for (i = 0; i < nslices; ++i)