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,2016,2017 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source 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.
47 #include "gromacs/fileio/confio.h"
48 #include "gromacs/fileio/trxio.h"
49 #include "gromacs/fileio/xvgr.h"
50 #include "gromacs/math/functions.h"
51 #include "gromacs/math/utilities.h"
52 #include "gromacs/math/vec.h"
53 #include "gromacs/mdtypes/md_enums.h"
54 #include "gromacs/topology/index.h"
55 #include "gromacs/topology/topology.h"
56 #include "gromacs/trajectory/trajectoryframe.h"
57 #include "gromacs/utility/arraysize.h"
58 #include "gromacs/utility/cstringutil.h"
59 #include "gromacs/utility/fatalerror.h"
60 #include "gromacs/utility/futil.h"
61 #include "gromacs/utility/smalloc.h"
62 #include "gromacs/utility/strdb.h"
65 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
85 typedef struct structure_factor
99 extern int* create_indexed_atom_type(reduced_atom_t* atm, int size)
102 * create an index of the atom types found in a group
103 * i.e.: for water index_atp[0]=type_number_of_O and
104 * index_atp[1]=type_number_of_H
106 * the last element is set to 0
108 int *index_atp, i, i_tmp, j;
110 reduced_atom* att = static_cast<reduced_atom*>(atm);
114 index_atp[0] = att[0].t;
115 for (i = 1; i < size; i++)
117 for (j = 0; j < i_tmp; j++)
119 if (att[i].t == index_atp[j])
124 if (j == i_tmp) /* i.e. no indexed atom type is == to atm[i].t */
127 srenew(index_atp, i_tmp * sizeof(int));
128 index_atp[i_tmp - 1] = att[i].t;
132 srenew(index_atp, i_tmp * sizeof(int));
133 index_atp[i_tmp - 1] = 0;
138 extern t_complex*** rc_tensor_allocation(int x, int y, int z)
145 snew(t[0][0], x * y * z);
147 for (j = 1; j < y; j++)
149 t[0][j] = t[0][j - 1] + z;
151 for (i = 1; i < x; i++)
154 t[i][0] = t[i - 1][0] + y * z;
155 for (j = 1; j < y; j++)
157 t[i][j] = t[i][j - 1] + z;
164 extern void compute_structure_factor(structure_factor_t* sft,
173 structure_factor* sf = static_cast<structure_factor*>(sft);
174 reduced_atom* redt = static_cast<reduced_atom*>(red);
178 real kdotx, asf, kx, ky, kz, krr;
179 int kr, maxkx, maxky, maxkz, i, j, k, p, *counter;
182 k_factor[XX] = 2 * M_PI / box[XX][XX];
183 k_factor[YY] = 2 * M_PI / box[YY][YY];
184 k_factor[ZZ] = 2 * M_PI / box[ZZ][ZZ];
186 maxkx = gmx::roundToInt(end_q / k_factor[XX]);
187 maxky = gmx::roundToInt(end_q / k_factor[YY]);
188 maxkz = gmx::roundToInt(end_q / k_factor[ZZ]);
190 snew(counter, sf->n_angles);
192 tmpSF = rc_tensor_allocation(maxkx, maxky, maxkz);
195 * compute real and imaginary part of the structure factor for every
198 fprintf(stderr, "\n");
199 for (i = 0; i < maxkx; i++)
201 fprintf(stderr, "\rdone %3.1f%% ", (100.0 * (i + 1)) / maxkx);
203 kx = i * k_factor[XX];
204 for (j = 0; j < maxky; j++)
206 ky = j * k_factor[YY];
207 for (k = 0; k < maxkz; k++)
209 if (i != 0 || j != 0 || k != 0)
211 kz = k * k_factor[ZZ];
212 krr = std::sqrt(gmx::square(kx) + gmx::square(ky) + gmx::square(kz));
213 if (krr >= start_q && krr <= end_q)
215 kr = gmx::roundToInt(krr / sf->ref_k);
216 if (kr < sf->n_angles)
218 counter[kr]++; /* will be used for the copmutation
220 for (p = 0; p < isize; p++)
222 asf = sf_table[redt[p].t][kr];
224 kdotx = kx * redt[p].x[XX] + ky * redt[p].x[YY] + kz * redt[p].x[ZZ];
226 tmpSF[i][j][k].re += std::cos(kdotx) * asf;
227 tmpSF[i][j][k].im += std::sin(kdotx) * asf;
234 } /* end loop on i */
236 * compute the square modulus of the structure factor, averaging on the surface
237 * kx*kx + ky*ky + kz*kz = krr*krr
238 * note that this is correct only for a (on the macroscopic scale)
241 for (i = 0; i < maxkx; i++)
243 kx = i * k_factor[XX];
244 for (j = 0; j < maxky; j++)
246 ky = j * k_factor[YY];
247 for (k = 0; k < maxkz; k++)
249 kz = k * k_factor[ZZ];
250 krr = std::sqrt(gmx::square(kx) + gmx::square(ky) + gmx::square(kz));
251 if (krr >= start_q && krr <= end_q)
253 kr = gmx::roundToInt(krr / sf->ref_k);
254 if (kr < sf->n_angles && counter[kr] != 0)
257 (gmx::square(tmpSF[i][j][k].re) + gmx::square(tmpSF[i][j][k].im))
271 extern gmx_structurefactors_t* gmx_structurefactors_init(const char* datfn)
274 /* Read the database for the structure factor of the different atoms */
277 gmx_structurefactors* gsf;
278 double a1, a2, a3, a4, b1, b2, b3, b4, c;
284 gmx::FilePtr fp = gmx::openLibraryFile(datfn);
288 snew(gsf->atomnm, nralloc);
289 snew(gsf->a, nralloc);
290 snew(gsf->b, nralloc);
291 snew(gsf->c, nralloc);
292 snew(gsf->p, nralloc);
294 gsf->nratoms = line_no;
295 while (get_a_line(fp.get(), line, STRLEN))
298 if (sscanf(line, "%s %d %lf %lf %lf %lf %lf %lf %lf %lf %lf", atomn, &p, &a1, &a2, &a3, &a4,
299 &b1, &b2, &b3, &b4, &c)
302 gsf->atomnm[i] = gmx_strdup(atomn);
316 gsf->nratoms = line_no;
317 if (line_no == nralloc)
320 srenew(gsf->atomnm, nralloc);
321 srenew(gsf->a, nralloc);
322 srenew(gsf->b, nralloc);
323 srenew(gsf->c, nralloc);
324 srenew(gsf->p, nralloc);
329 fprintf(stderr, "WARNING: Error in file %s at line %d ignored\n", datfn, line_no);
333 srenew(gsf->atomnm, gsf->nratoms);
334 srenew(gsf->a, gsf->nratoms);
335 srenew(gsf->b, gsf->nratoms);
336 srenew(gsf->c, gsf->nratoms);
337 srenew(gsf->p, gsf->nratoms);
339 return static_cast<gmx_structurefactors_t*>(gsf);
343 extern void rearrange_atoms(reduced_atom_t* positions,
347 const t_topology* top,
349 gmx_structurefactors_t* gsf)
350 /* given the group's index, return the (continuous) array of atoms */
354 reduced_atom* pos = static_cast<reduced_atom*>(positions);
358 for (i = 0; i < isize; i++)
360 pos[i].t = return_atom_type(*(top->atoms.atomname[index[i]]), gsf);
363 for (i = 0; i < isize; i++)
365 copy_rvec(fr->x[index[i]], pos[i].x);
370 extern int return_atom_type(const char* name, gmx_structurefactors_t* gsf)
377 t_united_h uh[] = { { "CH1", 1 }, { "CH2", 2 }, { "CH3", 3 }, { "CS1", 1 }, { "CS2", 2 },
378 { "CS3", 3 }, { "CP1", 1 }, { "CP2", 2 }, { "CP3", 3 } };
385 gmx_structurefactors* gsft = static_cast<gmx_structurefactors*>(gsf);
387 NCMT = gsft->nratoms;
391 for (i = 0; (i < asize(uh)); i++)
393 if (std::strcmp(name, uh[i].name) == 0)
395 return NCMT - 1 + uh[i].nh;
399 for (i = 0; (i < NCMT); i++)
401 if (std::strncmp(name, gsft->atomnm[i], std::strlen(gsft->atomnm[i])) == 0)
410 gmx_fatal(FARGS, "\nError: atom (%s) not in list (%d types checked)!\n", name, i);
415 for (i = 0; i < cnt; i++)
417 if (std::strlen(gsft->atomnm[tndx[i]]) > static_cast<size_t>(nrc))
419 nrc = std::strlen(gsft->atomnm[tndx[i]]);
428 extern int gmx_structurefactors_get_sf(gmx_structurefactors_t* gsf, int elem, real a[4], real b[4], real* c)
433 gmx_structurefactors* gsft = static_cast<gmx_structurefactors*>(gsf);
436 for (i = 0; i < 4; i++)
438 a[i] = gsft->a[elem][i];
439 b[i] = gsft->b[elem][i];
447 extern int do_scattering_intensity(const char* fnTPS,
456 const gmx_output_env_t* oenv)
458 int i, *isize, flags = TRX_READ_X, **index_atp;
465 reduced_atom_t** red;
466 structure_factor* sf;
472 gmx_structurefactors_t* gmx_sf;
479 gmx_sf = gmx_structurefactors_init(fnDAT);
481 gmx_structurefactors_get_sf(gmx_sf, 0, a, b, &c);
486 /* Read the topology informations */
487 read_tps_conf(fnTPS, &top, &ePBC, &xtop, nullptr, box, TRUE);
490 /* groups stuff... */
495 fprintf(stderr, "\nSelect %d group%s\n", ng, ng == 1 ? "" : "s");
498 get_index(&top.atoms, fnNDX, ng, isize, index, grpname);
502 rd_index(fnNDX, ng, isize, index, grpname);
505 /* The first time we read data is a little special */
506 read_first_frame(oenv, &status, fnTRX, &fr, flags);
508 sf->total_n_atoms = fr.natoms;
513 r_tmp = std::max(box[XX][XX], box[YY][YY]);
514 r_tmp = std::max(box[ZZ][ZZ], r_tmp);
516 sf->ref_k = (2.0 * M_PI) / (r_tmp);
517 /* ref_k will be the reference momentum unit */
518 sf->n_angles = gmx::roundToInt(end_q / sf->ref_k);
521 for (i = 0; i < ng; i++)
523 snew(sf->F[i], sf->n_angles);
525 for (i = 0; i < ng; i++)
527 snew(red[i], isize[i]);
528 rearrange_atoms(red[i], &fr, index[i], isize[i], &top, TRUE, gmx_sf);
529 index_atp[i] = create_indexed_atom_type(red[i], isize[i]);
532 sf_table = compute_scattering_factor_table(gmx_sf, static_cast<structure_factor_t*>(sf));
535 /* This is the main loop over frames */
540 for (i = 0; i < ng; i++)
542 rearrange_atoms(red[i], &fr, index[i], isize[i], &top, FALSE, gmx_sf);
544 compute_structure_factor(static_cast<structure_factor_t*>(sf), box, red[i], isize[i],
545 start_q, end_q, i, sf_table);
549 while (read_next_frame(oenv, status, &fr));
551 save_data(static_cast<structure_factor_t*>(sf), fnXVG, ng, start_q, end_q, oenv);
557 gmx_structurefactors_done(gmx_sf);
563 extern void save_data(structure_factor_t* sft,
568 const gmx_output_env_t* oenv)
573 double *tmp, polarization_factor, A;
575 structure_factor* sf = static_cast<structure_factor*>(sft);
577 fp = xvgropen(file, "Scattering Intensity", "q (1/nm)", "Intensity (a.u.)", oenv);
581 for (g = 0; g < ngrps; g++)
583 for (i = 0; i < sf->n_angles; i++)
587 * theta is half the angle between incoming and scattered vectors.
589 * polar. fact. = 0.5*(1+cos^2(2*theta)) = 1 - 0.5 * sin^2(2*theta)
591 * sin(theta) = q/(2k) := A -> sin^2(theta) = 4*A^2 (1-A^2) ->
592 * -> 0.5*(1+cos^2(2*theta)) = 1 - 2 A^2 (1-A^2)
594 A = static_cast<double>(i * sf->ref_k) / (2.0 * sf->momentum);
595 polarization_factor = 1 - 2.0 * gmx::square(A) * (1 - gmx::square(A));
596 sf->F[g][i] *= polarization_factor;
599 for (i = 0; i < sf->n_angles; i++)
601 if (i * sf->ref_k >= start_q && i * sf->ref_k <= end_q)
603 fprintf(fp, "%10.5f ", i * sf->ref_k);
604 for (g = 0; g < ngrps; g++)
606 fprintf(fp, " %10.5f ", (sf->F[g][i]) / (sf->total_n_atoms * sf->nSteps));
616 extern double CMSF(gmx_structurefactors_t* gsf, int type, int nh, double lambda, double sin_theta)
618 * return Cromer-Mann fit for the atomic scattering factor:
619 * sin_theta is the sine of half the angle between incoming and scattered
620 * vectors. See g_sq.h for a short description of CM fit.
624 double tmp = 0.0, k2;
633 * f0[k] = c + [SUM a_i*EXP(-b_i*(k^2)) ]
643 tmp = (CMSF(gsf, return_atom_type("C", gsf), 0, lambda, sin_theta)
644 + nh * CMSF(gsf, return_atom_type("H", gsf), 0, lambda, sin_theta));
649 k2 = (gmx::square(sin_theta) / gmx::square(10.0 * lambda));
650 gmx_structurefactors_get_sf(gsf, type, a, b, &c);
652 for (i = 0; (i < 4); i++)
654 tmp += a[i] * exp(-b[i] * k2);
661 extern real** gmx_structurefactors_table(gmx_structurefactors_t* gsf, real momentum, real ref_k, real lambda, int n_angles)
669 gmx_structurefactors* gsft = static_cast<gmx_structurefactors*>(gsf);
671 NCMT = gsft->nratoms;
674 snew(sf_table, nsftable);
675 for (i = 0; (i < nsftable); i++)
677 snew(sf_table[i], n_angles);
678 for (j = 0; j < n_angles; j++)
680 q = static_cast<double>(j * ref_k);
681 /* theta is half the angle between incoming
682 and scattered wavevectors. */
683 sin_theta = q / (2.0 * momentum);
686 sf_table[i][j] = CMSF(gsf, i, 0, lambda, sin_theta);
690 sf_table[i][j] = CMSF(gsf, i, i - NCMT + 1, lambda, sin_theta);
697 extern void gmx_structurefactors_done(gmx_structurefactors_t* gsf)
701 gmx_structurefactors* sf;
702 sf = static_cast<gmx_structurefactors*>(gsf);
704 for (i = 0; i < sf->nratoms; i++)
708 sfree(sf->atomnm[i]);
720 extern real** compute_scattering_factor_table(gmx_structurefactors_t* gsf, structure_factor_t* sft)
723 * this function build up a table of scattering factors for every atom
724 * type and for every scattering angle.
727 double hc = 1239.842;
730 structure_factor* sf = static_cast<structure_factor*>(sft);
733 /* \hbar \omega \lambda = hc = 1239.842 eV * nm */
734 sf->momentum = (static_cast<double>(2. * 1000.0 * M_PI * sf->energy) / hc);
735 sf->lambda = hc / (1000.0 * sf->energy);
736 fprintf(stderr, "\nwavelenght = %f nm\n", sf->lambda);
738 sf_table = gmx_structurefactors_table(gsf, sf->momentum, sf->ref_k, sf->lambda, sf->n_angles);