Analysis - Linear Density module#670
Conversation
There was a problem hiding this comment.
Generally, you don't want to use lists [a, b, c, d] for numerical stuff, use a numpy array np.array([elem.total_mass() for elem in group]). When you put this into weights= in numpy histogram, it's likely being converted into an array there (if not already an array), so with each loop iteration you're having to make this array many times.
|
Thanks @richardjgowers for the very useful inputs! I will look into them and commit once it's done |
|
I modified the code according to Richard's suggestions, I hope it is fine! Let me know if anything else can be improved. For example, when computing the density profiles on groups of atoms, the analysis is much slower. I wonder if it is only because of the additional computational load due to the calculation of the COMs or if there is any part of the code that has been written poorly/in a non efficient way. |
|
I think the groups are slow because they're done individually rather than all together. What would be really cool, (and is a separate issue), would be a subroutine which takes a series of groupings, and returns the center of mass/geometry for each. So for groups [1, 2, 4, 5], [6, 8, 10] and [11, 15] we could do # Array of all particles
[1, 2, 4, 5, 6, 8, 10, 11, 15]
# Pointers to starts and stops within the above array
starts = [0, 4, 8]
stops = [4, 8, -1]
do i = 1, nCOGs:
do j = starts(i), stops(i):
COG(i) = COG(i) + positions(j)Then we'd probably write it in C like these are: |
|
That was actually my main guess. I think it would be cool to be able to obtain, at once, an array of properties for all residues/segments in an atom group as follows: universe = mda.Universe(PSF,DCD)
selection = universe.select_atoms('all')
# arrays containing properties for each residue or each segment
residues_COMs = selection.residues.centers_of_mass
segments_principal_axes = selection.segments.principal_axes
# instead of the slow
residues_COMs = np.array([elem.center_of_mass for elem in selection.residues])
segments_principal_axes = np.array([elem.principalAxes() for elem in selection.segments])At least in our field of work this would come very handy as it often happens that we need to access these properties at once. Do you think this would be feasible and that the syntax above would be intuitive? |
|
Yeah, I think @dotsdl and I have talked about this before. In the upcoming topology revolution, we've done this for masses, so |
There was a problem hiding this comment.
You don't need to call np.array on these because they're always returned as arrays anyway
|
Docstrings everywhere please, following these rules: To make it look something like: |
|
Hi! I have fixed hopefully most of the smaller issues and added docstrings wherever necessary. Being my first time writing docstrings in numpy style, please let me know if they can be improved in any way.
|
There was a problem hiding this comment.
You need to copy the same thing from the top of all the other files here
|
Ok, so for the auto documentation to work, you need to go to http://pythonhosted.org/MDAnalysis/documentation_pages/analysis_modules.html |
|
Ok almost done,
If you could add a quick and easy regression test using the PSF/DCD system that would be good. Just run it once on that trajectory, and add the results as a test so if it stops working we know. Finally, if you want, you can add a few lines to the next release blog post, you can clone that repo and checkout the Then maybe have another @MDAnalysis/coredevs have a look at this too |
|
So, I am writing a test case but I found an inconsistency with the new parallel version of If there are no objections I guess it would be better to restore the previous behaviour ( |
There was a problem hiding this comment.
You don't need to ask for Universe here, as it's always available as selection.universe
|
Yeah, I think we needed the Universe in setup frames for parallel, can't remember why either. And yeah, for now just make this work in serial, then we can go back and fix it when parallel happens. |
…ort order, added save method, use of in-place sum of np arrays
…tible with AnalysisBase of MDAnalysis 0.13.0
mattiafelice-palermo
force-pushed
the
feature-LinearDensity
branch
from
3b103e0 to
fd44b72
Compare
…oved print statements from test
|
Ok so I think the option is now working as expected. If any other coredev wants to have a look at it please feel free :) Thanks! |
Analysis - Linear Density module
|
Hey Mattia. Looking good, congrats on your first PR! |
|
That's great, thank you Richard! :) |
2 participants
As discussed in here.
This is a quite simple analysis module that, given a NVT trajectory in an orthorombic simulation cell, outputs a column based text file with mass and charge densities along the three axes of the cell.
A minimal working example would be as follows:
where
groupingspecifies the group of atoms on which the density distribution will be computed (it takes the usual argumentsatoms,residues,fragmentsetc.) andbinsizeis the resolution of the computed histogram.The analysis module outputs a file named
dcdfilename.description_grouping.ldensfrom which it is easy to create plots of density profiles. I don't know if this module replicates the functionality of the already present density module (I could not understand from here) so I apologize in advance if this is a replica of already existing work!