4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
31 * For more info, check our website at http://www.gromacs.org
34 * Green Red Orange Magenta Azure Cyan Skyblue
55 #include "shift_util.h"
68 * f0[k] = c + [SUM a_i*EXP(-b_i*(k^2)) ]
72 const t_CM_table CM_t
[] =
75 { "H", { 0.489918, 0.262003, 0.196767, 0.049879 },
76 { 20.6593, 7.74039, 49.5519, 2.20159 },
78 { "HO", { 0.489918, 0.262003, 0.196767, 0.049879 },
79 { 20.6593, 7.74039, 49.5519, 2.20159 },
81 { "HW", { 0.489918, 0.262003, 0.196767, 0.049879 },
82 { 20.6593, 7.74039, 49.5519, 2.20159 },
84 { "C", { 2.26069, 1.56165, 1.05075, 0.839259 },
85 { 22.6907, 0.656665, 9.75618, 55.5949 },
87 { "CB", { 2.26069, 1.56165, 1.05075, 0.839259 },
88 { 22.6907, 0.656665, 9.75618, 55.5949 },
90 { "CS1", { 0., 0., 0., 0.},{ 0., 0., 0., 0.}, 0.},
91 { "CS2", { 0., 0., 0., 0.},{ 0., 0., 0., 0.}, 0.},
92 { "N", { 12.2126, 3.13220, 2.01250, 1.16630 },
93 { 0.005700, 9.89330, 28.9975, 0.582600 },
95 { "O", { 3.04850, 2.28680, 1.54630, 0.867000 },
96 { 13.2771, 5.70110, 0.323900, 32.9089 },
98 { "OW", { 3.04850, 2.28680, 1.54630, 0.867000 },
99 { 13.2771, 5.70110, 0.323900, 32.9089 },
101 { "OWT3", { 3.04850, 2.28680, 1.54630, 0.867000 },
102 { 13.2771, 5.70110, 0.323900, 32.9089 },
104 { "OA", { 3.04850, 2.28680, 1.54630, 0.867000 },
105 { 13.2771, 5.70110, 0.323900, 32.9089 },
107 { "OS", { 3.04850, 2.28680, 1.54630, 0.867000 },
108 { 13.2771, 5.70110, 0.323900, 32.9089 },
110 { "OSE", { 3.04850, 2.28680, 1.54630, 0.867000 },
111 { 13.2771, 5.70110, 0.323900, 32.9089 },
113 { "Na", { 3.25650, 3.93620, 1.39980, 1.00320 }, /* Na 1+ */
114 { 2.66710, 6.11530, 0.200100, 14.0390 },
116 { "CH3", { 0., 0., 0., 0.},{ 0., 0., 0., 0.}, 0.},
117 { "CH2", { 0., 0., 0., 0.},{ 0., 0., 0., 0.}, 0.},
118 { "CH1", { 0., 0., 0., 0.},{ 0., 0., 0., 0.}, 0.},
143 static void do_rdf(char *fnNDX
,char *fnTPS
,char *fnTRX
,
144 char *fnRDF
,char *fnCNRDF
, char *fnHQ
,
145 bool bCM
,real cutoff
,real binwidth
,real fade
)
149 char outf1
[STRLEN
],outf2
[STRLEN
];
151 int g
,ng
,natoms
,i
,j
,k
,nbin
,j0
,j1
,n
,nframes
;
154 int *isize
,isize_cm
=0,nrdf
=0,max_i
;
155 atom_id
**index
,*index_cm
=NULL
;
156 unsigned long int *sum
;
157 real t
,boxmin
,hbox
,hbox2
,cut2
,r
,r2
,invbinw
,normfac
;
158 real segvol
,spherevol
,prev_spherevol
,**rdf
;
159 rvec
*x
,xcom
,dx
,*x_i1
,xi
;
160 real
*inv_segvol
,vol
,vol_sum
,rho
;
161 bool *bExcl
,bTop
,bNonSelfExcl
;
164 atom_id ix
,jx
,***pairs
;
170 bTop
=read_tps_conf(fnTPS
,title
,&top
,&x
,NULL
,box
,TRUE
);
173 /* get exclusions from topology */
174 excl
=&(top
.atoms
.excl
);
176 fprintf(stderr
,"\nHow many groups do you want to calculate the RDF of?\n");
183 fprintf(stderr
,"\nSelect a reference group and %d group%s\n",
186 get_index(&top
.atoms
,fnNDX
,ng
+1,isize
,index
,grpname
);
188 rd_index(fnNDX
,ng
+1,isize
,index
,grpname
);
190 natoms
=read_first_x(&status
,fnTRX
,&t
,&x
,box
);
192 fatal_error(0,"Could not read coordinates from statusfile\n");
194 /* check with topology */
195 if ( natoms
> top
.atoms
.nr
)
196 fatal_error(0,"Trajectory (%d atoms) does not match topology (%d atoms)",
197 natoms
,top
.atoms
.nr
);
198 /* check with index groups */
199 for (i
=0; i
<=ng
; i
++)
200 for (j
=0; j
<isize
[i
]; j
++)
201 if ( index
[i
][j
] >= natoms
)
202 fatal_error(0,"Atom index (%d) in index group %s (%d atoms) larger "
203 "than number of atoms in trajectory (%d atoms)",
204 index
[i
][j
],grpname
[i
],isize
[i
],natoms
);
207 /* move index[0] to index_cm and make 'dummy' index[0] */
209 snew(index_cm
,isize_cm
);
210 for(i
=0; i
<isize
[0]; i
++)
211 index_cm
[i
]=index
[0][i
];
214 srenew(index
[0],isize
[0]);
215 /* make space for center of mass */
219 /* initialize some handy things */
220 boxmin
= min( norm(box
[XX
]), min( norm(box
[YY
]), norm(box
[ZZ
]) ) );
222 nbin
= (int)(hbox
/ binwidth
) + 1;
223 invbinw
= 1.0 / binwidth
;
233 for(g
=0; g
<ng
; g
++) {
234 if (isize
[g
+1] > max_i
)
237 /* this is THE array */
238 snew(count
[g
],nbin
+1);
240 /* make pairlist array for groups and exclusions */
241 snew(pairs
[g
],isize
[0]);
242 snew(npairs
[g
],isize
[0]);
243 for(i
=0; i
<isize
[0]; i
++) {
245 for(j
=0; j
< natoms
; j
++)
249 for( j
= excl
->index
[ix
]; j
< excl
->index
[ix
+1]; j
++)
250 bExcl
[excl
->a
[j
]]=TRUE
;
252 snew(pairs
[g
][i
], isize
[g
+1]);
253 bNonSelfExcl
= FALSE
;
254 for(j
=0; j
<isize
[g
+1]; j
++) {
259 /* Check if we have exclusions other than self exclusions */
260 bNonSelfExcl
= bNonSelfExcl
|| (ix
!= jx
);
264 srenew(pairs
[g
][i
],npairs
[g
][i
]);
266 /* Save a LOT of memory and some cpu cycles */
278 /* Must init pbc every step because of pressure coupling */
280 rm_pbc(&top
.idef
,natoms
,box
,x
,x
);
286 /* calculate centre of mass */
288 for(i
=0; (i
< isize_cm
); i
++) {
290 rvec_inc(xcom
,x
[ix
]);
292 /* store it in the first 'group' */
293 for(j
=0; (j
<DIM
); j
++)
294 x
[index
[0][0]][j
] = xcom
[j
] / isize_cm
;
297 for(g
=0; g
<ng
; g
++) {
298 /* Copy the indexed coordinates to a continuous array */
299 for(i
=0; i
<isize
[g
+1]; i
++)
300 copy_rvec(x
[index
[g
+1][i
]],x_i1
[i
]);
302 for(i
=0; i
<isize
[0]; i
++) {
303 copy_rvec(x
[index
[0][i
]],xi
);
304 if (npairs
[g
][i
] >= 0)
305 /* Expensive loop, because of indexing */
306 for(j
=0; j
<npairs
[g
][i
]; j
++) {
310 if (r2
>cut2
&& r2
<=hbox2
)
311 count
[g
][(int)(sqrt(r2
)*invbinw
)]++;
314 /* Cheaper loop, no exclusions */
315 for(j
=0; j
<isize
[g
+1]; j
++) {
316 pbc_dx(xi
,x_i1
[j
],dx
);
318 if (r2
>cut2
&& r2
<=hbox2
)
319 count
[g
][(int)(sqrt(r2
)*invbinw
)]++;
324 } while (read_next_x(status
,&t
,natoms
,x
,box
));
325 fprintf(stderr
,"\n");
332 vol
= vol_sum
/nframes
;
334 /* Calculate volume of sphere segments */
335 snew(inv_segvol
,nbin
);
337 for(i
=0; (i
<nbin
); i
++) {
339 spherevol
=(4.0/3.0)*M_PI
*r
*r
*r
;
340 segvol
=spherevol
-prev_spherevol
;
341 inv_segvol
[i
]=1.0/segvol
;
342 prev_spherevol
=spherevol
;
346 for(g
=0; g
<ng
; g
++) {
347 /* We have to normalize by dividing by the number of frames */
348 rho
= isize
[g
+1]/vol
;
349 normfac
= 1.0/((rho
*nframes
)*isize
[0]);
351 /* Do the normalization */
352 nrdf
= max(nbin
-1,1+(2*fade
/binwidth
));
354 for(i
=0; (i
<nbin
-1); i
++) {
355 r
= (i
+0.5)*binwidth
;
356 if ((fade
> 0) && (r
>= fade
))
357 rdf
[g
][i
] = 1+(count
[g
][i
]*inv_segvol
[i
]*normfac
-1)*exp(-16*sqr(r
/fade
-1));
359 rdf
[g
][i
] = count
[g
][i
]*inv_segvol
[i
]*normfac
;
361 for( ; (i
<nrdf
); i
++)
365 fp
=xvgropen(fnRDF
,"Radial Distribution","r","");
367 fprintf(fp
,"@ subtitle \"%s-%s\"\n",grpname
[0],grpname
[1]);
369 fprintf(fp
,"@ subtitle \"reference %s\"\n",grpname
[0]);
370 xvgr_legend(fp
,ng
,grpname
+1);
372 for(i
=0; (i
<nrdf
); i
++) {
373 fprintf(fp
,"%10g", (i
+0.5)*binwidth
);
375 fprintf(fp
," %10g",rdf
[g
][i
]);
382 /* h(Q) function: fourier transform of rdf */
385 real
*hq
,*integrand
,Q
;
387 /* Get a better number density later! */
390 snew(integrand
,nrdf
);
391 for(i
=0; (i
<nhq
); i
++) {
394 for(j
=1; (j
<nrdf
); j
++) {
395 r
= (j
+0.5)*binwidth
;
396 integrand
[j
] = (Q
== 0) ? 1.0 : sin(Q
*r
)/(Q
*r
);
397 integrand
[j
] *= 4.0*M_PI
*rho
*r
*r
*(rdf
[0][j
]-1.0);
399 hq
[i
] = print_and_integrate(debug
,nrdf
,binwidth
,integrand
,NULL
,0);
401 fp
=xvgropen(fnHQ
,"h(Q)","Q(/nm)","h(Q)");
402 for(i
=0; (i
<nhq
); i
++)
403 fprintf(fp
,"%10g %10g\n",i
*0.5,hq
[i
]);
411 normfac
= 1.0/(isize
[0]*nframes
);
412 fp
=xvgropen(fnCNRDF
,"Cumulative Number RDF","r","number");
414 fprintf(fp
,"@ subtitle \"%s-%s\"\n",grpname
[0],grpname
[1]);
416 fprintf(fp
,"@ subtitle \"reference %s\"\n",grpname
[0]);
417 xvgr_legend(fp
,ng
,grpname
+1);
420 for(i
=0; (i
<nbin
-1); i
++) {
421 fprintf(fp
,"%10g",i
*binwidth
);
422 for(g
=0; g
<ng
; g
++) {
423 fprintf(fp
," %10g",(real
)((double)sum
[g
]*normfac
));
424 sum
[g
] += count
[g
][i
];
431 do_view(fnCNRDF
,NULL
);
445 int comp_xdata(const void *a
,const void *b
)
455 else if (xb
->ndata
== 0)
458 tmp
= xa
->kkk
- xb
->kkk
;
467 static t_xdata
*init_xdata(int nx
,int ny
)
469 int ix
,iy
,i
,j
,maxkx
,maxky
;
475 for(i
=0; (i
<nx
); i
++) {
476 for(j
=0; (j
<ny
); j
++) {
477 ix
= abs((i
< maxkx
) ? i
: (i
- nx
));
478 iy
= abs((j
< maxky
) ? j
: (j
- ny
));
479 data
[ix
*ny
+iy
].kkk
= sqrt(ix
*ix
+iy
*iy
);
485 static void extract_sq(t_fftgrid
*fftgrid
,int nbin
,real k_max
,real lambda
,
486 real count
[],rvec box
,int npixel
,real
*map
[],
489 int nx
,ny
,nz
,nx2
,ny2
,nz2
,la2
,la12
;
491 int i
,j
,k
,maxkx
,maxky
,maxkz
,n
,ind
,ix
,iy
;
492 real k1
,kxy2
,kz2
,k2
,z
,kxy
,kxy_max
,cos_theta2
,ttt
,factor
;
497 kxy_max = k_max/sqrt(3);*/
498 unpack_fftgrid(fftgrid
,&nx
,&ny
,&nz
,&nx2
,&ny2
,&nz2
,
499 &la2
,&la12
,FALSE
,(t_fft_r
**)&ptr
);
500 /* This bit copied from pme.c */
505 for(i
=0; (i
<nx
); i
++) {
506 #define IDX(i,n) (i<=n/2) ? (i) : (i-n)
508 for(j
=0; (j
<ny
); j
++) {
512 for(k
=0; (k
<maxkz
); k
++,p0
++) {
513 if ((i
==0) && (j
==0) && (k
==0))
516 z
= sqrt(sqr(p0
->re
)+sqr(p0
->im
));
517 kxy2
= sqr(kk
[XX
]) + sqr(kk
[YY
]);
518 k2
= kxy2
+sqr(kk
[ZZ
]);
523 * R. W. James (1962),
524 * The Optical Principles of the Diffraction of X-Rays,
525 * Oxbow press, Woodbridge Connecticut
526 * the intensity is proportional to (eq. 9.10):
527 * I = C (1+cos^2 [2 theta])/2 FFT
529 * cos[2 theta] = cos^2[theta] - sin^2[theta] = 2 cos^2[theta] - 1
530 * we can compute the prefactor straight from cos[theta]
532 cos_theta2
= kxy2
/k2
;
533 /*ttt = z*0.5*(1+sqr(2*cos_theta2-1));*/
534 ttt
= z
*0.5*(1+cos_theta2
);
536 ix
= ((i
< maxkx
) ? i
: (i
- nx
));
537 iy
= ((j
< maxky
) ? j
: (j
- ny
));
538 map
[npixel
/2+ix
][npixel
/2+iy
] += ttt
;
539 data
[abs(ix
)*ny
+abs(iy
)].ndata
+= 1;
540 data
[abs(ix
)*ny
+abs(iy
)].intensity
+= ttt
;
543 fprintf(stderr
,"k (%g) > k_max (%g)\n",k1
,k_max
);
554 static void do_sq(char *fnNDX
,char *fnTPS
,char *fnTRX
,char *fnSQ
,
555 char *fnXPM
,real grid
,real lambda
,real distance
,
556 int npixel
,int nlevel
)
559 t_element elem
[] = { { "H", 1 }, { "C", 6 }, { "N", 7 }, { "O", 8 }, { "F", 9 }, { "S", 16 } };
560 #define NELEM asize(elem)
562 char title
[STRLEN
],*aname
;
563 int natoms
,i
,j
,k
,nbin
,j0
,j1
,n
,nframes
,pme_order
;
568 real I0
,C
,t
,k_max
,factor
,yfactor
,segvol
;
569 rvec
*x
,*xndx
,box_size
,kk
,lll
;
570 real fj0
,*fj
,max_spacing
,r
,lambda_1
;
573 int nx
,ny
,nz
,nelectron
;
574 atom_id ix
,jx
,**pairs
;
581 /* bTop=read_tps_conf(fnTPS,title,&top,&x,NULL,box,TRUE); */
583 fprintf(stderr
,"\nSelect group for structure factor computation:\n");
584 get_index(&top
.atoms
,fnNDX
,1,&isize
,&index
,&grpname
);
585 natoms
=read_first_x(&status
,fnTRX
,&t
,&x
,box
);
586 if (isize
< top
.atoms
.nr
)
590 fprintf(stderr
,"\n");
595 fatal_error(0,"Could not read coordinates from statusfile\n");
596 /* check with topology */
597 if ( natoms
> top
.atoms
.nr
)
598 fatal_error(0,"Trajectory (%d atoms) does not match topology (%d atoms)",
599 natoms
,top
.atoms
.nr
);
601 /* Atomic scattering factors */
605 for(i
=0; (i
<isize
); i
++) {
606 aname
= *(top
.atoms
.atomname
[index
[i
]]);
608 if (top
.atoms
.atom
[i
].ptype
== eptAtom
) {
609 for(j
=0; (j
<NELEM
); j
++)
610 if (aname
[0] == elem
[j
].name
[0]) {
615 fprintf(stderr
,"Warning: don't know number of electrons for atom %s\n",aname
);
617 /* Correct for partial charge */
618 fj
[i
] = fj0
- top
.atoms
.atom
[index
[i
]].q
;
625 /* Dump scattering factors */
626 for(i
=0; (i
<isize
); i
++)
627 fprintf(debug
,"Atom %3s-%5d q = %10.4f f = %10.4f\n",
628 *(top
.atoms
.atomname
[index
[i
]]),index
[i
],
629 top
.atoms
.atom
[index
[i
]].q
,fj
[i
]);
632 /* Constant for scattering */
633 C
= sqr(1.0/(ELECTRONMASS_keV
*KILO
*ELECTRONVOLT
*1e7
*distance
));
634 fprintf(stderr
,"C is %g\n",C
);
636 /* This bit is dimensionless */
638 max_spacing
= calc_grid(box
,grid
,&nx
,&ny
,&nz
,1);
639 pme_order
= max(4,1+(0.2/grid
));
641 data
= init_xdata(nx
,ny
);
643 fprintf(stderr
,"Largest grid spacing: %g nm, pme_order %d, %dx%d pixel on image\n",
644 max_spacing
,pme_order
,npixel
,npixel
);
645 init_pme(stdout
,NULL
,nx
,ny
,nz
,pme_order
,isize
,FALSE
,eewg3D
);
647 /* Determine largest k vector length. */
648 k_max
= 1+sqrt(sqr(1+nx
/2)+sqr(1+ny
/2)+sqr(1+nz
/2));
650 /* this is the S(q) array */
654 for(i
=0; (i
<npixel
); i
++)
659 /* Put the atoms with scattering factor on a grid. Misuses
660 * an old routine from the PPPM code.
662 for(j
=0; (j
<DIM
); j
++)
663 box_size
[j
] = box
[j
][j
];
665 /* Scale coordinates to the wavelength */
666 for(i
=0; (i
<isize
); i
++)
667 copy_rvec(x
[index
[i
]],xndx
[i
]);
669 /* put local atoms on grid. */
670 fftgrid
= spread_on_grid(stdout
,isize
,pme_order
,xndx
,fj
,box
,FALSE
);
672 /* FFT the density */
673 gmxfft3D(fftgrid
,FFTW_FORWARD
,NULL
);
675 /* Extract the Sq function and sum it into the average array */
676 extract_sq(fftgrid
,nbin
,k_max
,lambda
,count
,box_size
,npixel
,map
,data
);
679 } while (read_next_x(status
,&t
,natoms
,x
,box
));
680 fprintf(stderr
,"\n");
686 /* Normalize it ?? */
687 factor
= k_max
/(nbin
);
688 yfactor
= (1.0/nframes
)/*(1.0/fftgrid->nxyz)*/;
689 fp
=xvgropen(fnSQ
,"Structure Factor","q (1/nm)","S(q)");
690 fprintf(fp
,"@ subtitle \"Lambda = %g nm. Grid spacing = %g nm\"\n",
693 for(i
=0; i
<nbin
; i
++) {
694 r
= (i
+0.5)*factor
*2*M_PI
;
695 segvol
= 4*M_PI
*sqr(r
)*factor
;
696 fprintf(fp
,"%10g %10g\n",r
,count
[i
]*yfactor
/segvol
);
703 t_rgb rhi
= { 0,0,0 }, rlo
= { 1,1,1 };
704 real
*tx
,*ty
,hi
,inv_nframes
;
707 inv_nframes
= 1.0/nframes
;
710 for(i
=0; (i
<npixel
); i
++) {
714 for(j
=0; (j
<npixel
); j
++) {
715 map
[i
][j
] *= inv_nframes
;
716 hi
= max(hi
,map
[i
][j
]);
720 fp
= ffopen(fnXPM
,"w");
721 write_xpm(fp
,"Diffraction Image","Intensity","kx","ky",
722 nbin
,nbin
,tx
,ty
,map
,0,hi
,rlo
,rhi
,&nlevel
);
727 /* qsort(data,nx*ny,sizeof(data[0]),comp_xdata);
728 fp = ffopen("test.xvg","w");
729 for(i=0; (i<nx*ny); i++) {
730 if (data[i].ndata != 0) {
731 fprintf(fp,"%10.3f %10.3f\n",data[i].kkk,data[i].intensity/data[i].ndata);
739 t_complex
*** rc_tensor_allocation(int x
, int y
, int z
)
745 t
= (t_complex
***)calloc(x
,sizeof(t_complex
**));
746 if(!t
) exit(fprintf(stderr
,"\nallocation error"));
747 t
[0] = (t_complex
**)calloc(x
*y
,sizeof(t_complex
*));
748 if(!t
[0]) exit(fprintf(stderr
,"\nallocation error"));
749 t
[0][0] = (t_complex
*)calloc(x
*y
*z
,sizeof(t_complex
));
750 if(!t
[0][0]) exit(fprintf(stderr
,"\nallocation error"));
752 for( j
= 1 ; j
< y
; j
++)
753 t
[0][j
] = t
[0][j
-1] + z
;
754 for( i
= 1 ; i
< x
; i
++) {
756 t
[i
][0] = t
[i
-1][0] + y
*z
;
757 for( j
= 1 ; j
< y
; j
++)
758 t
[i
][j
] = t
[i
][j
-1] + z
;
763 int return_atom_type (char *type
)
767 for (i
= 0; (i
< asize(CM_t
)); i
++)
768 if (!strcmp (type
, CM_t
[i
].Label
))
770 fatal_error(0,"\nError: atom type (%s) not in list (%d types checked)!\n",
774 double CMSF (int type
, double lambda
, double sin_theta
)
776 * return Cromer-Mann fit for the atomic scattering factor:
777 * sin_theta is the sine of half the angle between incoming and scattered
778 * vectors. See g_sq.h for a short description of CM fit.
782 double tmp
= 0.0, k2
;
784 k2
= (sqr (sin_theta
) / sqr (10.0 * lambda
));
789 if (!strcmp (CM_t
[type
].Label
, "CS2") ||
790 !strcmp (CM_t
[type
].Label
, "CH2") ||
791 !strcmp (CM_t
[type
].Label
, "CH3") ||
792 !strcmp (CM_t
[type
].Label
, "CS3")) {
793 tmp
= CMSF (return_atom_type ("C"), lambda
, sin_theta
);
794 if (!strcmp (CM_t
[type
].Label
, "CH3") ||
795 !strcmp (CM_t
[type
].Label
, "CS3"))
796 tmp
+= 3.0 * CMSF (return_atom_type ("H"), lambda
, sin_theta
);
798 tmp
+= 2.0 * CMSF (return_atom_type ("H"), lambda
, sin_theta
);
803 for (i
= 0; i
< 4; i
++)
804 tmp
+= CM_t
[type
].a
[i
] * exp (-CM_t
[type
].b
[i
] * k2
);
809 void compute_scattering_factor_table (structure_factor
* sf
)
812 * this function build up a table of scattering factors for every atom
813 * type and for every scattering angle.
818 /* \hbar \omega \lambda = hc = 1239.842 eV * nm */
819 sf
->momentum
= ((double) (2. * 1000.0 * M_PI
* sf
->energy
) / 1239.842);
820 sf
->lambda
= 1239.842 / (1000.0 * sf
->energy
);
821 fprintf (stderr
, "\nwavelenght = %f nm\n", sf
->lambda
);
822 snew (sf_table
,asize (CM_t
));
823 for (i
= 0; (i
< asize(CM_t
)); i
++) {
824 snew (sf_table
[i
], sf
->n_angles
);
825 for (j
= 0; j
< sf
->n_angles
; j
++) {
826 q
= ((double) j
* sf
->ref_k
);
827 /* theta is half the angle between incoming and scattered wavevectors. */
828 sin_theta
= q
/ (2.0 * sf
->momentum
);
829 sf_table
[i
][j
] = CMSF (i
, sf
->lambda
, sin_theta
);
834 int * create_indexed_atom_type (reduced_atom
* atm
, int size
)
837 * create an index of the atom types found in a group
838 * i.e.: for water index_atp[0]=type_number_of_O and
839 * index_atp[1]=type_number_of_H
841 * the last element is set to 0
843 int *index_atp
, i
, i_tmp
, j
;
847 index_atp
[0] = atm
[0].t
;
848 for (i
= 1; i
< size
; i
++) {
849 for (j
= 0; j
< i_tmp
; j
++)
850 if (atm
[i
].t
== index_atp
[j
])
852 if (j
== i_tmp
) { /* i.e. no indexed atom type is == to atm[i].t */
854 srenew (index_atp
, i_tmp
* sizeof (int));
855 index_atp
[i_tmp
- 1] = atm
[i
].t
;
859 srenew (index_atp
, i_tmp
* sizeof (int));
860 index_atp
[i_tmp
- 1] = 0;
864 void rearrange_atoms (reduced_atom
* positions
, t_trxframe
* fr
, atom_id
* index
,
865 int isize
, t_topology
* top
, real
** sf_table
, bool flag
)
866 /* given the group's index, return the (continuous) array of atoms */
871 for (i
= 0; i
< isize
; i
++)
873 return_atom_type (*(top
->atoms
.atomtype
[index
[i
]]));
874 for (i
= 0; i
< isize
; i
++)
875 copy_rvec (fr
->x
[index
[i
]], positions
[i
].x
);
879 int atp_size (int *index_atp
)
888 void compute_structure_factor (structure_factor
* sf
, matrix box
,
889 reduced_atom
* red
, int isize
, real start_q
,
890 real end_q
, int group
)
894 real kdotx
, asf
, kx
, ky
, kz
, krr
;
895 int kr
, maxkx
, maxky
, maxkz
, i
, j
, k
, p
, *counter
;
898 k_factor
[XX
] = 2 * M_PI
/ box
[XX
][XX
];
899 k_factor
[YY
] = 2 * M_PI
/ box
[YY
][YY
];
900 k_factor
[ZZ
] = 2 * M_PI
/ box
[ZZ
][ZZ
];
902 maxkx
= (int) rint (end_q
/ k_factor
[XX
]);
903 maxky
= (int) rint (end_q
/ k_factor
[YY
]);
904 maxkz
= (int) rint (end_q
/ k_factor
[ZZ
]);
906 snew (counter
, sf
->n_angles
);
908 tmpSF
= rc_tensor_allocation(maxkx
,maxky
,maxkz
);
911 * compute real and imaginary part of the structure factor for every
914 fprintf(stderr
,"\n");
915 for (i
= 0; i
< maxkx
; i
++) {
916 fprintf (stderr
,"\rdone %3.1f%% ", (double)(100.0*(i
+1))/maxkx
);
918 kx
= i
* k_factor
[XX
];
919 for (j
= 0; j
< maxky
; j
++) {
920 ky
= j
* k_factor
[YY
];
921 for (k
= 0; k
< maxkz
; k
++)
922 if (i
!= 0 || j
!= 0 || k
!= 0) {
923 kz
= k
* k_factor
[ZZ
];
924 krr
= sqrt (sqr (kx
) + sqr (ky
) + sqr (kz
));
925 if (krr
>= start_q
&& krr
<= end_q
) {
926 kr
= (int) rint (krr
/sf
->ref_k
);
927 if (kr
< sf
->n_angles
) {
928 counter
[kr
]++; /* will be used for the copmutation
930 for (p
= 0; p
< isize
; p
++) {
932 asf
= sf_table
[red
[p
].t
][kr
];
934 kdotx
= kx
* red
[p
].x
[XX
] +
935 ky
* red
[p
].x
[YY
] + kz
* red
[p
].x
[ZZ
];
937 tmpSF
[i
][j
][k
].re
+= cos (kdotx
) * asf
;
938 tmpSF
[i
][j
][k
].im
+= sin (kdotx
) * asf
;
944 } /* end loop on i */
946 * compute the square modulus of the structure factor, averaging on the surface
947 * kx*kx + ky*ky + kz*kz = krr*krr
948 * note that this is correct only for a (on the macroscopic scale)
951 for (i
= 0; i
< maxkx
; i
++) {
952 kx
= i
* k_factor
[XX
]; for (j
= 0; j
< maxky
; j
++) {
953 ky
= j
* k_factor
[YY
]; for (k
= 0; k
< maxkz
; k
++) {
954 kz
= k
* k_factor
[ZZ
]; krr
= sqrt (sqr (kx
) + sqr (ky
)
955 + sqr (kz
)); if (krr
>= start_q
&& krr
<= end_q
) {
956 kr
= (int) rint (krr
/ sf
->ref_k
); if (kr
<
957 sf
->n_angles
&& counter
[kr
] != 0)
959 (sqr (tmpSF
[i
][j
][k
].re
) +
960 sqr (tmpSF
[i
][j
][k
].im
))/ counter
[kr
];
964 } sfree (counter
); free(tmpSF
[0][0]); free(tmpSF
[0]); free(tmpSF
);
967 void save_data (structure_factor
* sf
, char *file
, int ngrps
, real start_q
,
973 double *tmp
, polarization_factor
, A
;
975 fp
= xvgropen (file
, "Scattering Intensity", "q (1/nm)",
980 for (g
= 0; g
< ngrps
; g
++)
981 for (i
= 0; i
< sf
->n_angles
; i
++) {
983 * theta is half the angle between incoming and scattered vectors.
985 * polar. fact. = 0.5*(1+cos^2(2*theta)) = 1 - 0.5 * sin^2(2*theta)
987 * sin(theta) = q/(2k) := A -> sin^2(theta) = 4*A^2 (1-A^2) ->
988 * -> 0.5*(1+cos^2(2*theta)) = 1 - 2 A^2 (1-A^2)
990 A
= (double) (i
* sf
->ref_k
) / (2.0 * sf
->momentum
);
991 polarization_factor
= 1 - 2.0 * sqr (A
) * (1 - sqr (A
));
992 sf
->F
[g
][i
] *= polarization_factor
;
994 for (i
= 0; i
< sf
->n_angles
; i
++) {
995 if (i
* sf
->ref_k
>= start_q
&& i
* sf
->ref_k
<= end_q
) {
996 fprintf (fp
, "%10.5f ", i
* sf
->ref_k
);
997 for (g
= 0; g
< ngrps
; g
++)
998 fprintf (fp
, " %10.5f ", (sf
->F
[g
][i
]) /( sf
->total_n_atoms
*
1007 do_scattering_intensity (char* fnTPS
, char* fnNDX
, char* fnXVG
, char *fnTRX
,
1008 real start_q
,real end_q
, real energy
)
1010 int i
,ng
,*isize
,status
,flags
= TRX_READ_X
,**index_atp
;
1011 char **grpname
,title
[STRLEN
];
1016 structure_factor
*sf
;
1022 sf
->energy
= energy
;
1023 /* Read the topology informations */
1025 read_tps_conf (fnTPS
, title
, &top
, &xtop
, NULL
, box
, TRUE
);
1028 /* groups stuff... */
1030 "\nHow many groups do you want to calculate the I(q) of?\n");
1039 fprintf (stderr
, "\nSelect %d group%s\n", ng
,
1040 ng
== 1 ? "" : "s");
1042 get_index (&top
.atoms
, fnNDX
, ng
, isize
, index
, grpname
);
1044 rd_index (fnNDX
, ng
, isize
, index
, grpname
);
1046 /* The first time we read data is a little special */
1047 read_first_frame (&status
, fnTRX
, &fr
, flags
);
1049 sf
->total_n_atoms
= fr
.natoms
;
1052 snew (index_atp
, ng
);
1054 r_tmp
= max (box
[XX
][XX
], box
[YY
][YY
]);
1055 r_tmp
= (double) max (box
[ZZ
][ZZ
], r_tmp
);
1057 sf
->ref_k
= (2.0 * M_PI
) / (r_tmp
);
1058 /* ref_k will be the reference momentum unit */
1059 sf
->n_angles
= (int) rint (end_q
/ sf
->ref_k
);
1062 for (i
= 0; i
< ng
; i
++)
1063 snew (sf
->F
[i
], sf
->n_angles
);
1064 for (i
= 0; i
< ng
; i
++) {
1065 snew (red
[i
], isize
[i
]);
1066 rearrange_atoms (red
[i
], &fr
, index
[i
], isize
[i
], &top
, sf_table
, 1);
1067 index_atp
[i
] = create_indexed_atom_type (red
[i
], isize
[i
]);
1069 compute_scattering_factor_table (sf
);
1070 /* This is the main loop over frames */
1075 for (i
= 0; i
< ng
; i
++) {
1076 rearrange_atoms (red
[i
], &fr
, index
[i
], isize
[i
], &top
,
1079 compute_structure_factor (sf
, box
, red
[i
], isize
[i
],
1083 while (read_next_frame (status
, &fr
));
1085 save_data (sf
, fnXVG
, ng
, start_q
, end_q
);
1090 int gmx_rdf(int argc
,char *argv
[])
1092 static char *desc
[] = {
1093 "The structure of liquids can be studied by either neutron or X-ray",
1094 "scattering. The most common way to describe liquid structure is by a",
1095 "radial distribution function. However, this is not easy to obtain from",
1096 "a scattering experiment.[PAR]",
1097 "g_rdf calculates radial distribution functions in different ways.",
1098 "The normal method is around a (set of) particle(s), the other method",
1099 "is around the center of mass of a set of particles.[PAR]",
1100 "If a run input file is supplied ([TT]-s[tt]), exclusions defined",
1101 "in that file are taken into account when calculating the rdf.",
1102 "The option [TT]-cut[tt] is meant as an alternative way to avoid",
1103 "intramolecular peaks in the rdf plot.",
1104 "It is however better to supply a run input file with a higher number of",
1105 "exclusions. For eg. benzene a topology with nrexcl set to 5",
1106 "would eliminate all intramolecular contributions to the rdf.",
1107 "Note that all atoms in the selected groups are used, also the ones",
1108 "that don't have Lennard-Jones interactions.[PAR]",
1109 "Option [TT]-cn[tt] produces the cumulative number rdf.[PAR]"
1110 "To bridge the gap between theory and experiment structure factors can",
1111 "be computed (option [TT]-sq[tt]). The algorithm uses FFT, the grid"
1112 "spacing of which is determined by option [TT]-grid[tt]."
1114 static bool bCM
=FALSE
;
1115 static real cutoff
=0,binwidth
=0.001,grid
=0.05,fade
=0.0,lambda
=0.1,distance
=10;
1116 static int npixel
=256,nlevel
=20;
1117 static real start_q
=0.0, end_q
=60.0, energy
=12.0;
1119 { "-bin", FALSE
, etREAL
, {&binwidth
},
1121 { "-com", FALSE
, etBOOL
, {&bCM
},
1122 "RDF with respect to the center of mass of first group" },
1123 { "-cut", FALSE
, etREAL
, {&cutoff
},
1124 "Shortest distance (nm) to be considered"},
1125 { "-fade", FALSE
, etREAL
, {&fade
},
1126 "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." },
1127 /* { "-grid", FALSE, etREAL, {&grid},
1128 "Grid spacing (in nm) for FFTs when computing structure factors" },*/
1129 { "-npixel", FALSE
, etINT
, {&npixel
},
1130 "[HIDDEN]# pixels per edge of the square detector plate" },
1131 { "-nlevel", FALSE
, etINT
, {&nlevel
},
1132 "Number of different colors in the diffraction image" },
1133 { "-distance", FALSE
, etREAL
, {&distance
},
1134 "[HIDDEN]Distance (in cm) from the sample to the detector" },
1135 /* { "-wave", FALSE, etREAL, {&lambda},
1136 "Wavelength for X-rays/Neutrons for scattering. 0.1 nm corresponds to roughly 12 keV" },
1138 {"-startq", FALSE
, etREAL
, {&start_q
},
1139 "Starting q (1/nm) "},
1140 {"-endq", FALSE
, etREAL
, {&end_q
},
1142 {"-energy", FALSE
, etREAL
, {&energy
},
1143 "Energy of the incoming X-ray (keV) "}
1145 #define NPA asize(pa)
1150 { efTRX
, "-f", NULL
, ffREAD
},
1151 { efTPS
, NULL
, NULL
, ffOPTRD
},
1152 { efNDX
, NULL
, NULL
, ffOPTRD
},
1153 { efXVG
, "-o", "rdf", ffOPTWR
},
1154 { efXVG
, "-sq", "sq", ffOPTWR
},
1155 { efXVG
, "-cn", "rdf_cn", ffOPTWR
},
1156 { efXVG
, "-hq", "hq", ffOPTWR
},
1157 /* { efXPM, "-image", "sq", ffOPTWR }*/
1159 #define NFILE asize(fnm)
1161 CopyRight(stderr
,argv
[0]);
1162 parse_common_args(&argc
,argv
,PCA_CAN_VIEW
| PCA_CAN_TIME
| PCA_BE_NICE
,
1163 NFILE
,fnm
,NPA
,pa
,asize(desc
),desc
,0,NULL
);
1165 fnTPS
= ftp2fn_null(efTPS
,NFILE
,fnm
);
1166 fnNDX
= ftp2fn_null(efNDX
,NFILE
,fnm
);
1167 // bSQ = opt2bSet("-sq",NFILE,fnm) || opt2parg_bSet("-grid",NPA,pa);
1168 bSQ
= opt2bSet("-sq",NFILE
,fnm
);
1169 bRDF
= opt2bSet("-o",NFILE
,fnm
) || !bSQ
;
1173 fatal_error(0,"Need a tps file for calculating structure factors\n");
1176 if (!fnTPS
&& !fnNDX
)
1177 fatal_error(0,"Neither index file nor topology file specified\n"
1182 do_scattering_intensity(fnTPS
,fnNDX
,opt2fn("-sq",NFILE
,fnm
),ftp2fn(efTRX
,NFILE
,fnm
),
1183 start_q
, end_q
, energy
);
1184 /* old structure factor code */
1185 /* do_sq(fnNDX,fnTPS,ftp2fn(efTRX,NFILE,fnm),opt2fn("-sq",NFILE,fnm),
1186 ftp2fn(efXPM,NFILE,fnm),grid,lambda,distance,npixel,nlevel);
1189 do_rdf(fnNDX
,fnTPS
,ftp2fn(efTRX
,NFILE
,fnm
),
1190 opt2fn("-o",NFILE
,fnm
),opt2fn_null("-cn",NFILE
,fnm
),
1191 opt2fn_null("-hq",NFILE
,fnm
),
1192 bCM
,cutoff
,binwidth
,fade
);