Density profile in a pore

[Kami Kahzadi]

 Hello,

I am trying to generate density profile in a certain region (pore) of simulation box.

The density profile (molecules/nm3) needs to be obtained for 1.5<y<3.2 and 11<z<18 nm. Density vs y(pore distance). 

I want to generate it exactly like the attached picture, and I need help with my code attached.

My question is : how to plot and get the number of particles for a specific region vs distance? I have used atom selection based on the documentation.

Best regards,

Kamyab 

[Oliver Beckstein]

Hello Kamyab,

Welcome to the MDAnalysis mailing list!

On Sep 15, 2023, at 2:59 AM, Kami Kahzadi <kahza...@gmail.com> wrote:

 Hello,

I am trying to generate density profile in a certain region (pore) of simulation box.

The code that you’re using generates a density on a 3D grid. You can export it to a dx file and view it (see the docs).

The density profile (molecules/nm3) needs to be obtained for 1.5<y<3.2 and 11<z<18 nm. Density vs y(pore distance). 

I want to generate it exactly like the attached picture, and I need help with my code attached.

The 3D DensityAnalysis code does not contain tools to get 1D densities. We have code for linear density (LinearDensity) and for radial densities (RDF) but not for cylindrical geometries. You’ll have to write your own code. 

You can also look at MAICoS https://maicos.readthedocs.io/en/main/examples/radial-distribution-function.html, which uses MDAnalysis. Perhaps there’s something in there that works for you. GROMACS might also have a tool to do it.

My question is : how to plot and get the number of particles for a specific region vs distance? I have used atom selection based on the documentation.

I am assuming that you’re looking for cylindrical pore geometries. I would iterate over frames in the trajectory https://userguide.mdanalysis.org/stable/examples/analysis/custom_trajectory_analysis.html and for each particle compute the distance from the pore axis. Then use numpy’s histogram function to histogram all data on a radial grid and properly normalize the density with the Jacobian for cylindrical coordinates.

Best,

Oliver

Best regards,

Kamyab 

import MDAnalysis

import matplotlib.pyplot as plt

from MDAnalysis.analysis.density import *

from MDAnalysis import Universe

u = MDAnalysis.Universe('.tpr','.xtc')

binding_selection = u.select_atoms("name (binding atom) and prop z>111 and prop z<188 and prop y<32 and prop y>15")

binding_COM = ligand_selection.center_of_mass()

density_atom = u.select_atoms('resname (target molecule) and around 10 name (binding atom) and prop z>111 and prop z<188', updating= True)

density = DensityAnalysis(density_atom,gridcenter=binding_COM, xdim=41.44 , ydim=10,

    zdim=78.32, delta=1).run(start=3000, step=5000)

density.results.density.convert_density('nm^{-3}')

--
You received this message because you are subscribed to the Google Groups "MDnalysis discussion" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mdnalysis-discus...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/CAKXcTVV7CvmTZPCiAcjTbwpvkZc2JUCT5gr1aePNtRWvCzjhUw%40mail.gmail.com.
<density.py><IMG_20230909_182418.jpg>

--

Oliver Beckstein (he/his/him)

email: orbe...@mdanalysis.org
twitter: @orbeckst

GitHub: @orbeckst

MDAnalysis – a NumFOCUS fiscally sponsored project

https://www.mdanalysis.org/