2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team,
6 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
48 #include "gmx_fatal.h"
66 typedef struct gmx_structurefactors {
68 int *p; /* proton number */
69 int *n; /* neutron number */
70 /* Parameters for the Cromer Mann fit */
71 real **a; /* parameter a */
72 real **b; /* parameter b */
73 real *c; /* parameter c */
74 char **atomnm; /* atomname */
76 } gmx_structurefactors;
78 typedef struct reduced_atom{
84 typedef struct structure_factor
98 extern int * create_indexed_atom_type (reduced_atom_t * atm, int size)
101 * create an index of the atom types found in a group
102 * i.e.: for water index_atp[0]=type_number_of_O and
103 * index_atp[1]=type_number_of_H
105 * the last element is set to 0
107 int *index_atp, i, i_tmp, j;
109 reduced_atom *att=(reduced_atom *)atm;
113 index_atp[0] = att[0].t;
114 for (i = 1; i < size; i++) {
115 for (j = 0; j < i_tmp; j++)
116 if (att[i].t == index_atp[j])
118 if (j == i_tmp) { /* i.e. no indexed atom type is == to atm[i].t */
120 srenew (index_atp, i_tmp * sizeof (int));
121 index_atp[i_tmp - 1] = att[i].t;
125 srenew (index_atp, i_tmp * sizeof (int));
126 index_atp[i_tmp - 1] = 0;
132 extern t_complex *** rc_tensor_allocation(int x, int y, int z)
137 t = (t_complex ***)calloc(x,sizeof(t_complex**));
138 if(!t) exit(fprintf(stderr,"\nallocation error"));
139 t[0] = (t_complex **)calloc(x*y,sizeof(t_complex*));
140 if(!t[0]) exit(fprintf(stderr,"\nallocation error"));
141 t[0][0] = (t_complex *)calloc(x*y*z,sizeof(t_complex));
142 if(!t[0][0]) exit(fprintf(stderr,"\nallocation error"));
144 for( j = 1 ; j < y ; j++)
145 t[0][j] = t[0][j-1] + z;
146 for( i = 1 ; i < x ; i++) {
148 t[i][0] = t[i-1][0] + y*z;
149 for( j = 1 ; j < y ; j++)
150 t[i][j] = t[i][j-1] + z;
156 extern void compute_structure_factor (structure_factor_t * sft, matrix box,
157 reduced_atom_t * red, int isize, real start_q,
158 real end_q, int group,real **sf_table)
160 structure_factor *sf=(structure_factor *)sft;
161 reduced_atom *redt=(reduced_atom *)red;
165 real kdotx, asf, kx, ky, kz, krr;
166 int kr, maxkx, maxky, maxkz, i, j, k, p, *counter;
169 k_factor[XX] = 2 * M_PI / box[XX][XX];
170 k_factor[YY] = 2 * M_PI / box[YY][YY];
171 k_factor[ZZ] = 2 * M_PI / box[ZZ][ZZ];
173 maxkx = (int) (end_q / k_factor[XX] + 0.5);
174 maxky = (int) (end_q / k_factor[YY] + 0.5);
175 maxkz = (int) (end_q / k_factor[ZZ] + 0.5);
177 snew (counter, sf->n_angles);
179 tmpSF = rc_tensor_allocation(maxkx,maxky,maxkz);
182 * compute real and imaginary part of the structure factor for every
185 fprintf(stderr,"\n");
186 for (i = 0; i < maxkx; i++) {
187 fprintf (stderr,"\rdone %3.1f%% ", (double)(100.0*(i+1))/maxkx);
188 kx = i * k_factor[XX];
189 for (j = 0; j < maxky; j++) {
190 ky = j * k_factor[YY];
191 for (k = 0; k < maxkz; k++)
192 if (i != 0 || j != 0 || k != 0) {
193 kz = k * k_factor[ZZ];
194 krr = sqrt (sqr (kx) + sqr (ky) + sqr (kz));
195 if (krr >= start_q && krr <= end_q) {
196 kr = (int) (krr/sf->ref_k + 0.5);
197 if (kr < sf->n_angles) {
198 counter[kr]++; /* will be used for the copmutation
200 for (p = 0; p < isize; p++) {
201 asf = sf_table[redt[p].t][kr];
203 kdotx = kx * redt[p].x[XX] +
204 ky * redt[p].x[YY] + kz * redt[p].x[ZZ];
206 tmpSF[i][j][k].re += cos (kdotx) * asf;
207 tmpSF[i][j][k].im += sin (kdotx) * asf;
213 } /* end loop on i */
215 * compute the square modulus of the structure factor, averaging on the surface
216 * kx*kx + ky*ky + kz*kz = krr*krr
217 * note that this is correct only for a (on the macroscopic scale)
220 for (i = 0; i < maxkx; i++) {
221 kx = i * k_factor[XX]; for (j = 0; j < maxky; j++) {
222 ky = j * k_factor[YY]; for (k = 0; k < maxkz; k++) {
223 kz = k * k_factor[ZZ]; krr = sqrt (sqr (kx) + sqr (ky)
224 + sqr (kz)); if (krr >= start_q && krr <= end_q) {
225 kr = (int) (krr / sf->ref_k + 0.5);
226 if (kr < sf->n_angles && counter[kr] != 0)
228 (sqr (tmpSF[i][j][k].re) +
229 sqr (tmpSF[i][j][k].im))/ counter[kr];
233 } sfree (counter); free(tmpSF[0][0]); free(tmpSF[0]); free(tmpSF);
237 extern gmx_structurefactors_t *gmx_structurefactors_init(const char *datfn) {
239 /* Read the database for the structure factor of the different atoms */
243 gmx_structurefactors *gsf;
244 double a1,a2,a3,a4,b1,b2,b3,b4,c;
254 snew(gsf->atomnm,nralloc);
255 snew(gsf->a,nralloc);
256 snew(gsf->b,nralloc);
257 snew(gsf->c,nralloc);
258 snew(gsf->p,nralloc);
260 gsf->nratoms=line_no;
261 while(get_a_line(fp,line,STRLEN)) {
263 if (sscanf(line,"%s %d %lf %lf %lf %lf %lf %lf %lf %lf %lf",
264 atomn,&p,&a1,&a2,&a3,&a4,&b1,&b2,&b3,&b4,&c) == 11) {
265 gsf->atomnm[i]=strdup(atomn);
279 gsf->nratoms=line_no;
280 if (line_no==nralloc){
282 srenew(gsf->atomnm,nralloc);
283 srenew(gsf->a,nralloc);
284 srenew(gsf->b,nralloc);
285 srenew(gsf->c,nralloc);
286 srenew(gsf->p,nralloc);
290 fprintf(stderr,"WARNING: Error in file %s at line %d ignored\n",
294 srenew(gsf->atomnm,gsf->nratoms);
295 srenew(gsf->a,gsf->nratoms);
296 srenew(gsf->b,gsf->nratoms);
297 srenew(gsf->c,gsf->nratoms);
298 srenew(gsf->p,gsf->nratoms);
302 return (gmx_structurefactors_t *) gsf;
307 extern void rearrange_atoms (reduced_atom_t * positions, t_trxframe *fr, atom_id * index,
308 int isize, t_topology * top, gmx_bool flag,gmx_structurefactors_t *gsf)
309 /* given the group's index, return the (continuous) array of atoms */
313 reduced_atom *pos=(reduced_atom *)positions;
316 for (i = 0; i < isize; i++)
318 return_atom_type (*(top->atoms.atomname[index[i]]),gsf);
319 for (i = 0; i < isize; i++)
320 copy_rvec (fr->x[index[i]], pos[i].x);
322 positions=(reduced_atom_t *)pos;
326 extern int return_atom_type (const char *name,gmx_structurefactors_t *gsf)
333 { "CH1", 1 }, { "CH2", 2 }, { "CH3", 3 },
334 { "CS1", 1 }, { "CS2", 2 }, { "CS3", 3 },
335 { "CP1", 1 }, { "CP2", 2 }, { "CP3", 3 }
343 gmx_structurefactors *gsft=(gmx_structurefactors *)gsf;
349 for(i=0; (i<asize(uh)); i++)
350 if (strcmp(name,uh[i].name) == 0)
351 return NCMT-1+uh[i].nh;
353 for(i=0; (i<NCMT); i++){
354 if (strncmp (name, gsft->atomnm[i],strlen(gsft->atomnm[i])) == 0){
361 gmx_fatal(FARGS,"\nError: atom (%s) not in list (%d types checked)!\n",
366 if(strlen(gsft->atomnm[tndx[i]])>(size_t)nrc){
367 nrc=strlen(gsft->atomnm[tndx[i]]);
378 extern int gmx_structurefactors_get_sf(gmx_structurefactors_t *gsf, int elem, real a[4], real b[4], real *c){
382 gmx_structurefactors *gsft=(gmx_structurefactors *)gsf;
386 a[i]=gsft->a[elem][i];
387 b[i]=gsft->b[elem][i];
395 extern int do_scattering_intensity (const char* fnTPS, const char* fnNDX,
396 const char* fnXVG, const char *fnTRX,
398 real start_q,real end_q,
399 real energy,int ng,const output_env_t oenv)
401 int i,*isize,flags = TRX_READ_X,**index_atp;
403 char **grpname,title[STRLEN];
408 reduced_atom_t **red;
409 structure_factor *sf;
416 gmx_structurefactors_t *gmx_sf;
424 gmx_sf=gmx_structurefactors_init(fnDAT);
426 success=gmx_structurefactors_get_sf(gmx_sf,0, a, b, &c);
431 /* Read the topology informations */
432 read_tps_conf (fnTPS, title, &top, &ePBC, &xtop, NULL, box, TRUE);
435 /* groups stuff... */
440 fprintf (stderr, "\nSelect %d group%s\n", ng,
443 get_index (&top.atoms, fnNDX, ng, isize, index, grpname);
445 rd_index (fnNDX, ng, isize, index, grpname);
447 /* The first time we read data is a little special */
448 read_first_frame (oenv,&status, fnTRX, &fr, flags);
450 sf->total_n_atoms = fr.natoms;
453 snew (index_atp, ng);
455 r_tmp = max (box[XX][XX], box[YY][YY]);
456 r_tmp = (double) max (box[ZZ][ZZ], r_tmp);
458 sf->ref_k = (2.0 * M_PI) / (r_tmp);
459 /* ref_k will be the reference momentum unit */
460 sf->n_angles = (int) (end_q / sf->ref_k + 0.5);
463 for (i = 0; i < ng; i++)
464 snew (sf->F[i], sf->n_angles);
465 for (i = 0; i < ng; i++) {
466 snew (red[i], isize[i]);
467 rearrange_atoms (red[i], &fr, index[i], isize[i], &top, TRUE,gmx_sf);
468 index_atp[i] = create_indexed_atom_type (red[i], isize[i]);
471 sf_table = compute_scattering_factor_table (gmx_sf,(structure_factor_t *)sf,&nsftable);
474 /* This is the main loop over frames */
478 for (i = 0; i < ng; i++) {
479 rearrange_atoms (red[i], &fr, index[i], isize[i], &top,FALSE,gmx_sf);
481 compute_structure_factor ((structure_factor_t *)sf, box, red[i], isize[i],
482 start_q, end_q, i, sf_table);
486 while (read_next_frame (oenv,status, &fr));
488 save_data ((structure_factor_t *)sf, fnXVG, ng, start_q, end_q,oenv);
494 gmx_structurefactors_done(gmx_sf);
500 extern void save_data (structure_factor_t *sft, const char *file, int ngrps,
501 real start_q, real end_q, const output_env_t oenv)
506 double *tmp, polarization_factor, A;
508 structure_factor *sf=(structure_factor *)sft;
510 fp = xvgropen (file, "Scattering Intensity", "q (1/nm)",
511 "Intensity (a.u.)",oenv);
515 for (g = 0; g < ngrps; g++)
516 for (i = 0; i < sf->n_angles; i++) {
519 * theta is half the angle between incoming and scattered vectors.
521 * polar. fact. = 0.5*(1+cos^2(2*theta)) = 1 - 0.5 * sin^2(2*theta)
523 * sin(theta) = q/(2k) := A -> sin^2(theta) = 4*A^2 (1-A^2) ->
524 * -> 0.5*(1+cos^2(2*theta)) = 1 - 2 A^2 (1-A^2)
526 A = (double) (i * sf->ref_k) / (2.0 * sf->momentum);
527 polarization_factor = 1 - 2.0 * sqr (A) * (1 - sqr (A));
528 sf->F[g][i] *= polarization_factor;
530 for (i = 0; i < sf->n_angles; i++) {
531 if (i * sf->ref_k >= start_q && i * sf->ref_k <= end_q) {
532 fprintf (fp, "%10.5f ", i * sf->ref_k);
533 for (g = 0; g < ngrps; g++)
534 fprintf (fp, " %10.5f ", (sf->F[g][i]) /( sf->total_n_atoms*
544 extern double CMSF (gmx_structurefactors_t *gsf,int type,int nh,double lambda, double sin_theta)
546 * return Cromer-Mann fit for the atomic scattering factor:
547 * sin_theta is the sine of half the angle between incoming and scattered
548 * vectors. See g_sq.h for a short description of CM fit.
552 double tmp = 0.0, k2;
561 * f0[k] = c + [SUM a_i*EXP(-b_i*(k^2)) ]
570 tmp = (CMSF (gsf,return_atom_type ("C",gsf),0,lambda, sin_theta) +
571 nh*CMSF (gsf,return_atom_type ("H",gsf),0,lambda, sin_theta));
575 k2 = (sqr (sin_theta) / sqr (10.0 * lambda));
576 success=gmx_structurefactors_get_sf(gsf,type,a,b,&c);
578 for (i = 0; (i < 4); i++)
579 tmp += a[i] * exp (-b[i] * k2);
586 extern real **gmx_structurefactors_table(gmx_structurefactors_t *gsf,real momentum, real ref_k, real lambda, int n_angles){
593 gmx_structurefactors *gsft=(gmx_structurefactors *)gsf;
598 snew (sf_table,nsftable);
599 for (i = 0; (i < nsftable); i++) {
600 snew (sf_table[i], n_angles);
601 for (j = 0; j < n_angles; j++) {
602 q = ((double) j * ref_k);
603 /* theta is half the angle between incoming
604 and scattered wavevectors. */
605 sin_theta = q / (2.0 * momentum);
607 sf_table[i][j] = CMSF (gsf,i,0,lambda, sin_theta);
610 sf_table[i][j] = CMSF (gsf,i,i-NCMT+1,lambda, sin_theta);
616 extern void gmx_structurefactors_done(gmx_structurefactors_t *gsf){
619 gmx_structurefactors *sf;
620 sf=(gmx_structurefactors *) gsf;
622 for(i=0;i<sf->nratoms;i++){
625 sfree(sf->atomnm[i]);
638 extern real **compute_scattering_factor_table (gmx_structurefactors_t *gsf,structure_factor_t *sft,int *nsftable)
641 * this function build up a table of scattering factors for every atom
642 * type and for every scattering angle.
648 structure_factor *sf=(structure_factor *)sft;
651 /* \hbar \omega \lambda = hc = 1239.842 eV * nm */
652 sf->momentum = ((double) (2. * 1000.0 * M_PI * sf->energy) / hc);
653 sf->lambda = hc / (1000.0 * sf->energy);
654 fprintf (stderr, "\nwavelenght = %f nm\n", sf->lambda);
656 sf_table=gmx_structurefactors_table(gsf,sf->momentum,sf->ref_k,sf->lambda,sf->n_angles);