Structure factor and X-ray scattering from radial distribution functions
Project description
Structure factor and scattering from radial distribution functions
A package to calculate structure factors and x-ray (solution) scattering signals from Radial Distribution Functions (RDFs), sampled from Molecular Dynamics (MD) simulations.
Also includes 3 different finite-size RDF-corrections, as well as a set of Fourier truncation window functions.
Please see (and cite) this work for the necessary background. There is also a pre-print freely available here.
Please read the docs for tutorials on how to use the tools included.
If you want even more details and inspiration, the the example notebook contains more usage-examples, and the this notebook contains all the scripts used to make the plots in the paper.
Example: X-ray Scattering of Water:
from grsq.grsq import rdfset_from_dir
V = 122900.85207774633 # volume of MD box
stoich = {'H_v': 8190, 'O_v': 4095} # number of atoms in MD box
# Helper function to create RDFSet
rdfs = rdfset_from_dir('tests/data/xray_water/',
volume=V, stoich=stoich)
qvec = np.arange(0, 10, 0.05) # create q vector
ig_u = rdfs.get_iq(qvec) # Calculate scattering
rdfs.vdv_correct() # apply finite-size corrections
ig_c = rdfs.get_iq(qvec) # Calculate scattering again
fig, ax = plt.subplots(1, 1, figsize=(9, 5))
ax.plot(qvec, ig_u, label='Raw')
ax.plot(qvec, ig_c, label='Corrected')
ax.set_xlim([0, 10])
ax.set_xlabel('Q (Å$^{-1})$')
ax.set_ylabel('I(Q)')
ax.legend(loc='best');
fig.tight_layout()
{width=85%}
Installation
The package is available from PyPi:
pip install grsq
Or clone this repository and add it to your $PYTHONPATH
.
Introduction
For each pairwise RDF you have sampled from your system, you should create an RDF
object, in this example between the Oxygens and Hydrogens of the solvent-region:
dat = np.genfromtxt('some_data_file.dat') # has (r,g(r)) as columns in this example
rdf = RDF(dat[:, 0], dat[:, 1], # r and g(r)
'O', 'H', # atom types
'solvent', 'solvent', # atom type regions
n1=4095, n2=8190, # number of atoms of each type
volume=122900.85, # volume in Å of simulation cell
qvec=qvec)
There are 2 possible regions for your atom types to exist in: 'solute' and 'solvent'.
From a single RDF
object, you can calculate the structure factor:
rdf.structure_factor()
as well as the first (atomic), and 2nd term of the x-ray scattering signal:
rdf.term1() + rdf.term2()
You most likely have a large amount of pairwise RDFs from a system of a solute comprised of various elements in a solvent. To make things simple, we can put them all into an RDFSet
instance. How to do so most easily depends on how you store your RDF data, naming, data format, etc...
If you follow the naming convention in the documentation of the helper function rdfset_from_dir
, this function can help you do so, as in the water example above.
The RDFSet Class
An OrderedDict of RDF objects, which you have to create yourself, you might e.g. already have created a list_of_rdf_objects
:
from grsq.grsq import RDFSet
rdfs = RDFSet()
for rdf in list_of_rdf_objects:
rdfs.add_rdf(rdf)
Or have a look at:
from grsq.grsq import rdfset_from_dir
?rdfset_from_dir
Having an RDFSet
is very convenient, as you get access to a bunch of methods:
rdfs.show()
: Returns information of the setrdfs.get_iq(qvec=None, cross=False, damping=None)
Get the total coherent x-ray scattering signal from entire set. You can specifyqvec
andDamping
to overwrite whatever is already contained in each RDF object.rdfs.get_solute()
: Returns the scattering of only the solute-solute atom type RDFsrdfs.get_solvent()
: Returns the scattering of only the solvent-solvent atom type RDFsrdfs.get_cross()
: Returns the scattering of only the solute-solvent and solvent-solute atom type RDFsrdfs.get_dv()
: Returns the excluded volume-scattering
All methods take qvec
and damping
as arguments. To use fourier damping functions, parse a damping
object to the method:
from grsq.damping import Damping
damping = Damping('lorch', L=r_max)
rdfs.get_solvent(qvec=qvec, damping=damping)
Finite Size corrections
The following finite size corrections are available (see this work for details):
-
rdf.volume_correct(Rl)
: Apply the volume correction to the RDF, using a spherical $V_l$ volume of radiusRl
: $$ g_{lm}^\mathrm{\infty}(r) = \frac{\rho_m}{\rho_\mathrm{eff}} g_{lm}^N(r) = g_{lm}^N(r)\rho_m \frac{V_\mathrm{cell} - V_l}{N_m - \delta_{lm}}, $$ -
rdf.perera_correct(Rl, r_avg)
: Apply the correction from Perera et al. $$ g^\mathrm{\infty}{lm}(r) = g^N{lm}(r) \left [ 1 + \frac{1 - g^{N,0}}{2} \left( 1 + \tanh \left( \frac{r - \kappa_{lm}}{\alpha_{lm}} \right) \right) \right ]. $$ Where $g^{N,0}$ is obtained by settingr_avg
to get the average value of g(r > r_avg), andRl
= $\kappa_{lm} / 2$, to be consistent with the volume correction, as well as the recommendations from the original paper. -
rdf.vdv_correct()
: Apply the Ganguly / van der Vegt correction: $$ g^\mathrm{\infty}\mathrm{lm}(r) = g^N\mathrm{lm}(r) \frac{N_m}{N_m - \left[ (\Delta N(r) + \delta_{lm}) \left( 1 - \frac{(4/3) \pi r^3}{V_\mathrm{cell}} \right)^{-1} \right]} $$
Support methods:
rdf.fit(Rl_guess, fit_start=25, fit_stop=50)
: Fit the excluded volume radius.
Can be used together withrdf.volume_correct(Rl)
to estimateRl
. It minimizes g(r) - 1 fromfit_start
tofit_stop
Try ?RDF
or help(RDF)
for more info.
The RDFSet
also has the method rdfs.vdv_correct()
that apply the van der Vegt correction to all relevant RDFs automatically.
Copyright
Copyright (c) 2022, Asmus Ougaard Dohn
Project details
Release history Release notifications | RSS feed
Download files
Download the file for your platform. If you're not sure which to choose, learn more about installing packages.