Computing mass and charge density on each axis
Here we compute the mass and charge density of water along the three cartesian axes of a fixed-volume unit cell (i.e. from a simulation in the NVT ensemble).
Last updated: December 2022 with MDAnalysis 2.4.0-dev0
Minimum version of MDAnalysis: 0.17.0
Packages required:
[1]:
import MDAnalysis as mda
from MDAnalysis.tests.datafiles import waterPSF, waterDCD
from MDAnalysis.analysis import lineardensity as lin
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
Loading files
The test files we are working with are a cube of water.
[2]:
u = mda.Universe(waterPSF, waterDCD)
/home/pbarletta/mambaforge/envs/guide/lib/python3.9/site-packages/MDAnalysis/coordinates/DCD.py:165: DeprecationWarning: DCDReader currently makes independent timesteps by copying self.ts while other readers update self.ts inplace. This behavior will be changed in 3.0 to be the same as other readers. Read more at https://github.com/MDAnalysis/mdanalysis/issues/3889 to learn if this change in behavior might affect you.
warnings.warn("DCDReader currently makes independent timesteps"
MDAnalysis.analysis.lineardensity.LinearDensity (API docs) will partition each of your axes into bins of user-specified binsize (in angstrom), and give the average mass density and average charge density of your atom group selection.
This analysis is only suitable for a trajectory with a fixed box size. While passing a trajectory with a variable box size will not raise an error, LinearDensity will not account for changing dimensions. It will only evaluate the density of your atoms in the bins created from the trajectory frame when the class is first initialised.
Below, we iterate through the trajectory to verify that its box dimensions remain constant.
[3]:
for ts in u.trajectory:
print(ts.dimensions)
[50. 50. 50. 90. 90. 90.]
[50. 50. 50. 90. 90. 90.]
[50. 50. 50. 90. 90. 90.]
[50. 50. 50. 90. 90. 90.]
[50. 50. 50. 90. 90. 90.]
[50. 50. 50. 90. 90. 90.]
[50. 50. 50. 90. 90. 90.]
[50. 50. 50. 90. 90. 90.]
[50. 50. 50. 90. 90. 90.]
[50. 50. 50. 90. 90. 90.]
You can choose to compute the density of individual atoms, residues, segments, or fragments (groups of bonded atoms with no bonds to any atom outside the group). By default, the grouping is for atoms.
[4]:
density = lin.LinearDensity(u.atoms,
grouping='atoms').run()
The results of the analysis are in density.results.
[5]:
density.nbins
[5]:
200
[6]:
density.results['x']['mass_density']
[6]:
array([0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0.00053562, 0.00080344, 0.00876945, 0.03507781, 0.00107125,
0.00348155, 0.00241031, 0.02791523, 0.04277601, 0.0175389 ,
0.00160687, 0.00133906, 0.00026781, 0. , 0.00107125,
0.00107125, 0.00053562, 0. , 0.03400656, 0.0196814 ,
0.02339659, 0.0135559 , 0.00026781, 0.00107125, 0.00107125,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ])
[7]:
density.results['x'].keys()
[7]:
KeysView({'dim': 0, 'slice_volume': 625.0, 'mass_density': array([0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0.00053562, 0.00080344, 0.00876945, 0.03507781, 0.00107125,
0.00348155, 0.00241031, 0.02791523, 0.04277601, 0.0175389 ,
0.00160687, 0.00133906, 0.00026781, 0. , 0.00107125,
0.00107125, 0.00053562, 0. , 0.03400656, 0.0196814 ,
0.02339659, 0.0135559 , 0.00026781, 0.00107125, 0.00107125,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ]), 'mass_density_stddev': array([0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0.00107125, 0.00122727, 0.01688797, 0.01691979, 0.00177646,
0.00241031, 0.00279604, 0.02179554, 0.02689655, 0.02096112,
0.001312 , 0.00133906, 0.00080344, 0. , 0.001312 ,
0.001312 , 0.00107125, 0. , 0.01700328, 0.03402765,
0.02131476, 0.01957657, 0.00080344, 0.001312 , 0.001312 ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ]), 'charge_density': array([ 0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0.00022158, 0.00033237, -0.00033237, -0.00132949, 0.00044316,
0.00144029, 0.00099712, -0.00033237, -0.00210503, -0.00066475,
0.00066475, 0.00055396, 0.00011079, 0. , 0.00044316,
0.00044316, 0.00022158, 0. , -0.00177266, 0.00022158,
-0.00022158, -0.00033237, 0.00011079, 0.00044316, 0.00044316,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ]), 'charge_density_stddev': array([0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0.00044316, 0.00050771, 0.00099712, 0.00108553, 0.00073491,
0.00099712, 0.00115669, 0.00111344, 0.00144029, 0.00112985,
0.00054276, 0.00055396, 0.00033237, 0. , 0.00054276,
0.00054276, 0.00044316, 0. , 0.00088633, 0.0018406 ,
0.00129204, 0.00111344, 0.00033237, 0.00054276, 0.00054276,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ]), 'hist_bin_edges': array([ 0. , 0.25, 0.5 , 0.75, 1. , 1.25, 1.5 , 1.75, 2. ,
2.25, 2.5 , 2.75, 3. , 3.25, 3.5 , 3.75, 4. , 4.25,
4.5 , 4.75, 5. , 5.25, 5.5 , 5.75, 6. , 6.25, 6.5 ,
6.75, 7. , 7.25, 7.5 , 7.75, 8. , 8.25, 8.5 , 8.75,
9. , 9.25, 9.5 , 9.75, 10. , 10.25, 10.5 , 10.75, 11. ,
11.25, 11.5 , 11.75, 12. , 12.25, 12.5 , 12.75, 13. , 13.25,
13.5 , 13.75, 14. , 14.25, 14.5 , 14.75, 15. , 15.25, 15.5 ,
15.75, 16. , 16.25, 16.5 , 16.75, 17. , 17.25, 17.5 , 17.75,
18. , 18.25, 18.5 , 18.75, 19. , 19.25, 19.5 , 19.75, 20. ,
20.25, 20.5 , 20.75, 21. , 21.25, 21.5 , 21.75, 22. , 22.25,
22.5 , 22.75, 23. , 23.25, 23.5 , 23.75, 24. , 24.25, 24.5 ,
24.75, 25. , 25.25, 25.5 , 25.75, 26. , 26.25, 26.5 , 26.75,
27. , 27.25, 27.5 , 27.75, 28. , 28.25, 28.5 , 28.75, 29. ,
29.25, 29.5 , 29.75, 30. , 30.25, 30.5 , 30.75, 31. , 31.25,
31.5 , 31.75, 32. , 32.25, 32.5 , 32.75, 33. , 33.25, 33.5 ,
33.75, 34. , 34.25, 34.5 , 34.75, 35. , 35.25, 35.5 , 35.75,
36. , 36.25, 36.5 , 36.75, 37. , 37.25, 37.5 , 37.75, 38. ,
38.25, 38.5 , 38.75, 39. , 39.25, 39.5 , 39.75, 40. , 40.25,
40.5 , 40.75, 41. , 41.25, 41.5 , 41.75, 42. , 42.25, 42.5 ,
42.75, 43. , 43.25, 43.5 , 43.75, 44. , 44.25, 44.5 , 44.75,
45. , 45.25, 45.5 , 45.75, 46. , 46.25, 46.5 , 46.75, 47. ,
47.25, 47.5 , 47.75, 48. , 48.25, 48.5 , 48.75, 49. , 49.25,
49.5 , 49.75, 50. ], dtype=float32)})
[8]:
density.results['y']['dim']
[8]:
1
[9]:
plt.plot(np.linspace(0, 50, 200), density.results['x']['mass_density'])
[9]:
[<matplotlib.lines.Line2D at 0x7f5d682b3b20>]
References
[1] Richard J. Gowers, Max Linke, Jonathan Barnoud, Tyler J. E. Reddy, Manuel N. Melo, Sean L. Seyler, Jan Domański, David L. Dotson, Sébastien Buchoux, Ian M. Kenney, and Oliver Beckstein. MDAnalysis: A Python Package for the Rapid Analysis of Molecular Dynamics Simulations. Proceedings of the 15th Python in Science Conference, pages 98–105, 2016. 00152. URL: https://conference.scipy.org/proceedings/scipy2016/oliver_beckstein.html, doi:10.25080/Majora-629e541a-00e.
[2] Naveen Michaud-Agrawal, Elizabeth J. Denning, Thomas B. Woolf, and Oliver Beckstein. MDAnalysis: A toolkit for the analysis of molecular dynamics simulations. Journal of Computational Chemistry, 32(10):2319–2327, July 2011. 00778. URL: http://doi.wiley.com/10.1002/jcc.21787, doi:10.1002/jcc.21787.