gmx_order_4.5.4_modifiedbyChrisNeale.c
1 
/*


2 
*

3 
* This source code is part of

4 
*

5 
* G R O M A C S

6 
*

7 
* GROningen MAchine for Chemical Simulations

8 
*

9 
* VERSION 3.2.0

10 
* Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.

11 
* Copyright (c) 19912000, University of Groningen, The Netherlands.

12 
* Copyright (c) 20012004, The GROMACS development team,

13 
* check out http://www.gromacs.org for more information.

14 

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.

19 
*

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.

26 
*

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.

29 
*

30 
* For more info, check our website at http://www.gromacs.org

31 
*

32 
* And Hey:

33 
* Green Red Orange Magenta Azure Cyan Skyblue

34 
*/

35 
#ifdef HAVE_CONFIG_H

36 
#include <config.h> 
37 
#endif

38  
39 
#include <math.h> 
40 
#include <ctype.h> 
41  
42 
#include "sysstuff.h" 
43 
#include "string.h" 
44 
#include "typedefs.h" 
45 
#include "smalloc.h" 
46 
#include "macros.h" 
47 
#include "gstat.h" 
48 
#include "vec.h" 
49 
#include "xvgr.h" 
50 
#include "pbc.h" 
51 
#include "copyrite.h" 
52 
#include "futil.h" 
53 
#include "statutil.h" 
54 
#include "index.h" 
55 
#include "tpxio.h" 
56 
#include "confio.h" 
57 
#include "cmat.h" 
58 
#include "gmx_ana.h" 
59  
60  
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 boxaxis */

66 
/* is calculated. In that case i = j = axis, and delta(ij) = 1. */

67 
/* */

68 
/* Peter Tieleman, April 1995 */

69 
/* P.J. van Maaren, November 2005 Added tetrahedral stuff */

70 
/****************************************************************************/

71  
72 
static void find_nearest_neighbours(t_topology top, int ePBC, 
73 
int natoms, matrix box,

74 
rvec x[],int maxidx,atom_id index[],

75 
real time, 
76 
real *sgmean, real *skmean, 
77 
int nslice,int slice_dim, 
78 
real sgslice[],real skslice[], 
79 
gmx_rmpbc_t gpbc) 
80 
{ 
81 
FILE *fpoutdist; 
82 
char fnsgdist[32]; 
83 
int ix,jx,nsgbin, *sgbin;

84 
int i1,i2,i,ibin,j,k,l,n,*nn[4]; 
85 
rvec dx,dx1,dx2,rj,rk,urk,urj; 
86 
real cost,cost2,*sgmol,*skmol,rmean,rmean2,r2,box2,*r_nn[4];

87 
t_pbc pbc; 
88 
t_mat *dmat; 
89 
t_dist *d; 
90 
int m1,mm,sl_index;

91 
int **nnb,*sl_count;

92 
real onethird=1.0/3.0; 
93 
/* dmat = init_mat(maxidx, FALSE); */

94 
box2 = box[XX][XX] * box[XX][XX]; 
95 
snew(sl_count,nslice); 
96 
for (i=0; (i<4); i++) { 
97 
snew(r_nn[i],natoms); 
98 
snew(nn[i],natoms); 
99 

100 
for (j=0; (j<natoms); j++) { 
101 
r_nn[i][j] = box2; 
102 
} 
103 
} 
104 

105 
snew(sgmol,maxidx); 
106 
snew(skmol,maxidx); 
107  
108 
/* Must init pbc every step because of pressure coupling */

109 
set_pbc(&pbc,ePBC,box); 
110 

111 
gmx_rmpbc(gpbc,natoms,box,x); 
112  
113 
nsgbin = 1 + 1/0.0005; 
114 
snew(sgbin,nsgbin); 
115  
116 
*sgmean = 0.0; 
117 
*skmean = 0.0; 
118 
l=0;

119 
for (i=0; (i<maxidx); i++) { /* loop over index file */ 
120 
ix = index[i]; 
121 
for (j=0; (j<maxidx); j++) { 
122 
if (i == j) continue; 
123  
124 
jx = index[j]; 
125 

126 
pbc_dx(&pbc,x[ix],x[jx],dx); 
127 
r2=iprod(dx,dx); 
128  
129 
/* set_mat_entry(dmat,i,j,r2); */

130  
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; 
146 
} 
147 
} 
148  
149  
150 
/* calculate mean distance between nearest neighbours */

151 
rmean = 0;

152 
for (j=0; (j<4); j++) { 
153 
r_nn[j][i] = sqrt(r_nn[j][i]); 
154 
rmean += r_nn[j][i]; 
155 
} 
156 
rmean /= 4;

157 

158 
n = 0;

159 
sgmol[i] = 0.0; 
160 
skmol[i] = 0.0; 
161  
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); 
168  
169 
unitv(rk,urk); 
170 
unitv(rj,urj); 
171 

172 
cost = iprod(urk,urj) + onethird; 
173 
cost2 = cost * cost; 
174  
175 
/* sgmol[i] += 3*cost2/32; */

176 
sgmol[i] += cost2; 
177  
178 
/* determine distribution */

179 
ibin = nsgbin * cost2; 
180 
if (ibin < nsgbin)

181 
sgbin[ibin]++; 
182 
/* printf("%d %d %f %d %d\n", j, k, cost * cost, ibin, sgbin[ibin]);*/

183 
l++; 
184 
n++; 
185 
} 
186 
} 
187  
188 
/* normalize sgmol between 0.0 and 1.0 */

189 
sgmol[i] = 3*sgmol[i]/32; 
190 
*sgmean += sgmol[i]; 
191  
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]) );

198 
*/

199 
} 
200 

201 
*skmean += skmol[i]; 
202 

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 */

209 

210 
*sgmean /= maxidx; 
211 
*skmean /= maxidx; 
212 

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]; 
217 
} 
218 
} 
219 
sfree(sl_count); 
220 
sfree(sgbin); 
221 
sfree(sgmol); 
222 
sfree(skmol); 
223 
for (i=0; (i<4); i++) { 
224 
sfree(r_nn[i]); 
225 
sfree(nn[i]); 
226 
} 
227 
} 
228  
229  
230 
static void calc_tetra_order_parm(const char *fnNDX,const char *fnTPS, 
231 
const char *fnTRX, const char *sgfn, 
232 
const char *skfn, 
233 
int nslice,int slice_dim, 
234 
const char *sgslfn,const char *skslfn, 
235 
const output_env_t oenv)

236 
{ 
237 
FILE *fpsg=NULL,*fpsk=NULL; 
238 
t_topology top; 
239 
int ePBC;

240 
char title[STRLEN],fn[STRLEN],subtitle[STRLEN];

241 
t_trxstatus *status; 
242 
int natoms;

243 
real t; 
244 
rvec *xtop,*x; 
245 
matrix box; 
246 
real sg,sk; 
247 
atom_id **index; 
248 
char **grpname;

249 
int i,*isize,ng,nframes;

250 
real *sg_slice,*sg_slice_tot,*sk_slice,*sk_slice_tot; 
251 
gmx_rmpbc_t gpbc=NULL;

252  
253  
254 
read_tps_conf(fnTPS,title,&top,&ePBC,&xtop,NULL,box,FALSE);

255  
256 
snew(sg_slice,nslice); 
257 
snew(sk_slice,nslice); 
258 
snew(sg_slice_tot,nslice); 
259 
snew(sk_slice_tot,nslice); 
260 
ng = 1;

261 
/* get index groups */

262 
printf("Select the group that contains the atoms you want to use for the tetrahedrality order parameter calculation:\n");

263 
snew(grpname,ng); 
264 
snew(index,ng); 
265 
snew(isize,ng); 
266 
get_index(&top.atoms,fnNDX,ng,isize,index,grpname); 
267  
268 
/* Analyze trajectory */

269 
natoms=read_first_x(oenv,&status,fnTRX,&t,&x,box); 
270 
if ( natoms > top.atoms.nr )

271 
gmx_fatal(FARGS,"Topology (%d atoms) does not match trajectory (%d atoms)",

272 
top.atoms.nr,natoms); 
273 
check_index(NULL,ng,index[0],NULL,natoms); 
274  
275 
fpsg=xvgropen(sgfn,"S\\sg\\N Angle Order Parameter","Time (ps)","S\\sg\\N", 
276 
oenv); 
277 
fpsk=xvgropen(skfn,"S\\sk\\N Distance Order Parameter","Time (ps)","S\\sk\\N", 
278 
oenv); 
279  
280 
/* loop over frames */

281 
gpbc = gmx_rmpbc_init(&top.idef,ePBC,natoms,box); 
282 
nframes = 0;

283 
do {

284 
find_nearest_neighbours(top,ePBC,natoms,box,x,isize[0],index[0],t, 
285 
&sg,&sk,nslice,slice_dim,sg_slice,sk_slice,gpbc); 
286 
for(i=0; (i<nslice); i++) { 
287 
sg_slice_tot[i] += sg_slice[i]; 
288 
sk_slice_tot[i] += sk_slice[i]; 
289 
} 
290 
fprintf(fpsg,"%f %f\n", t, sg);

291 
fprintf(fpsk,"%f %f\n", t, sk);

292 
nframes++; 
293 
} while (read_next_x(oenv,status,&t,natoms,x,box));

294 
close_trj(status); 
295 
gmx_rmpbc_done(gpbc); 
296  
297 
sfree(grpname); 
298 
sfree(index); 
299 
sfree(isize); 
300  
301 
ffclose(fpsg); 
302 
ffclose(fpsk); 
303 

304 
fpsg = xvgropen(sgslfn, 
305 
"S\\sg\\N Angle Order Parameter / Slab","(nm)","S\\sg\\N", 
306 
oenv); 
307 
fpsk = xvgropen(skslfn, 
308 
"S\\sk\\N Distance Order Parameter / Slab","(nm)","S\\sk\\N", 
309 
oenv); 
310 
for(i=0; (i<nslice); i++) { 
311 
fprintf(fpsg,"%10g %10g\n",(i+0.5)*box[slice_dim][slice_dim]/nslice, 
312 
sg_slice_tot[i]/nframes); 
313 
fprintf(fpsk,"%10g %10g\n",(i+0.5)*box[slice_dim][slice_dim]/nslice, 
314 
sk_slice_tot[i]/nframes); 
315 
} 
316 
ffclose(fpsg); 
317 
ffclose(fpsk); 
318 
} 
319  
320  
321 
/* Print name of first atom in all groups in index file */

322 
static void print_types(atom_id index[], atom_id a[], int ngrps, 
323 
char *groups[], t_topology *top)

324 
{ 
325 
int i;

326  
327 
fprintf(stderr,"Using following groups: \n");

328 
for(i = 0; i < ngrps; i++) 
329 
fprintf(stderr,"Groupname: %s First atomname: %s First atomnr %u\n",

330 
groups[i], *(top>atoms.atomname[a[index[i]]]), a[index[i]]); 
331 
fprintf(stderr,"\n");

332 
} 
333  
334 
static void check_length(real length, int a, int b) 
335 
{ 
336 
if (length > 0.3) 
337 
fprintf(stderr,"WARNING: distance between atoms %d and "

338 
"%d > 0.3 nm (%f). Index file might be corrupt.\n",

339 
a, b, length); 
340 
} 
341  
342 
void calc_order(const char *fn, atom_id *index, atom_id *a, rvec **order, 
343 
real ***slOrder, real *slWidth, int nslices, gmx_bool bSliced,

344 
gmx_bool bUnsat, t_topology *top, int ePBC, int ngrps, int axis, 
345 
gmx_bool permolecule, gmx_bool radial, gmx_bool distcalc, const char *radfn, 
346 
real ***distvals, 
347 
const output_env_t oenv)

348 
{ 
349 
/* if permolecule = TRUE, order parameters will be calculed per molecule

350 
* and stored in slOrder with #slices = # molecules */

351 
rvec *x0, /* coordinates with pbc */

352 
*x1, /* coordinates without pbc */

353 
dist; /* vector between two atoms */

354 
matrix box; /* box (3x3) */

355 
t_trxstatus *status; 
356 
rvec cossum, /* sum of vector angles for three axes */

357 
Sx, Sy, Sz, /* the three molecular axes */

358 
tmp1, tmp2, /* temp. rvecs for calculating dot products */

359 
frameorder; /* order parameters for one frame */

360 
real *slFrameorder; /* order parameter for one frame, per slice */

361 
real length, /* total distance between two atoms */

362 
t, /* time from trajectory */

363 
z_ave,z1,z2; /* average z, used to det. which slice atom is in */

364 
int natoms, /* nr. atoms in trj */ 
365 
nr_tails, /* nr tails, to check if index file is correct */

366 
size=0, /* nr. of atoms in group. same as nr_tails */ 
367 
i,j,m,k,l,teller = 0,

368 
slice, /* current slice number */

369 
nr_frames = 0,

370 
*slCount; /* nr. of atoms in one slice */

371 
real dbangle = 0, /* angle between double bond and axis */ 
372 
sdbangle = 0;/* sum of these angles */ 
373 
gmx_bool use_unitvector = FALSE; /* use a specified unit vector instead of axis to specify unit normal*/

374 
rvec direction, com,dref,dvec; 
375 
int comsize, distsize;

376 
atom_id *comidx=NULL, *distidx=NULL; 
377 
char *grpname=NULL; 
378 
t_pbc pbc; 
379 
real arcdist; 
380 
gmx_rmpbc_t gpbc=NULL;

381  
382 
// BEGIN CN MOD 

383 
// variables for modification by CN

384 
rvec ab,ac,bc,e,ce,eg,ef,g,f,cg,cf; 
385 
float dab,def,deg,gdot,fdot,edot;

386 
// END CN MOD 

387  
388  
389 
/* PBC added for centerofmass vector*/

390 
/* Initiate the pbc structure */

391 
memset(&pbc,0,sizeof(pbc)); 
392  
393 
if ((natoms = read_first_x(oenv,&status,fn,&t,&x0,box)) == 0) 
394 
gmx_fatal(FARGS,"Could not read coordinates from statusfile\n");

395  
396 
nr_tails = index[1]  index[0]; 
397 
fprintf(stderr,"Number of elements in first group: %d\n",nr_tails);

398 
/* take first group as standard. Not rocksolid, but might catch error in index*/

399  
400 
if (permolecule)

401 
{ 
402 
nslices=nr_tails; 
403 
bSliced=FALSE; /*force slices off */

404 
fprintf(stderr,"Calculating order parameters for each of %d molecules\n",

405 
nslices); 
406 
} 
407 

408 
if (radial)

409 
{ 
410 
use_unitvector=TRUE; 
411 
fprintf(stderr,"Select an index group to calculate the radial membrane normal\n");

412 
get_index(&top>atoms,radfn,1,&comsize,&comidx,&grpname);

413 
if (distcalc)

414 
{ 
415 
if (grpname!=NULL) 
416 
sfree(grpname); 
417 
fprintf(stderr,"Select an index group to use as distance reference\n");

418 
get_index(&top>atoms,radfn,1,&distsize,&distidx,&grpname);

419 
bSliced=FALSE; /*force slices off*/

420 
} 
421 
} 
422  
423 
if (use_unitvector && bSliced)

424 
fprintf(stderr,"Warning: slicing and specified unit vectors are not currently compatible\n");

425  
426 
snew(slCount, nslices); 
427 
snew(*slOrder, nslices); 
428 
for(i = 0; i < nslices; i++) 
429 
snew((*slOrder)[i],ngrps); 
430 
if (distcalc)

431 
{ 
432 
snew(*distvals, nslices); 
433 
for(i = 0; i < nslices; i++) 
434 
snew((*distvals)[i],ngrps); 
435 
} 
436 
snew(*order,ngrps); 
437 
snew(slFrameorder, nslices); 
438 
snew(x1, natoms); 
439 

440 
if (bSliced) {

441 
*slWidth = box[axis][axis]/nslices; 
442 
fprintf(stderr,"Box divided in %d slices. Initial width of slice: %f\n",

443 
nslices, *slWidth); 
444 
} 
445  
446 
#if 0

447 
nr_tails = index[1]  index[0];

448 
fprintf(stderr,"Number of elements in first group: %d\n",nr_tails);

449 
/* take first group as standard. Not rocksolid, but might catch error

450 
in index*/

451 
#endif

452  
453 
teller = 0;

454  
455 
gpbc = gmx_rmpbc_init(&top>idef,ePBC,natoms,box); 
456 
/*********** Start processing trajectory ***********/

457 
do {

458 
if (bSliced)

459 
*slWidth = box[axis][axis]/nslices; 
460 
teller++; 
461 

462 
set_pbc(&pbc,ePBC,box); 
463 
gmx_rmpbc_copy(gpbc,natoms,box,x0,x1); 
464  
465 
/* Now loop over all groups. There are ngrps groups, the order parameter can

466 
be calculated for grp 1 to grp ngrps  1. For each group, loop over all

467 
atoms in group, which is index[i] to (index[i+1]  1) See block.h. Of

468 
course, in this case index[i+1] index[i] has to be the same for all

469 
groups, namely the number of tails. i just runs over all atoms in a tail,

470 
so for DPPC ngrps = 16 and i runs from 1 to 14, including 14

471 
*/

472 

473 
if (radial)

474 
{ 
475 
/*centerofmass determination*/

476 
com[XX]=0.0; com[YY]=0.0; com[ZZ]=0.0; 
477 
for (j=0;j<comsize;j++) 
478 
rvec_inc(com,x1[comidx[j]]); 
479 
svmul(1.0/comsize,com,com); 
480 
} 
481 
if (distcalc)

482 
{ 
483 
dref[XX]=0.0; dref[YY]=0.0;dref[ZZ]=0.0; 
484 
for (j=0;j<distsize;j++) 
485 
rvec_inc(dist,x1[distidx[j]]); 
486 
svmul(1.0/distsize,dref,dref); 
487 
pbc_dx(&pbc,dref,com,dvec); 
488 
unitv(dvec,dvec); 
489 
} 
490 

491 
for (i = 1; i < ngrps  1; i++) { 
492 
clear_rvec(frameorder); 
493 

494 
size = index[i+1]  index[i];

495 
if (size != nr_tails)

496 
gmx_fatal(FARGS,"grp %d does not have same number of"

497 
" elements as grp 1\n",i);

498 

499 
for (j = 0; j < size; j++) { 
500 
if (radial)

501 
/*create unit vector*/

502 
{ 
503 
pbc_dx(&pbc,x1[a[index[i]+j]],com,direction); 
504 
unitv(direction,direction); 
505 
/*DEBUG*/

506 
/*if (j==0)

507 
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],

508 
direction[0],direction[1],direction[2]);*/

509 
} 
510  
511 
// BEGIN CN MOD 

512 
if (bUnsat) {

513 
//exactly as for the case without unsat, but here we select the ce vector to represent the hydrogen

514 
// place the 2 imaginary hydrogen atoms by enscribing a tetrahedron in a cube

515 
// a,b,c are the 3 backbone atoms, where a and b are in two cube corners and c is at the center

516 
rvec_sub(x1[a[index[i1]+j]], x1[a[index[i+1]+j]], ab); 
517 
rvec_sub(x1[a[index[i]+j]], x1[a[index[i+1]+j]], ac);

518 
rvec_sub(x1[a[index[i]+j]], x1[a[index[i1]+j]], bc);

519 
// e is a point in the middle of the face opposite from the ab face

520 
for(m=0;m<DIM;m++){ 
521 
ce[m]=(ac[m]+bc[m])/2;

522 
} 
523 
rvec_add(x1[a[index[i]+j]],ce,e); 
524 
//UNCOMMENT THE NEXT LINE TO PRINT OUT THE POSITION OF THE PSEUDOATOM FOR TESTING <<<<<<<<<< LOOK LOOK LOOK LOOK

525 
//fprintf(stdout,"ATOM E %f %f %f\n",e[0],e[1],e[2]);

526 
//now take the angle between ce and z and this gets plugged into S_CD=1/2<3cos^2(ang)1>

527 
edot=ce[2]/norm(ce);

528  
529 
// using frameorder[ZZ] to store the result

530 
frameorder[ZZ] += 0.5 * (3 * edot*edot  1); 
531  
532 
} else {

533 
// place the 2 imaginary hydrogen atoms by enscribing a tetrahedron in a cube

534 
// a,b,c are the 3 backbone atoms, where a and b are in two cube corners and c is at the center

535 
rvec_sub(x1[a[index[i1]+j]], x1[a[index[i+1]+j]], ab); 
536 
rvec_sub(x1[a[index[i]+j]], x1[a[index[i+1]+j]], ac);

537 
rvec_sub(x1[a[index[i]+j]], x1[a[index[i1]+j]], bc);

538 
// e is a point in the middle of the face opposite from the ab face

539 
for(m=0;m<DIM;m++){ 
540 
ce[m]=(ac[m]+bc[m])/2;

541 
} 
542 
rvec_add(x1[a[index[i]+j]],ce,e); 
543 
dab=norm(ab); 
544 
//the eg vector is the cross product of the ac and ce vectors (and ef is cross prod of bc and ce)

545 
cprod(ac,ce,eg); 
546 
cprod(bc,ce,ef); 
547 
deg=norm(eg); 
548 
def=norm(ef); 
549 
//progress along eg by the correct distance from e to find g (and along ef from e to find f)

550 
for(m=0;m<DIM;m++){ 
551 
g[m]=e[m]eg[m]*(dab/deg)/2;

552 
f[m]=e[m]ef[m]*(dab/def)/2;

553 
} 
554 
//now find the cg and cf vectors, which are the CD vectors for S_CD

555 
rvec_sub(g,x1[a[index[i]+j]],cg); 
556 
rvec_sub(f,x1[a[index[i]+j]],cf); 
557 
//UNCOMMENT THE NEXT TWO LINES TO PRINT OUT THE POSITION OF THE PSEUDOATOM FOR TESTING <<<<<<<<<< LOOK LOOK LOOK LOOK

558 
//fprintf(stdout,"ATOM G %f %f %f\n",g[0],g[1],g[2]);

559 
//fprintf(stdout,"ATOM F %f %f %f\n",f[0],f[1],f[2]);

560  
561 
//now take the angle between cg and z  and also for cf and z  and this gets plugged into S_CD=1/2<3cos^2(ang)1>

562 
gdot=cg[2]/norm(cg);

563 
fdot=cf[2]/norm(cf);

564  
565 
// using frameorder[ZZ] to store the result  mult by 0.25 instead of 0.5 because we are adding two values

566 
frameorder[ZZ] += 0.25 * (3 * gdot*gdot  1); 
567 
frameorder[ZZ] += 0.25 * (3 * fdot*fdot  1); 
568 
} 
569  
570 
// END CN MOD 

571 
} /* end loop j, over all atoms in group */

572 

573 
for (m = 0; m < DIM; m++) 
574 
(*order)[i][m] += (frameorder[m]/size); 
575 

576 
if (!permolecule)

577 
{ /*Skip following if doing permolecule*/

578 
for (k = 0; k < nslices; k++) { 
579 
if (slCount[k]) { /* if no elements, nothing has to be added */ 
580 
(*slOrder)[k][i] += slFrameorder[k]/slCount[k]; 
581 
slFrameorder[k] = 0; slCount[k] = 0; 
582 
} 
583 
} 
584 
} /* end loop i, over all groups in indexfile */

585 
} 
586 
nr_frames++; 
587 

588 
} while (read_next_x(oenv,status,&t,natoms,x0,box));

589 
/*********** done with status file **********/

590 

591 
fprintf(stderr,"\nRead trajectory. Printing parameters to file\n");

592 
gmx_rmpbc_done(gpbc); 
593  
594 
/* average over frames */

595 
for (i = 1; i < ngrps  1; i++) { 
596 
svmul(1.0/nr_frames, (*order)[i], (*order)[i]); 
597 
fprintf(stderr,"Atom %d Tensor: x=%g , y=%g, z=%g\n",i,(*order)[i][XX],

598 
(*order)[i][YY], (*order)[i][ZZ]); 
599 
if (bSliced  permolecule) {

600 
for (k = 0; k < nslices; k++) 
601 
(*slOrder)[k][i] /= nr_frames; 
602 
} 
603 
if (distcalc)

604 
for (k = 0; k < nslices; k++) 
605 
(*distvals)[k][i] /= nr_frames; 
606 
} 
607  
608 
if (bUnsat)

609 
fprintf(stderr,"Average angle between double bond and normal: %f\n",

610 
180*sdbangle/(nr_frames * size*M_PI));

611  
612 
sfree(x0); /* free memory used by coordinate arrays */

613 
sfree(x1); 
614 
if (comidx!=NULL) 
615 
sfree(comidx); 
616 
if (distidx!=NULL) 
617 
sfree(distidx); 
618 
if (grpname!=NULL) 
619 
sfree(grpname); 
620 
} 
621  
622  
623 
void order_plot(rvec order[], real *slOrder[], const char *afile, const char *bfile, 
624 
const char *cfile, int ngrps, int nslices, real slWidth, gmx_bool bSzonly, 
625 
gmx_bool permolecule, real **distvals, const output_env_t oenv)

626 
{ 
627 
FILE *ord, *slOrd; /* xvgr files with order parameters */

628 
int atom, slice; /* atom corresponding to order para.*/ 
629 
char buf[256]; /* for xvgr title */ 
630 
real S; /* order parameter averaged over all atoms */

631  
632 
if (permolecule)

633 
{ 
634 
sprintf(buf,"Scd order parameters");

635 
ord = xvgropen(afile,buf,"Atom","S",oenv); 
636 
sprintf(buf, "Orderparameters per atom per slice");

637 
slOrd = xvgropen(bfile, buf, "Molecule", "S",oenv); 
638 
for (atom = 1; atom < ngrps  1; atom++) { 
639 
fprintf(ord,"%12d %12g\n", atom, 1 * (0.6667 * order[atom][XX] + 
640 
0.333 * order[atom][YY])); 
641 
} 
642  
643 
for (slice = 0; slice < nslices; slice++) { 
644 
fprintf(slOrd,"%12d\t",slice);

645 
if (distvals)

646 
fprintf(slOrd,"%12g\t", distvals[slice][1]); /*use distance value at second carbon*/ 
647 
for (atom = 1; atom < ngrps  1; atom++) 
648 
fprintf(slOrd,"%12g\t", slOrder[slice][atom]);

649 
fprintf(slOrd,"\n");

650 
} 
651  
652 
} 
653 
else if (bSzonly) { 
654 
sprintf(buf,"Orderparameters Sz per atom");

655 
ord = xvgropen(afile,buf,"Atom","S",oenv); 
656 
fprintf(stderr,"ngrps = %d, nslices = %d",ngrps, nslices);

657  
658 
sprintf(buf, "Orderparameters per atom per slice");

659 
slOrd = xvgropen(bfile, buf, "Slice", "S",oenv); 
660 

661 
for (atom = 1; atom < ngrps  1; atom++) 
662 
fprintf(ord,"%12d %12g\n", atom, order[atom][ZZ]);

663  
664 
for (slice = 0; slice < nslices; slice++) { 
665 
S = 0;

666 
for (atom = 1; atom < ngrps  1; atom++) 
667 
S += slOrder[slice][atom]; 
668 
fprintf(slOrd,"%12g %12g\n", slice*slWidth, S/atom);

669 
} 
670  
671 
} else {

672 
sprintf(buf,"Order tensor diagonal elements");

673 
ord = xvgropen(afile,buf,"Atom","S",oenv); 
674 
sprintf(buf,"Deuterium order parameters");

675 
slOrd = xvgropen(cfile,buf, "Atom", "Scd",oenv); 
676  
677 
for (atom = 1; atom < ngrps  1; atom++) { 
678 
fprintf(ord,"%12d %12g %12g %12g\n", atom, order[atom][XX],

679 
order[atom][YY], order[atom][ZZ]); 
680 
// BEGIN CN MOD 

681 
// will place the order parameter value directly in the ZZ index.

682 
// not sure why I need to take the negative value....

683 
fprintf(slOrd,"%12d %12g\n", atom, 1 * order[atom][ZZ]); 
684 
// END CN MOD 

685 
} 
686 

687 
ffclose(ord); 
688 
ffclose(slOrd); 
689 
} 
690 
} 
691  
692 
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) 
693 
{ 
694 
/*function to write order parameters as B factors in PDB file using

695 
first frame of trajectory*/

696 
t_trxstatus *status; 
697 
int natoms;

698 
t_trxframe fr, frout; 
699 
t_atoms useatoms; 
700 
int i,j,ctr,nout;

701  
702 
ngrps=2; /*we don't have an order parameter for the first or 
703 
last atom in each chain*/

704 
nout=nslices*ngrps; 
705 
natoms=read_first_frame(oenv,&status,ftp2fn(efTRX,nfile,fnm),&fr, 
706 
TRX_NEED_X); 
707 
close_trj(status); 
708 
frout = fr; 
709 
frout.natoms=nout; 
710 
frout.bF=FALSE; 
711 
frout.bV=FALSE; 
712 
frout.x=0;

713 
snew(frout.x,nout); 
714 

715 
init_t_atoms(&useatoms,nout,TRUE); 
716 
useatoms.nr=nout; 
717  
718 
/*initialize PDBinfo*/

719 
for (i=0;i<useatoms.nr;++i) 
720 
{ 
721 
useatoms.pdbinfo[i].type=0;

722 
useatoms.pdbinfo[i].occup=0.0; 
723 
useatoms.pdbinfo[i].bfac=0.0; 
724 
useatoms.pdbinfo[i].bAnisotropic=FALSE; 
725 
} 
726  
727 
for (j=0,ctr=0;j<nslices;j++) 
728 
for (i=0;i<ngrps;i++,ctr++) 
729 
{ 
730 
/*iterate along each chain*/

731 
useatoms.pdbinfo[ctr].bfac=order[j][i+1];

732 
if (distvals)

733 
useatoms.pdbinfo[ctr].occup=distvals[j][i+1];

734 
copy_rvec(fr.x[a[index[i+1]+j]],frout.x[ctr]);

735 
useatoms.atomname[ctr]=top>atoms.atomname[a[index[i+1]+j]];

736 
useatoms.atom[ctr]=top>atoms.atom[a[index[i+1]+j]];

737 
useatoms.nres=max(useatoms.nres,useatoms.atom[ctr].resind+1);

738 
useatoms.resinfo[useatoms.atom[ctr].resind]=top>atoms.resinfo[useatoms.atom[ctr].resind]; /*copy resinfo*/

739 
} 
740  
741 
write_sto_conf(opt2fn("ob",nfile,fnm),"Order parameters",&useatoms,frout.x,NULL,frout.ePBC,frout.box); 
742 

743 
sfree(frout.x); 
744 
free_t_atoms(&useatoms,FALSE); 
745 
} 
746  
747 
int gmx_order(int argc,char *argv[]) 
748 
{ 
749 
const char *desc[] = { 
750 
"Compute the order parameter per atom for carbon tails. For atom i the",

751 
"vector i1, i+1 is used together with an axis. ",

752 
"The index file should contain only the groups to be used for calculations,",

753 
"with each group of equivalent carbons along the relevant acyl chain in its own",

754 
"group. There should not be any generic groups (like System, Protein) in the index",

755 
"file to avoid confusing the program (this is not relevant to tetrahedral order",

756 
"parameters however, which only work for water anyway).[PAR]",

757 
"The program can also give all",

758 
"diagonal elements of the order tensor and even calculate the deuterium",

759 
"order parameter Scd (default). If the option [TT]szonly[tt] is given, only one",

760 
"order tensor component (specified by the [TT]d[tt] option) is given and the",

761 
"order parameter per slice is calculated as well. If [TT]szonly[tt] is not",

762 
"selected, all diagonal elements and the deuterium order parameter is",

763 
"given.[PAR]"

764 
"The tetrahedrality order parameters can be determined",

765 
"around an atom. Both angle an distance order parameters are calculated. See",

766 
"P.L. Chau and A.J. Hardwick, Mol. Phys., 93, (1998), 511518.",

767 
"for more details.[BR]",

768 
""

769 
}; 
770  
771 
static int nslices = 1; /* nr of slices defined */ 
772 
static gmx_bool bSzonly = FALSE; /* True if only Sz is wanted */ 
773 
static gmx_bool bUnsat = FALSE; /* True if carbons are unsat. */ 
774 
static const char *normal_axis[] = { NULL, "z", "x", "y", NULL }; 
775 
static gmx_bool permolecule = FALSE; /*compute on a permolecule basis */ 
776 
static gmx_bool radial = FALSE; /*compute a radial membrane normal */ 
777 
static gmx_bool distcalc = FALSE; /*calculate distance from a reference group */ 
778 
t_pargs pa[] = { 
779 
{ "d", FALSE, etENUM, {normal_axis},

780 
"Direction of the normal on the membrane" },

781 
{ "sl", FALSE, etINT, {&nslices},

782 
"Calculate order parameter as function of boxlength, dividing the box"

783 
" in #nr slices." },

784 
{ "szonly", FALSE, etBOOL,{&bSzonly},

785 
"Only give Sz element of order tensor. (axis can be specified with [TT]d[tt])" },

786 
{ "unsat", FALSE, etBOOL,{&bUnsat},

787 
"Calculate order parameters for unsaturated carbons. Note that this can"

788 
"not be mixed with normal order parameters." },

789 
{ "permolecule", FALSE, etBOOL,{&permolecule},

790 
"Compute permolecule Scd order parameters" },

791 
{ "radial", FALSE, etBOOL,{&radial},

792 
"Compute a radial membrane normal" },

793 
{ "calcdist", FALSE, etBOOL,{&distcalc},

794 
"Compute distance from a reference (currently defined only for radial and permolecule)" },

795 
}; 
796  
797 
rvec *order; /* order par. for each atom */

798 
real **slOrder; /* same, per slice */

799 
real slWidth = 0.0; /* width of a slice */ 
800 
char **grpname; /* groupnames */ 
801 
int ngrps, /* nr. of groups */ 
802 
i, 
803 
axis=0; /* normal axis */ 
804 
t_topology *top; /* topology */

805 
int ePBC;

806 
atom_id *index, /* indices for a */

807 
*a; /* atom numbers in each group */

808 
t_blocka *block; /* data from index file */

809 
t_filenm fnm[] = { /* files for g_order */

810 
{ efTRX, "f", NULL, ffREAD }, /* trajectory file */ 
811 
{ efNDX, "n", NULL, ffREAD }, /* index file */ 
812 
{ efNDX, "nr", NULL, ffREAD }, /* index for radial axis calculation */ 
813 
{ efTPX, NULL, NULL, ffREAD }, /* topology file */ 
814 
{ efXVG,"o","order", ffWRITE }, /* xvgr output file */ 
815 
{ efXVG,"od","deuter", ffWRITE }, /* xvgr output file */ 
816 
{ efPDB,"ob",NULL, ffWRITE }, /* write Scd as B factors to PDB if permolecule */ 
817 
{ efXVG,"os","sliced", ffWRITE }, /* xvgr output file */ 
818 
{ efXVG,"Sg","sgang", ffOPTWR }, /* xvgr output file */ 
819 
{ efXVG,"Sk","skdist", ffOPTWR }, /* xvgr output file */ 
820 
{ efXVG,"Sgsl","sgangslice", ffOPTWR }, /* xvgr output file */ 
821 
{ efXVG,"Sksl","skdistslice", ffOPTWR }, /* xvgr output file */ 
822 
}; 
823 
gmx_bool bSliced = FALSE; /* True if box is sliced */

824 
#define NFILE asize(fnm)

825 
real **distvals=NULL;

826 
const char *sgfnm,*skfnm,*ndxfnm,*tpsfnm,*trxfnm; 
827 
output_env_t oenv; 
828  
829 
CopyRight(stderr,argv[0]);

830 

831 
parse_common_args(&argc,argv,PCA_CAN_VIEW  PCA_CAN_TIME  PCA_BE_NICE, 
832 
NFILE,fnm,asize(pa),pa,asize(desc),desc,0, NULL,&oenv); 
833 
if (nslices < 1) 
834 
gmx_fatal(FARGS,"Can not have nslices < 1");

835 
sgfnm = opt2fn_null("Sg",NFILE,fnm);

836 
skfnm = opt2fn_null("Sk",NFILE,fnm);

837 
ndxfnm = opt2fn_null("n",NFILE,fnm);

838 
tpsfnm = ftp2fn(efTPX,NFILE,fnm); 
839 
trxfnm = ftp2fn(efTRX,NFILE,fnm); 
840 

841 
/* Calculate axis */

842 
if (strcmp(normal_axis[0],"x") == 0) axis = XX; 
843 
else if (strcmp(normal_axis[0],"y") == 0) axis = YY; 
844 
else if (strcmp(normal_axis[0],"z") == 0) axis = ZZ; 
845 
else gmx_fatal(FARGS,"Invalid axis, use x, y or z"); 
846 

847 
switch (axis) {

848 
case 0: 
849 
fprintf(stderr,"Taking x axis as normal to the membrane\n");

850 
break;

851 
case 1: 
852 
fprintf(stderr,"Taking y axis as normal to the membrane\n");

853 
break;

854 
case 2: 
855 
fprintf(stderr,"Taking z axis as normal to the membrane\n");

856 
break;

857 
} 
858 

859 
/* tetraheder order parameter */

860 
if (skfnm  sgfnm) {

861 
/* If either of theoptions is set we compute both */

862 
sgfnm = opt2fn("Sg",NFILE,fnm);

863 
skfnm = opt2fn("Sk",NFILE,fnm);

864 
calc_tetra_order_parm(ndxfnm,tpsfnm,trxfnm,sgfnm,skfnm,nslices,axis, 
865 
opt2fn("Sgsl",NFILE,fnm),opt2fn("Sksl",NFILE,fnm), 
866 
oenv); 
867 
/* view xvgr files */

868 
do_view(oenv,opt2fn("Sg",NFILE,fnm), NULL); 
869 
do_view(oenv,opt2fn("Sk",NFILE,fnm), NULL); 
870 
if (nslices > 1) { 
871 
do_view(oenv,opt2fn("Sgsl",NFILE,fnm), NULL); 
872 
do_view(oenv,opt2fn("Sksl",NFILE,fnm), NULL); 
873 
} 
874 
} 
875 
else {

876 
/* tail order parameter */

877 

878 
if (nslices > 1) { 
879 
bSliced = TRUE; 
880 
fprintf(stderr,"Dividing box in %d slices.\n\n", nslices);

881 
} 
882 

883 
if (bSzonly)

884 
fprintf(stderr,"Only calculating Sz\n");

885 
if (bUnsat)

886 
fprintf(stderr,"Taking carbons as unsaturated!\n");

887 

888 
top = read_top(ftp2fn(efTPX,NFILE,fnm),&ePBC); /* read topology file */

889 

890 
block = init_index(ftp2fn(efNDX,NFILE,fnm),&grpname); 
891 
index = block>index; /* get indices from t_block block */

892 
a = block>a; /* see block.h */

893 
ngrps = block>nr; 
894 

895 
if (permolecule)

896 
{ 
897 
nslices = index[1]  index[0]; /*I think this assumes contiguous lipids in topology*/ 
898 
fprintf(stderr,"Calculating Scd order parameters for each of %d molecules\n",nslices);

899 
} 
900 

901 
if (distcalc)

902 
{ 
903 
radial=TRUE; 
904 
fprintf(stderr,"Calculating radial distances\n");

905 
if (!permolecule)

906 
gmx_fatal(FARGS,"Cannot yet output radial distances without permolecule\n");

907 
} 
908 

909 
/* show atomtypes, to check if index file is correct */

910 
print_types(index, a, ngrps, grpname, top); 
911  
912 
calc_order(ftp2fn(efTRX,NFILE,fnm), index, a, &order, 
913 
&slOrder, &slWidth, nslices, bSliced, bUnsat, 
914 
top, ePBC, ngrps, axis,permolecule,radial,distcalc,opt2fn_null("nr",NFILE,fnm),&distvals, oenv);

915 

916 
if (radial)

917 
ngrps; /*don't print the last groupwas used for

918 
centerofmass determination*/

919 

920 
order_plot(order, slOrder, opt2fn("o",NFILE,fnm), opt2fn("os",NFILE,fnm), 
921 
opt2fn("od",NFILE,fnm), ngrps, nslices, slWidth, bSzonly,permolecule,distvals,oenv);

922  
923 
if (opt2bSet("ob",NFILE,fnm)) 
924 
{ 
925 
if (!permolecule)

926 
fprintf(stderr, 
927 
"Won't write Bfactors with averaged order parameters; use permolecule\n");

928 
else

929 
write_bfactors(fnm,NFILE,index,a,nslices,ngrps,slOrder,top,distvals,oenv); 
930 
} 
931  
932 

933 
do_view(oenv,opt2fn("o",NFILE,fnm), NULL); /* view xvgr file */ 
934 
do_view(oenv,opt2fn("os",NFILE,fnm), NULL); /* view xvgr file */ 
935 
do_view(oenv,opt2fn("od",NFILE,fnm), NULL); /* view xvgr file */ 
936 
} 
937 

938 
thanx(stderr); 
939 
if (distvals!=NULL) 
940 
{ 
941 
for (i=0;i<nslices;++i) 
942 
sfree(distvals[i]); 
943 
sfree(distvals); 
944 
} 
945 

946 
return 0; 
947 
} 