grsq package

grsq.grsq

Contains the main classes to calculate structure factors and X-ray Scattering signals from pair RDFs, as well as correct for finite-size-sampled RDFs.

Each RDF is treated as an object RDF, which contains all the neccessary attributes of an RDF sampled from MD simulations.

All RDFs of a given system can then be collected into an RDFSet, which can calculate the full scattering signal, from the entire subsystem, the scattering from the solute-atoms only, the solvent-atoms only, or the cross term comprised of all solute-solvent pair-RDFs. The class can also do finite-size corrections accross all (solvent-solvent and solute-solvent) RDFs as well.

class grsq.grsq.RDF(r, g, name1, name2, region1, region2, n1=None, n2=None, damp=None, volume=None, qvec=None, rho_norm='KD', r_max=None, r_avg=None)

Bases: object

Class for representing a pair-RDF sampled between two atom types.

r: np.ndarray (1, ): r vector: The r vector defining the bins in which the RDF was sampled, in Å.
g: np.ndarray, same dim as r: g(r) values of each r bin
name1: str: Name of atomtype/element of ‘left’ atoms
name2: str: Name of atomtype/element of ‘right’ atoms
region1: str, ‘solute’ or ‘solvent’: Which region does the ‘left’ atom belong to
region2: str, ‘solute’, or ‘solvent’: Which region does the ‘right’ atom belong to.
n1: int: Number of atoms of this type
n2: int: Number of atoms of this type
damp: grsg.damping.Damping: An instance of the grsq.damping.Damping class that provides the damping function used in the RDF calculations. See the documentation for Damping for more information on the methods and attributes available in the Damping class.
qvec: np.ndarray (1, ): Scattering vector, in Å^-1
rho_norm: str, ‘KD’ or ‘N’: Specify which bulk density was used in in obtaining g(r)

‘(N - KD_ij) / V’: As in e.g. VMD

‘N/V’: As in e.g. MDAnalysis
r_max: float: Specify at which r to stop the integration
r_avg: float: Use avg. g(r > r_avg) for g0

atomic_f0(custom_cm=False): Calculate atomic form factors and update RDF object in place

coordination_number(r_min, r_max): Calculate the coordination number in the region: r_min >= r > r_max

correct(method='volume', **kwargs)

Correct for finite size sampling.

See 10.1063/5.0164365 for an overview, but please cite the original papers of the method(s) you use also.

method: str: ‘volume’, ‘perera’, or ‘vegt’: volume: Uses the volume-correction. See: 10.1063/5.0164365

perera: See 10.1016/j.molliq.2010.05.006

vegt: 10.1021/ct301017q

fit(Ri_guess, fit_start=25, fit_stop=50): Fit the exluded volume by minimizing g(r) - 1 in the range from fit_start to fit_stop

get_info(): Show contents of RDF object

i_term1(custom_cm=False): Coherent X-Ray Scattering, atomic term

i_term2(custom_cm=False): Coherent X-Ray Scattering, interatomic term

perera_correct(Ri, r_avg)

See 10.1016/j.molliq.2010.05.006

Ri: float: Radius of exluced volume - paper recommends kappa is 2 x this.
r_avg: float: Use avg g(r > r_avg)

structure_factor(): Calculate the structure factor for a given RDF

vdv_correct(): See 10.1080/08927022.2017.1416114 For best explanation. Only implemented for g(r)’s sampled to r_max <= 1/2 x boxsize

volume_correct(Ri)

class grsq.grsq.RDFSet(*arg, **kwargs)

Bases: OrderedDict

Container object for a set of RDFs.

For more information, see the Tutorials, Main Tools page in the project documentation.

Can be used to calculate the solute-solute, cross, and solvent-solvent scattering terms. Given an rdfs = RDFSet(...)

>>> qvec = np.arange(0, 10, 0.01)  # define your q range
>>> damp = Damping('lorch', L=25)  # use a damping window
>>> I_v = rdfs.get_solvent(qvec=qvec, damp=damp)
>>> I_c = rdfs.get_cross(qvec=qvec, damp=damp)
>>> I_u = rdfs.get_solute(qvec=qvec)  # no damping for solute

You can also apply the van der Vegt finite size correction to all (cross and solvent) RDFs:

>>> rdfs.vdv_correct()

And calculate the displaced-volume scattering:

>>> rdfs.get_dv(qvec=qvec, damp=damp)

add_flipped(rdf)

add_rdf(rdf)

get_cross(qvec=None, damping=None, custom_cm=False): Get cross-term from subset of RDFSet. It will automatically add Y-X to a set that only has X-Y. Which is what you normally want.

get_dv(qvec=None, alpha=1, damping=None, custom_cm=False)

Displaced Volume Scattering

Make f_dv form factors for all solute types.
Calculate total S_dv

get_iq(qvec=None, cross=False, damping=None, custom_cm=False): Calculate Coherent Scattering from set

get_solute(qvec=None, damping=None, custom_cm=False)

Get solute-solute scattering from subset of RDFSet It will automatically add Y-X to a set that only has X-Y. Which is what you normally want.

The option to use damping is included for consistency, but be aware that it does not make much sense to use it for the solute-solute term, as all RDFs should go to zero at long r values.

get_solvent(qvec=None, damping=None, custom_cm=False): Get solvent-only from subset of RDFSet It will automatically add Y-X to a set that only has X-Y. Which is what you normally want.

show(): Show contents of RDFSet.

vdv_correct(): Apply van der Vegt correction to all solute-solvent and solvent-solvent rdfs

grsq.grsq.rdfset_from_dir(dirname, prmtop=None, volume=None, stoich=None)

Read in a directory of rdf files with the naming convention:: gEl1_w1-El2_w2.dat
where:: El1: Element name 1 (also used for the IAM)

El2: Element name 2 (also used for the IAM)

w1: Region of element 1 (u, v) (i.e. solUte, solVent)

w2: Region of element 2 (u, v)

So e.g. water would have the RDFs: gO_v-O_v.dat, gO_v-H_v.dat, etc..
prmtop (str): path to structure file for a single frame:: The method can read the stoichometry from either a prmtop or an xyzfile.

If an xyz file is used, you need to specify the volume in which the RDFs

where sampled manually.

volume (float): RDF sampling volume. Only use if xyz is specified in prmtop

stoich (dict):: {‘El1_w1’:N1, ‘El2_w2:N2 …}

Dictionary containing number of atoms of each element/atom type.

Use if no prmtop/xyz is used at all

grsq.damping module

class grsq.damping.Damping(style, L=None, r_max=None, r_cut=None)

Bases: object

Class that applies a damping window function to ( g(r) - 1 ) in the generalized debye integral.

Instantiate the object and add it to your RDF(), e.g.

>>> L = 25  # for a 50 Å cell
>>> damp = Damping('lorch', L=L)
>>> rdf = RDF(..., damp=damp)

You can also instruct an RDFSet to use the same damping for each RDF to calculate the full scattering signal

>>> rdfs = RDFSet(...)
>>> damp = Damping('lorch', L=25)
>>> rdfs.get_iq(damping=damp)

Parameters:

style (str): can be either of these:: ‘lorch’: sin(pi*L) / pi*L

‘dhabal’: as in 10.1039/c6cp07599a

‘zederkof’: as above but with an r_cut as well. 10.1021/jacs.2c04505

‘panman’: Exponential fade function, e.g. PhysRevLett.125.226001
L: Float (for Lorch damping): Damping factor for the lorch damping, usually 1/2 of

the box length.
r_max: Float (for Dhabal, Zederkof, Panman): The r-value where the rdf should be 0
r_cut: Float (only for zederkof damping): When to start the smooth cutoff. Everything before

is untouched.

damp(R): Apply the chosen damping

dhabal(R): Damping like in 10.1039/c6cp07599a

lorch(R): Lorch-like damping

panman(R): Damping like in PhysRevLett.125.226001

zederkof(R): Damping like in 10.1021/jacs.2c04505

grsq.debye module

class grsq.debye.Debye(qvec=None, custom_cm=None)

Bases: object

Discrete Interatomic-Distance Debye Scattering.

      ___N  ___N
      \    \                       sin(q · r_ij)
I(q) = )    )   f_i(q) · f_j(q) · ---------------
      /__  /__                       q · r_ij
      i=1  j=1

where:

I(q)    : Scattering intensity as a function of the momentum transfer q.
N       : Total number of atoms.
f_i, f_j: Atomic form actors for atoms i and j, respectively.
r_ij    : Distance between atoms i and j.
q       : Momentum transfer, or scattering vector.

The atomic form factors are generated from the Cromer Mann parameters from Waasmaier, A. Kirfel Acta Cryst. A51, 416-431 (1995), using the python library periodictable NB: This assumes neutral atoms. See the tutorials for how to use parameters for charged atoms.

Usage:

>>> from grsq.debye import Debye
>>> from ase.io import read
>>> atoms = read('./tests/data/xray_2particle/test.xyz')  # read e.g. xyz file as ASE atoms object
>>> deb = Debye(qvec=np.arange(0, 20, 0.01))
>>> i_deb = deb.debye(atoms)

Parameters:

qvec: np.ndarray (1, )

your Q vector

custom_cm: None or dict

Will use either the standard Waasmeier & Kirfel Cromer-Mann parameters from periodictable or you can enter your own dictionary of CM parameters in the following format:

>>> cm = {'Pt': {'a': np.array([27.00590, 17.76390, 15.71310, 5.783700]),
>>>              'b': np.array([1.512930, 8.811740, 0.4245930, 38.61030]),
>>>              'c': 11.68830}}

atomic_f0(atom)

Calculate atomic form factor for a single atom type.

atom: str: Name of the atom type

debye(atoms, custom_elements=None)

Debye Scattering from ASE atoms object

atoms: ase.Atoms object: atoms to calculate the scattering from
custom_elements: list: list of strings defining the atom types in your Atoms object. If None, the elements from the Atoms object will be used.

debye_histogrammed(atoms, dr=0.05, use_numba=False)

Histogrammed Debye Approximation.

Parameters:

atoms: (ASE Atoms object): The atoms to calculate the scattering from
dr: (float): Histogram bin size. Smaller -> slower, more accurate
numba: (bool): True: Calculate distances with parallal numba implementation, or False: with simple numpy broadcasting. Numba might speed up _very_ large calculations slightly.

debye_numba(atoms, custom_elements=None): Faster implementation, using numba. Useful if you need to do a _lot_ of evaluations.

debye_selective(atoms, idx1, idx2, custom_elements=None): Only terms between atoms in lists idx1 and idx2 (but not between atoms in either list)

update_f0(atoms, custom_elements=None): Update atomic form factors in the object

grsq.debye.get_cm(species)

grsq.debye.populate_cm(): Helper method to get all Cromer Mann parameters from the periodictable library

grsq.accel module

grsq.accel.get_dists_numba(pos_l, pos_m): Get interatomic distances between l and m atoms Using numba’s parallelisation.

grsq.accel.numbye(qvec, natoms, f0, rs, s): Used by Debye().debye_numba()

Module contents

Structure factor and scattering from radial distribution functions