a coarse-grained model that shields two adjacent beads (1-2 and 1-3 interactions) within the same chain

[yingyingwukim]

• I am currently working on calculating the RDF for polymer chains and need to exclude the effects of bond angle interactions. Let me explain the purpose of my code, and I would appreciate your help in implementing a coarse-grained model that shields two adjacent beads (1-2 and 1-3 interactions) within the same chain.
• (1)First, the all-atom system consists of 40 chains, and the all-atom molecular chains are mapped to centroids based on a mapping rule, with these centroids serving as the new beads in the coarse-grained model.
• (2)After mapping, the coarse-grained model consists of linear molecular chains arranged as S-S-S-S-E-P-P-E-S-S-S-S.
• (3)The system now contains three types of beads: S, E, and P. I need to calculate the RDF for each bead type while excluding adjacent beads within the same chain. I would like your guidance on how to achieve this calculation.
• (4)I am using the InterRDF(atom1, atom2, nbins=200, range=(0.0, 18), exclusion_block=(3,3)) method to exclude neighboring interactions. However, the radial distribution functions for S-E and P-P show significant discrepancies compared to literature values. I suspect that the method for excluding neighbors may not be implemented correctly. The documentation on the MDAnalysis website is somewhat limited, and I hope you can assist me with this.

`# -- coding: utf-8 --
"""
Modified script to construct a coarse-grained model using the center of mass of residues from an all-atom model,
establishing the coarse-grained positions and bond connections, and excluding 1-2 and 1-3 interactions within the same chain for RDF calculations.

Calculating the RDF of all atoms—3 types of beads: E, P, S—6 types of RDF:

g(r)EE, g(r)EP, g(r)ES,
g(r)PP, g(r)PS,
g(r)SS.

Function 1 - Calculate center of mass.
Function 2 - Select atom groups and calculate center of mass.
Function 3 - Calculate RDF.
@author: Yingying Wu
"""

import numpy as np
import matplotlib.pyplot as plt
from MDAnalysis import Universe
from MDAnalysis.analysis.rdf import InterRDF
from tqdm import tqdm
import os

Function: Calculate the center of mass for a given atom group

def calculate_com(atom_group):
return atom_group.center_of_mass(unwrap=True)

Function: Select residues, calculate their center of mass, and categorize by residue type

def select_and_calculate_com(u, num_chains, num_residues):
"""
Select residues for each chain, calculate their center of mass, and categorize them by residue type.
"""
com_list_all = []

for chain_id in range(num_chains):
    start_residue = chain_id * num_residues
    end_residue = start_residue + num_residues

    n = chain_id * 151  # Atom ID offset for each chain

    # Calculate the center of mass for each residue in the current chain
    com_list = []
    for i, residue in enumerate(u.residues[start_residue:end_residue]):
        if i == 6:  # For the 7th residue (index 6), include specific atoms to calculate the center of mass
            cor_atoms = u.select_atoms(f"bynum {n+75} or bynum {n+76} or bynum {n+77}")
            combined_group = residue.atoms.union(cor_atoms)
            com = calculate_com(combined_group)
        else:
            com = calculate_com(residue.atoms)
        com_list.append(com)

    com_list_all.extend(com_list)  # Add the current chain's centers of mass to the total list

return com_list_all

Function: Calculate RDF between centers of mass while excluding 1-2 and 1-3 interactions within the same chain

def compute_rdf(atom1, atom2):
"""
Exclude 1-2 and 1-3 interactions within the same chain (i.e., bond angle effects on g(r))
SSSSEPPESSSS
"""
rdf = InterRDF(atom1, atom2, nbins=200, range=(0.0, 18), exclusion_block=(3,3))
rdf.run()
return rdf.results.bins, rdf.results.rdf

Main program: Initialize Universe and perform coarse-grained mapping

u = Universe("G:/01ibi-pu/PU-IBI/AA/TPU/TPU40chains-eq-03.data", "G:/01ibi-pu/PU-IBI/AA/TPU/unwrap.dcd", format="LAMMPS")
num_chains = 40
num_residues = 12
start = 1
end = 1001
n_frames = end - start # Number of frames
residue_types = "S" * 4 + "EPPE" + "S" * 4 # Define residue type pattern

Create an empty coarse-grained model Universe

new_u = Universe.empty(num_chains * num_residues, trajectory=True, n_frames=len(u.trajectory))
new_u.add_TopologyAttr('name')
new_u.add_TopologyAttr('type')

bonds = []
for ts_idx, ts in enumerate(tqdm(u.trajectory, desc="Calculating Centers of Mass")):
com_list_all = select_and_calculate_com(u, num_chains, num_residues)
new_positions = np.array(com_list_all)
new_u.trajectory[ts_idx].positions = new_positions
new_u.trajectory[ts_idx].dimensions = u.dimensions
new_u.atoms.names = list(residue_types) * num_chains
new_u.atoms.types = list(residue_types) * num_chains

Select each type of coarse-grained atom from the new Universe

S_atoms = new_u.select_atoms('name S')
E_atoms = new_u.select_atoms('name E')
P_atoms = new_u.select_atoms('name P')

rdf_pairs = [
(S_atoms, S_atoms, 'SS'),
(S_atoms, E_atoms, 'SE'),
(S_atoms, P_atoms, 'SP'),
(E_atoms, E_atoms, 'EE'),
(E_atoms, P_atoms, 'EP'),
(P_atoms, P_atoms, 'PP')
]

rdf_results = {}

for (atom1, atom2, label) in rdf_pairs:
bins, rdf_data = compute_rdf(atom1, atom2)
rdf_results[label] = rdf_data

Set the output directory for RDF results

output_dir = 'G:/01ibi-pu/PU-IBI/AA/TPU/rdf-TPU'
os.makedirs(output_dir, exist_ok=True)

Plot and save RDF graphs for each pair

for key, rdf_data in rdf_results.items():
plt.figure(figsize=(7, 6))
y_max = 2 if key in ['SS', 'SE', 'SP'] else 3
plt.plot(bins, rdf_data, label=f'g(r) {key}')
plt.xlim(0, 18)
plt.ylim(0, y_max)
plt.xlabel('Distance (Å)', fontdict={'family': 'Arial', 'size': 22})
plt.ylabel('g(r)', fontdict={'family': 'Arial', 'size': 22})
plt.axhline(1, c='blue', linestyle='--')
plt.title(f'RDF {key}')
plt.legend()
plt.grid(True)
plt.savefig(os.path.join(output_dir, f'rdf_{key}.png'), bbox_inches='tight', dpi=800)
plt.show()
np.savetxt(os.path.join(output_dir, f'rdf_{key}.txt'), np.column_stack((bins, rdf_data)), fmt="%.4f,%.6f", delimiter=",", header="Distance, g(r)", comments='')
`