3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
66 static void check_box_c(matrix box
)
68 if (fabs(box
[ZZ
][XX
]) > GMX_REAL_EPS
*box
[ZZ
][ZZ
] ||
69 fabs(box
[ZZ
][YY
]) > GMX_REAL_EPS
*box
[ZZ
][ZZ
])
71 "The last box vector is not parallel to the z-axis: %f %f %f",
72 box
[ZZ
][XX
],box
[ZZ
][YY
],box
[ZZ
][ZZ
]);
75 static void calc_comg(int is
,int *coi
,int *index
,gmx_bool bMass
,t_atom
*atom
,
82 if (bMass
&& atom
==NULL
)
83 gmx_fatal(FARGS
,"No masses available while mass weighting was requested");
88 for(i
=coi
[c
]; i
<coi
[c
+1]; i
++) {
92 xc
[d
] += m
*x
[index
[i
]][d
];
95 rvec_inc(xc
,x
[index
[i
]]);
99 svmul(1/mtot
,xc
,x_comg
[c
]);
103 static void split_group(int isize
,int *index
,char *grpname
,
104 t_topology
*top
,char type
,
105 int *is_out
,int **coi_out
)
110 int cur
,mol
,res
,i
,a
,i1
;
112 /* Split up the group in molecules or residues */
118 atom
= top
->atoms
.atom
;
121 gmx_fatal(FARGS
,"Unknown rdf option '%s'",type
);
127 for(i
=0; i
<isize
; i
++) {
130 /* Check if the molecule number has changed */
131 i1
= mols
->index
[mol
+1];
134 i1
= mols
->index
[mol
+1];
140 } else if (type
== 'r') {
141 /* Check if the residue index has changed */
142 res
= atom
[a
].resind
;
151 printf("Group '%s' of %d atoms consists of %d %s\n",
153 (type
=='m' ? "molecules" : "residues"));
159 static void do_rdf(const char *fnNDX
,const char *fnTPS
,const char *fnTRX
,
160 const char *fnRDF
,const char *fnCNRDF
, const char *fnHQ
,
161 gmx_bool bCM
,const char *close
,
162 const char **rdft
,gmx_bool bXY
,gmx_bool bPBC
,gmx_bool bNormalize
,
163 real cutoff
,real binwidth
,real fade
,int ng
,
164 const output_env_t oenv
)
168 char outf1
[STRLEN
],outf2
[STRLEN
];
169 char title
[STRLEN
],gtitle
[STRLEN
],refgt
[30];
170 int g
,natoms
,i
,ii
,j
,k
,nbin
,j0
,j1
,n
,nframes
;
173 int *isize
,isize_cm
=0,nrdf
=0,max_i
,isize0
,isize_g
;
174 atom_id
**index
,*index_cm
=NULL
;
175 #if (defined SIZEOF_LONG_LONG_INT) && (SIZEOF_LONG_LONG_INT >= 8)
180 real t
,rmax2
,cut2
,r
,r2
,r2ii
,invhbinw
,normfac
;
181 real segvol
,spherevol
,prev_spherevol
,**rdf
;
182 rvec
*x
,dx
,*x0
=NULL
,*x_i1
,xi
;
183 real
*inv_segvol
,invvol
,invvol_sum
,rho
;
184 gmx_bool bClose
,*bExcl
,bTop
,bNonSelfExcl
;
187 atom_id ix
,jx
,***pairs
;
188 t_topology
*top
=NULL
;
189 int ePBC
=-1,ePBCrdf
=-1;
194 gmx_rmpbc_t gpbc
=NULL
;
195 int *is
=NULL
,**coi
=NULL
,cur
,mol
,i1
,res
,a
;
199 bClose
= (close
[0] != 'n');
203 bTop
=read_tps_conf(fnTPS
,title
,top
,&ePBC
,&x
,NULL
,box
,TRUE
);
204 if (bTop
&& !(bCM
|| bClose
))
205 /* get exclusions from topology */
206 excl
= &(top
->excls
);
211 fprintf(stderr
,"\nSelect a reference group and %d group%s\n",
214 get_index(&(top
->atoms
),fnNDX
,ng
+1,isize
,index
,grpname
);
215 atom
= top
->atoms
.atom
;
217 rd_index(fnNDX
,ng
+1,isize
,index
,grpname
);
224 split_group(isize
[0],index
[0],grpname
[0],top
,close
[0],&is
[0],&coi
[0]);
227 if (rdft
[0][0] != 'a') {
228 /* Split up all the groups in molecules or residues */
231 for(g
=((bCM
|| bClose
) ? 1 : 0); g
<ng
+1; g
++) {
232 split_group(isize
[g
],index
[g
],grpname
[g
],top
,rdft
[0][0],&is
[g
],&coi
[g
]);
238 snew(coi
[0],is
[0]+1);
240 coi
[0][1] = isize
[0];
243 } else if (bClose
|| rdft
[0][0] != 'a') {
250 natoms
=read_first_x(oenv
,&status
,fnTRX
,&t
,&x
,box
);
252 gmx_fatal(FARGS
,"Could not read coordinates from statusfile\n");
254 /* check with topology */
255 if ( natoms
> top
->atoms
.nr
)
256 gmx_fatal(FARGS
,"Trajectory (%d atoms) does not match topology (%d atoms)",
257 natoms
,top
->atoms
.nr
);
258 /* check with index groups */
259 for (i
=0; i
<ng
+1; i
++)
260 for (j
=0; j
<isize
[i
]; j
++)
261 if ( index
[i
][j
] >= natoms
)
262 gmx_fatal(FARGS
,"Atom index (%d) in index group %s (%d atoms) larger "
263 "than number of atoms in trajectory (%d atoms)",
264 index
[i
][j
],grpname
[i
],isize
[i
],natoms
);
266 /* initialize some handy things */
268 ePBC
= guess_ePBC(box
);
270 copy_mat(box
,box_pbc
);
275 case epbcXY
: ePBCrdf
= epbcXY
; break;
276 case epbcNONE
: ePBCrdf
= epbcNONE
; break;
278 gmx_fatal(FARGS
,"xy-rdf's are not supported for pbc type'%s'",
282 /* Make sure the z-height does not influence the cut-off */
283 box_pbc
[ZZ
][ZZ
] = 2*max(box
[XX
][XX
],box
[YY
][YY
]);
288 rmax2
= 0.99*0.99*max_cutoff2(bXY
? epbcXY
: epbcXYZ
,box_pbc
);
290 rmax2
= sqr(3*max(box
[XX
][XX
],max(box
[YY
][YY
],box
[ZZ
][ZZ
])));
292 fprintf(debug
,"rmax2 = %g\n",rmax2
);
294 /* We use the double amount of bins, so we can correctly
295 * write the rdf and rdf_cn output at i*binwidth values.
297 nbin
= (int)(sqrt(rmax2
) * 2 / binwidth
);
298 invhbinw
= 2.0 / binwidth
;
307 for(g
=0; g
<ng
; g
++) {
308 if (isize
[g
+1] > max_i
)
311 /* this is THE array */
312 snew(count
[g
],nbin
+1);
314 /* make pairlist array for groups and exclusions */
315 snew(pairs
[g
],isize
[0]);
316 snew(npairs
[g
],isize
[0]);
317 for(i
=0; i
<isize
[0]; i
++) {
318 /* We can only have exclusions with atomic rdfs */
319 if (!(bCM
|| bClose
|| rdft
[0][0] != 'a')) {
321 for(j
=0; j
< natoms
; j
++)
325 for( j
= excl
->index
[ix
]; j
< excl
->index
[ix
+1]; j
++)
326 bExcl
[excl
->a
[j
]]=TRUE
;
328 snew(pairs
[g
][i
], isize
[g
+1]);
329 bNonSelfExcl
= FALSE
;
330 for(j
=0; j
<isize
[g
+1]; j
++) {
335 /* Check if we have exclusions other than self exclusions */
340 srenew(pairs
[g
][i
],npairs
[g
][i
]);
342 /* Save a LOT of memory and some cpu cycles */
356 if (bPBC
&& (NULL
!= top
))
357 gpbc
= gmx_rmpbc_init(&top
->idef
,ePBC
,natoms
,box
);
360 /* Must init pbc every step because of pressure coupling */
361 copy_mat(box
,box_pbc
);
364 gmx_rmpbc(gpbc
,natoms
,box
,x
);
367 clear_rvec(box_pbc
[ZZ
]);
369 set_pbc(&pbc
,ePBCrdf
,box_pbc
);
372 /* Set z-size to 1 so we get the surface iso the volume */
375 invvol
= 1/det(box_pbc
);
376 invvol_sum
+= invvol
;
379 /* Calculate center of mass of the whole group */
380 calc_comg(is
[0],coi
[0],index
[0],TRUE
,atom
,x
,x0
);
381 } else if (!bClose
&& rdft
[0][0] != 'a') {
382 calc_comg(is
[0],coi
[0],index
[0],rdft
[0][6]=='m',atom
,x
,x0
);
385 for(g
=0; g
<ng
; g
++) {
386 if (rdft
[0][0] == 'a') {
387 /* Copy the indexed coordinates to a continuous array */
388 for(i
=0; i
<isize
[g
+1]; i
++)
389 copy_rvec(x
[index
[g
+1][i
]],x_i1
[i
]);
391 /* Calculate the COMs/COGs and store in x_i1 */
392 calc_comg(is
[g
+1],coi
[g
+1],index
[g
+1],rdft
[0][6]=='m',atom
,x
,x_i1
);
395 for(i
=0; i
<isize0
; i
++) {
397 /* Special loop, since we need to determine the minimum distance
398 * over all selected atoms in the reference molecule/residue.
400 if (rdft
[0][0] == 'a')
401 isize_g
= isize
[g
+1];
404 for(j
=0; j
<isize_g
; j
++) {
406 /* Loop over the selected atoms in the reference molecule */
407 for(ii
=coi
[0][i
]; ii
<coi
[0][i
+1]; ii
++) {
409 pbc_dx(&pbc
,x
[index
[0][ii
]],x_i1
[j
],dx
);
411 rvec_sub(x
[index
[0][ii
]],x_i1
[j
],dx
);
413 r2ii
= dx
[XX
]*dx
[XX
] + dx
[YY
]*dx
[YY
];
419 if (r2
>cut2
&& r2
<=rmax2
)
420 count
[g
][(int)(sqrt(r2
)*invhbinw
)]++;
423 /* Real rdf between points in space */
424 if (bCM
|| rdft
[0][0] != 'a') {
427 copy_rvec(x
[index
[0][i
]],xi
);
429 if (rdft
[0][0] == 'a' && npairs
[g
][i
] >= 0) {
430 /* Expensive loop, because of indexing */
431 for(j
=0; j
<npairs
[g
][i
]; j
++) {
434 pbc_dx(&pbc
,xi
,x
[jx
],dx
);
436 rvec_sub(xi
,x
[jx
],dx
);
439 r2
= dx
[XX
]*dx
[XX
] + dx
[YY
]*dx
[YY
];
442 if (r2
>cut2
&& r2
<=rmax2
)
443 count
[g
][(int)(sqrt(r2
)*invhbinw
)]++;
446 /* Cheaper loop, no exclusions */
447 if (rdft
[0][0] == 'a')
448 isize_g
= isize
[g
+1];
451 for(j
=0; j
<isize_g
; j
++) {
453 pbc_dx(&pbc
,xi
,x_i1
[j
],dx
);
455 rvec_sub(xi
,x_i1
[j
],dx
);
457 r2
= dx
[XX
]*dx
[XX
] + dx
[YY
]*dx
[YY
];
460 if (r2
>cut2
&& r2
<=rmax2
)
461 count
[g
][(int)(sqrt(r2
)*invhbinw
)]++;
468 } while (read_next_x(oenv
,status
,&t
,natoms
,x
,box
));
469 fprintf(stderr
,"\n");
471 if (bPBC
&& (NULL
!= top
))
472 gmx_rmpbc_done(gpbc
);
479 invvol
= invvol_sum
/nframes
;
481 /* Calculate volume of sphere segments or length of circle segments */
482 snew(inv_segvol
,(nbin
+1)/2);
484 for(i
=0; (i
<(nbin
+1)/2); i
++) {
485 r
= (i
+ 0.5)*binwidth
;
489 spherevol
=(4.0/3.0)*M_PI
*r
*r
*r
;
491 segvol
=spherevol
-prev_spherevol
;
492 inv_segvol
[i
]=1.0/segvol
;
493 prev_spherevol
=spherevol
;
497 for(g
=0; g
<ng
; g
++) {
498 /* We have to normalize by dividing by the number of frames */
499 if (rdft
[0][0] == 'a')
500 normfac
= 1.0/(nframes
*invvol
*isize0
*isize
[g
+1]);
502 normfac
= 1.0/(nframes
*invvol
*isize0
*is
[g
+1]);
504 /* Do the normalization */
505 nrdf
= max((nbin
+1)/2,1+2*fade
/binwidth
);
507 for(i
=0; i
<(nbin
+1)/2; i
++) {
512 j
= count
[g
][i
*2-1] + count
[g
][i
*2];
513 if ((fade
> 0) && (r
>= fade
))
514 rdf
[g
][i
] = 1 + (j
*inv_segvol
[i
]*normfac
-1)*exp(-16*sqr(r
/fade
-1));
517 rdf
[g
][i
] = j
*inv_segvol
[i
]*normfac
;
519 rdf
[g
][i
] = j
/(binwidth
*isize0
*nframes
);
522 for( ; (i
<nrdf
); i
++)
526 if (rdft
[0][0] == 'a') {
527 sprintf(gtitle
,"Radial distribution");
529 sprintf(gtitle
,"Radial distribution of %s %s",
530 rdft
[0][0]=='m' ? "molecule" : "residue",
531 rdft
[0][6]=='m' ? "COM" : "COG");
533 fp
=xvgropen(fnRDF
,gtitle
,"r","",oenv
);
535 sprintf(refgt
," %s","COM");
537 sprintf(refgt
," closest atom in %s.",close
);
539 sprintf(refgt
,"%s","");
542 if (output_env_get_print_xvgr_codes(oenv
))
543 fprintf(fp
,"@ subtitle \"%s%s - %s\"\n",grpname
[0],refgt
,grpname
[1]);
546 if (output_env_get_print_xvgr_codes(oenv
))
547 fprintf(fp
,"@ subtitle \"reference %s%s\"\n",grpname
[0],refgt
);
548 xvgr_legend(fp
,ng
,(const char**)(grpname
+1),oenv
);
550 for(i
=0; (i
<nrdf
); i
++) {
551 fprintf(fp
,"%10g",i
*binwidth
);
553 fprintf(fp
," %10g",rdf
[g
][i
]);
558 do_view(oenv
,fnRDF
,NULL
);
560 /* h(Q) function: fourier transform of rdf */
563 real
*hq
,*integrand
,Q
;
565 /* Get a better number density later! */
566 rho
= isize
[1]*invvol
;
568 snew(integrand
,nrdf
);
569 for(i
=0; (i
<nhq
); i
++) {
572 for(j
=1; (j
<nrdf
); j
++) {
574 integrand
[j
] = (Q
== 0) ? 1.0 : sin(Q
*r
)/(Q
*r
);
575 integrand
[j
] *= 4.0*M_PI
*rho
*r
*r
*(rdf
[0][j
]-1.0);
577 hq
[i
] = print_and_integrate(debug
,nrdf
,binwidth
,integrand
,NULL
,0);
579 fp
=xvgropen(fnHQ
,"h(Q)","Q(/nm)","h(Q)",oenv
);
580 for(i
=0; (i
<nhq
); i
++)
581 fprintf(fp
,"%10g %10g\n",i
*0.5,hq
[i
]);
583 do_view(oenv
,fnHQ
,NULL
);
589 normfac
= 1.0/(isize0
*nframes
);
590 fp
=xvgropen(fnCNRDF
,"Cumulative Number RDF","r","number",oenv
);
592 if (output_env_get_print_xvgr_codes(oenv
))
593 fprintf(fp
,"@ subtitle \"%s-%s\"\n",grpname
[0],grpname
[1]);
596 if (output_env_get_print_xvgr_codes(oenv
))
597 fprintf(fp
,"@ subtitle \"reference %s\"\n",grpname
[0]);
598 xvgr_legend(fp
,ng
,(const char**)(grpname
+1),oenv
);
601 for(i
=0; (i
<=nbin
/2); i
++) {
602 fprintf(fp
,"%10g",i
*binwidth
);
603 for(g
=0; g
<ng
; g
++) {
604 fprintf(fp
," %10g",(real
)((double)sum
[g
]*normfac
));
606 sum
[g
] += count
[g
][i
*2] + count
[g
][i
*2+1];
613 do_view(oenv
,fnCNRDF
,NULL
);
622 int gmx_rdf(int argc
,char *argv
[])
624 const char *desc
[] = {
625 "The structure of liquids can be studied by either neutron or X-ray",
626 "scattering. The most common way to describe liquid structure is by a",
627 "radial distribution function. However, this is not easy to obtain from",
628 "a scattering experiment.[PAR]",
629 "g_rdf calculates radial distribution functions in different ways.",
630 "The normal method is around a (set of) particle(s), the other methods",
631 "are around the center of mass of a set of particles ([TT]-com[tt])",
632 "or to the closest particle in a set ([TT]-surf[tt]).",
633 "With all methods rdf's can also be calculated around axes parallel",
634 "to the z-axis with option [TT]-xy[tt].",
635 "With option [TT]-surf[tt] normalization can not be used.[PAR]",
636 "The option [TT]-rdf[tt] sets the type of rdf to be computed.",
637 "Default is for atoms or particles, but one can also select center",
638 "of mass or geometry of molecules or residues. In all cases only",
639 "the atoms in the index groups are taken into account.",
640 "For molecules and/or the center of mass option a run input file",
642 "Other weighting than COM or COG can currently only be achieved",
643 "by providing a run input file with different masses.",
644 "Options [TT]-com[tt] and [TT]-surf[tt] also work in conjunction",
645 "with [TT]-rdf[tt].[PAR]",
646 "If a run input file is supplied ([TT]-s[tt]) and [TT]-rdf[tt] is set",
647 "to [TT]atom[tt], exclusions defined",
648 "in that file are taken into account when calculating the rdf.",
649 "The option [TT]-cut[tt] is meant as an alternative way to avoid",
650 "intramolecular peaks in the rdf plot.",
651 "It is however better to supply a run input file with a higher number of",
652 "exclusions. For eg. benzene a topology with nrexcl set to 5",
653 "would eliminate all intramolecular contributions to the rdf.",
654 "Note that all atoms in the selected groups are used, also the ones",
655 "that don't have Lennard-Jones interactions.[PAR]",
656 "Option [TT]-cn[tt] produces the cumulative number rdf,",
657 "i.e. the average number of particles within a distance r.[PAR]",
658 "To bridge the gap between theory and experiment structure factors can",
659 "be computed (option [TT]-sq[tt]). The algorithm uses FFT, the grid",
660 "spacing of which is determined by option [TT]-grid[tt]."
662 static gmx_bool bCM
=FALSE
,bXY
=FALSE
,bPBC
=TRUE
,bNormalize
=TRUE
;
663 static real cutoff
=0,binwidth
=0.002,grid
=0.05,fade
=0.0,lambda
=0.1,distance
=10;
664 static int npixel
=256,nlevel
=20,ngroups
=1;
665 static real start_q
=0.0, end_q
=60.0, energy
=12.0;
667 static const char *closet
[]= { NULL
, "no", "mol", "res", NULL
};
668 static const char *rdft
[]={ NULL
, "atom", "mol_com", "mol_cog", "res_com", "res_cog", NULL
};
671 { "-bin", FALSE
, etREAL
, {&binwidth
},
673 { "-com", FALSE
, etBOOL
, {&bCM
},
674 "RDF with respect to the center of mass of first group" },
675 { "-surf", FALSE
, etENUM
, {closet
},
676 "RDF with respect to the surface of the first group" },
677 { "-rdf", FALSE
, etENUM
, {rdft
},
679 { "-pbc", FALSE
, etBOOL
, {&bPBC
},
680 "Use periodic boundary conditions for computing distances. Without PBC the maximum range will be three times the largest box edge." },
681 { "-norm", FALSE
, etBOOL
, {&bNormalize
},
682 "Normalize for volume and density" },
683 { "-xy", FALSE
, etBOOL
, {&bXY
},
684 "Use only the x and y components of the distance" },
685 { "-cut", FALSE
, etREAL
, {&cutoff
},
686 "Shortest distance (nm) to be considered"},
687 { "-ng", FALSE
, etINT
, {&ngroups
},
688 "Number of secondary groups to compute RDFs around a central group" },
689 { "-fade", FALSE
, etREAL
, {&fade
},
690 "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." },
691 { "-grid", FALSE
, etREAL
, {&grid
},
692 "[HIDDEN]Grid spacing (in nm) for FFTs when computing structure factors" },
693 { "-npixel", FALSE
, etINT
, {&npixel
},
694 "[HIDDEN]# pixels per edge of the square detector plate" },
695 { "-nlevel", FALSE
, etINT
, {&nlevel
},
696 "Number of different colors in the diffraction image" },
697 { "-distance", FALSE
, etREAL
, {&distance
},
698 "[HIDDEN]Distance (in cm) from the sample to the detector" },
699 { "-wave", FALSE
, etREAL
, {&lambda
},
700 "[HIDDEN]Wavelength for X-rays/Neutrons for scattering. 0.1 nm corresponds to roughly 12 keV" },
702 {"-startq", FALSE
, etREAL
, {&start_q
},
703 "Starting q (1/nm) "},
704 {"-endq", FALSE
, etREAL
, {&end_q
},
706 {"-energy", FALSE
, etREAL
, {&energy
},
707 "Energy of the incoming X-ray (keV) "}
709 #define NPA asize(pa)
710 const char *fnTPS
,*fnNDX
,*fnDAT
=NULL
;
715 { efTRX
, "-f", NULL
, ffREAD
},
716 { efTPS
, NULL
, NULL
, ffOPTRD
},
717 { efNDX
, NULL
, NULL
, ffOPTRD
},
718 { efDAT
, "-d", "sfactor", ffOPTRD
},
719 { efXVG
, "-o", "rdf", ffOPTWR
},
720 { efXVG
, "-sq", "sq", ffOPTWR
},
721 { efXVG
, "-cn", "rdf_cn", ffOPTWR
},
722 { efXVG
, "-hq", "hq", ffOPTWR
},
723 /* { efXPM, "-image", "sq", ffOPTWR }*/
725 #define NFILE asize(fnm)
727 CopyRight(stderr
,argv
[0]);
728 parse_common_args(&argc
,argv
,PCA_CAN_VIEW
| PCA_CAN_TIME
| PCA_BE_NICE
,
729 NFILE
,fnm
,NPA
,pa
,asize(desc
),desc
,0,NULL
,&oenv
);
731 bSQ
= opt2bSet("-sq",NFILE
,fnm
);
733 please_cite(stdout
,"Cromer1968a");
735 bRDF
= opt2bSet("-o",NFILE
,fnm
) || !bSQ
;
736 if (bSQ
|| bCM
|| closet
[0][0]!='n' || rdft
[0][0]=='m' || rdft
[0][6]=='m') {
737 fnTPS
= ftp2fn(efTPS
,NFILE
,fnm
);
738 fnDAT
= ftp2fn(efDAT
,NFILE
,fnm
);
740 fnTPS
= ftp2fn_null(efTPS
,NFILE
,fnm
);
742 fnNDX
= ftp2fn_null(efNDX
,NFILE
,fnm
);
744 if (closet
[0][0] != 'n') {
746 gmx_fatal(FARGS
,"Can not have both -com and -surf");
749 fprintf(stderr
,"Turning of normalization because of option -surf\n");
754 if (!bSQ
&& (!fnTPS
&& !fnNDX
))
755 gmx_fatal(FARGS
,"Neither index file nor topology file specified\n"
759 do_scattering_intensity(fnTPS
,fnNDX
,opt2fn("-sq",NFILE
,fnm
),
760 ftp2fn(efTRX
,NFILE
,fnm
),fnDAT
,
761 start_q
, end_q
, energy
, ngroups
,oenv
);
764 do_rdf(fnNDX
,fnTPS
,ftp2fn(efTRX
,NFILE
,fnm
),
765 opt2fn("-o",NFILE
,fnm
),opt2fn_null("-cn",NFILE
,fnm
),
766 opt2fn_null("-hq",NFILE
,fnm
),
767 bCM
,closet
[0],rdft
,bXY
,bPBC
,bNormalize
,cutoff
,binwidth
,fade
,ngroups
,