Average structure used in a reference structure

[royyraa]

Hello everyone,

I am trying to create an average structure from a set of DCD files and then use that structure as a reference to calculate the RMSD. However, there seems to be an error in calculating the RMSD with reference to the average structure. It is likely that the code is not properly adding the average structure to the reference. I am attaching the code snippet here. Any kind of help or suggestions would be appreciated.

 universes = [mda.Universe("test.psf", dcd) for dcd in dcd_files]
 n_atoms = universes[0].atoms.n_atoms
 avg_positions = np.zeros((n_atoms, 3))
 total_frames = 0
 for u in universes:
        for ts in u.trajectory:
            avg_positions += u.atoms.positions
            total_frames += 1
    if total_frames == 0:
        print(f"No frames found in {folder_name}")
        continue

    avg_positions /= total_frames
    reference = mda.Universe("test.psf", dcd_files[0])
    reference.atoms.positions = avg_positions
    ref_atoms = reference.select_atoms('backbone')
    rmsd_results = []

    for dcd in dcd_files:
        u = mda.Universe("test.psf", dcd)
        mobile_atoms = u.select_atoms('backbone')
        R = rms.RMSD(mobile_atoms, ref_atoms, select='backbone')
        R.run()
        rmsd_values = R.rmsd[:, 2]

[RMeli]

Hi @royyraa. What is the error you are getting?

Not sure if this is a copy-paste issue, but the indentation of your code is not what I would expect. Would you mind fixing it so that we can try to understand what the real issue is?

[royyraa]

    universes = [mda.Universe("test.psf", dcd) for dcd in dcd_files]
    n_atoms = universes[0].atoms.n_atoms
    avg_positions = np.zeros((n_atoms, 3))

    total_frames = 0

    for u in universes:
        for ts in u.trajectory:
            avg_positions += u.atoms.positions
            total_frames += 1

    if total_frames == 0:
        print(f"No frames found in {folder_name}")
        continue

    avg_positions /= total_frames

    # Use one of the existing DCD files to create the reference universe
    reference = mda.Universe("test.psf", dcd_files[0])
    reference.atoms.positions = avg_positions

    # Select the atom group for the reference universe
    ref_atoms = reference.select_atoms('backbone')

    if len(ref_atoms) == 0:
        print(f"No backbone atoms found in the reference universe for {folder_name}")
        continue

    # Iterate over all DCD files and perform RMSD calculations
    rmsd_results = []
    for dcd in dcd_files:
        u = mda.Universe("test.psf", dcd)
        mobile_atoms = u.select_atoms('backbone')

        if len(mobile_atoms) == 0:
            print(f"No backbone atoms found in the mobile universe for {dcd}")
            continue

        # Perform the RMSD calculation
        R = rms.RMSD(mobile_atoms, ref_atoms, select='backbone')
        R.run()

        # Store the RMSD values
        rmsd_values = R.rmsd[:, 2]  # Extract the RMSD values

        if len(rmsd_values) == 0:
            print(f"RMSD calculation resulted in empty data for {dcd}")
            continue

        rmsd_results.append(rmsd_values)

[RMeli]

Thanks @royyraa, the indentation looks correct now. What is the error you are getting?

[royyraa]

Sorry for the issue with the indentation. Here is the properly indented code. I believe that the average structure is not being used as a reference in the RMSD calculation. Instead, it seems to be using the first DCD file in the array. I hope this helps clarify my issue. Please let me know if you have any other questions

[RMeli]

Your code can be run by using the PSF, DCD and DCD2 files provided by MDAnalysis tests:

    import MDAnalysis as mda
    from MDAnalysisTests.datafiles import PSF, DCD, DCD2
    from MDAnalysis.analysis import rsm

    dcd_files = [DCD, DCD2]

    universes = [mda.Universe(PSF, dcd, frames=10) for dcd in dcd_files]
    n_atoms = universes[0].atoms.n_atoms
    avg_positions = np.zeros((n_atoms, 3))

    total_frames = 0

    for u in universes:
        for ts in u.trajectory:
            avg_positions += u.atoms.positions
            total_frames += 1

    avg_positions /= total_frames

    # Use one of the existing DCD files to create the reference universe
    reference = mda.Universe(PSF, dcd_files[0])
    reference.atoms.positions = avg_positions

    # Select the atom group for the reference universe
    ref_atoms = reference.select_atoms('backbone')

    # Iterate over all DCD files and perform RMSD calculations
    rmsd_results = []
    for dcd in dcd_files:
        u = mda.Universe(PSF, dcd)
        mobile_atoms = u.select_atoms('backbone')

        if len(mobile_atoms) == 0:
            print(f"No backbone atoms found in the mobile universe for {dcd}")
            continue

        # Perform the RMSD calculation
        R = rms.RMSD(mobile_atoms, ref_atoms, select='backbone')
        R.run()

        # Store the RMSD values
        rmsd_values = R.rmsd[:, 2]  # Extract the RMSD values

        rmsd_results.append(rmsd_values)

I run out of time, but I'll investigate in the coming days.