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, 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
,
79 int natoms
, matrix box
,
80 rvec x
[], int maxidx
, int index
[],
81 real
*sgmean
, real
*skmean
,
82 int nslice
, int slice_dim
,
83 real sgslice
[], real skslice
[],
86 int ix
, jx
, nsgbin
, *sgbin
;
87 int i
, ibin
, j
, k
, l
, n
, *nn
[4];
88 rvec dx
, rj
, rk
, urk
, urj
;
89 real cost
, cost2
, *sgmol
, *skmol
, rmean
, rmean2
, r2
, box2
, *r_nn
[4];
93 real onethird
= 1.0/3.0;
94 /* dmat = init_mat(maxidx, FALSE); */
95 box2
= box
[XX
][XX
] * box
[XX
][XX
];
96 snew(sl_count
, nslice
);
97 for (i
= 0; (i
< 4); i
++)
99 snew(r_nn
[i
], natoms
);
102 for (j
= 0; (j
< natoms
); j
++)
111 /* Must init pbc every step because of pressure coupling */
112 set_pbc(&pbc
, ePBC
, box
);
114 gmx_rmpbc(gpbc
, natoms
, box
, x
);
116 nsgbin
= 2001; // Calculated as (1 + 1/0.0005)
122 for (i
= 0; (i
< maxidx
); i
++) /* loop over index file */
125 for (j
= 0; (j
< maxidx
); j
++)
134 pbc_dx(&pbc
, x
[ix
], x
[jx
], dx
);
137 /* set_mat_entry(dmat,i,j,r2); */
139 /* determine the nearest neighbours */
142 r_nn
[3][i
] = r_nn
[2][i
]; nn
[3][i
] = nn
[2][i
];
143 r_nn
[2][i
] = r_nn
[1][i
]; nn
[2][i
] = nn
[1][i
];
144 r_nn
[1][i
] = r_nn
[0][i
]; nn
[1][i
] = nn
[0][i
];
145 r_nn
[0][i
] = r2
; nn
[0][i
] = j
;
147 else if (r2
< r_nn
[1][i
])
149 r_nn
[3][i
] = r_nn
[2][i
]; nn
[3][i
] = nn
[2][i
];
150 r_nn
[2][i
] = r_nn
[1][i
]; nn
[2][i
] = nn
[1][i
];
151 r_nn
[1][i
] = r2
; nn
[1][i
] = j
;
153 else if (r2
< r_nn
[2][i
])
155 r_nn
[3][i
] = r_nn
[2][i
]; nn
[3][i
] = nn
[2][i
];
156 r_nn
[2][i
] = r2
; nn
[2][i
] = j
;
158 else if (r2
< r_nn
[3][i
])
160 r_nn
[3][i
] = r2
; nn
[3][i
] = j
;
165 /* calculate mean distance between nearest neighbours */
167 for (j
= 0; (j
< 4); j
++)
169 r_nn
[j
][i
] = std::sqrt(r_nn
[j
][i
]);
178 /* Chau1998a eqn 3 */
179 /* angular part tetrahedrality order parameter per atom */
180 for (j
= 0; (j
< 3); j
++)
182 for (k
= j
+1; (k
< 4); k
++)
184 pbc_dx(&pbc
, x
[ix
], x
[index
[nn
[k
][i
]]], rk
);
185 pbc_dx(&pbc
, x
[ix
], x
[index
[nn
[j
][i
]]], rj
);
190 cost
= iprod(urk
, urj
) + onethird
;
193 /* sgmol[i] += 3*cost2/32; */
196 /* determine distribution */
197 ibin
= static_cast<int>(nsgbin
* cost2
);
202 /* printf("%d %d %f %d %d\n", j, k, cost * cost, ibin, sgbin[ibin]);*/
208 /* normalize sgmol between 0.0 and 1.0 */
209 sgmol
[i
] = 3*sgmol
[i
]/32;
212 /* distance part tetrahedrality order parameter per atom */
213 rmean2
= 4 * 3 * rmean
* rmean
;
214 for (j
= 0; (j
< 4); j
++)
216 skmol
[i
] += (rmean
- r_nn
[j
][i
]) * (rmean
- r_nn
[j
][i
]) / rmean2
;
217 /* printf("%d %f (%f %f %f %f) \n",
218 i, skmol[i], rmean, rmean2, r_nn[j][i], (rmean - r_nn[j][i]) );
224 /* Compute sliced stuff */
225 sl_index
= static_cast<int>(std::round((1+x
[i
][slice_dim
]/box
[slice_dim
][slice_dim
])*nslice
)) % nslice
;
226 sgslice
[sl_index
] += sgmol
[i
];
227 skslice
[sl_index
] += skmol
[i
];
228 sl_count
[sl_index
]++;
229 } /* loop over entries in index file */
234 for (i
= 0; (i
< nslice
); i
++)
238 sgslice
[i
] /= sl_count
[i
];
239 skslice
[i
] /= sl_count
[i
];
246 for (i
= 0; (i
< 4); i
++)
254 static void calc_tetra_order_parm(const char *fnNDX
, const char *fnTPS
,
255 const char *fnTRX
, const char *sgfn
,
257 int nslice
, int slice_dim
,
258 const char *sgslfn
, const char *skslfn
,
259 const gmx_output_env_t
*oenv
)
261 FILE *fpsg
= NULL
, *fpsk
= NULL
;
272 int i
, *isize
, ng
, nframes
;
273 real
*sg_slice
, *sg_slice_tot
, *sk_slice
, *sk_slice_tot
;
274 gmx_rmpbc_t gpbc
= NULL
;
277 read_tps_conf(fnTPS
, &top
, &ePBC
, &xtop
, NULL
, box
, FALSE
);
279 snew(sg_slice
, nslice
);
280 snew(sk_slice
, nslice
);
281 snew(sg_slice_tot
, nslice
);
282 snew(sk_slice_tot
, nslice
);
284 /* get index groups */
285 printf("Select the group that contains the atoms you want to use for the tetrahedrality order parameter calculation:\n");
289 get_index(&top
.atoms
, fnNDX
, ng
, isize
, index
, grpname
);
291 /* Analyze trajectory */
292 natoms
= read_first_x(oenv
, &status
, fnTRX
, &t
, &x
, box
);
293 if (natoms
> top
.atoms
.nr
)
295 gmx_fatal(FARGS
, "Topology (%d atoms) does not match trajectory (%d atoms)",
296 top
.atoms
.nr
, natoms
);
298 check_index(NULL
, ng
, index
[0], NULL
, natoms
);
300 fpsg
= xvgropen(sgfn
, "S\\sg\\N Angle Order Parameter", "Time (ps)", "S\\sg\\N",
302 fpsk
= xvgropen(skfn
, "S\\sk\\N Distance Order Parameter", "Time (ps)", "S\\sk\\N",
305 /* loop over frames */
306 gpbc
= gmx_rmpbc_init(&top
.idef
, ePBC
, natoms
);
310 find_nearest_neighbours(ePBC
, natoms
, box
, x
, isize
[0], index
[0],
311 &sg
, &sk
, nslice
, slice_dim
, sg_slice
, sk_slice
, gpbc
);
312 for (i
= 0; (i
< nslice
); i
++)
314 sg_slice_tot
[i
] += sg_slice
[i
];
315 sk_slice_tot
[i
] += sk_slice
[i
];
317 fprintf(fpsg
, "%f %f\n", t
, sg
);
318 fprintf(fpsk
, "%f %f\n", t
, sk
);
321 while (read_next_x(oenv
, status
, &t
, x
, box
));
323 gmx_rmpbc_done(gpbc
);
332 fpsg
= xvgropen(sgslfn
,
333 "S\\sg\\N Angle Order Parameter / Slab", "(nm)", "S\\sg\\N",
335 fpsk
= xvgropen(skslfn
,
336 "S\\sk\\N Distance Order Parameter / Slab", "(nm)", "S\\sk\\N",
338 for (i
= 0; (i
< nslice
); i
++)
340 fprintf(fpsg
, "%10g %10g\n", (i
+0.5)*box
[slice_dim
][slice_dim
]/nslice
,
341 sg_slice_tot
[i
]/nframes
);
342 fprintf(fpsk
, "%10g %10g\n", (i
+0.5)*box
[slice_dim
][slice_dim
]/nslice
,
343 sk_slice_tot
[i
]/nframes
);
350 /* Print name of first atom in all groups in index file */
351 static void print_types(int index
[], int a
[], int ngrps
,
352 char *groups
[], const t_topology
*top
)
356 fprintf(stderr
, "Using following groups: \n");
357 for (i
= 0; i
< ngrps
; i
++)
359 fprintf(stderr
, "Groupname: %s First atomname: %s First atomnr %d\n",
360 groups
[i
], *(top
->atoms
.atomname
[a
[index
[i
]]]), a
[index
[i
]]);
362 fprintf(stderr
, "\n");
365 static void check_length(real length
, int a
, int b
)
369 fprintf(stderr
, "WARNING: distance between atoms %d and "
370 "%d > 0.3 nm (%f). Index file might be corrupt.\n",
375 void calc_order(const char *fn
, int *index
, int *a
, rvec
**order
,
376 real
***slOrder
, real
*slWidth
, int nslices
, gmx_bool bSliced
,
377 gmx_bool bUnsat
, const t_topology
*top
, int ePBC
, int ngrps
, int axis
,
378 gmx_bool permolecule
, gmx_bool radial
, gmx_bool distcalc
, const char *radfn
,
380 const gmx_output_env_t
*oenv
)
382 /* if permolecule = TRUE, order parameters will be calculed per molecule
383 * and stored in slOrder with #slices = # molecules */
384 rvec
*x0
, /* coordinates with pbc */
385 *x1
, /* coordinates without pbc */
386 dist
; /* vector between two atoms */
387 matrix box
; /* box (3x3) */
389 rvec cossum
, /* sum of vector angles for three axes */
390 Sx
, Sy
, Sz
, /* the three molecular axes */
391 tmp1
, tmp2
, /* temp. rvecs for calculating dot products */
392 frameorder
; /* order parameters for one frame */
393 real
*slFrameorder
; /* order parameter for one frame, per slice */
394 real length
, /* total distance between two atoms */
395 t
, /* time from trajectory */
396 z_ave
, z1
, z2
; /* average z, used to det. which slice atom is in */
397 int natoms
, /* nr. atoms in trj */
398 nr_tails
, /* nr tails, to check if index file is correct */
399 size
= 0, /* nr. of atoms in group. same as nr_tails */
400 i
, j
, m
, k
, teller
= 0,
401 slice
, /* current slice number */
403 int *slCount
; /* nr. of atoms in one slice */
404 real sdbangle
= 0; /* sum of these angles */
405 gmx_bool use_unitvector
= FALSE
; /* use a specified unit vector instead of axis to specify unit normal*/
406 rvec direction
, com
, dref
, dvec
;
407 int comsize
, distsize
;
408 int *comidx
= NULL
, *distidx
= NULL
;
409 char *grpname
= NULL
;
411 real arcdist
, tmpdist
;
412 gmx_rmpbc_t gpbc
= NULL
;
414 /* PBC added for center-of-mass vector*/
415 /* Initiate the pbc structure */
416 std::memset(&pbc
, 0, sizeof(pbc
));
418 if ((natoms
= read_first_x(oenv
, &status
, fn
, &t
, &x0
, box
)) == 0)
420 gmx_fatal(FARGS
, "Could not read coordinates from statusfile\n");
423 nr_tails
= index
[1] - index
[0];
424 fprintf(stderr
, "Number of elements in first group: %d\n", nr_tails
);
425 /* take first group as standard. Not rocksolid, but might catch error in index*/
430 bSliced
= FALSE
; /*force slices off */
431 fprintf(stderr
, "Calculating order parameters for each of %d molecules\n",
437 use_unitvector
= TRUE
;
438 fprintf(stderr
, "Select an index group to calculate the radial membrane normal\n");
439 get_index(&top
->atoms
, radfn
, 1, &comsize
, &comidx
, &grpname
);
447 fprintf(stderr
, "Select an index group to use as distance reference\n");
448 get_index(&top
->atoms
, radfn
, 1, &distsize
, &distidx
, &grpname
);
449 bSliced
= FALSE
; /*force slices off*/
452 if (use_unitvector
&& bSliced
)
454 fprintf(stderr
, "Warning: slicing and specified unit vectors are not currently compatible\n");
457 snew(slCount
, nslices
);
458 snew(*slOrder
, nslices
);
459 for (i
= 0; i
< nslices
; i
++)
461 snew((*slOrder
)[i
], ngrps
);
465 snew(*distvals
, nslices
);
466 for (i
= 0; i
< nslices
; i
++)
468 snew((*distvals
)[i
], ngrps
);
472 snew(slFrameorder
, nslices
);
477 *slWidth
= box
[axis
][axis
]/nslices
;
478 fprintf(stderr
, "Box divided in %d slices. Initial width of slice: %f\n",
484 nr_tails
= index
[1] - index
[0];
485 fprintf(stderr
, "Number of elements in first group: %d\n", nr_tails
);
486 /* take first group as standard. Not rocksolid, but might catch error
492 gpbc
= gmx_rmpbc_init(&top
->idef
, ePBC
, natoms
);
493 /*********** Start processing trajectory ***********/
498 *slWidth
= box
[axis
][axis
]/nslices
;
502 set_pbc(&pbc
, ePBC
, box
);
503 gmx_rmpbc_copy(gpbc
, natoms
, box
, x0
, x1
);
505 /* Now loop over all groups. There are ngrps groups, the order parameter can
506 be calculated for grp 1 to grp ngrps - 1. For each group, loop over all
507 atoms in group, which is index[i] to (index[i+1] - 1) See block.h. Of
508 course, in this case index[i+1] -index[i] has to be the same for all
509 groups, namely the number of tails. i just runs over all atoms in a tail,
510 so for DPPC ngrps = 16 and i runs from 1 to 14, including 14
516 /*center-of-mass determination*/
517 com
[XX
] = 0.0; com
[YY
] = 0.0; com
[ZZ
] = 0.0;
518 for (j
= 0; j
< comsize
; j
++)
520 rvec_inc(com
, x1
[comidx
[j
]]);
522 svmul(1.0/comsize
, com
, com
);
526 dref
[XX
] = 0.0; dref
[YY
] = 0.0; dref
[ZZ
] = 0.0;
527 for (j
= 0; j
< distsize
; j
++)
529 rvec_inc(dist
, x1
[distidx
[j
]]);
531 svmul(1.0/distsize
, dref
, dref
);
534 pbc_dx(&pbc
, dref
, com
, dvec
);
539 for (i
= 1; i
< ngrps
- 1; i
++)
541 clear_rvec(frameorder
);
543 size
= index
[i
+1] - index
[i
];
544 if (size
!= nr_tails
)
546 gmx_fatal(FARGS
, "grp %d does not have same number of"
547 " elements as grp 1\n", i
);
550 for (j
= 0; j
< size
; j
++)
553 /*create unit vector*/
555 pbc_dx(&pbc
, x1
[a
[index
[i
]+j
]], com
, direction
);
556 unitv(direction
, direction
);
559 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],
560 direction[0],direction[1],direction[2]);*/
565 /* Using convention for unsaturated carbons */
566 /* first get Sz, the vector from Cn to Cn+1 */
567 rvec_sub(x1
[a
[index
[i
+1]+j
]], x1
[a
[index
[i
]+j
]], dist
);
569 check_length(length
, a
[index
[i
]+j
], a
[index
[i
+1]+j
]);
570 svmul(1.0/length
, dist
, Sz
);
572 /* this is actually the cosine of the angle between the double bond
573 and axis, because Sz is normalized and the two other components of
574 the axis on the bilayer are zero */
577 sdbangle
+= gmx_angle(direction
, Sz
); /*this can probably be optimized*/
581 sdbangle
+= std::acos(Sz
[axis
]);
586 /* get vector dist(Cn-1,Cn+1) for tail atoms */
587 rvec_sub(x1
[a
[index
[i
+1]+j
]], x1
[a
[index
[i
-1]+j
]], dist
);
588 length
= norm(dist
); /* determine distance between two atoms */
589 check_length(length
, a
[index
[i
-1]+j
], a
[index
[i
+1]+j
]);
591 svmul(1.0/length
, dist
, Sz
);
592 /* Sz is now the molecular axis Sz, normalized and all that */
595 /* now get Sx. Sx is normal to the plane of Cn-1, Cn and Cn+1 so
596 we can use the outer product of Cn-1->Cn and Cn+1->Cn, I hope */
597 rvec_sub(x1
[a
[index
[i
+1]+j
]], x1
[a
[index
[i
]+j
]], tmp1
);
598 rvec_sub(x1
[a
[index
[i
-1]+j
]], x1
[a
[index
[i
]+j
]], tmp2
);
599 cprod(tmp1
, tmp2
, Sx
);
600 svmul(1.0/norm(Sx
), Sx
, Sx
);
602 /* now we can get Sy from the outer product of Sx and Sz */
604 svmul(1.0/norm(Sy
), Sy
, Sy
);
606 /* the square of cosine of the angle between dist and the axis.
607 Using the innerproduct, but two of the three elements are zero
608 Determine the sum of the orderparameter of all atoms in group
612 cossum
[XX
] = gmx::square(iprod(Sx
, direction
)); /* this is allowed, since Sa is normalized */
613 cossum
[YY
] = gmx::square(iprod(Sy
, direction
));
614 cossum
[ZZ
] = gmx::square(iprod(Sz
, direction
));
618 cossum
[XX
] = gmx::square(Sx
[axis
]); /* this is allowed, since Sa is normalized */
619 cossum
[YY
] = gmx::square(Sy
[axis
]);
620 cossum
[ZZ
] = gmx::square(Sz
[axis
]);
623 for (m
= 0; m
< DIM
; m
++)
625 frameorder
[m
] += 0.5 * (3.0 * cossum
[m
] - 1.0);
630 /* get average coordinate in box length for slicing,
631 determine which slice atom is in, increase count for that
632 slice. slFrameorder and slOrder are reals, not
633 rvecs. Only the component [axis] of the order tensor is
634 kept, until I find it necessary to know the others too
637 z1
= x1
[a
[index
[i
-1]+j
]][axis
];
638 z2
= x1
[a
[index
[i
+1]+j
]][axis
];
639 z_ave
= 0.5 * (z1
+ z2
);
642 z_ave
+= box
[axis
][axis
];
644 if (z_ave
> box
[axis
][axis
])
646 z_ave
-= box
[axis
][axis
];
649 slice
= static_cast<int>((0.5 + (z_ave
/ (*slWidth
))) - 1);
650 slCount
[slice
]++; /* determine slice, increase count */
652 slFrameorder
[slice
] += 0.5 * (3 * cossum
[axis
] - 1);
654 else if (permolecule
)
656 /* store per-molecule order parameter
657 * To just track single-axis order: (*slOrder)[j][i] += 0.5 * (3 * iprod(cossum,direction) - 1);
658 * following is for Scd order: */
659 (*slOrder
)[j
][i
] += -1* (1.0/3.0 * (3 * cossum
[XX
] - 1) + 1.0/3.0 * 0.5 * (3.0 * cossum
[YY
] - 1));
665 /* bin order parameter by arc distance from reference group*/
666 arcdist
= gmx_angle(dvec
, direction
);
667 (*distvals
)[j
][i
] += arcdist
;
671 /* Want minimum lateral distance to first group calculated */
672 tmpdist
= trace(box
); /* should be max value */
673 for (k
= 0; k
< distsize
; k
++)
675 pbc_dx(&pbc
, x1
[distidx
[k
]], x1
[a
[index
[i
]+j
]], dvec
);
676 /* at the moment, just remove dvec[axis] */
678 tmpdist
= std::min(tmpdist
, norm2(dvec
));
680 //fprintf(stderr, "Min dist %f; trace %f\n", tmpdist, trace(box));
681 (*distvals
)[j
][i
] += std::sqrt(tmpdist
);
684 } /* end loop j, over all atoms in group */
686 for (m
= 0; m
< DIM
; m
++)
688 (*order
)[i
][m
] += (frameorder
[m
]/size
);
692 { /*Skip following if doing per-molecule*/
693 for (k
= 0; k
< nslices
; k
++)
695 if (slCount
[k
]) /* if no elements, nothing has to be added */
697 (*slOrder
)[k
][i
] += slFrameorder
[k
]/slCount
[k
];
698 slFrameorder
[k
] = 0; slCount
[k
] = 0;
701 } /* end loop i, over all groups in indexfile */
706 while (read_next_x(oenv
, status
, &t
, x0
, box
));
707 /*********** done with status file **********/
709 fprintf(stderr
, "\nRead trajectory. Printing parameters to file\n");
710 gmx_rmpbc_done(gpbc
);
712 /* average over frames */
713 for (i
= 1; i
< ngrps
- 1; i
++)
715 svmul(1.0/nr_frames
, (*order
)[i
], (*order
)[i
]);
716 fprintf(stderr
, "Atom %d Tensor: x=%g , y=%g, z=%g\n", i
, (*order
)[i
][XX
],
717 (*order
)[i
][YY
], (*order
)[i
][ZZ
]);
718 if (bSliced
|| permolecule
)
720 for (k
= 0; k
< nslices
; k
++)
722 (*slOrder
)[k
][i
] /= nr_frames
;
727 for (k
= 0; k
< nslices
; k
++)
729 (*distvals
)[k
][i
] /= nr_frames
;
736 fprintf(stderr
, "Average angle between double bond and normal: %f\n",
737 180*sdbangle
/(nr_frames
* size
*M_PI
));
740 sfree(x0
); /* free memory used by coordinate arrays */
757 void order_plot(rvec order
[], real
*slOrder
[], const char *afile
, const char *bfile
,
758 const char *cfile
, int ngrps
, int nslices
, real slWidth
, gmx_bool bSzonly
,
759 gmx_bool permolecule
, real
**distvals
, const gmx_output_env_t
*oenv
)
761 FILE *ord
, *slOrd
; /* xvgr files with order parameters */
762 int atom
, slice
; /* atom corresponding to order para.*/
763 char buf
[256]; /* for xvgr title */
764 real S
; /* order parameter averaged over all atoms */
768 sprintf(buf
, "Scd order parameters");
769 ord
= xvgropen(afile
, buf
, "Atom", "S", oenv
);
770 sprintf(buf
, "Orderparameters per atom per slice");
771 slOrd
= xvgropen(bfile
, buf
, "Molecule", "S", oenv
);
772 for (atom
= 1; atom
< ngrps
- 1; atom
++)
774 fprintf(ord
, "%12d %12g\n", atom
, -1.0 * (2.0/3.0 * order
[atom
][XX
] +
775 1.0/3.0 * order
[atom
][YY
]));
778 for (slice
= 0; slice
< nslices
; slice
++)
780 fprintf(slOrd
, "%12d\t", slice
);
783 fprintf(slOrd
, "%12g\t", distvals
[slice
][1]); /*use distance value at second carbon*/
785 for (atom
= 1; atom
< ngrps
- 1; atom
++)
787 fprintf(slOrd
, "%12g\t", slOrder
[slice
][atom
]);
789 fprintf(slOrd
, "\n");
795 sprintf(buf
, "Orderparameters Sz per atom");
796 ord
= xvgropen(afile
, buf
, "Atom", "S", oenv
);
797 fprintf(stderr
, "ngrps = %d, nslices = %d", ngrps
, nslices
);
799 sprintf(buf
, "Orderparameters per atom per slice");
800 slOrd
= xvgropen(bfile
, buf
, "Slice", "S", oenv
);
802 for (atom
= 1; atom
< ngrps
- 1; atom
++)
804 fprintf(ord
, "%12d %12g\n", atom
, order
[atom
][ZZ
]);
807 for (slice
= 0; slice
< nslices
; slice
++)
810 for (atom
= 1; atom
< ngrps
- 1; atom
++)
812 S
+= slOrder
[slice
][atom
];
814 fprintf(slOrd
, "%12g %12g\n", slice
*slWidth
, S
/atom
);
820 sprintf(buf
, "Order tensor diagonal elements");
821 ord
= xvgropen(afile
, buf
, "Atom", "S", oenv
);
822 sprintf(buf
, "Deuterium order parameters");
823 slOrd
= xvgropen(cfile
, buf
, "Atom", "Scd", oenv
);
825 for (atom
= 1; atom
< ngrps
- 1; atom
++)
827 fprintf(ord
, "%12d %12g %12g %12g\n", atom
, order
[atom
][XX
],
828 order
[atom
][YY
], order
[atom
][ZZ
]);
829 fprintf(slOrd
, "%12d %12g\n", atom
, -1.0 * (2.0/3.0 * order
[atom
][XX
] +
830 1.0/3.0 * order
[atom
][YY
]));
838 void write_bfactors(t_filenm
*fnm
, int nfile
, int *index
, int *a
, int nslices
, int ngrps
, real
**order
, const t_topology
*top
, real
**distvals
, gmx_output_env_t
*oenv
)
840 /*function to write order parameters as B factors in PDB file using
841 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 read_first_frame(oenv
, &status
, ftp2fn(efTRX
, nfile
, fnm
), &fr
, TRX_NEED_X
);
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
= std::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 done_atom(&useatoms
);
896 int gmx_order(int argc
, char *argv
[])
898 const char *desc
[] = {
899 "[THISMODULE] computes 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 "[THISMODULE] 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.",
919 static int nslices
= 1; /* nr of slices defined */
920 static gmx_bool bSzonly
= FALSE
; /* True if only Sz is wanted */
921 static gmx_bool bUnsat
= FALSE
; /* True if carbons are unsat. */
922 static const char *normal_axis
[] = { NULL
, "z", "x", "y", NULL
};
923 static gmx_bool permolecule
= FALSE
; /*compute on a per-molecule basis */
924 static gmx_bool radial
= FALSE
; /*compute a radial membrane normal */
925 static gmx_bool distcalc
= FALSE
; /*calculate distance from a reference group */
927 { "-d", FALSE
, etENUM
, {normal_axis
},
928 "Direction of the normal on the membrane" },
929 { "-sl", FALSE
, etINT
, {&nslices
},
930 "Calculate order parameter as function of box length, dividing the box"
931 " into this number of slices." },
932 { "-szonly", FALSE
, etBOOL
, {&bSzonly
},
933 "Only give Sz element of order tensor. (axis can be specified with [TT]-d[tt])" },
934 { "-unsat", FALSE
, etBOOL
, {&bUnsat
},
935 "Calculate order parameters for unsaturated carbons. Note that this can"
936 "not be mixed with normal order parameters." },
937 { "-permolecule", FALSE
, etBOOL
, {&permolecule
},
938 "Compute per-molecule Scd order parameters" },
939 { "-radial", FALSE
, etBOOL
, {&radial
},
940 "Compute a radial membrane normal" },
941 { "-calcdist", FALSE
, etBOOL
, {&distcalc
},
942 "Compute distance from a reference" },
945 rvec
*order
; /* order par. for each atom */
946 real
**slOrder
; /* same, per slice */
947 real slWidth
= 0.0; /* width of a slice */
948 char **grpname
; /* groupnames */
949 int ngrps
, /* nr. of groups */
951 axis
= 0; /* normal axis */
952 t_topology
*top
; /* topology */
954 int *index
, /* indices for a */
955 *a
; /* atom numbers in each group */
956 t_blocka
*block
; /* data from index file */
957 t_filenm fnm
[] = { /* files for g_order */
958 { efTRX
, "-f", NULL
, ffREAD
}, /* trajectory file */
959 { efNDX
, "-n", NULL
, ffREAD
}, /* index file */
960 { efNDX
, "-nr", NULL
, ffOPTRD
}, /* index for radial axis calculation */
961 { efTPR
, NULL
, NULL
, ffREAD
}, /* topology file */
962 { efXVG
, "-o", "order", ffWRITE
}, /* xvgr output file */
963 { efXVG
, "-od", "deuter", ffWRITE
}, /* xvgr output file */
964 { efPDB
, "-ob", NULL
, ffOPTWR
}, /* write Scd as B factors to PDB if permolecule */
965 { efXVG
, "-os", "sliced", ffWRITE
}, /* xvgr output file */
966 { efXVG
, "-Sg", "sg-ang", ffOPTWR
}, /* xvgr output file */
967 { efXVG
, "-Sk", "sk-dist", ffOPTWR
}, /* xvgr output file */
968 { efXVG
, "-Sgsl", "sg-ang-slice", ffOPTWR
}, /* xvgr output file */
969 { efXVG
, "-Sksl", "sk-dist-slice", ffOPTWR
}, /* xvgr output file */
971 gmx_bool bSliced
= FALSE
; /* True if box is sliced */
972 #define NFILE asize(fnm)
973 real
**distvals
= NULL
;
974 const char *sgfnm
, *skfnm
, *ndxfnm
, *tpsfnm
, *trxfnm
;
975 gmx_output_env_t
*oenv
;
977 if (!parse_common_args(&argc
, argv
, PCA_CAN_VIEW
| PCA_CAN_TIME
,
978 NFILE
, fnm
, asize(pa
), pa
, asize(desc
), desc
, 0, NULL
, &oenv
))
984 gmx_fatal(FARGS
, "Can not have nslices < 1");
986 sgfnm
= opt2fn_null("-Sg", NFILE
, fnm
);
987 skfnm
= opt2fn_null("-Sk", NFILE
, fnm
);
988 ndxfnm
= opt2fn_null("-n", NFILE
, fnm
);
989 tpsfnm
= ftp2fn(efTPR
, NFILE
, fnm
);
990 trxfnm
= ftp2fn(efTRX
, NFILE
, fnm
);
993 GMX_RELEASE_ASSERT(normal_axis
[0] != NULL
, "Options inconsistency; normal_axis[0] is NULL");
994 if (std::strcmp(normal_axis
[0], "x") == 0)
998 else if (std::strcmp(normal_axis
[0], "y") == 0)
1002 else if (std::strcmp(normal_axis
[0], "z") == 0)
1008 gmx_fatal(FARGS
, "Invalid axis, use x, y or z");
1014 fprintf(stderr
, "Taking x axis as normal to the membrane\n");
1017 fprintf(stderr
, "Taking y axis as normal to the membrane\n");
1020 fprintf(stderr
, "Taking z axis as normal to the membrane\n");
1024 /* tetraheder order parameter */
1027 /* If either of theoptions is set we compute both */
1028 sgfnm
= opt2fn("-Sg", NFILE
, fnm
);
1029 skfnm
= opt2fn("-Sk", NFILE
, fnm
);
1030 calc_tetra_order_parm(ndxfnm
, tpsfnm
, trxfnm
, sgfnm
, skfnm
, nslices
, axis
,
1031 opt2fn("-Sgsl", NFILE
, fnm
), opt2fn("-Sksl", NFILE
, fnm
),
1033 /* view xvgr files */
1034 do_view(oenv
, opt2fn("-Sg", NFILE
, fnm
), NULL
);
1035 do_view(oenv
, opt2fn("-Sk", NFILE
, fnm
), NULL
);
1038 do_view(oenv
, opt2fn("-Sgsl", NFILE
, fnm
), NULL
);
1039 do_view(oenv
, opt2fn("-Sksl", NFILE
, fnm
), NULL
);
1044 /* tail order parameter */
1049 fprintf(stderr
, "Dividing box in %d slices.\n\n", nslices
);
1054 fprintf(stderr
, "Only calculating Sz\n");
1058 fprintf(stderr
, "Taking carbons as unsaturated!\n");
1061 top
= read_top(ftp2fn(efTPR
, NFILE
, fnm
), &ePBC
); /* read topology file */
1063 block
= init_index(ftp2fn(efNDX
, NFILE
, fnm
), &grpname
);
1064 index
= block
->index
; /* get indices from t_block block */
1065 a
= block
->a
; /* see block.h */
1070 nslices
= index
[1] - index
[0]; /*I think this assumes contiguous lipids in topology*/
1071 fprintf(stderr
, "Calculating Scd order parameters for each of %d molecules\n", nslices
);
1076 fprintf(stderr
, "Calculating radial distances\n");
1079 gmx_fatal(FARGS
, "Cannot yet output radial distances without permolecule\n");
1083 /* show atomtypes, to check if index file is correct */
1084 print_types(index
, a
, ngrps
, grpname
, top
);
1086 calc_order(ftp2fn(efTRX
, NFILE
, fnm
), index
, a
, &order
,
1087 &slOrder
, &slWidth
, nslices
, bSliced
, bUnsat
,
1088 top
, ePBC
, ngrps
, axis
, permolecule
, radial
, distcalc
, opt2fn_null("-nr", NFILE
, fnm
), &distvals
, oenv
);
1092 ngrps
--; /*don't print the last group--was used for
1093 center-of-mass determination*/
1096 order_plot(order
, slOrder
, opt2fn("-o", NFILE
, fnm
), opt2fn("-os", NFILE
, fnm
),
1097 opt2fn("-od", NFILE
, fnm
), ngrps
, nslices
, slWidth
, bSzonly
, permolecule
, distvals
, oenv
);
1099 if (opt2bSet("-ob", NFILE
, fnm
))
1104 "Won't write B-factors with averaged order parameters; use -permolecule\n");
1108 write_bfactors(fnm
, NFILE
, index
, a
, nslices
, ngrps
, slOrder
, top
, distvals
, oenv
);
1113 do_view(oenv
, opt2fn("-o", NFILE
, fnm
), NULL
); /* view xvgr file */
1114 do_view(oenv
, opt2fn("-os", NFILE
, fnm
), NULL
); /* view xvgr file */
1115 do_view(oenv
, opt2fn("-od", NFILE
, fnm
), NULL
); /* view xvgr file */
1118 if (distvals
!= NULL
)
1120 for (i
= 0; i
< nslices
; ++i
)