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
61 /****************************************************************************/
62 /* This program calculates the order parameter per atom for an interface or */
63 /* bilayer, averaged over time. */
64 /* S = 1/2 * (3 * cos(i)cos(j) - delta(ij)) */
65 /* It is assumed that the order parameter with respect to a box-axis */
66 /* is calculated. In that case i = j = axis, and delta(ij) = 1. */
68 /* Peter Tieleman, April 1995 */
69 /* P.J. van Maaren, November 2005 Added tetrahedral stuff */
70 /****************************************************************************/
72 static void find_nearest_neighbours(t_topology top
, int ePBC
,
73 int natoms
, matrix box
,
74 rvec x
[],int maxidx
,atom_id index
[],
76 real
*sgmean
, real
*skmean
,
77 int nslice
,int slice_dim
,
78 real sgslice
[],real skslice
[])
82 int ix
,jx
,nsgbin
, *sgbin
;
83 int i1
,i2
,i
,ibin
,j
,k
,l
,n
,*nn
[4];
84 rvec dx
,dx1
,dx2
,rj
,rk
,urk
,urj
;
85 real cost
,cost2
,*sgmol
,*skmol
,rmean
,rmean2
,r2
,box2
,*r_nn
[4];
91 real onethird
=1.0/3.0;
93 /* 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
++) {
101 for (j
=0; (j
<natoms
); j
++) {
109 /* Must init pbc every step because of pressure coupling */
110 set_pbc(&pbc
,ePBC
,box
);
111 rm_pbc(&(top
.idef
),ePBC
,natoms
,box
,x
,x
);
113 nsgbin
= 1 + 1/0.0005;
119 for (i
=0; (i
<maxidx
); i
++) { /* loop over index file */
121 for (j
=0; (j
<maxidx
); j
++) {
122 if (i
== j
) continue;
126 pbc_dx(&pbc
,x
[ix
],x
[jx
],dx
);
129 /* set_mat_entry(dmat,i,j,r2); */
131 /* determine the nearest neighbours */
132 if (r2
< r_nn
[0][i
]) {
133 r_nn
[3][i
] = r_nn
[2][i
]; nn
[3][i
] = nn
[2][i
];
134 r_nn
[2][i
] = r_nn
[1][i
]; nn
[2][i
] = nn
[1][i
];
135 r_nn
[1][i
] = r_nn
[0][i
]; nn
[1][i
] = nn
[0][i
];
136 r_nn
[0][i
] = r2
; nn
[0][i
] = j
;
137 } else if (r2
< r_nn
[1][i
]) {
138 r_nn
[3][i
] = r_nn
[2][i
]; nn
[3][i
] = nn
[2][i
];
139 r_nn
[2][i
] = r_nn
[1][i
]; nn
[2][i
] = nn
[1][i
];
140 r_nn
[1][i
] = r2
; nn
[1][i
] = j
;
141 } else if (r2
< r_nn
[2][i
]) {
142 r_nn
[3][i
] = r_nn
[2][i
]; nn
[3][i
] = nn
[2][i
];
143 r_nn
[2][i
] = r2
; nn
[2][i
] = j
;
144 } else if (r2
< r_nn
[3][i
]) {
145 r_nn
[3][i
] = r2
; nn
[3][i
] = j
;
150 /* calculate mean distance between nearest neighbours */
152 for (j
=0; (j
<4); j
++) {
153 r_nn
[j
][i
] = sqrt(r_nn
[j
][i
]);
162 /* Chau1998a eqn 3 */
163 /* angular part tetrahedrality order parameter per atom */
164 for (j
=0; (j
<3); j
++) {
165 for (k
=j
+1; (k
<4); k
++) {
166 pbc_dx(&pbc
,x
[ix
],x
[index
[nn
[k
][i
]]],rk
);
167 pbc_dx(&pbc
,x
[ix
],x
[index
[nn
[j
][i
]]],rj
);
172 cost
= iprod(urk
,urj
) + onethird
;
175 /* sgmol[i] += 3*cost2/32; */
178 /* determine distribution */
179 ibin
= nsgbin
* cost2
;
182 /* printf("%d %d %f %d %d\n", j, k, cost * cost, ibin, sgbin[ibin]);*/
188 /* normalize sgmol between 0.0 and 1.0 */
189 sgmol
[i
] = 3*sgmol
[i
]/32;
192 /* distance part tetrahedrality order parameter per atom */
193 rmean2
= 4 * 3 * rmean
* rmean
;
194 for (j
=0; (j
<4); j
++) {
195 skmol
[i
] += (rmean
- r_nn
[j
][i
]) * (rmean
- r_nn
[j
][i
]) / rmean2
;
196 /* printf("%d %f (%f %f %f %f) \n",
197 i, skmol[i], rmean, rmean2, r_nn[j][i], (rmean - r_nn[j][i]) );
203 /* Compute sliced stuff */
204 sl_index
= gmx_nint((1+x
[i
][slice_dim
]/box
[slice_dim
][slice_dim
])*nslice
) % nslice
;
205 sgslice
[sl_index
] += sgmol
[i
];
206 skslice
[sl_index
] += skmol
[i
];
207 sl_count
[sl_index
]++;
208 } /* loop over entries in index file */
213 for(i
=0; (i
<nslice
); i
++) {
214 if (sl_count
[i
] > 0) {
215 sgslice
[i
] /= sl_count
[i
];
216 skslice
[i
] /= sl_count
[i
];
223 for (i
=0; (i
<4); i
++) {
230 static void calc_tetra_order_parm(const char *fnNDX
,const char *fnTPS
,
231 const char *fnTRX
, const char *sgfn
,
233 int nslice
,int slice_dim
,
234 const char *sgslfn
,const char *skslfn
,
235 const output_env_t oenv
)
237 FILE *fpsg
=NULL
,*fpsk
=NULL
;
240 char title
[STRLEN
],fn
[STRLEN
],subtitle
[STRLEN
];
249 int i
,*isize
,ng
,nframes
;
250 real
*sg_slice
,*sg_slice_tot
,*sk_slice
,*sk_slice_tot
;
252 read_tps_conf(fnTPS
,title
,&top
,&ePBC
,&xtop
,NULL
,box
,FALSE
);
254 snew(sg_slice
,nslice
);
255 snew(sk_slice
,nslice
);
256 snew(sg_slice_tot
,nslice
);
257 snew(sk_slice_tot
,nslice
);
259 /* get index groups */
260 printf("Select the group that contains the atoms you want to use for the tetrahedrality order parameter calculation:\n");
264 get_index(&top
.atoms
,fnNDX
,ng
,isize
,index
,grpname
);
266 /* Analyze trajectory */
267 natoms
=read_first_x(oenv
,&status
,fnTRX
,&t
,&x
,box
);
268 if ( natoms
> top
.atoms
.nr
)
269 gmx_fatal(FARGS
,"Topology (%d atoms) does not match trajectory (%d atoms)",
270 top
.atoms
.nr
,natoms
);
271 check_index(NULL
,ng
,index
[0],NULL
,natoms
);
273 fpsg
=xvgropen(sgfn
,"S\\sg\\N Angle Order Parameter","Time (ps)","S\\sg\\N",
275 fpsk
=xvgropen(skfn
,"S\\sk\\N Distance Order Parameter","Time (ps)","S\\sk\\N",
278 /* loop over frames */
281 find_nearest_neighbours(top
,ePBC
,natoms
,box
,x
,isize
[0],index
[0],t
,
282 &sg
,&sk
,nslice
,slice_dim
,sg_slice
,sk_slice
);
283 for(i
=0; (i
<nslice
); i
++) {
284 sg_slice_tot
[i
] += sg_slice
[i
];
285 sk_slice_tot
[i
] += sk_slice
[i
];
287 fprintf(fpsg
,"%f %f\n", t
, sg
);
288 fprintf(fpsk
,"%f %f\n", t
, sk
);
290 } while (read_next_x(oenv
,status
,&t
,natoms
,x
,box
));
300 fpsg
= xvgropen(sgslfn
,
301 "S\\sg\\N Angle Order Parameter / Slab","(nm)","S\\sg\\N",
303 fpsk
= xvgropen(skslfn
,
304 "S\\sk\\N Distance Order Parameter / Slab","(nm)","S\\sk\\N",
306 for(i
=0; (i
<nslice
); i
++) {
307 fprintf(fpsg
,"%10g %10g\n",(i
+0.5)*box
[slice_dim
][slice_dim
]/nslice
,
308 sg_slice_tot
[i
]/nframes
);
309 fprintf(fpsk
,"%10g %10g\n",(i
+0.5)*box
[slice_dim
][slice_dim
]/nslice
,
310 sk_slice_tot
[i
]/nframes
);
317 /* Print name of first atom in all groups in index file */
318 static void print_types(atom_id index
[], atom_id a
[], int ngrps
,
319 char *groups
[], t_topology
*top
)
323 fprintf(stderr
,"Using following groups: \n");
324 for(i
= 0; i
< ngrps
; i
++)
325 fprintf(stderr
,"Groupname: %s First atomname: %s First atomnr %u\n",
326 groups
[i
], *(top
->atoms
.atomname
[a
[index
[i
]]]), a
[index
[i
]]);
327 fprintf(stderr
,"\n");
330 static void check_length(real length
, int a
, int b
)
333 fprintf(stderr
,"WARNING: distance between atoms %d and "
334 "%d > 0.3 nm (%f). Index file might be corrupt.\n",
338 void calc_order(const char *fn
, atom_id
*index
, atom_id
*a
, rvec
**order
,
339 real
***slOrder
, real
*slWidth
, int nslices
, bool bSliced
,
340 bool bUnsat
, t_topology
*top
, int ePBC
, int ngrps
, int axis
,
341 bool permolecule
, bool radial
, bool distcalc
, const char *radfn
,
343 const output_env_t oenv
)
345 /* if permolecule = TRUE, order parameters will be calculed per molecule
346 * and stored in slOrder with #slices = # molecules */
347 rvec
*x0
, /* coordinates with pbc */
348 *x1
, /* coordinates without pbc */
349 dist
; /* vector between two atoms */
350 matrix box
; /* box (3x3) */
352 rvec cossum
, /* sum of vector angles for three axes */
353 Sx
, Sy
, Sz
, /* the three molecular axes */
354 tmp1
, tmp2
, /* temp. rvecs for calculating dot products */
355 frameorder
; /* order parameters for one frame */
356 real
*slFrameorder
; /* order parameter for one frame, per slice */
357 real length
, /* total distance between two atoms */
358 t
, /* time from trajectory */
359 z_ave
,z1
,z2
; /* average z, used to det. which slice atom is in */
360 int natoms
, /* nr. atoms in trj */
361 nr_tails
, /* nr tails, to check if index file is correct */
362 size
=0, /* nr. of atoms in group. same as nr_tails */
363 i
,j
,m
,k
,l
,teller
= 0,
364 slice
, /* current slice number */
366 *slCount
; /* nr. of atoms in one slice */
367 real dbangle
= 0, /* angle between double bond and axis */
368 sdbangle
= 0;/* sum of these angles */
369 bool use_unitvector
= FALSE
; /* use a specified unit vector instead of axis to specify unit normal*/
370 rvec direction
, com
,dref
,dvec
;
371 int comsize
, distsize
;
372 atom_id
*comidx
=NULL
, *distidx
=NULL
;
377 /* PBC added for center-of-mass vector*/
378 /* Initiate the pbc structure */
379 memset(&pbc
,0,sizeof(pbc
));
381 if ((natoms
= read_first_x(oenv
,&status
,fn
,&t
,&x0
,box
)) == 0)
382 gmx_fatal(FARGS
,"Could not read coordinates from statusfile\n");
384 nr_tails
= index
[1] - index
[0];
385 fprintf(stderr
,"Number of elements in first group: %d\n",nr_tails
);
386 /* take first group as standard. Not rocksolid, but might catch error in index*/
391 bSliced
=FALSE
; /*force slices off */
392 fprintf(stderr
,"Calculating order parameters for each of %d molecules\n",
399 fprintf(stderr
,"Select an index group to calculate the radial membrane normal\n");
400 get_index(&top
->atoms
,radfn
,1,&comsize
,&comidx
,&grpname
);
405 fprintf(stderr
,"Select an index group to use as distance reference\n");
406 get_index(&top
->atoms
,radfn
,1,&distsize
,&distidx
,&grpname
);
407 bSliced
=FALSE
; /*force slices off*/
411 if (use_unitvector
&& bSliced
)
412 fprintf(stderr
,"Warning: slicing and specified unit vectors are not currently compatible\n");
414 snew(slCount
, nslices
);
415 snew(*slOrder
, nslices
);
416 for(i
= 0; i
< nslices
; i
++)
417 snew((*slOrder
)[i
],ngrps
);
420 snew(*distvals
, nslices
);
421 for(i
= 0; i
< nslices
; i
++)
422 snew((*distvals
)[i
],ngrps
);
425 snew(slFrameorder
, nslices
);
429 *slWidth
= box
[axis
][axis
]/nslices
;
430 fprintf(stderr
,"Box divided in %d slices. Initial width of slice: %f\n",
435 nr_tails
= index
[1] - index
[0];
436 fprintf(stderr
,"Number of elements in first group: %d\n",nr_tails
);
437 /* take first group as standard. Not rocksolid, but might catch error
443 /*********** Start processing trajectory ***********/
446 *slWidth
= box
[axis
][axis
]/nslices
;
449 set_pbc(&pbc
,ePBC
,box
);
450 rm_pbc(&(top
->idef
),ePBC
,top
->atoms
.nr
,box
,x0
,x1
);
452 /* Now loop over all groups. There are ngrps groups, the order parameter can
453 be calculated for grp 1 to grp ngrps - 1. For each group, loop over all
454 atoms in group, which is index[i] to (index[i+1] - 1) See block.h. Of
455 course, in this case index[i+1] -index[i] has to be the same for all
456 groups, namely the number of tails. i just runs over all atoms in a tail,
457 so for DPPC ngrps = 16 and i runs from 1 to 14, including 14
462 /*center-of-mass determination*/
463 com
[XX
]=0.0; com
[YY
]=0.0; com
[ZZ
]=0.0;
464 for (j
=0;j
<comsize
;j
++)
465 rvec_inc(com
,x1
[comidx
[j
]]);
466 svmul(1.0/comsize
,com
,com
);
470 dref
[XX
]=0.0; dref
[YY
]=0.0;dref
[ZZ
]=0.0;
471 for (j
=0;j
<distsize
;j
++)
472 rvec_inc(dist
,x1
[distidx
[j
]]);
473 svmul(1.0/distsize
,dref
,dref
);
474 pbc_dx(&pbc
,dref
,com
,dvec
);
478 for (i
= 1; i
< ngrps
- 1; i
++) {
479 clear_rvec(frameorder
);
481 size
= index
[i
+1] - index
[i
];
482 if (size
!= nr_tails
)
483 gmx_fatal(FARGS
,"grp %d does not have same number of"
484 " elements as grp 1\n",i
);
486 for (j
= 0; j
< size
; j
++) {
488 /*create unit vector*/
490 pbc_dx(&pbc
,x1
[a
[index
[i
]+j
]],com
,direction
);
491 unitv(direction
,direction
);
494 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],
495 direction[0],direction[1],direction[2]);*/
499 /* Using convention for unsaturated carbons */
500 /* first get Sz, the vector from Cn to Cn+1 */
501 rvec_sub(x1
[a
[index
[i
+1]+j
]], x1
[a
[index
[i
]+j
]], dist
);
503 check_length(length
, a
[index
[i
]+j
], a
[index
[i
+1]+j
]);
504 svmul(1/length
, dist
, Sz
);
506 /* this is actually the cosine of the angle between the double bond
507 and axis, because Sz is normalized and the two other components of
508 the axis on the bilayer are zero */
511 sdbangle
+= acos(iprod(direction
,Sz
)); /*this can probably be optimized*/
514 sdbangle
+= acos(Sz
[axis
]);
516 /* get vector dist(Cn-1,Cn+1) for tail atoms */
517 rvec_sub(x1
[a
[index
[i
+1]+j
]], x1
[a
[index
[i
-1]+j
]], dist
);
518 length
= norm(dist
); /* determine distance between two atoms */
519 check_length(length
, a
[index
[i
-1]+j
], a
[index
[i
+1]+j
]);
521 svmul(1/length
, dist
, Sz
);
522 /* Sz is now the molecular axis Sz, normalized and all that */
525 /* now get Sx. Sx is normal to the plane of Cn-1, Cn and Cn+1 so
526 we can use the outer product of Cn-1->Cn and Cn+1->Cn, I hope */
527 rvec_sub(x1
[a
[index
[i
+1]+j
]], x1
[a
[index
[i
]+j
]], tmp1
);
528 rvec_sub(x1
[a
[index
[i
-1]+j
]], x1
[a
[index
[i
]+j
]], tmp2
);
529 cprod(tmp1
, tmp2
, Sx
);
530 svmul(1/norm(Sx
), Sx
, Sx
);
532 /* now we can get Sy from the outer product of Sx and Sz */
534 svmul(1/norm(Sy
), Sy
, Sy
);
536 /* the square of cosine of the angle between dist and the axis.
537 Using the innerproduct, but two of the three elements are zero
538 Determine the sum of the orderparameter of all atoms in group
542 cossum
[XX
] = sqr(iprod(Sx
,direction
)); /* this is allowed, since Sa is normalized */
543 cossum
[YY
] = sqr(iprod(Sy
,direction
));
544 cossum
[ZZ
] = sqr(iprod(Sz
,direction
));
548 cossum
[XX
] = sqr(Sx
[axis
]); /* this is allowed, since Sa is normalized */
549 cossum
[YY
] = sqr(Sy
[axis
]);
550 cossum
[ZZ
] = sqr(Sz
[axis
]);
553 for (m
= 0; m
< DIM
; m
++)
554 frameorder
[m
] += 0.5 * (3 * cossum
[m
] - 1);
557 /* get average coordinate in box length for slicing,
558 determine which slice atom is in, increase count for that
559 slice. slFrameorder and slOrder are reals, not
560 rvecs. Only the component [axis] of the order tensor is
561 kept, until I find it necessary to know the others too
564 z1
= x1
[a
[index
[i
-1]+j
]][axis
];
565 z2
= x1
[a
[index
[i
+1]+j
]][axis
];
566 z_ave
= 0.5 * (z1
+ z2
);
568 z_ave
+= box
[axis
][axis
];
569 if (z_ave
> box
[axis
][axis
])
570 z_ave
-= box
[axis
][axis
];
572 slice
= (int)(0.5 + (z_ave
/ (*slWidth
))) - 1;
573 slCount
[slice
]++; /* determine slice, increase count */
575 slFrameorder
[slice
] += 0.5 * (3 * cossum
[axis
] - 1);
577 else if (permolecule
)
579 /* store per-molecule order parameter
580 * To just track single-axis order: (*slOrder)[j][i] += 0.5 * (3 * iprod(cossum,direction) - 1);
581 * following is for Scd order: */
582 (*slOrder
)[j
][i
] += -1* (0.3333 * (3 * cossum
[XX
] - 1) + 0.3333 * 0.5 * (3 * cossum
[YY
] - 1));
586 /* bin order parameter by arc distance from reference group*/
587 arcdist
=acos(iprod(dvec
,direction
));
588 (*distvals
)[j
][i
]+=arcdist
;
590 } /* end loop j, over all atoms in group */
592 for (m
= 0; m
< DIM
; m
++)
593 (*order
)[i
][m
] += (frameorder
[m
]/size
);
596 { /*Skip following if doing per-molecule*/
597 for (k
= 0; k
< nslices
; k
++) {
598 if (slCount
[k
]) { /* if no elements, nothing has to be added */
599 (*slOrder
)[k
][i
] += slFrameorder
[k
]/slCount
[k
];
600 slFrameorder
[k
] = 0; slCount
[k
] = 0;
603 } /* end loop i, over all groups in indexfile */
607 } while (read_next_x(oenv
,status
,&t
,natoms
,x0
,box
));
608 /*********** done with status file **********/
610 fprintf(stderr
,"\nRead trajectory. Printing parameters to file\n");
612 /* average over frames */
613 for (i
= 1; i
< ngrps
- 1; i
++) {
614 svmul(1.0/nr_frames
, (*order
)[i
], (*order
)[i
]);
615 fprintf(stderr
,"Atom %d Tensor: x=%g , y=%g, z=%g\n",i
,(*order
)[i
][XX
],
616 (*order
)[i
][YY
], (*order
)[i
][ZZ
]);
617 if (bSliced
|| permolecule
) {
618 for (k
= 0; k
< nslices
; k
++)
619 (*slOrder
)[k
][i
] /= nr_frames
;
622 for (k
= 0; k
< nslices
; k
++)
623 (*distvals
)[k
][i
] /= nr_frames
;
627 fprintf(stderr
,"Average angle between double bond and normal: %f\n",
628 180*sdbangle
/(nr_frames
* size
*M_PI
));
630 sfree(x0
); /* free memory used by coordinate arrays */
641 void order_plot(rvec order
[], real
*slOrder
[], const char *afile
, const char *bfile
,
642 const char *cfile
, int ngrps
, int nslices
, real slWidth
, bool bSzonly
,
643 bool permolecule
, real
**distvals
, const output_env_t oenv
)
645 FILE *ord
, *slOrd
; /* xvgr files with order parameters */
646 int atom
, slice
; /* atom corresponding to order para.*/
647 char buf
[256]; /* for xvgr title */
648 real S
; /* order parameter averaged over all atoms */
652 sprintf(buf
,"Scd order parameters");
653 ord
= xvgropen(afile
,buf
,"Atom","S",oenv
);
654 sprintf(buf
, "Orderparameters per atom per slice");
655 slOrd
= xvgropen(bfile
, buf
, "Molecule", "S",oenv
);
656 for (atom
= 1; atom
< ngrps
- 1; atom
++) {
657 fprintf(ord
,"%12d %12g\n", atom
, -1 * (0.6667 * order
[atom
][XX
] +
658 0.333 * order
[atom
][YY
]));
661 for (slice
= 0; slice
< nslices
; slice
++) {
662 fprintf(slOrd
,"%12d\t",slice
);
664 fprintf(slOrd
,"%12g\t", distvals
[slice
][1]); /*use distance value at second carbon*/
665 for (atom
= 1; atom
< ngrps
- 1; atom
++)
666 fprintf(slOrd
,"%12g\t", slOrder
[slice
][atom
]);
672 sprintf(buf
,"Orderparameters Sz per atom");
673 ord
= xvgropen(afile
,buf
,"Atom","S",oenv
);
674 fprintf(stderr
,"ngrps = %d, nslices = %d",ngrps
, nslices
);
676 sprintf(buf
, "Orderparameters per atom per slice");
677 slOrd
= xvgropen(bfile
, buf
, "Slice", "S",oenv
);
679 for (atom
= 1; atom
< ngrps
- 1; atom
++)
680 fprintf(ord
,"%12d %12g\n", atom
, order
[atom
][ZZ
]);
682 for (slice
= 0; slice
< nslices
; slice
++) {
684 for (atom
= 1; atom
< ngrps
- 1; atom
++)
685 S
+= slOrder
[slice
][atom
];
686 fprintf(slOrd
,"%12g %12g\n", slice
*slWidth
, S
/atom
);
690 sprintf(buf
,"Order tensor diagonal elements");
691 ord
= xvgropen(afile
,buf
,"Atom","S",oenv
);
692 sprintf(buf
,"Deuterium order parameters");
693 slOrd
= xvgropen(cfile
,buf
, "Atom", "Scd",oenv
);
695 for (atom
= 1; atom
< ngrps
- 1; atom
++) {
696 fprintf(ord
,"%12d %12g %12g %12g\n", atom
, order
[atom
][XX
],
697 order
[atom
][YY
], order
[atom
][ZZ
]);
698 fprintf(slOrd
,"%12d %12g\n", atom
, -1 * (0.6667 * order
[atom
][XX
] +
699 0.333 * order
[atom
][YY
]));
707 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
)
709 /*function to write order parameters as B factors in PDB file using
710 first frame of trajectory*/
713 t_trxframe fr
, frout
;
717 ngrps
-=2; /*we don't have an order parameter for the first or
718 last atom in each chain*/
720 natoms
=read_first_frame(oenv
,&status
,ftp2fn(efTRX
,nfile
,fnm
),&fr
,
730 init_t_atoms(&useatoms
,nout
,TRUE
);
733 /*initialize PDBinfo*/
734 for (i
=0;i
<useatoms
.nr
;++i
)
736 useatoms
.pdbinfo
[i
].type
=0;
737 useatoms
.pdbinfo
[i
].occup
=0.0;
738 useatoms
.pdbinfo
[i
].bfac
=0.0;
739 useatoms
.pdbinfo
[i
].bAnisotropic
=FALSE
;
742 for (j
=0,ctr
=0;j
<nslices
;j
++)
743 for (i
=0;i
<ngrps
;i
++,ctr
++)
745 /*iterate along each chain*/
746 useatoms
.pdbinfo
[ctr
].bfac
=order
[j
][i
+1];
748 useatoms
.pdbinfo
[ctr
].occup
=distvals
[j
][i
+1];
749 copy_rvec(fr
.x
[a
[index
[i
+1]+j
]],frout
.x
[ctr
]);
750 useatoms
.atomname
[ctr
]=top
->atoms
.atomname
[a
[index
[i
+1]+j
]];
751 useatoms
.atom
[ctr
]=top
->atoms
.atom
[a
[index
[i
+1]+j
]];
752 useatoms
.nres
=max(useatoms
.nres
,useatoms
.atom
[ctr
].resind
+1);
753 useatoms
.resinfo
[useatoms
.atom
[ctr
].resind
]=top
->atoms
.resinfo
[useatoms
.atom
[ctr
].resind
]; /*copy resinfo*/
756 write_sto_conf(opt2fn("-ob",nfile
,fnm
),"Order parameters",&useatoms
,frout
.x
,NULL
,frout
.ePBC
,frout
.box
);
759 free_t_atoms(&useatoms
,FALSE
);
762 int gmx_order(int argc
,char *argv
[])
764 const char *desc
[] = {
765 "Compute the order parameter per atom for carbon tails. For atom i the",
766 "vector i-1, i+1 is used together with an axis. ",
767 "The index file should contain only the groups to be used for calculations,",
768 "with each group of equivalent carbons along the relevant acyl chain in its own",
769 "group. There should not be any generic groups (like System, Protein) in the index",
770 "file to avoid confusing the program (this is not relevant to tetrahedral order",
771 "parameters however, which only work for water anyway).[PAR]",
772 "The program can also give all",
773 "diagonal elements of the order tensor and even calculate the deuterium",
774 "order parameter Scd (default). If the option -szonly is given, only one",
775 "order tensor component (specified by the -d option) is given and the",
776 "order parameter per slice is calculated as well. If -szonly is not",
777 "selected, all diagonal elements and the deuterium order parameter is",
779 "The tetrahedrality order parameters can be determined",
780 "around an atom. Both angle an distance order parameters are calculated. See",
781 "P.-L. Chau and A.J. Hardwick, Mol. Phys., 93, (1998), 511-518.",
782 "for more details.[BR]",
786 static int nslices
= 1; /* nr of slices defined */
787 static bool bSzonly
= FALSE
; /* True if only Sz is wanted */
788 static bool bUnsat
= FALSE
; /* True if carbons are unsat. */
789 static const char *normal_axis
[] = { NULL
, "z", "x", "y", NULL
};
790 static bool permolecule
= FALSE
; /*compute on a per-molecule basis */
791 static bool radial
= FALSE
; /*compute a radial membrane normal */
792 static bool distcalc
= FALSE
; /*calculate distance from a reference group */
794 { "-d", FALSE
, etENUM
, {normal_axis
},
795 "Direction of the normal on the membrane" },
796 { "-sl", FALSE
, etINT
, {&nslices
},
797 "Calculate order parameter as function of boxlength, dividing the box"
799 { "-szonly", FALSE
, etBOOL
,{&bSzonly
},
800 "Only give Sz element of order tensor. (axis can be specified with -d)" },
801 { "-unsat", FALSE
, etBOOL
,{&bUnsat
},
802 "Calculate order parameters for unsaturated carbons. Note that this can"
803 "not be mixed with normal order parameters." },
804 { "-permolecule", FALSE
, etBOOL
,{&permolecule
},
805 "Compute per-molecule Scd order parameters" },
806 { "-radial", FALSE
, etBOOL
,{&radial
},
807 "Compute a radial membrane normal" },
808 { "-calcdist", FALSE
, etBOOL
,{&distcalc
},
809 "Compute distance from a reference (currently defined only for radial and permolecule)" },
812 rvec
*order
; /* order par. for each atom */
813 real
**slOrder
; /* same, per slice */
814 real slWidth
= 0.0; /* width of a slice */
815 char **grpname
; /* groupnames */
816 int ngrps
, /* nr. of groups */
818 axis
=0; /* normal axis */
819 t_topology
*top
; /* topology */
821 atom_id
*index
, /* indices for a */
822 *a
; /* atom numbers in each group */
823 t_blocka
*block
; /* data from index file */
824 t_filenm fnm
[] = { /* files for g_order */
825 { efTRX
, "-f", NULL
, ffREAD
}, /* trajectory file */
826 { efNDX
, "-n", NULL
, ffREAD
}, /* index file */
827 { efNDX
, "-nr", NULL
, ffREAD
}, /* index for radial axis calculation */
828 { efTPX
, NULL
, NULL
, ffREAD
}, /* topology file */
829 { efXVG
,"-o","order", ffWRITE
}, /* xvgr output file */
830 { efXVG
,"-od","deuter", ffWRITE
}, /* xvgr output file */
831 { efPDB
,"-ob",NULL
, ffWRITE
}, /* write Scd as B factors to PDB if permolecule */
832 { efXVG
,"-os","sliced", ffWRITE
}, /* xvgr output file */
833 { efXVG
,"-Sg","sg-ang", ffOPTWR
}, /* xvgr output file */
834 { efXVG
,"-Sk","sk-dist", ffOPTWR
}, /* xvgr output file */
835 { efXVG
,"-Sgsl","sg-ang-slice", ffOPTWR
}, /* xvgr output file */
836 { efXVG
,"-Sksl","sk-dist-slice", ffOPTWR
}, /* xvgr output file */
838 bool bSliced
= FALSE
; /* True if box is sliced */
839 #define NFILE asize(fnm)
840 real
**distvals
=NULL
;
841 const char *sgfnm
,*skfnm
,*ndxfnm
,*tpsfnm
,*trxfnm
;
844 CopyRight(stderr
,argv
[0]);
846 parse_common_args(&argc
,argv
,PCA_CAN_VIEW
| PCA_CAN_TIME
| PCA_BE_NICE
,
847 NFILE
,fnm
,asize(pa
),pa
,asize(desc
),desc
,0, NULL
,&oenv
);
849 gmx_fatal(FARGS
,"Can not have nslices < 1");
850 sgfnm
= opt2fn_null("-Sg",NFILE
,fnm
);
851 skfnm
= opt2fn_null("-Sk",NFILE
,fnm
);
852 ndxfnm
= opt2fn_null("-n",NFILE
,fnm
);
853 tpsfnm
= ftp2fn(efTPX
,NFILE
,fnm
);
854 trxfnm
= ftp2fn(efTRX
,NFILE
,fnm
);
857 if (strcmp(normal_axis
[0],"x") == 0) axis
= XX
;
858 else if (strcmp(normal_axis
[0],"y") == 0) axis
= YY
;
859 else if (strcmp(normal_axis
[0],"z") == 0) axis
= ZZ
;
860 else gmx_fatal(FARGS
,"Invalid axis, use x, y or z");
864 fprintf(stderr
,"Taking x axis as normal to the membrane\n");
867 fprintf(stderr
,"Taking y axis as normal to the membrane\n");
870 fprintf(stderr
,"Taking z axis as normal to the membrane\n");
874 /* tetraheder order parameter */
875 if (skfnm
|| sgfnm
) {
876 /* If either of theoptions is set we compute both */
877 sgfnm
= opt2fn("-Sg",NFILE
,fnm
);
878 skfnm
= opt2fn("-Sk",NFILE
,fnm
);
879 calc_tetra_order_parm(ndxfnm
,tpsfnm
,trxfnm
,sgfnm
,skfnm
,nslices
,axis
,
880 opt2fn("-Sgsl",NFILE
,fnm
),opt2fn("-Sksl",NFILE
,fnm
),
882 /* view xvgr files */
883 do_view(oenv
,opt2fn("-Sg",NFILE
,fnm
), NULL
);
884 do_view(oenv
,opt2fn("-Sk",NFILE
,fnm
), NULL
);
886 do_view(oenv
,opt2fn("-Sgsl",NFILE
,fnm
), NULL
);
887 do_view(oenv
,opt2fn("-Sksl",NFILE
,fnm
), NULL
);
891 /* tail order parameter */
895 fprintf(stderr
,"Dividing box in %d slices.\n\n", nslices
);
899 fprintf(stderr
,"Only calculating Sz\n");
901 fprintf(stderr
,"Taking carbons as unsaturated!\n");
903 top
= read_top(ftp2fn(efTPX
,NFILE
,fnm
),&ePBC
); /* read topology file */
905 block
= init_index(ftp2fn(efNDX
,NFILE
,fnm
),&grpname
);
906 index
= block
->index
; /* get indices from t_block block */
907 a
= block
->a
; /* see block.h */
912 nslices
= index
[1] - index
[0]; /*I think this assumes contiguous lipids in topology*/
913 fprintf(stderr
,"Calculating Scd order parameters for each of %d molecules\n",nslices
);
919 fprintf(stderr
,"Calculating radial distances\n");
921 gmx_fatal(FARGS
,"Cannot yet output radial distances without permolecule\n");
924 /* show atomtypes, to check if index file is correct */
925 print_types(index
, a
, ngrps
, grpname
, top
);
927 calc_order(ftp2fn(efTRX
,NFILE
,fnm
), index
, a
, &order
,
928 &slOrder
, &slWidth
, nslices
, bSliced
, bUnsat
,
929 top
, ePBC
, ngrps
, axis
,permolecule
,radial
,distcalc
,opt2fn_null("-nr",NFILE
,fnm
),&distvals
, oenv
);
932 ngrps
--; /*don't print the last group--was used for
933 center-of-mass determination*/
935 order_plot(order
, slOrder
, opt2fn("-o",NFILE
,fnm
), opt2fn("-os",NFILE
,fnm
),
936 opt2fn("-od",NFILE
,fnm
), ngrps
, nslices
, slWidth
, bSzonly
,permolecule
,distvals
,oenv
);
938 if (opt2bSet("-ob",NFILE
,fnm
))
942 "Won't write B-factors with averaged order parameters; use -permolecule\n");
944 write_bfactors(fnm
,NFILE
,index
,a
,nslices
,ngrps
,slOrder
,top
,distvals
,oenv
);
948 do_view(oenv
,opt2fn("-o",NFILE
,fnm
), NULL
); /* view xvgr file */
949 do_view(oenv
,opt2fn("-os",NFILE
,fnm
), NULL
); /* view xvgr file */
950 do_view(oenv
,opt2fn("-od",NFILE
,fnm
), NULL
); /* view xvgr file */
956 for (i
=0;i
<nslices
;++i
)