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='')

[orbeckst]

Sorry that you haven't received an answer.

Perhaps you could simplify your question and ask something more specific? I assume that most people on the forum don't have the time to read through your code in detail. If you can reduce a question to a few lines of code with a concise description then this might make it easier. (see also #4558 )

It is good that you're showing your code but

However, the radial distribution functions for S-E and P-P show significant discrepancies compared to literature values.

unless someone is an expert on this very specific system, nobody will be able to help unless you give more evidence, e.g., what your results are and what the literature values are (citation and comparison to your data).

I suspect that the method for excluding neighbors may not be implemented correctly.

It's possible but then you should be able to produce a test case where you 100% know the answer and where you can validate the answer yourself. We'd be very happy to see a clear test case that demonstrates a bug. Once we have such an example it's relatively easy to fix code.