Project

General

Profile

gmx_rdf.c

Patched version of gmd_rdf.c - David van der Spoel, 04/18/2006 08:28 AM

 
1
/*
2
 * $Id: gmx_rdf.c,v 1.12.2.3 2006/03/23 07:43:05 spoel Exp $
3
 * 
4
 *                This source code is part of
5
 * 
6
 *                 G   R   O   M   A   C   S
7
 * 
8
 *          GROningen MAchine for Chemical Simulations
9
 * 
10
 *                        VERSION 3.2.0
11
 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12
 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13
 * Copyright (c) 2001-2004, The GROMACS development team,
14
 * check out http://www.gromacs.org for more information.
15

16
 * This program is free software; you can redistribute it and/or
17
 * modify it under the terms of the GNU General Public License
18
 * as published by the Free Software Foundation; either version 2
19
 * of the License, or (at your option) any later version.
20
 * 
21
 * If you want to redistribute modifications, please consider that
22
 * scientific software is very special. Version control is crucial -
23
 * bugs must be traceable. We will be happy to consider code for
24
 * inclusion in the official distribution, but derived work must not
25
 * be called official GROMACS. Details are found in the README & COPYING
26
 * files - if they are missing, get the official version at www.gromacs.org.
27
 * 
28
 * To help us fund GROMACS development, we humbly ask that you cite
29
 * the papers on the package - you can find them in the top README file.
30
 * 
31
 * For more info, check our website at http://www.gromacs.org
32
 * 
33
 * And Hey:
34
 * Green Red Orange Magenta Azure Cyan Skyblue
35
 */
36

    
37
#ifdef HAVE_CONFIG_H
38
#include <config.h>
39
#endif
40

    
41
#include <math.h>
42
#include <ctype.h>
43
#include "string2.h"
44
#include "sysstuff.h"
45
#include "typedefs.h"
46
#include "macros.h"
47
#include "vec.h"
48
#include "pbc.h"
49
#include "rmpbc.h"
50
#include "xvgr.h"
51
#include "copyrite.h"
52
#include "futil.h"
53
#include "statutil.h"
54
#include "tpxio.h"
55
#include "index.h"
56
#include "smalloc.h"
57
#include "fftgrid.h"
58
#include "calcgrid.h"
59
#include "nrnb.h"
60
#include "shift_util.h"
61
#include "pme.h"
62
#include "gstat.h"
63
#include "matio.h"
64

    
65
typedef struct
66
{
67
  char *label;
68
  int  elem,mass;
69
  real a[4], b[4], c;
70
} t_CM_table;
71

    
72
/*
73
 * 
74
 * f0[k] = c + [SUM a_i*EXP(-b_i*(k^2)) ]
75
 *             i=1,4
76
 */
77

    
78
const t_CM_table CM_t[] =
79
{
80

    
81
  { "H", 1,  1,   { 0.489918, 0.262003, 0.196767, 0.049879 },
82
    { 20.6593, 7.74039, 49.5519, 2.20159 },
83
    0.001305 },
84
  { "C", 6,  12, { 2.26069, 1.56165, 1.05075, 0.839259 },
85
    { 22.6907, 0.656665, 9.75618, 55.5949 },
86
    0.286977 },
87
  { "N", 7,  14,   { 12.2126, 3.13220, 2.01250, 1.16630 },     
88
    { 0.005700, 9.89330, 28.9975, 0.582600 },
89
    -11.529 },
90
  { "O", 8,  16, { 3.04850, 2.28680, 1.54630, 0.867000 },  
91
    { 13.2771, 5.70110, 0.323900, 32.9089 },
92
    0.250800 },
93
  { "Na", 11, 23,  { 3.25650, 3.93620, 1.39980, 1.00320 },       /*  Na 1+ */
94
    { 2.66710, 6.11530, 0.200100, 14.0390 }, 
95
    0.404000 }
96
};
97

    
98
#define NCMT asize(CM_t)
99

    
100
typedef struct
101
{
102
  int     n_angles;
103
  int     n_groups;
104
  double  lambda;
105
  double  energy;
106
  double  momentum;
107
  double  ref_k;
108
  double  **F;
109
  int     nSteps;
110
  int     total_n_atoms;
111
} structure_factor;
112

    
113
typedef struct
114
{
115
  rvec x;
116
  int  t;
117
} reduced_atom;
118

    
119
static void check_box_c(matrix box)
120
{
121
  if (fabs(box[ZZ][XX]) > GMX_REAL_EPS*box[ZZ][ZZ] ||
122
      fabs(box[ZZ][YY]) > GMX_REAL_EPS*box[ZZ][ZZ])
123
    gmx_fatal(FARGS,
124
              "The last box vector is not parallel to the z-axis: %f %f %f",
125
              box[ZZ][XX],box[ZZ][YY],box[ZZ][ZZ]);
126
}
127

    
128
static void do_rdf(char *fnNDX,char *fnTPS,char *fnTRX,
129
                   char *fnRDF,char *fnCNRDF, char *fnHQ,
130
                   bool bCM,bool bXY,bool bPBC,
131
                   real cutoff,real binwidth,real fade,int ng)
132
{
133
  FILE       *fp;
134
  int        status;
135
  char       outf1[STRLEN],outf2[STRLEN];
136
  char       title[STRLEN];
137
  int        g,natoms,i,j,k,nbin,j0,j1,n,nframes;
138
  int        **count;
139
  char       **grpname;
140
  int        *isize,isize_cm=0,nrdf=0,max_i;
141
  atom_id    **index,*index_cm=NULL;
142
#if (defined SIZEOF_LONG_LONG_INT) && (SIZEOF_LONG_LONG_INT >= 8)    
143
  long long int *sum;
144
#else
145
  double     *sum;
146
#endif
147
  real       t,rmax2,cut2,r,r2,invbinw,normfac;
148
  real       segvol,spherevol,prev_spherevol,**rdf;
149
  rvec       *x,xcom,dx,*x_i1,xi;
150
  real       *inv_segvol,invvol,invvol_sum,rho;
151
  bool       *bExcl,bTop,bNonSelfExcl;
152
  matrix     box,box_pbc;
153
  int        **npairs,*nself;
154
  atom_id    ix,jx,***pairs;
155
  t_topology top;
156
  t_block    *excl;
157
  t_pbc      pbc;
158

    
159
  excl=NULL;
160
  
161
  if (fnTPS) {
162
    bTop=read_tps_conf(fnTPS,title,&top,&x,NULL,box,TRUE);
163
    mk_single_top(&top);
164
    if (bTop && !bCM)
165
      /* get exclusions from topology */
166
      excl=&(top.atoms.excl);
167
  }
168
  snew(grpname,ng+1);
169
  snew(isize,ng+1);
170
  snew(index,ng+1);
171
  fprintf(stderr,"\nSelect a reference group and %d group%s\n",
172
          ng,ng==1?"":"s");
173
  if (fnTPS)
174
    get_index(&top.atoms,fnNDX,ng+1,isize,index,grpname);
175
  else
176
    rd_index(fnNDX,ng+1,isize,index,grpname);
177
  
178
  natoms=read_first_x(&status,fnTRX,&t,&x,box);
179
  if ( !natoms )
180
    gmx_fatal(FARGS,"Could not read coordinates from statusfile\n");
181
  if (fnTPS)
182
    /* check with topology */
183
    if ( natoms > top.atoms.nr ) 
184
      gmx_fatal(FARGS,"Trajectory (%d atoms) does not match topology (%d atoms)",
185
                  natoms,top.atoms.nr);
186
  /* check with index groups */
187
  for (i=0; i<=ng; i++)
188
    for (j=0; j<isize[i]; j++)
189
      if ( index[i][j] >= natoms )
190
        gmx_fatal(FARGS,"Atom index (%d) in index group %s (%d atoms) larger "
191
                    "than number of atoms in trajectory (%d atoms)",
192
                    index[i][j],grpname[i],isize[i],natoms);
193
  
194
  if (bCM) {
195
    /* move index[0] to index_cm and make 'dummy' index[0] */
196
    isize_cm=isize[0];
197
    snew(index_cm,isize_cm);
198
    for(i=0; i<isize[0]; i++)
199
      index_cm[i]=index[0][i];
200
    isize[0]=1;
201
    index[0][0]=natoms;
202
    srenew(index[0],isize[0]);
203
    /* make space for center of mass */
204
    srenew(x,natoms+1);
205
  }
206
  
207
  /* initialize some handy things */
208
  copy_mat(box,box_pbc);
209
  if (bXY) {
210
    check_box_c(box);
211
    /* Make sure the z-height does not influence the cut-off */
212
    box_pbc[ZZ][ZZ] = 2*max(box[XX][XX],box[YY][YY]);
213
  }
214
  rmax2   = 0.99*0.99*max_cutoff2(box_pbc);
215
  nbin    = (int)(sqrt(rmax2) / binwidth) + 1;
216
  invbinw = 1.0 / binwidth;
217
  cut2   = sqr(cutoff);
218

    
219
  snew(count,ng);
220
  snew(pairs,ng);
221
  snew(npairs,ng);
222
  snew(nself,ng);
223

    
224
  snew(bExcl,natoms);
225
  max_i = 0;
226
  for(g=0; g<ng; g++) {
227
    if (isize[g+1] > max_i)
228
      max_i = isize[g+1];
229

    
230
    /* this is THE array */
231
    snew(count[g],nbin+1);
232
  
233
    /* make pairlist array for groups and exclusions */
234
    snew(pairs[g],isize[0]);
235
    snew(npairs[g],isize[0]);
236
    for(i=0; i<isize[0]; i++) {
237
      ix = index[0][i];
238
      for(j=0; j < natoms; j++)
239
        bExcl[j] = FALSE;
240
      /* exclusions? */
241
      if (excl)
242
        for( j = excl->index[ix]; j < excl->index[ix+1]; j++)
243
          bExcl[excl->a[j]]=TRUE;
244
      k = 0;
245
      snew(pairs[g][i], isize[g+1]);
246
      bNonSelfExcl = FALSE;
247
      for(j=0; j<isize[g+1]; j++) {
248
        jx = index[g+1][j];
249
        if (!bExcl[jx])
250
          pairs[g][i][k++]=jx;
251
        else if (ix != jx)
252
          /* Check if we have exclusions other than self exclusions */
253
          bNonSelfExcl = TRUE;
254
        if (ix == jx)
255
          nself[g]++;
256
      }
257
      if (bNonSelfExcl) {
258
        npairs[g][i]=k;
259
        srenew(pairs[g][i],npairs[g][i]);
260
      } else {
261
        /* Save a LOT of memory and some cpu cycles */
262
        npairs[g][i]=-1;
263
        sfree(pairs[g][i]);
264
      }
265
    }
266
  }
267
  sfree(bExcl);
268

    
269
  snew(x_i1,max_i);
270
  nframes = 0;
271
  invvol_sum = 0;
272
  do {
273
    /* Must init pbc every step because of pressure coupling */
274
    copy_mat(box,box_pbc);
275
    if (bPBC) {
276
      rm_pbc(&top.idef,natoms,box,x,x);
277
      if (bXY) {
278
        check_box_c(box);
279
        box_pbc[ZZ][ZZ] = 2*max(box[XX][XX],box[YY][YY]);
280
      }
281
      set_pbc(&pbc,box_pbc);
282

    
283
      if (bXY)
284
        /* Set z-size to 1 so we get the surface iso the volume */
285
        box_pbc[ZZ][ZZ] = 1;
286
    }
287
    invvol = 1/det(box_pbc);
288
    invvol_sum += invvol;
289

    
290
    if (bCM) {
291
      /* calculate centre of mass */
292
      clear_rvec(xcom);
293
      for(i=0; (i < isize_cm); i++) {
294
        ix = index_cm[i];
295
        rvec_inc(xcom,x[ix]);
296
      }
297
      /* store it in the first 'group' */
298
      for(j=0; (j<DIM); j++)
299
        x[index[0][0]][j] = xcom[j] / isize_cm;
300
    }
301

    
302
    for(g=0; g<ng; g++) {
303
      /* Copy the indexed coordinates to a continuous array */
304
      for(i=0; i<isize[g+1]; i++)
305
        copy_rvec(x[index[g+1][i]],x_i1[i]);
306
    
307
      for(i=0; i<isize[0]; i++) {
308
        copy_rvec(x[index[0][i]],xi);
309
        if (npairs[g][i] >= 0)
310
          /* Expensive loop, because of indexing */
311
          for(j=0; j<npairs[g][i]; j++) {
312
            jx=pairs[g][i][j];
313
            if (bPBC)
314
              pbc_dx(&pbc,xi,x[jx],dx);
315
            else
316
              rvec_sub(xi,x[jx],dx);
317
              
318
            if (bXY)
319
              r2 = dx[XX]*dx[XX] + dx[YY]*dx[YY];
320
            else 
321
              r2=iprod(dx,dx);
322
            if (r2>cut2 && r2<=rmax2)
323
              count[g][(int)(sqrt(r2)*invbinw)]++;
324
          }
325
        else
326
          /* Cheaper loop, no exclusions */
327
          for(j=0; j<isize[g+1]; j++) {
328
            if (bPBC)
329
              pbc_dx(&pbc,xi,x_i1[j],dx);
330
            else
331
              rvec_sub(xi,x_i1[j],dx);
332
            if (bXY)
333
              r2 = dx[XX]*dx[XX] + dx[YY]*dx[YY];
334
            else 
335
              r2=iprod(dx,dx);
336
            if (r2>cut2 && r2<=rmax2)
337
              count[g][(int)(sqrt(r2)*invbinw)]++;
338
          }
339
      }
340
    }
341
    nframes++;
342
  } while (read_next_x(status,&t,natoms,x,box));
343
  fprintf(stderr,"\n");
344
  
345
  close_trj(status);
346
  
347
  sfree(x);
348
  
349
  /* Average volume */
350
  invvol = invvol_sum/nframes;
351

    
352
  /* Calculate volume of sphere segments or length of circle segments */
353
  snew(inv_segvol,nbin);
354
  prev_spherevol=0;
355
  for(i=0; (i<nbin); i++) {
356
    r = (i+1)*binwidth;
357
    if (bXY) {
358
      spherevol=M_PI*r*r;
359
    } else {
360
      spherevol=(4.0/3.0)*M_PI*r*r*r;
361
    }
362
    segvol=spherevol-prev_spherevol;
363
    inv_segvol[i]=1.0/segvol;
364
    prev_spherevol=spherevol;
365
  }
366
  
367
  snew(rdf,ng);
368
  for(g=0; g<ng; g++) {
369
    /* We have to normalize by dividing by the number of frames */
370
    normfac = 1.0/(nframes*invvol*(isize[0]*isize[g+1] - nself[g]));
371
      
372
    /* Do the normalization */
373
    nrdf = max(nbin-1,1+(2*fade/binwidth));
374
    snew(rdf[g],nrdf);
375
    for(i=0; (i<nbin-1); i++) {
376
      r = (i+0.5)*binwidth;
377
      if ((fade > 0) && (r >= fade))
378
        rdf[g][i] = 1+(count[g][i]*inv_segvol[i]*normfac-1)*exp(-16*sqr(r/fade-1));
379
      else
380
        rdf[g][i] = count[g][i]*inv_segvol[i]*normfac;
381
    }
382
    for( ; (i<nrdf); i++)
383
      rdf[g][i] = 1.0;
384
  }
385
    
386
  fp=xvgropen(fnRDF,"Radial Distribution","r","");
387
  if (ng==1)
388
    fprintf(fp,"@ subtitle \"%s-%s\"\n",grpname[0],grpname[1]);
389
  else {
390
    fprintf(fp,"@ subtitle \"reference %s\"\n",grpname[0]);
391
    xvgr_legend(fp,ng,grpname+1);
392
  }
393
  for(i=0; (i<nrdf); i++) {
394
    fprintf(fp,"%10g", (i+0.5)*binwidth);
395
    for(g=0; g<ng; g++)
396
      fprintf(fp," %10g",rdf[g][i]);
397
    fprintf(fp,"\n");
398
  }
399
  ffclose(fp);
400
  
401
  do_view(fnRDF,NULL);
402

    
403
  /* h(Q) function: fourier transform of rdf */  
404
  if (fnHQ) {
405
    int nhq = 401;
406
    real *hq,*integrand,Q;
407
    
408
    /* Get a better number density later! */
409
    rho = isize[1]*invvol;
410
    snew(hq,nhq);
411
    snew(integrand,nrdf);
412
    for(i=0; (i<nhq); i++) {
413
      Q = i*0.5;
414
      integrand[0] = 0;
415
      for(j=1; (j<nrdf); j++) {
416
        r = (j+0.5)*binwidth;
417
        integrand[j]  = (Q == 0) ? 1.0 : sin(Q*r)/(Q*r);
418
        integrand[j] *= 4.0*M_PI*rho*r*r*(rdf[0][j]-1.0);
419
      }
420
      hq[i] = print_and_integrate(debug,nrdf,binwidth,integrand,NULL,0);
421
    }
422
    fp=xvgropen(fnHQ,"h(Q)","Q(/nm)","h(Q)");
423
    for(i=0; (i<nhq); i++) 
424
      fprintf(fp,"%10g %10g\n",i*0.5,hq[i]);
425
    ffclose(fp);
426
    do_view(fnHQ,NULL);
427
    sfree(hq);
428
    sfree(integrand);
429
  }
430
  
431
  if (fnCNRDF) {  
432
    normfac = 1.0/(isize[0]*nframes);
433
    fp=xvgropen(fnCNRDF,"Cumulative Number RDF","r","number");
434
    if (ng==1)
435
      fprintf(fp,"@ subtitle \"%s-%s\"\n",grpname[0],grpname[1]);
436
    else {
437
      fprintf(fp,"@ subtitle \"reference %s\"\n",grpname[0]);
438
      xvgr_legend(fp,ng,grpname+1);
439
    }
440
    snew(sum,ng);
441
    for(i=0; (i<nbin-1); i++) {
442
      fprintf(fp,"%10g",i*binwidth);
443
      for(g=0; g<ng; g++) {
444
        fprintf(fp," %10g",(real)((double)sum[g]*normfac));
445
        sum[g] += count[g][i];
446
      }
447
      fprintf(fp,"\n");
448
    }
449
    ffclose(fp);
450
    sfree(sum);
451
    
452
    do_view(fnCNRDF,NULL);
453
  }
454

    
455
  for(g=0; g<ng; g++)
456
    sfree(rdf[g]);
457
  sfree(rdf);
458
}
459

    
460
typedef struct {
461
  int  ndata;
462
  real kkk;
463
  real intensity;
464
} t_xdata;
465

    
466
int comp_xdata(const void *a,const void *b)
467
{
468
  t_xdata *xa,*xb;
469
  real tmp;
470
  
471
  xa = (t_xdata *)a;
472
  xb = (t_xdata *)b;
473
  
474
  if (xa->ndata == 0)
475
    return 1;
476
  else if (xb->ndata == 0)
477
    return -1;
478
  else {
479
    tmp = xa->kkk - xb->kkk;
480
    if (tmp < 0)
481
      return -1;
482
    else if (tmp > 0)
483
      return 1;
484
    else return 0;
485
  }
486
}
487

    
488
static t_xdata *init_xdata(int nx,int ny)
489
{
490
  int     ix,iy,i,j,maxkx,maxky;
491
  t_xdata *data;
492
  
493
  maxkx = (nx+1)/2;
494
  maxky = (ny+1)/2;
495
  snew(data,nx*ny);
496
  for(i=0; (i<nx); i++) {
497
    for(j=0; (j<ny); j++) {
498
      ix = abs((i < maxkx) ? i : (i - nx)); 
499
      iy = abs((j < maxky) ? j : (j - ny)); 
500
      data[ix*ny+iy].kkk = sqrt(ix*ix+iy*iy);
501
    }
502
  }
503
  return data;
504
}
505

    
506
static void extract_sq(t_fftgrid *fftgrid,int nbin,real k_max,real lambda,
507
                       real count[],rvec box,int npixel,real *map[],
508
                       t_xdata data[])
509
{
510
  int     nx,ny,nz,nx2,ny2,nz2,la2,la12;
511
  t_complex *ptr,*p0;
512
  int     i,j,k,maxkx,maxky,maxkz,n,ind,ix,iy;
513
  real    k1,kxy2,kz2,k2,z,kxy,kxy_max,cos_theta2,ttt,factor;
514
  rvec    lll,kk;
515
  
516
  /*calc_lll(box,lll);
517
    k_max   = nbin/factor;
518
    kxy_max = k_max/sqrt(3);*/
519
  unpack_fftgrid(fftgrid,&nx,&ny,&nz,&nx2,&ny2,&nz2,
520
                 &la2,&la12,FALSE,(real **)&ptr);
521
  /* This bit copied from pme.c */
522
  maxkx = (nx+1)/2;
523
  maxky = (ny+1)/2;
524
  maxkz = nz/2+1;
525
  factor = nbin/k_max;
526
  for(i=0; (i<nx); i++) {
527
#define IDX(i,n)  (i<=n/2) ? (i) : (i-n)
528
    kk[XX] = IDX(i,nx);
529
    for(j=0; (j<ny); j++) {
530
      kk[YY] = IDX(j,ny);
531
      ind = INDEX(i,j,0);
532
      p0  = ptr + ind;
533
      for(k=0; (k<maxkz); k++,p0++) {
534
        if ((i==0) && (j==0) && (k==0))
535
          continue;
536
        kk[ZZ] = IDX(k,nz);
537
        z   = sqrt(sqr(p0->re)+sqr(p0->im));
538
        kxy2 = sqr(kk[XX]) + sqr(kk[YY]);
539
        k2   = kxy2+sqr(kk[ZZ]);
540
        k1   = sqrt(k2);
541
        ind  = k1*factor;
542
        if (ind < nbin) {
543
          /* According to:
544
           * R. W. James (1962), 
545
           * The Optical Principles of the Diffraction of X-Rays,
546
           * Oxbow press, Woodbridge Connecticut
547
           * the intensity is proportional to (eq. 9.10):
548
           * I = C (1+cos^2 [2 theta])/2 FFT
549
           * And since
550
           * cos[2 theta] = cos^2[theta] - sin^2[theta] = 2 cos^2[theta] - 1 
551
           * we can compute the prefactor straight from cos[theta]
552
           */
553
          cos_theta2  = kxy2/k2;
554
          /*ttt         = z*0.5*(1+sqr(2*cos_theta2-1));*/
555
          ttt         = z*0.5*(1+cos_theta2);
556
          count[ind] += ttt;
557
          ix = ((i < maxkx) ? i : (i - nx)); 
558
          iy = ((j < maxky) ? j : (j - ny));
559
          map[npixel/2+ix][npixel/2+iy] += ttt; 
560
          data[abs(ix)*ny+abs(iy)].ndata     += 1;
561
          data[abs(ix)*ny+abs(iy)].intensity += ttt;
562
        }
563
        else
564
          fprintf(stderr,"k (%g) > k_max (%g)\n",k1,k_max);
565
      }
566
    }
567
  }
568
}
569

    
570
typedef struct {
571
  char *name;
572
  int  nelec;
573
} t_element;
574

    
575
static void do_sq(char *fnNDX,char *fnTPS,char *fnTRX,char *fnSQ,
576
                  char *fnXPM,real grid,real lambda,real distance,
577
                  int npixel,int nlevel)
578
{
579
  FILE       *fp;
580
  t_element  elem[] = { { "H", 1 }, { "C", 6 }, { "N", 7 }, { "O", 8 }, { "F", 9 }, { "S", 16 } };
581
#define NELEM asize(elem)
582
  int        status;
583
  char       title[STRLEN],*aname;
584
  int        natoms,i,j,k,nbin,j0,j1,n,nframes,pme_order;
585
  real       *count,**map;
586
  char       *grpname;
587
  int        isize;
588
  atom_id    *index;
589
  real       I0,C,t,k_max,factor,yfactor,segvol;
590
  rvec       *x,*xndx,box_size,kk,lll;
591
  real       fj0,*fj,max_spacing,r,lambda_1;
592
  bool       *bExcl;
593
  matrix     box;
594
  int        nx,ny,nz,nelectron;
595
  atom_id    ix,jx,**pairs;
596
  splinevec  *theta;
597
  t_topology top;
598
  t_fftgrid  *fftgrid;
599
  t_nrnb     nrnb;
600
  t_xdata    *data;
601
    
602
  /*  bTop=read_tps_conf(fnTPS,title,&top,&x,NULL,box,TRUE); */
603

    
604
  fprintf(stderr,"\nSelect group for structure factor computation:\n");
605
  get_index(&top.atoms,fnNDX,1,&isize,&index,&grpname);
606
  natoms=read_first_x(&status,fnTRX,&t,&x,box);
607
  if (isize < top.atoms.nr)
608
    snew(xndx,isize);
609
  else
610
    xndx = x;
611
  fprintf(stderr,"\n");
612
  
613
  init_nrnb(&nrnb);
614
    
615
  if ( !natoms )
616
    gmx_fatal(FARGS,"Could not read coordinates from statusfile\n");
617
  /* check with topology */
618
  if ( natoms > top.atoms.nr ) 
619
    gmx_fatal(FARGS,"Trajectory (%d atoms) does not match topology (%d atoms)",
620
                natoms,top.atoms.nr);
621
        
622
  /* Atomic scattering factors */
623
  snew(fj,isize);
624
  I0 = 0;
625
  nelectron = 0;
626
  for(i=0; (i<isize); i++) {
627
    aname = *(top.atoms.atomname[index[i]]);
628
    fj0 = 1;
629
    if (top.atoms.atom[i].ptype == eptAtom) {
630
      for(j=0; (j<NELEM); j++)
631
        if (aname[0] == elem[j].name[0]) {
632
          fj0 = elem[j].nelec;
633
          break;
634
        }
635
      if (j == NELEM)
636
        fprintf(stderr,"Warning: don't know number of electrons for atom %s\n",aname);
637
    }
638
    /* Correct for partial charge */
639
    fj[i] = fj0 - top.atoms.atom[index[i]].q;
640
    
641
    nelectron += fj[i];
642
    
643
    I0 += sqr(fj[i]);
644
  }
645
  if (debug) {
646
    /* Dump scattering factors */
647
    for(i=0; (i<isize); i++)
648
      fprintf(debug,"Atom %3s-%5d q = %10.4f  f = %10.4f\n",
649
              *(top.atoms.atomname[index[i]]),index[i],
650
              top.atoms.atom[index[i]].q,fj[i]);
651
  }
652

    
653
  /* Constant for scattering */
654
  C = sqr(1.0/(ELECTRONMASS_keV*KILO*ELECTRONVOLT*1e7*distance));
655
  fprintf(stderr,"C is %g\n",C);
656
  
657
  /* This bit is dimensionless */
658
  nx = ny = nz = 0;
659
  max_spacing = calc_grid(box,grid,&nx,&ny,&nz,1);        
660
  pme_order   = max(4,1+(0.2/grid));
661
  npixel      = max(nx,ny);
662
  data        = init_xdata(nx,ny);
663
  
664
  fprintf(stderr,"Largest grid spacing: %g nm, pme_order %d, %dx%d pixel on image\n",
665
          max_spacing,pme_order,npixel,npixel);
666
  fftgrid = init_pme(stdout,NULL,nx,ny,nz,pme_order,isize,FALSE,FALSE,eewg3D);
667
    
668
  /* Determine largest k vector length. */
669
  k_max = 1+sqrt(sqr(1+nx/2)+sqr(1+ny/2)+sqr(1+nz/2));
670

    
671
  /* this is the S(q) array */
672
  nbin = npixel;
673
  snew(count,nbin+1);
674
  snew(map,npixel);
675
  for(i=0; (i<npixel); i++)
676
    snew(map[i],npixel);
677
  
678
  nframes = 0;
679
  do {
680
    /* Put the atoms with scattering factor on a grid. Misuses
681
     * an old routine from the PPPM code.
682
     */
683
    for(j=0; (j<DIM); j++)
684
      box_size[j] = box[j][j];
685
    
686
    /* Scale coordinates to the wavelength */
687
    for(i=0; (i<isize); i++)
688
      copy_rvec(x[index[i]],xndx[i]);
689
      
690
    /* put local atoms on grid. */
691
    spread_on_grid(stdout,fftgrid,isize,pme_order,xndx,fj,box,FALSE,TRUE);
692

    
693
    /* FFT the density */
694
    gmxfft3D(fftgrid,GMX_FFT_REAL_TO_COMPLEX,NULL);  
695
    
696
    /* Extract the Sq function and sum it into the average array */
697
    extract_sq(fftgrid,nbin,k_max,lambda,count,box_size,npixel,map,data);
698
    
699
    nframes++;
700
  } while (read_next_x(status,&t,natoms,x,box));
701
  fprintf(stderr,"\n");
702
  
703
  close_trj(status);
704
  
705
  sfree(x);
706

    
707
  /* Normalize it ?? */  
708
  factor  = k_max/(nbin);
709
  yfactor = (1.0/nframes)/*(1.0/fftgrid->nxyz)*/;
710
  fp=xvgropen(fnSQ,"Structure Factor","q (1/nm)","S(q)");
711
  fprintf(fp,"@ subtitle \"Lambda = %g nm. Grid spacing = %g nm\"\n",
712
          lambda,grid);
713
  factor *= lambda;
714
  for(i=0; i<nbin; i++) {
715
    r      = (i+0.5)*factor*2*M_PI;
716
    segvol = 4*M_PI*sqr(r)*factor;
717
    fprintf(fp,"%10g %10g\n",r,count[i]*yfactor/segvol);
718
  }
719
  ffclose(fp);
720
  
721
  do_view(fnSQ,NULL);
722

    
723
  if (fnXPM) {
724
    t_rgb rhi = { 0,0,0 }, rlo = { 1,1,1 };
725
    real *tx,*ty,hi,inv_nframes;
726
    
727
    hi = 0;
728
    inv_nframes = 1.0/nframes;
729
    snew(tx,npixel);
730
    snew(ty,npixel);
731
    for(i=0; (i<npixel); i++) {
732
      tx[i] = i-npixel/2;
733
      ty[i] = i-npixel/2;
734

    
735
      for(j=0; (j<npixel); j++) { 
736
        map[i][j] *= inv_nframes;
737
        hi         = max(hi,map[i][j]);
738
      }
739
    }
740
      
741
    fp = ffopen(fnXPM,"w");
742
    write_xpm(fp,0,"Diffraction Image","Intensity","kx","ky",
743
              nbin,nbin,tx,ty,map,0,hi,rlo,rhi,&nlevel);
744
    fclose(fp);
745
    sfree(tx);
746
    sfree(ty);
747

    
748
    /* qsort(data,nx*ny,sizeof(data[0]),comp_xdata);    
749
       fp = ffopen("test.xvg","w");
750
       for(i=0; (i<nx*ny); i++) {
751
       if (data[i].ndata != 0) {
752
       fprintf(fp,"%10.3f  %10.3f\n",data[i].kkk,data[i].intensity/data[i].ndata);
753
       }
754
       }
755
       fclose(fp);
756
    */
757
  }
758
}
759

    
760
t_complex *** rc_tensor_allocation(int x, int y, int z)
761
{
762
  t_complex ***t;
763
  int i,j;
764
  
765
  snew(t,x);
766
  t = (t_complex ***)calloc(x,sizeof(t_complex**));
767
  if(!t) exit(fprintf(stderr,"\nallocation error"));
768
  t[0] = (t_complex **)calloc(x*y,sizeof(t_complex*));
769
  if(!t[0]) exit(fprintf(stderr,"\nallocation error"));
770
  t[0][0] = (t_complex *)calloc(x*y*z,sizeof(t_complex));
771
  if(!t[0][0]) exit(fprintf(stderr,"\nallocation error"));
772
  
773
  for( j = 1 ; j < y ; j++) 
774
    t[0][j] = t[0][j-1] + z;
775
  for( i = 1 ; i < x ; i++) {
776
    t[i] = t[i-1] + y;
777
    t[i][0] = t[i-1][0] + y*z;
778
    for( j = 1 ; j < y ; j++) 
779
      t[i][j] = t[i][j-1] + z;
780
  }
781
  return t;
782
}
783
    
784
int return_atom_type (char *name)
785
{
786
  typedef struct {
787
    char *name;
788
    int  nh;
789
  } t_united_h;
790
  t_united_h uh[] = {
791
    { "CH1", 1 }, { "CH2", 2 }, { "CH3", 3 }, 
792
    { "CS1", 1 }, { "CS2", 2 }, { "CS3", 3 }, 
793
    { "CP1", 1 }, { "CP2", 2 }, { "CP3", 3 }
794
  };
795
  int i;
796

    
797
  for(i-0; (i<asize(uh)); i++) 
798
    if (strcmp(name,uh[i].name) == 0)
799
      return NCMT-1+uh[i].nh;
800
      
801
  for(i=0; (i<NCMT); i++)
802
    if (strncmp (name, CM_t[i].label,strlen(CM_t[i].label)) == 0)
803
      return i;
804
  gmx_fatal(FARGS,"\nError: atom (%s) not in list (%d types checked)!\n", 
805
            name,i);
806
  
807
  return 0;
808
}
809

    
810
double CMSF (int type,int nh,double lambda, double sin_theta)
811
/* 
812
 * return Cromer-Mann fit for the atomic scattering factor:
813
 * sin_theta is the sine of half the angle between incoming and scattered
814
 * vectors. See g_sq.h for a short description of CM fit.
815
 */
816
{
817
  int i;
818
  double tmp = 0.0, k2;
819
  
820
  /*
821
   *  united atoms case
822
   *  CH2 / CH3 groups  
823
   */
824
  if (nh > 0) {
825
    tmp = (CMSF (return_atom_type ("C"),0,lambda, sin_theta) +
826
           nh*CMSF (return_atom_type ("H"),0,lambda, sin_theta));
827
  }
828
  /* all atom case */
829
  else {
830
    k2 = (sqr (sin_theta) / sqr (10.0 * lambda));
831
    tmp = CM_t[type].c;
832
    for (i = 0; (i < 4); i++)
833
      tmp += CM_t[type].a[i] * exp (-CM_t[type].b[i] * k2);
834
  }
835
  return tmp;
836
}
837

    
838
real **compute_scattering_factor_table (structure_factor * sf,int *nsftable)
839
{
840
/*
841
 *  this function build up a table of scattering factors for every atom
842
 *  type and for every scattering angle.
843
 */
844
    int i, j;
845
    double sin_theta,q,hc=1239.842;
846
    real ** sf_table;
847

    
848
    /* \hbar \omega \lambda = hc = 1239.842 eV * nm */
849
    sf->momentum = ((double) (2. * 1000.0 * M_PI * sf->energy) / hc);
850
    sf->lambda = hc / (1000.0 * sf->energy);
851
    fprintf (stderr, "\nwavelenght = %f nm\n", sf->lambda);
852
    *nsftable = NCMT+3;
853
    snew (sf_table,*nsftable);
854
    for (i = 0; (i < *nsftable); i++) {
855
        snew (sf_table[i], sf->n_angles);
856
        for (j = 0; j < sf->n_angles; j++) {
857
            q = ((double) j * sf->ref_k);
858
            /* theta is half the angle between incoming 
859
               and scattered wavevectors. */
860
            sin_theta = q / (2.0 * sf->momentum);
861
            if (i < NCMT)
862
              sf_table[i][j] = CMSF (i,0,sf->lambda, sin_theta);
863
            else
864
              sf_table[i][j] = CMSF (i,i-NCMT+1,sf->lambda, sin_theta);
865
        }
866
    }
867
    return sf_table;
868
}
869

    
870
int * create_indexed_atom_type (reduced_atom * atm, int size)
871
{
872
/* 
873
 * create an index of the atom types found in a  group
874
 * i.e.: for water index_atp[0]=type_number_of_O and 
875
 *                 index_atp[1]=type_number_of_H
876
 * 
877
 * the last element is set to 0 
878
 */
879
    int *index_atp, i, i_tmp, j;
880

    
881
    snew (index_atp, 1);
882
    i_tmp = 1;
883
    index_atp[0] = atm[0].t;
884
    for (i = 1; i < size; i++) {
885
        for (j = 0; j < i_tmp; j++)
886
            if (atm[i].t == index_atp[j])
887
                break;
888
        if (j == i_tmp) {        /* i.e. no indexed atom type is  == to atm[i].t */
889
            i_tmp++;
890
            srenew (index_atp, i_tmp * sizeof (int));
891
            index_atp[i_tmp - 1] = atm[i].t;
892
        }
893
    }
894
    i_tmp++;
895
    srenew (index_atp, i_tmp * sizeof (int));
896
    index_atp[i_tmp - 1] = 0;
897
    return index_atp;
898
}
899

    
900
void rearrange_atoms (reduced_atom * positions, t_trxframe *fr, atom_id * index,
901
                      int isize, t_topology * top, bool flag)
902
/* given the group's index, return the (continuous) array of atoms */
903
{
904
  int i;
905
  
906
  if (flag)
907
    for (i = 0; i < isize; i++)
908
      positions[i].t =
909
        return_atom_type (*(top->atoms.atomname[index[i]]));
910
  for (i = 0; i < isize; i++)
911
    copy_rvec (fr->x[index[i]], positions[i].x);
912
}
913

    
914

    
915
int atp_size (int *index_atp)
916
{
917
    int i = 0;
918

    
919
    while (index_atp[i])
920
        i++;
921
    return i;
922
}
923

    
924
void compute_structure_factor (structure_factor * sf, matrix box,
925
                               reduced_atom * red, int isize, real start_q,
926
                               real end_q, int group,real **sf_table)
927
{
928
    t_complex ***tmpSF;
929
    rvec k_factor;
930
    real kdotx, asf, kx, ky, kz, krr;
931
    int kr, maxkx, maxky, maxkz, i, j, k, p, *counter;
932

    
933

    
934
    k_factor[XX] = 2 * M_PI / box[XX][XX];
935
    k_factor[YY] = 2 * M_PI / box[YY][YY];
936
    k_factor[ZZ] = 2 * M_PI / box[ZZ][ZZ];
937

    
938
    maxkx = (int) rint (end_q / k_factor[XX]);
939
    maxky = (int) rint (end_q / k_factor[YY]);
940
    maxkz = (int) rint (end_q / k_factor[ZZ]);
941

    
942
    snew (counter, sf->n_angles);
943

    
944
    tmpSF = rc_tensor_allocation(maxkx,maxky,maxkz);
945
/*
946
 * The big loop...
947
 * compute real and imaginary part of the structure factor for every
948
 * (kx,ky,kz)) 
949
 */
950
    fprintf(stderr,"\n");
951
    for (i = 0; i < maxkx; i++) {
952
        fprintf (stderr,"\rdone %3.1f%%     ", (double)(100.0*(i+1))/maxkx);
953
        kx = i * k_factor[XX];
954
        for (j = 0; j < maxky; j++) {
955
            ky = j * k_factor[YY];
956
            for (k = 0; k < maxkz; k++)
957
                if (i != 0 || j != 0 || k != 0) {
958
                    kz = k * k_factor[ZZ];
959
                    krr = sqrt (sqr (kx) + sqr (ky) + sqr (kz));
960
                    if (krr >= start_q && krr <= end_q) {
961
                        kr = (int) rint (krr/sf->ref_k);
962
                        if (kr < sf->n_angles) {
963
                            counter[kr]++;  /* will be used for the copmutation 
964
                                               of the average*/
965
                            for (p = 0; p < isize; p++) {
966
                                    
967
                                asf = sf_table[red[p].t][kr];
968

    
969
                                kdotx = kx * red[p].x[XX] +
970
                                    ky * red[p].x[YY] + kz * red[p].x[ZZ];
971
                                
972
                                tmpSF[i][j][k].re += cos (kdotx) * asf;
973
                                tmpSF[i][j][k].im += sin (kdotx) * asf;
974
                            }
975
                        }
976
                    }
977
                }
978
        }
979
    }                                /* end loop on i */
980
/*
981
 *  compute the square modulus of the structure factor, averaging on the surface
982
 *  kx*kx + ky*ky + kz*kz = krr*krr 
983
 *  note that this is correct only for a (on the macroscopic scale)
984
 *  isotropic system. 
985
 */
986
    for (i = 0; i < maxkx; i++) {
987
        kx = i * k_factor[XX]; for (j = 0; j < maxky; j++) {
988
            ky = j * k_factor[YY]; for (k = 0; k < maxkz; k++) {
989
                kz = k * k_factor[ZZ]; krr = sqrt (sqr (kx) + sqr (ky)
990
                + sqr (kz)); if (krr >= start_q && krr <= end_q) {
991
                    kr = (int) rint (krr / sf->ref_k); if (kr <
992
                    sf->n_angles && counter[kr] != 0)
993
                        sf->F[group][kr] +=
994
                            (sqr (tmpSF[i][j][k].re) +
995
                             sqr (tmpSF[i][j][k].im))/ counter[kr];
996
                }
997
            }
998
        }
999
    } sfree (counter); free(tmpSF[0][0]); free(tmpSF[0]); free(tmpSF);
1000
}
1001

    
1002
void save_data (structure_factor * sf, char *file, int ngrps, real start_q,
1003
           real end_q)
1004
{
1005

    
1006
    FILE *fp;
1007
    int i, g = 0;
1008
    double *tmp, polarization_factor, A;
1009

    
1010
    fp = xvgropen (file, "Scattering Intensity", "q (1/nm)",
1011
                   "Intensity (a.u.)");
1012

    
1013
    snew (tmp, ngrps);
1014

    
1015
    for (g = 0; g < ngrps; g++)
1016
        for (i = 0; i < sf->n_angles; i++) {
1017
/*
1018
 *          theta is half the angle between incoming and scattered vectors.
1019
 *          
1020
 *          polar. fact. = 0.5*(1+cos^2(2*theta)) = 1 - 0.5 * sin^2(2*theta)
1021
 *          
1022
 *          sin(theta) = q/(2k) := A  ->  sin^2(theta) = 4*A^2 (1-A^2) ->
1023
 *          -> 0.5*(1+cos^2(2*theta)) = 1 - 2 A^2 (1-A^2)
1024
 */
1025
            A = (double) (i * sf->ref_k) / (2.0 * sf->momentum);
1026
            polarization_factor = 1 - 2.0 * sqr (A) * (1 - sqr (A));
1027
            sf->F[g][i] *= polarization_factor;
1028
        }
1029
    for (i = 0; i < sf->n_angles; i++) {
1030
        if (i * sf->ref_k >= start_q && i * sf->ref_k <= end_q) {
1031
            fprintf (fp, "%10.5f  ", i * sf->ref_k);
1032
            for (g = 0; g < ngrps; g++)
1033
               fprintf (fp, "  %10.5f ", (sf->F[g][i]) /( sf->total_n_atoms*
1034
                                                          sf->nSteps));   
1035
            fprintf (fp, "\n");
1036
        }
1037
    }
1038
    ffclose (fp);
1039
}
1040

    
1041
int
1042
do_scattering_intensity (char* fnTPS, char* fnNDX, char* fnXVG, char *fnTRX,
1043
                         real start_q,real end_q, real energy,int ng)
1044
{
1045
    int i,*isize,status,flags = TRX_READ_X,**index_atp;
1046
    char **grpname,title[STRLEN];
1047
    atom_id **index;
1048
    t_topology top;
1049
    t_trxframe fr;
1050
    reduced_atom **red;
1051
    structure_factor *sf;
1052
    rvec *xtop;
1053
    real **sf_table;
1054
    int nsftable;
1055
    matrix box;
1056
    double r_tmp;
1057

    
1058
    snew (sf, 1);
1059
    sf->energy = energy;
1060
    /* Read the topology informations */
1061
    
1062
    read_tps_conf (fnTPS, title, &top, &xtop, NULL, box, TRUE);
1063
    sfree (xtop);
1064

    
1065
    /* groups stuff... */
1066
    snew (isize, ng);
1067
    snew (index, ng);
1068
    snew (grpname, ng);
1069

    
1070
    fprintf (stderr, "\nSelect %d group%s\n", ng,
1071
             ng == 1 ? "" : "s");
1072
    if (fnTPS)
1073
        get_index (&top.atoms, fnNDX, ng, isize, index, grpname);
1074
    else
1075
        rd_index (fnNDX, ng, isize, index, grpname);
1076

    
1077
    /* The first time we read data is a little special */
1078
    read_first_frame (&status, fnTRX, &fr, flags);
1079

    
1080
    sf->total_n_atoms = fr.natoms;
1081
    
1082
    snew (red, ng);
1083
    snew (index_atp, ng);
1084

    
1085
    r_tmp = max (box[XX][XX], box[YY][YY]);
1086
    r_tmp = (double) max (box[ZZ][ZZ], r_tmp);
1087

    
1088
    sf->ref_k = (2.0 * M_PI) / (r_tmp);
1089
    /* ref_k will be the reference momentum unit */
1090
    sf->n_angles = (int) rint (end_q / sf->ref_k);
1091

    
1092
    snew (sf->F, ng);
1093
    for (i = 0; i < ng; i++)
1094
        snew (sf->F[i], sf->n_angles);
1095
    for (i = 0; i < ng; i++) {
1096
        snew (red[i], isize[i]);
1097
        rearrange_atoms (red[i], &fr, index[i], isize[i], &top, TRUE);
1098
        index_atp[i] = create_indexed_atom_type (red[i], isize[i]);
1099
    }
1100
    sf_table = compute_scattering_factor_table (sf,&nsftable);
1101
    /* This is the main loop over frames */
1102

    
1103
    do {
1104
        sf->nSteps++;
1105
        for (i = 0; i < ng; i++) {
1106
            rearrange_atoms (red[i], &fr, index[i], isize[i], &top,FALSE);
1107

    
1108
            compute_structure_factor (sf, box, red[i], isize[i],
1109
                                      start_q, end_q, i, sf_table);
1110
        }
1111
    }
1112
    while (read_next_frame (status, &fr));
1113

    
1114
    save_data (sf, fnXVG, ng, start_q, end_q);
1115

    
1116
    return 0;
1117
}
1118

    
1119
int gmx_rdf(int argc,char *argv[])
1120
{
1121
  static char *desc[] = {
1122
    "The structure of liquids can be studied by either neutron or X-ray",
1123
    "scattering. The most common way to describe liquid structure is by a",
1124
    "radial distribution function. However, this is not easy to obtain from",
1125
    "a scattering experiment.[PAR]",
1126
    "g_rdf calculates radial distribution functions in different ways.",
1127
    "The normal method is around a (set of) particle(s), the other method",
1128
    "is around the center of mass of a set of particles.",
1129
    "With both methods rdf's can also be calculated around axes parallel",
1130
    "to the z-axis with option [TT]-xy[tt].[PAR]",
1131
    "If a run input file is supplied ([TT]-s[tt]), exclusions defined",
1132
    "in that file are taken into account when calculating the rdf.",
1133
    "The option [TT]-cut[tt] is meant as an alternative way to avoid",
1134
    "intramolecular peaks in the rdf plot.",
1135
    "It is however better to supply a run input file with a higher number of",
1136
    "exclusions. For eg. benzene a topology with nrexcl set to 5",
1137
    "would eliminate all intramolecular contributions to the rdf.",
1138
    "Note that all atoms in the selected groups are used, also the ones",
1139
    "that don't have Lennard-Jones interactions.[PAR]",
1140
    "Option [TT]-cn[tt] produces the cumulative number rdf.[PAR]"
1141
    "To bridge the gap between theory and experiment structure factors can",
1142
    "be computed (option [TT]-sq[tt]). The algorithm uses FFT, the grid"
1143
    "spacing of which is determined by option [TT]-grid[tt]."
1144
  };
1145
  static bool bCM=FALSE,bXY=FALSE,bPBC=TRUE;
1146
  static real cutoff=0,binwidth=0.002,grid=0.05,fade=0.0,lambda=0.1,distance=10;
1147
  static int  npixel=256,nlevel=20,ngroups=1;
1148
  static real start_q=0.0, end_q=60.0, energy=12.0;
1149
  t_pargs pa[] = {
1150
    { "-bin",      FALSE, etREAL, {&binwidth},
1151
      "Binwidth (nm)" },
1152
    { "-com",      FALSE, etBOOL, {&bCM},
1153
      "RDF with respect to the center of mass of first group" },
1154
    { "-pbc",      FALSE, etBOOL, {&bPBC},
1155
      "Use periodic boundary conditions for computing distances" },
1156
    { "-xy",       FALSE, etBOOL, {&bXY},
1157
      "Use only the x and y components of the distance" },
1158
    { "-cut",      FALSE, etREAL, {&cutoff},
1159
      "Shortest distance (nm) to be considered"},
1160
    { "-ng",       FALSE, etINT, {&ngroups},
1161
      "Number of secondary groups to compute RDFs around a central group" },
1162
    { "-fade",     FALSE, etREAL, {&fade},
1163
      "From this distance onwards the RDF is tranformed by g'(r) = 1 + [g(r)-1] exp(-(r/fade-1)^2 to make it go to 1 smoothly. If fade is 0.0 nothing is done." },
1164
    { "-grid",     FALSE, etREAL, {&grid},
1165
      "[HIDDEN]Grid spacing (in nm) for FFTs when computing structure factors" },
1166
    { "-npixel",   FALSE, etINT,  {&npixel},
1167
      "[HIDDEN]# pixels per edge of the square detector plate" },
1168
    { "-nlevel",   FALSE, etINT,  {&nlevel},
1169
      "Number of different colors in the diffraction image" },
1170
    { "-distance", FALSE, etREAL, {&distance},
1171
      "[HIDDEN]Distance (in cm) from the sample to the detector" },
1172
    { "-wave",     FALSE, etREAL, {&lambda},
1173
      "[HIDDEN]Wavelength for X-rays/Neutrons for scattering. 0.1 nm corresponds to roughly 12 keV" },
1174
    
1175
    {"-startq", FALSE, etREAL, {&start_q},
1176
     "Starting q (1/nm) "},
1177
    {"-endq", FALSE, etREAL, {&end_q},
1178
     "Ending q (1/nm)"},
1179
    {"-energy", FALSE, etREAL, {&energy},
1180
     "Energy of the incoming X-ray (keV) "}
1181
  };
1182
#define NPA asize(pa)
1183
  char       *fnTPS,*fnNDX;
1184
  bool       bSQ,bRDF;
1185
  
1186
  t_filenm   fnm[] = {
1187
    { efTRX, "-f",  NULL,     ffREAD },
1188
    { efTPS, NULL,  NULL,     ffOPTRD },
1189
    { efNDX, NULL,  NULL,     ffOPTRD },
1190
    { efXVG, "-o",  "rdf",    ffOPTWR },
1191
    { efXVG, "-sq", "sq",     ffOPTWR },
1192
    { efXVG, "-cn", "rdf_cn", ffOPTWR },
1193
    { efXVG, "-hq", "hq",     ffOPTWR },
1194
/*    { efXPM, "-image", "sq",  ffOPTWR }*/
1195
  };
1196
#define NFILE asize(fnm)
1197
  
1198
  CopyRight(stderr,argv[0]);
1199
  parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
1200
                    NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL);
1201

    
1202
  fnTPS = ftp2fn_null(efTPS,NFILE,fnm);
1203
  fnNDX = ftp2fn_null(efNDX,NFILE,fnm);
1204
  /*  bSQ   = opt2bSet("-sq",NFILE,fnm) || opt2parg_bSet("-grid",NPA,pa); */
1205
  bSQ   = opt2bSet("-sq",NFILE,fnm);
1206
  bRDF  = opt2bSet("-o",NFILE,fnm) || !bSQ;
1207
  
1208
  if (bSQ) {
1209
    if (!fnTPS)
1210
      gmx_fatal(FARGS,"Need a tps file for calculating structure factors\n");
1211
  }
1212
  else {
1213
    if (!fnTPS && !fnNDX)
1214
      gmx_fatal(FARGS,"Neither index file nor topology file specified\n"
1215
                  "             Nothing to do!");
1216
  }
1217
 
1218
  if  (bSQ) 
1219
   do_scattering_intensity(fnTPS,fnNDX,opt2fn("-sq",NFILE,fnm),ftp2fn(efTRX,NFILE,fnm),
1220
                           start_q, end_q, energy, ngroups  );
1221
/* old structure factor code */
1222
/*    do_sq(fnNDX,fnTPS,ftp2fn(efTRX,NFILE,fnm),opt2fn("-sq",NFILE,fnm),
1223
          ftp2fn(efXPM,NFILE,fnm),grid,lambda,distance,npixel,nlevel);
1224
*/
1225
  if (bRDF) 
1226
    do_rdf(fnNDX,fnTPS,ftp2fn(efTRX,NFILE,fnm),
1227
           opt2fn("-o",NFILE,fnm),opt2fn_null("-cn",NFILE,fnm),
1228
           opt2fn_null("-hq",NFILE,fnm),
1229
           bCM,bXY,bPBC,cutoff,binwidth,fade,ngroups);
1230

    
1231
  thanx(stderr);
1232
  
1233
  return 0;
1234
}