Project

General

Profile

gmx_order_4.5.4_modifiedbyChrisNeale.c

Chris Neale, 06/30/2015 03:46 PM

 
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) 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.
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 box-axis        */
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 center-of-mass 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
                /*center-of-mass 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[i-1]+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[i-1]+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 PSEUDO-ATOM 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[i-1]+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[i-1]+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 PSEUDO-ATOM 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 per-molecule*/
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 i-1, 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), 511-518.",
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 per-molecule 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 per-molecule 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","sg-ang", ffOPTWR },      /* xvgr output file           */
819
    { efXVG,"-Sk","sk-dist", ffOPTWR },     /* xvgr output file           */
820
    { efXVG,"-Sgsl","sg-ang-slice", ffOPTWR },      /* xvgr output file           */
821
    { efXVG,"-Sksl","sk-dist-slice", 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 group--was used for 
918
                           center-of-mass 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 B-factors 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
}