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 * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
46 #include "gromacs/fileio/confio.h"
47 #include "gromacs/fileio/strdb.h"
48 #include "gromacs/fileio/trxio.h"
49 #include "gromacs/fileio/xvgr.h"
50 #include "gromacs/legacyheaders/macros.h"
51 #include "gromacs/legacyheaders/names.h"
52 #include "gromacs/legacyheaders/oenv.h"
53 #include "gromacs/legacyheaders/typedefs.h"
54 #include "gromacs/math/utilities.h"
55 #include "gromacs/math/vec.h"
56 #include "gromacs/topology/index.h"
57 #include "gromacs/utility/cstringutil.h"
58 #include "gromacs/utility/fatalerror.h"
59 #include "gromacs/utility/futil.h"
60 #include "gromacs/utility/smalloc.h"
63 typedef struct gmx_structurefactors {
65 int *p; /* proton number */
66 int *n; /* neutron number */
67 /* Parameters for the Cromer Mann fit */
68 real **a; /* parameter a */
69 real **b; /* parameter b */
70 real *c; /* parameter c */
71 char **atomnm; /* atomname */
73 } gmx_structurefactors;
75 typedef struct reduced_atom{
81 typedef struct structure_factor
95 extern int * create_indexed_atom_type (reduced_atom_t * atm, int size)
98 * create an index of the atom types found in a group
99 * i.e.: for water index_atp[0]=type_number_of_O and
100 * index_atp[1]=type_number_of_H
102 * the last element is set to 0
104 int *index_atp, i, i_tmp, j;
106 reduced_atom *att = (reduced_atom *)atm;
110 index_atp[0] = att[0].t;
111 for (i = 1; i < size; i++)
113 for (j = 0; j < i_tmp; j++)
115 if (att[i].t == index_atp[j])
120 if (j == i_tmp) /* i.e. no indexed atom type is == to atm[i].t */
123 srenew (index_atp, i_tmp * sizeof (int));
124 index_atp[i_tmp - 1] = att[i].t;
128 srenew (index_atp, i_tmp * sizeof (int));
129 index_atp[i_tmp - 1] = 0;
135 extern t_complex *** rc_tensor_allocation(int x, int y, int z)
142 snew(t[0][0], x*y*z);
144 for (j = 1; j < y; j++)
146 t[0][j] = t[0][j-1] + z;
148 for (i = 1; i < x; i++)
151 t[i][0] = t[i-1][0] + y*z;
152 for (j = 1; j < y; j++)
154 t[i][j] = t[i][j-1] + z;
161 extern void compute_structure_factor (structure_factor_t * sft, matrix box,
162 reduced_atom_t * red, int isize, real start_q,
163 real end_q, int group, real **sf_table)
165 structure_factor *sf = (structure_factor *)sft;
166 reduced_atom *redt = (reduced_atom *)red;
170 real kdotx, asf, kx, ky, kz, krr;
171 int kr, maxkx, maxky, maxkz, i, j, k, p, *counter;
174 k_factor[XX] = 2 * M_PI / box[XX][XX];
175 k_factor[YY] = 2 * M_PI / box[YY][YY];
176 k_factor[ZZ] = 2 * M_PI / box[ZZ][ZZ];
178 maxkx = static_cast<int>(end_q / k_factor[XX] + 0.5);
179 maxky = static_cast<int>(end_q / k_factor[YY] + 0.5);
180 maxkz = static_cast<int>(end_q / k_factor[ZZ] + 0.5);
182 snew (counter, sf->n_angles);
184 tmpSF = rc_tensor_allocation(maxkx, maxky, maxkz);
187 * compute real and imaginary part of the structure factor for every
190 fprintf(stderr, "\n");
191 for (i = 0; i < maxkx; i++)
193 fprintf (stderr, "\rdone %3.1f%% ", (100.0*(i+1))/maxkx);
194 kx = i * k_factor[XX];
195 for (j = 0; j < maxky; j++)
197 ky = j * k_factor[YY];
198 for (k = 0; k < maxkz; k++)
200 if (i != 0 || j != 0 || k != 0)
202 kz = k * k_factor[ZZ];
203 krr = std::sqrt (sqr (kx) + sqr (ky) + sqr (kz));
204 if (krr >= start_q && krr <= end_q)
206 kr = static_cast<int>(krr/sf->ref_k + 0.5);
207 if (kr < sf->n_angles)
209 counter[kr]++; /* will be used for the copmutation
211 for (p = 0; p < isize; p++)
213 asf = sf_table[redt[p].t][kr];
215 kdotx = kx * redt[p].x[XX] +
216 ky * redt[p].x[YY] + kz * redt[p].x[ZZ];
218 tmpSF[i][j][k].re += std::cos(kdotx) * asf;
219 tmpSF[i][j][k].im += std::sin(kdotx) * asf;
226 } /* end loop on i */
228 * compute the square modulus of the structure factor, averaging on the surface
229 * kx*kx + ky*ky + kz*kz = krr*krr
230 * note that this is correct only for a (on the macroscopic scale)
233 for (i = 0; i < maxkx; i++)
235 kx = i * k_factor[XX]; for (j = 0; j < maxky; j++)
237 ky = j * k_factor[YY]; for (k = 0; k < maxkz; k++)
239 kz = k * k_factor[ZZ]; krr = std::sqrt (sqr (kx) + sqr (ky)
240 + sqr (kz)); if (krr >= start_q && krr <= end_q)
242 kr = static_cast<int>(krr / sf->ref_k + 0.5);
243 if (kr < sf->n_angles && counter[kr] != 0)
246 (sqr (tmpSF[i][j][k].re) +
247 sqr (tmpSF[i][j][k].im))/ counter[kr];
260 extern gmx_structurefactors_t *gmx_structurefactors_init(const char *datfn)
263 /* Read the database for the structure factor of the different atoms */
267 gmx_structurefactors *gsf;
268 double a1, a2, a3, a4, b1, b2, b3, b4, c;
278 snew(gsf->atomnm, nralloc);
279 snew(gsf->a, nralloc);
280 snew(gsf->b, nralloc);
281 snew(gsf->c, nralloc);
282 snew(gsf->p, nralloc);
284 gsf->nratoms = line_no;
285 while (get_a_line(fp, line, STRLEN))
288 if (sscanf(line, "%s %d %lf %lf %lf %lf %lf %lf %lf %lf %lf",
289 atomn, &p, &a1, &a2, &a3, &a4, &b1, &b2, &b3, &b4, &c) == 11)
291 gsf->atomnm[i] = gmx_strdup(atomn);
305 gsf->nratoms = line_no;
306 if (line_no == nralloc)
309 srenew(gsf->atomnm, nralloc);
310 srenew(gsf->a, nralloc);
311 srenew(gsf->b, nralloc);
312 srenew(gsf->c, nralloc);
313 srenew(gsf->p, nralloc);
318 fprintf(stderr, "WARNING: Error in file %s at line %d ignored\n",
323 srenew(gsf->atomnm, gsf->nratoms);
324 srenew(gsf->a, gsf->nratoms);
325 srenew(gsf->b, gsf->nratoms);
326 srenew(gsf->c, gsf->nratoms);
327 srenew(gsf->p, gsf->nratoms);
331 return (gmx_structurefactors_t *) gsf;
336 extern void rearrange_atoms (reduced_atom_t * positions, t_trxframe *fr, atom_id * index,
337 int isize, t_topology * top, gmx_bool flag, gmx_structurefactors_t *gsf)
338 /* given the group's index, return the (continuous) array of atoms */
342 reduced_atom *pos = (reduced_atom *)positions;
346 for (i = 0; i < isize; i++)
349 return_atom_type (*(top->atoms.atomname[index[i]]), gsf);
352 for (i = 0; i < isize; i++)
354 copy_rvec (fr->x[index[i]], pos[i].x);
359 extern int return_atom_type (const char *name, gmx_structurefactors_t *gsf)
366 { "CH1", 1 }, { "CH2", 2 }, { "CH3", 3 },
367 { "CS1", 1 }, { "CS2", 2 }, { "CS3", 3 },
368 { "CP1", 1 }, { "CP2", 2 }, { "CP3", 3 }
376 gmx_structurefactors *gsft = (gmx_structurefactors *)gsf;
378 NCMT = gsft->nratoms;
382 for (i = 0; (i < asize(uh)); i++)
384 if (std::strcmp(name, uh[i].name) == 0)
386 return NCMT-1+uh[i].nh;
390 for (i = 0; (i < NCMT); i++)
392 if (std::strncmp (name, gsft->atomnm[i], std::strlen(gsft->atomnm[i])) == 0)
401 gmx_fatal(FARGS, "\nError: atom (%s) not in list (%d types checked)!\n",
407 for (i = 0; i < cnt; i++)
409 if (std::strlen(gsft->atomnm[tndx[i]]) > (size_t)nrc)
411 nrc = std::strlen(gsft->atomnm[tndx[i]]);
422 extern int gmx_structurefactors_get_sf(gmx_structurefactors_t *gsf, int elem, real a[4], real b[4], real *c)
427 gmx_structurefactors *gsft = (gmx_structurefactors *)gsf;
430 for (i = 0; i < 4; i++)
432 a[i] = gsft->a[elem][i];
433 b[i] = gsft->b[elem][i];
441 extern int do_scattering_intensity (const char* fnTPS, const char* fnNDX,
442 const char* fnXVG, const char *fnTRX,
444 real start_q, real end_q,
445 real energy, int ng, const output_env_t oenv)
447 int i, *isize, flags = TRX_READ_X, **index_atp;
449 char **grpname, title[STRLEN];
454 reduced_atom_t **red;
455 structure_factor *sf;
461 gmx_structurefactors_t *gmx_sf;
468 gmx_sf = gmx_structurefactors_init(fnDAT);
470 gmx_structurefactors_get_sf(gmx_sf, 0, a, b, &c);
475 /* Read the topology informations */
476 read_tps_conf (fnTPS, title, &top, &ePBC, &xtop, NULL, box, TRUE);
479 /* groups stuff... */
484 fprintf (stderr, "\nSelect %d group%s\n", ng,
488 get_index (&top.atoms, fnNDX, ng, isize, index, grpname);
492 rd_index (fnNDX, ng, isize, index, grpname);
495 /* The first time we read data is a little special */
496 read_first_frame (oenv, &status, fnTRX, &fr, flags);
498 sf->total_n_atoms = fr.natoms;
501 snew (index_atp, ng);
503 r_tmp = std::max(box[XX][XX], box[YY][YY]);
504 r_tmp = std::max(box[ZZ][ZZ], r_tmp);
506 sf->ref_k = (2.0 * M_PI) / (r_tmp);
507 /* ref_k will be the reference momentum unit */
508 sf->n_angles = static_cast<int>(end_q / sf->ref_k + 0.5);
511 for (i = 0; i < ng; i++)
513 snew (sf->F[i], sf->n_angles);
515 for (i = 0; i < ng; i++)
517 snew (red[i], isize[i]);
518 rearrange_atoms (red[i], &fr, index[i], isize[i], &top, TRUE, gmx_sf);
519 index_atp[i] = create_indexed_atom_type (red[i], isize[i]);
522 sf_table = compute_scattering_factor_table (gmx_sf, (structure_factor_t *)sf);
525 /* This is the main loop over frames */
530 for (i = 0; i < ng; i++)
532 rearrange_atoms (red[i], &fr, index[i], isize[i], &top, FALSE, gmx_sf);
534 compute_structure_factor ((structure_factor_t *)sf, box, red[i], isize[i],
535 start_q, end_q, i, sf_table);
539 while (read_next_frame (oenv, status, &fr));
541 save_data ((structure_factor_t *)sf, fnXVG, ng, start_q, end_q, oenv);
547 gmx_structurefactors_done(gmx_sf);
553 extern void save_data (structure_factor_t *sft, const char *file, int ngrps,
554 real start_q, real end_q, const output_env_t oenv)
559 double *tmp, polarization_factor, A;
561 structure_factor *sf = (structure_factor *)sft;
563 fp = xvgropen (file, "Scattering Intensity", "q (1/nm)",
564 "Intensity (a.u.)", oenv);
568 for (g = 0; g < ngrps; g++)
570 for (i = 0; i < sf->n_angles; i++)
574 * theta is half the angle between incoming and scattered vectors.
576 * polar. fact. = 0.5*(1+cos^2(2*theta)) = 1 - 0.5 * sin^2(2*theta)
578 * sin(theta) = q/(2k) := A -> sin^2(theta) = 4*A^2 (1-A^2) ->
579 * -> 0.5*(1+cos^2(2*theta)) = 1 - 2 A^2 (1-A^2)
581 A = static_cast<double>(i * sf->ref_k) / (2.0 * sf->momentum);
582 polarization_factor = 1 - 2.0 * sqr (A) * (1 - sqr (A));
583 sf->F[g][i] *= polarization_factor;
586 for (i = 0; i < sf->n_angles; i++)
588 if (i * sf->ref_k >= start_q && i * sf->ref_k <= end_q)
590 fprintf (fp, "%10.5f ", i * sf->ref_k);
591 for (g = 0; g < ngrps; g++)
593 fprintf (fp, " %10.5f ", (sf->F[g][i]) /( sf->total_n_atoms*
604 extern double CMSF (gmx_structurefactors_t *gsf, int type, int nh, double lambda, double sin_theta)
606 * return Cromer-Mann fit for the atomic scattering factor:
607 * sin_theta is the sine of half the angle between incoming and scattered
608 * vectors. See g_sq.h for a short description of CM fit.
612 double tmp = 0.0, k2;
621 * f0[k] = c + [SUM a_i*EXP(-b_i*(k^2)) ]
631 tmp = (CMSF (gsf, return_atom_type ("C", gsf), 0, lambda, sin_theta) +
632 nh*CMSF (gsf, return_atom_type ("H", gsf), 0, lambda, sin_theta));
637 k2 = (sqr (sin_theta) / sqr (10.0 * lambda));
638 gmx_structurefactors_get_sf(gsf, type, a, b, &c);
640 for (i = 0; (i < 4); i++)
642 tmp += a[i] * exp (-b[i] * k2);
650 extern real **gmx_structurefactors_table(gmx_structurefactors_t *gsf, real momentum, real ref_k, real lambda, int n_angles)
658 gmx_structurefactors *gsft = (gmx_structurefactors *)gsf;
660 NCMT = gsft->nratoms;
663 snew (sf_table, nsftable);
664 for (i = 0; (i < nsftable); i++)
666 snew (sf_table[i], n_angles);
667 for (j = 0; j < n_angles; j++)
669 q = static_cast<double>(j * ref_k);
670 /* theta is half the angle between incoming
671 and scattered wavevectors. */
672 sin_theta = q / (2.0 * momentum);
675 sf_table[i][j] = CMSF (gsf, i, 0, lambda, sin_theta);
679 sf_table[i][j] = CMSF (gsf, i, i-NCMT+1, lambda, sin_theta);
686 extern void gmx_structurefactors_done(gmx_structurefactors_t *gsf)
690 gmx_structurefactors *sf;
691 sf = (gmx_structurefactors *) gsf;
693 for (i = 0; i < sf->nratoms; i++)
697 sfree(sf->atomnm[i]);
710 extern real **compute_scattering_factor_table (gmx_structurefactors_t *gsf, structure_factor_t *sft)
713 * this function build up a table of scattering factors for every atom
714 * type and for every scattering angle.
717 double hc = 1239.842;
720 structure_factor *sf = (structure_factor *)sft;
723 /* \hbar \omega \lambda = hc = 1239.842 eV * nm */
724 sf->momentum = (static_cast<double>(2. * 1000.0 * M_PI * sf->energy) / hc);
725 sf->lambda = hc / (1000.0 * sf->energy);
726 fprintf (stderr, "\nwavelenght = %f nm\n", sf->lambda);
728 sf_table = gmx_structurefactors_table(gsf, sf->momentum, sf->ref_k, sf->lambda, sf->n_angles);