Density calculation of globular polymer system?

LAMMPS
1

Dear All,
I am working on a project in which we have formed a spherical polymer system of 100 chains with each chain having 140 atoms. I want to compute the density of the sphere bining with respect to radius of the sphere from the lammpstrj file obtained from the lammps trajectory. I tried with my own code but ended with error. Could anyone please help me out in this regard?

the code that I used
import numpy as np

def compute_density_from_lammps(trajectory_file, num_chains, chain_length, bin_size):
# Read LAMMPS trajectory file
with open(trajectory_file, ‘r’) as file:
lines = file.readlines()

# Extract atom coordinates at each time step
coordinates = []
atom_start = False
for line in lines:
    if 'ITEM: TIMESTEP' in line:
        atom_start = False
        if len(coordinates) > 0:
            break
    if atom_start:
        coords = line.split()[2:5]
        coordinates.append([float(coord) for coord in coords])
    if 'ITEM: ATOMS' in line:
        atom_start = True

coordinates = np.array(coordinates)

# Calculate the center coordinates of the system
center = np.mean(coordinates, axis=0)

# Calculate the distances of particles from the center
distances = np.linalg.norm(coordinates - center, axis=1)

# Define radial bins
max_distance = np.max(distances)
num_bins = int(max_distance / bin_size) + 1
bins = np.linspace(0, max_distance, num=num_bins)

# Count particles in each bin
atom_counts = np.histogram(distances, bins=bins)[0]

# Calculate bin volumes
bin_volumes = (4/3) * np.pi * (bins[1:]**3 - bins[:-1]**3)

# Compute densities
densities = atom_counts / bin_volumes

# Average densities over chains
average_density = np.mean(densities) * chain_length / num_chains

return average_density

Example usage

trajectory_file = ‘lammps_trajectory.txt’
num_chains = 100
chain_length = 140
bin_size = 1.0

density = compute_density_from_lammps(trajectory_file, num_chains, chain_length, bin_size)
print(“Average density:”, density)

Thank you all
With regards,
Kesavan

2

It’s hard to help without knowing what error messages you’re getting. Also, some parts of your script would be helped by you using an established module like MDAnalysis to handle the trajectory reading, instead of “rolling your own” and risking some error.

3

Thank you, I will try with Mdanalysis and post the results if things go right.

4

If you import your trajectory using MDAnalysis, then you can easily extract spherical density profiles using MAICoS.

5

Thank you Simon