How to compute the density of the spherical polymer using radius as a bin from the lammps trajectory?

Question

I am working on condensate studies and I did form a globular condensate of the polymer system using lammps. I wanted to find out the number density of condensate with respect to the radius as bins to find out the core density and surface density of the polymer. I tried using this python script but couldn't able to fetch any results. If anyone in the forum could help me out with this?

the code

import numpy as np

# Define the bin size
dr = 10  # adjust as needed

# Read in the coordinates from the Lammps trajectory file
# assuming x, y, z coordinates are in columns 1, 2, and 3
coords = np.loadtxt('test2.xyz', usecols=(1, 2, 3))

# Compute the distance of each atom from the center of the sphere
distances = np.linalg.norm(coords, axis=1)

# Bin the distances into r bins of equal size
bins = np.arange(0, np.max(distances) + dr, dr)
print(len(bins))
digitized = np.digitize(distances, bins)

print (len(digitized))
# Count the number of atoms in each bin
counts = np.bincount(digitized)
print(len(counts))

# Divide the count by the volume of each bin to get the number density
volumes = 4/3 * np.pi * (bins[1:]**3 - bins[:-1]*3)
print(len(volumes))
number_density = counts / volumes

# Compute the average number density over the trajectory
average_number_density = np.mean(number_density)

The xyz file looks like this

14000
 generated by VMD
  11       186.343002      194.738007      171.914001
  3        184.606995      195.947006      168.707001
  18       185.141006      197.257996      165.209000
  5        187.419006      199.007996      162.889008
  11       189.791000      198.328003      159.753006
  9        192.240005      195.429993      159.931000

I need to get the density of the molecules from the inner core to the outer surface by binning the spherical polymer number density vs R (bins).

I got the error operands could not be broadcast together with shapes (51,) (50,)

Answer 1

Understanding the Error:

Reason you get the error is because you are trying to divide two numpy arrays of unequal shape.

len(volumes)=32
len(counts)=33

Therefore, I changed the code for the volume calculation such that volume also has the same shape as the bins. By this I was able to get average_number_density.

Below are the lines I changed in the volume calculation.

offset = bins.tolist()[-1]*3
volumes = 4/3 * np.pi * np.array([ i**3- offset for i in bins])

This gives the average_number_density = 1.3246815571536801e-09

Full Code for completeness:

# Define the bin size
dr = 10 # adjust as needed

# Read in the coordinates from the Lammps trajectory file
# assuming x, y, z coordinates are in columns 1, 2, and 3
coords = np.loadtxt('test2.xyz', usecols=(1, 2, 3))

# Compute the distance of each atom from the center of the sphere
distances = np.linalg.norm(coords, axis=1)

# Bin the distances into r bins of equal size
bins = np.arange(0, np.max(distances) + dr, dr)
digitized = np.digitize(distances, bins)

# Count the number of atoms in each bin
counts = np.bincount(digitized)

# Divide the count by the volume of each bin to get the number density
offset = bins.tolist()[-1]*3
# volumes = 4/3 * np.pi * (bins[1:]**3 - bins[:-1]*3) 
volumes = 4/3 * np.pi * np.array([ i**3- offset for i in bins]) 

number_density = counts / volumes

# Compute the average number density over the trajectory
average_number_density = np.mean(number_density)

print(average_number_density)

And the XYZ file I used: (Notice the first two lines are taken out, otherwise numpy will complain of not having data in all columns)

  11       186.343002      194.738007      171.914001
  3        184.606995      195.947006      168.707001
  18       185.141006      197.257996      165.209000
  5        187.419006      199.007996      162.889008
  11       189.791000      198.328003      159.753006
  9        192.240005      195.429993      159.931000