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