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

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,)

• Welcome to Mater Modeling ! Please paste your entire input XYZ file as a code block. Also please elaborate what you mean by "I tried using this python script but couldn't able to fetch any results. " by adding what is your expected result, and what is the result you are currently getting. May 11, 2023 at 9:22
• Please add the entire stack trace of the error as a code block May 11, 2023 at 11:57

### 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

• Thank you for your kind help. May 16, 2023 at 11:19
• @kesavan if my answer solves your question then please mark the answer as accepted by clicking on the tick mark next to it. You also get +2 points for doing it. May 16, 2023 at 11:33