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
);