Constructing empty universe and adding residue attributes gives error

[nambiarnik]

Expected behavior

When I try to read a .gro file and create a new environment out of the data, It throws out errors whenever I use add_TopologyAttr function to specify residues or resids. I want to be able to use this data from .gro file to then write a duplicate .gro file, which will give me the opportunity to modify and add residues.

Actual behavior

After assigning resid or resname attributes of the new environment, throws out an index error whenever I try to call the residues or resids atomgroup of the new environment.
IndexError: index 400 is out of bounds for axis 0 with size 400index 400 is out of bounds for axis 0 with size 400

Code to reproduce the behavior

import MDAnalysis as mda
from MDAnalysis.coordinates.GRO import GROWriter
import numpy as np
from scipy.spatial.transform import Rotation as R


# Example usage
input_gro = "C:/Users/nambi/Downloads/replace bilayer/step9.gro"
output_gro = "C:/Users/nambi/Downloads/replace bilayer/modified_membrane.gro"

# Load membrane and replacement lipid
#gro file has 400 lipids only
membrane = mda.Universe(input_gro)

  
# Store modified atoms
new_positions = []
new_names = []
new_resids = []
new_residue_map = {}

# Copy over existing atoms except for the lipids being replaced
for atom in membrane.atoms:
    new_positions.append(atom.position)
    new_names.append(atom.name)
    new_resids.append(atom.resid)
    new_residue_map[atom.resid] = atom.resname  # Store original residue names

# Create a new Universe with modified topology and coordinates
unique_resids = sorted(set(new_resids))  # Get unique residue IDs
new_resnames = [new_residue_map[resid] for resid in unique_resids]  # Ensure correct resname mapping    
new_universe = mda.Universe.empty(n_atoms=len(new_positions), n_residues=len(unique_resids), atom_resindex=new_resids,trajectory=True)
new_universe.dimensions = membrane.dimensions
new_universe.add_TopologyAttr("names", new_names)
new_universe.add_TopologyAttr("resids", list(membrane.residues.resids)) 
new_universe.add_TopologyAttr("resnames", list(membrane.residues.resnames))  
new_universe.atoms.positions = np.array(new_positions)

# Write the modified structure
with GROWriter(output_gro) as writer:
    writer.write(new_universe.atoms)

Current version of MDAnalysis

• Which version are you using? 2.7.0
• Which version of Python (python -V)? 3.11
• Which operating system? Windows

[orbeckst]

Do you need to build a Universe with empty? Couldn't you create an AtomGroup ag and then create a new universe with Merge():

ag = u.select_atoms(...)
new_u = mda.Merge(ag)
new_u.atoms.write("new.gro")

[orbeckst]

@nambiarnik did using Merge solve your problem?

If you think that there's a bug in Universe.empty() then could you create a minimal working example with a file from out test suite (e.g., MDAnalysisTests.datafiles.GRO that clearly shows the issue?

[nambiarnik]

@orbeckst Merge() did work for my use case. I am not sure about Universe.empty()

[orbeckst]

Given that Merge solves your problem, it will be useful to have your use case archived in our discussion forum.

If any specific issues with Universe.empty() arise, you (or anyone else) should feel free to open a new issue. Thanks!

[orbeckst]

Ok, thanks for the feedback. I'll convert the issue to a discussion.