Error in RDF calculation

[LipikaBaidya]

I am using this to calculate RDF:

import MDAnalysis as mda
import MDAnalysis.analysis.rdf as mdana
import matplotlib.pyplot as plt
u = mda.Universe("pos_cm_0_001.pdb", "output.dcd")
u.dimensions = [72,72,72, 90, 90, 90]

s1 = u.select_atoms('name C')
s2 = u.select_atoms('name MG')

ags   = [s1,s2]
rdf = mdana.InterRDF(s1,s2,periodic=True)#,norm="rdf")
rdf.run()
plt.plot(rdf.bins, rdf.rdf)

But I am getting 0 values.

How to fix this??

[flyfir248]

Make sure that the selections for s1 and s2 are actually selecting atoms from your system. You can do this by printing the number of atoms selected.

print(f'Number of atoms in s1: {len(s1)}')
print(f'Number of atoms in s2: {len(s2)}')

try this :

import MDAnalysis as mda
import MDAnalysis.analysis.rdf as mdana
import matplotlib.pyplot as plt

# Load the universe
u = mda.Universe("pos_cm_0_001.pdb", "output.dcd")
u.dimensions = [72, 72, 72, 90, 90, 90]

# Select atoms
s1 = u.select_atoms('name C')
s2 = u.select_atoms('name MG')

# Check the number of atoms selected
print(f'Number of atoms in s1: {len(s1)}')
print(f'Number of atoms in s2: {len(s2)}')

# Check universe and trajectory
print(u)
print(u.trajectory)

# Check dimensions
print(u.dimensions)

# Calculate RDF
ags = [s1, s2]
rdf = mdana.InterRDF(s1, s2, periodic=True)
rdf.run()

# Check RDF values
print(rdf.rdf)

# Plot RDF
plt.plot(rdf.bins, rdf.rdf)
plt.xlabel('Distance (Å)')
plt.ylabel('RDF')
plt.title('Radial Distribution Function')
plt.show()

[LipikaBaidya]

The code is giving error: /Users/lipika/anaconda3/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 #3889 to learn if this change in behavior might affect you. warnings.warn("DCDReader currently makes independent timesteps" Number of atoms in s1: 30 Number of atoms in s2: 224 <Universe with 254 atoms> <DCDReader output.dcd with 100000 frames of 254 atoms> [72. 72. 72. 90. 90. 90.] Traceback (most recent call last): File "/Users/lipika/HM_rna/SIS_ION/rdf.py", line 27, in <module> rdf.run() File "/Users/lipika/anaconda3/lib/python3.9/site-packages/MDAnalysis/analysis/base.py", line 450, in run self._conclude() File "/Users/lipika/anaconda3/lib/python3.9/site-packages/MDAnalysis/analysis/rdf.py", line 314, in _conclude norm *= N / box_vol ZeroDivisionError: float division by zero

…

On Sun, Jul 7, 2024 at 4:04 AM Anoop Johny ***@***.***> wrote: Make sure that the selections for s1 and s2 are actually selecting atoms from your system. You can do this by printing the number of atoms selected. print(f'Number of atoms in s1: {len(s1)}') print(f'Number of atoms in s2: {len(s2)}') try this : import MDAnalysis as mda import MDAnalysis.analysis.rdf as mdana import matplotlib.pyplot as plt # Load the universe u = mda.Universe("pos_cm_0_001.pdb", "output.dcd") u.dimensions = [72, 72, 72, 90, 90, 90] # Select atoms s1 = u.select_atoms('name C') s2 = u.select_atoms('name MG') # Check the number of atoms selected print(f'Number of atoms in s1: {len(s1)}') print(f'Number of atoms in s2: {len(s2)}') # Check universe and trajectory print(u) print(u.trajectory) # Check dimensions print(u.dimensions) # Calculate RDF ags = [s1, s2] rdf = mdana.InterRDF(s1, s2, periodic=True) rdf.run() # Check RDF values print(rdf.rdf) # Plot RDF plt.plot(rdf.bins, rdf.rdf) plt.xlabel('Distance (Å)') plt.ylabel('RDF') plt.title('Radial Distribution Function') plt.show() — Reply to this email directly, view it on GitHub <#4631 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AQO2GLWXJLZAHOBUHMXTU73ZLDY7LAVCNFSM6AAAAABKKPFQ46VHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM4TSNZYGIZDK> . You are receiving this because you authored the thread.Message ID: ***@***.***>

[orbeckst]

Does your trajectory record box information?

If not then you have to add it manually but it’s not enough to set dimensions because like coordinates, it gets changed at each Timestep (to None if you don’t have box information). Instead you can use the transformation https://docs.mdanalysis.org/stable/documentation_pages/transformations/boxdimensions.html