.. _tutorials-label: ===================== Tutorials, Main Tools ===================== This section is meant to help you familiarize yourself with the workflow of the main grsq tools. The central hurdle is usually to make sure your RDF-data is correctly read in by ``grsq``. This might seem a bit cumbersome, but once you get the hang of it, it becomes really fast and easy to test out all the various corrections, calculate scattering between sub-sets of RDFs, scan damping parameters and so on. It's worth it, promise! For solute-solvent MD simulations, it is often useful to sub-partition all of your sampled pair-RDFs into whether: 1. Both atom types sampled are in the solute (the 'solute' term, giving the scattering :math:`I_{uu}(q)`) 2. Both atom types sampled are in the solvent (the 'solvent' term, :math:`I_{vv}(q)`) 3. One of the atom types is in the solute, and the other is in the solvent (the 'cross' terms :math:`I_{uv}(q)` and :math:`I_{vu}(q)`) Where the 'atom type' is defined as 'an atom with the same form factor' (often simply an element). The total scattering signal is thus a sum of these terms: .. math:: I(q) = \sum_{w_l}^{\{u,v\}} \sum_{w_m}^{\{u,v\}} I_{w_lw_m}(q) = I_{uu}(q) + I_{uv}(q) + I_{vu}(q) + I_{vv}(q). If we look at the central equation that we need to evaluate to get to an x-ray scattering signal from a set of RDFs (i.e. the generalized Debye equation, see this paper_ for details): .. _paper: https://doi.org/10.1063/5.0164365 .. _eq-grsq: .. math:: :label: eq-grsq :nowrap: \begin{align*} I_{w_lw_m}(q) = & \delta_{w_lw_m} \sum_{l} N_l f_l(q)^2 \\ & + \sum_{l} \sum_{m} f_l(q) f_m(q) \frac{N_l (N_m - \delta_{lm})}{V} \\ & 4 \pi \int_0^R (g_{lm}(r) - g^0_{w_lw_m}) \frac{\sin(qr)}{qr} r^2 \, dr, \end{align*} where :math:`f_l(q)` and :math:`f_m(q)` are the atomic form factors, :math:`\delta_{w_lw_m}` is the Kronecker delta, :math:`N_l` and :math:`N_m` are the numbers of atoms of type :math:`l` and :math:`m` respectively, :math:`g_{lm}(r)` is the RDF beween atom types :math:`l` and :math:`m`, and :math:`V` is the volume of your MD simulation cell. :math:`g^0_{w_lw_m} = 0` when :math:`w_l= u` and :math:`w_m = u`, and 1 otherwise. Lastly, :math:`R` is the :math:`r`-distance at which you want to truncate your integral. The exact equation goes to infinity, but your sampled RDFs dont. By default, `grsq` will integrate until the end of your sampled RDFs. We will return to this in one of the tutorials. Based on this, we can gather that the information you need to calculate the scattering from your set of RDFs is: 1. The g(r) sampled from your MD trajectory(/ies). This can be done efficiently with e.g. VMD_, MDTraj_, or MDAnalysis_ 2. The numbers of each atom type 3. The volume of your simulation cell The atomic form factors are most often tabulated within the Independent Atom Model (IAM), which grsg uses the python package periodictable_ to retrieve. So unless you need to specifically modify them, you don't need to think about getting those. That's it! .. _VMD: https://www.ks.uiuc.edu/Research/vmd/ .. _MDTraj: https://www.mdtraj.org .. _MDAnalysis: https://www.mdanalysis.org/ .. _periodictable: https://periodictable.readthedocs.io/en/latest/ .. note:: GRSQ expects Å for lengths, Å\ :sup:`-3`\ for volumes, and Å\ :sup:`-1`\ for \q\-vectors. Let's start simple: The Equivalence of The Discrete- and General Debye Equation =========================================================== In this example, we manually enter the information we need to check that the generalized Debye equation used on an RDF between two single atoms gives the same as simply using the discrete Debye using the discrete distance between the two atoms. ``grsq`` also includes various implementations of the discrete Debye function, See the :ref:`tutorials-other-label` for more information. .. literalinclude:: tutorials/grsqvsdebye.py .. figure:: tutorials/gfx/tut1_gr.png :alt: The sampled RDF used to calculate I(Q) :align: center :figclass: align-center .. figure:: tutorials/gfx/tut1_iq.png :alt: The calculated signals :align: center :figclass: align-center This was perhaps an excessive amount of manual work - but you don't need to. This only serves the purpose of looking a bit under the hood to see what's going on. Using RDFSets ============= The same :math:`I_{uu}(Q)` can be obtained much more quickly: .. literalinclude:: tutorials/simple2particle.py .. figure:: tutorials/gfx/tut1_iq_simple.png :alt: The calculated signals :align: center :figclass: align-center Here, we've used the helper function :py:meth:`grsq.grsq.rdfset_from_dir` to load all the rdf data files from a subdirectory - in this case there is only one. This requires that the files have a certain name-format, check its docstring for details. Really Using RDFSets ==================== An `RDFSet` is an `OrderedDict` with some extra functionalities. You can apply finite-size corrections across all (solvent-solvent and solute-solvent) RDFs, calculate the X-ray scattering terms: 1. Total signal (as in the previous tutorial): ``rdfs.get_iq(qvec)`` 2. Solute-Solvent (cross): ``rdfs.get_cross(qvec)`` 3. Solvent-Solvent: ``rdfs.get_solvent(qvec)`` 4. Displaced (Excluded) Volume: ``rdfs.get_dv(qvec)`` Define subsets of RDFs to analyse further, and add ``RDFSet``\s together as you do python ``list``\s. To test this out, we need some more RDFs. The gitlab repo of `grsq` contains some we can use, in the `data/rdfsets` directory here_ .. _here: https://gitlab.com/asod/grsq/-/tree/main/doc/source/tutorials/ :target: _blank Since you can have formatted your RDF files in a myriad of ways yourself, we refrain from using the helper function from before, and simply write a small piece of code that makes the RDFSet from scratch: .. literalinclude:: tutorials/generate_rdfset.py Note that we are manually adding g_BA(r) from the g_AB(r) data. The directory only contained g_AB(r) datafiles and not also g_BA(r), since they are identical, and it thus makes no sense to sample both of them. However, the `RDFSet` still needs both to get the right amplitude of the individual signals in the end. See our paper_ for details. Now, you should have the following: .. code-block:: none rdfs.show(): | # | KEY | STOICHOMETRY | REGION1, REGION2 | VOLUME | 000 | C-solute--C-solute | N1: 30, N2: 30 | solute , solute | 277088.4332 | | 001 | C-solute--Fe-solute | N1: 30, N2: 1 | solute , solute | 277088.4332 | | 002 | Fe-solute--C-solute | N1: 1, N2: 30 | solute , solute | 277088.4332 | | 003 | C-solute--H-solute | N1: 30, N2: 24 | solute , solute | 277088.4332 | | 004 | H-solute--C-solute | N1: 24, N2: 30 | solute , solute | 277088.4332 | | 005 | C-solute--H-solvent | N1: 30, N2: 14682 | solute , solvent | 277088.4332 | | 006 | H-solvent--C-solute | N1: 14682, N2: 30 | solvent, solute | 277088.4332 | | 007 | C-solute--N-solute | N1: 30, N2: 6 | solute , solute | 277088.4332 | | 008 | N-solute--C-solute | N1: 6, N2: 30 | solute , solute | 277088.4332 | | 009 | C-solute--O-solvent | N1: 30, N2: 7329 | solute , solvent | 277088.4332 | | 010 | O-solvent--C-solute | N1: 7329, N2: 30 | solvent, solute | 277088.4332 | | 011 | Fe-solute--H-solute | N1: 1, N2: 24 | solute , solute | 277088.4332 | | 012 | H-solute--Fe-solute | N1: 24, N2: 1 | solute , solute | 277088.4332 | | 013 | Fe-solute--H-solvent | N1: 1, N2: 14682 | solute , solvent | 277088.4332 | | 014 | H-solvent--Fe-solute | N1: 14682, N2: 1 | solvent, solute | 277088.4332 | | 015 | Fe-solute--N-solute | N1: 1, N2: 6 | solute , solute | 277088.4332 | | 016 | N-solute--Fe-solute | N1: 6, N2: 1 | solute , solute | 277088.4332 | | 017 | Fe-solute--O-solvent | N1: 1, N2: 7329 | solute , solvent | 277088.4332 | | 018 | O-solvent--Fe-solute | N1: 7329, N2: 1 | solvent, solute | 277088.4332 | | 019 | H-solute--H-solute | N1: 24, N2: 24 | solute , solute | 277088.4332 | | 020 | H-solute--H-solvent | N1: 24, N2: 14682 | solute , solvent | 277088.4332 | | 021 | H-solvent--H-solute | N1: 14682, N2: 24 | solvent, solute | 277088.4332 | | 022 | H-solute--N-solute | N1: 24, N2: 6 | solute , solute | 277088.4332 | | 023 | N-solute--H-solute | N1: 6, N2: 24 | solute , solute | 277088.4332 | | 024 | H-solute--O-solvent | N1: 24, N2: 7329 | solute , solvent | 277088.4332 | | 025 | O-solvent--H-solute | N1: 7329, N2: 24 | solvent, solute | 277088.4332 | | 026 | H-solvent--H-solvent | N1: 14682, N2: 14682 | solvent, solvent | 277088.4332 | | 027 | N-solute--H-solvent | N1: 6, N2: 14682 | solute , solvent | 277088.4332 | | 028 | H-solvent--N-solute | N1: 14682, N2: 6 | solvent, solute | 277088.4332 | | 029 | N-solute--N-solute | N1: 6, N2: 6 | solute , solute | 277088.4332 | | 030 | N-solute--O-solvent | N1: 6, N2: 7329 | solute , solvent | 277088.4332 | | 031 | O-solvent--N-solute | N1: 7329, N2: 6 | solvent, solute | 277088.4332 | | 032 | O-solvent--H-solvent | N1: 7329, N2: 14682 | solvent, solvent | 277088.4332 | | 033 | H-solvent--O-solvent | N1: 14682, N2: 7329 | solvent, solvent | 277088.4332 | | 034 | O-solvent--O-solvent | N1: 7329, N2: 7329 | solvent, solvent | 277088.4332 | That looks correct. We can quickly check how many combinations we should have: .. code-block:: python import itertools atom_types = [f'{atom}-{region}' for region, elements in stoich.items() for atom in elements] len(list(itertools.product(atom_types, repeat=2))) >>> 36 We only have 35, but there is also only one Fe atom, so there is no Fe-Fe RDFs. That adds up! Now we can apply low-q corrections (see the :ref:`next tutorial `), damping windows, and simulate the X-ray scattering signal: .. code-block:: python import matplotlib.pyplot as plt from grsq.damping import Damping # Applies van der Vegt correction to all relevant RDFs. # This overwrites the g(r) values with the corrected values # So if you do it twice, you have overcorrected! rdfs.vdv_correct() # Init a Lorch damping object with L = 1/2 box (to which the RDFs are sampled) damp = Damping('lorch', L=rdfs[0].r[-1]) # Define a q vector qvec = np.arange(0, 6, 0.01) # Calculate the cross term (for example): I_c = rdfs.get_cross(qvec=qvec, damping=damp) plt.plot(qvec, I_c) plt.xlim([0, 6]) plt.xlabel('$q$ (Å$^{-1}$)') plt.ylabel('$I_c(q)$ (a.u.)') plt.tight_layout() plt.show() .. figure:: tutorials/gfx/cross.png :alt: The resulting cross term :align: center :figclass: align-center Try yourself to calculate and plot the solute term, the solvent term, and the displaced-volume term. Sub-selections -------------- With the ``rdfs`` in hand, you can quickly analyse subsets of various pairs you want to investigate the scattering from, to get an overview of which pairs contribute what to the total signal(s). An example: .. code-block:: python # Fe-Solvent only fe_solv = RDFSet() for key, rdf in rdfs.items(): # Find whichever RDFs you're interested in, for example: if (('Fe' in key[0]) and (rdf.region2 =='solvent')): fe_solv[key] = rdf fe_solv.add_flipped(rdf) fe_solv.show() # N-Solvent only n_solv = RDFSet() for key, rdf in rdfs.items(): # Find whichever RDFs you're interested in, for example: if (('N' in key[0]) and (rdf.region2 =='solvent')): n_solv[key] = rdf n_solv.add_flipped(rdf) n_solv.show() # You can also add RDFSets: combined = fe_solv + n_solv I_combined = combined.get_cross(qvec=qvec, damping=damp) I_fe_solv = fe_solv.get_cross(qvec=qvec, damping=damp) I_n_solv = n_solv.get_cross(qvec=qvec, damping=damp) plt.plot(qvec, I_fe_solv, label='Fe-Solvent') plt.plot(qvec, I_n_solv, label='N-Solvent') plt.plot(qvec, I_n_solv + I_fe_solv, label='Summed (V1)') plt.plot(qvec, I_combined, '--', label='Summed (V2)') plt.xlim([0, 6]) plt.legend(loc='best') plt.xlabel('$q$ (Å$^{-1}$)') plt.ylabel('$I_c(q)$ (a.u.)') plt.tight_layout() plt.show() .. figure:: tutorials/gfx/rdfset_juggling.png :alt: More RDFSet uggling :align: center :figclass: align-center .. _fscor: Using and analysing finite-size corrections =========================================== As described in detail in our paper_, a useful method for checking how well your chosen finite-size correction performed is to calculate the structure factor from a corrected RDF at :math:`q=0` Å :math:`^{-1}` at increasing truncation-lengths :math:`R` of the integral in :ref:`the generalized Debye equation `. To illustrate this, let us recreate Fig. 3(B) in the paper_. .. literalinclude:: tutorials/s0r.py .. figure:: tutorials/gfx/tut3.png :alt: The calculated signals :align: center :figclass: align-center