Calculating Center for a subunit on a residue

[Onyema Augustine Chimezie]

I ran a simulation on a protein-DNA system and I want to calculate the center of mass of the DNA fragment alone. So i used the command as shown below. It had no error but I can print or list the center of mass of the selected residue. I want to get the center of mass from first frame to last frame. 

time = [ ]

cm = [ ]

dna1=u.select_atoms("resid 1 to 147")

for ts in u.trajectory:

    time = u.trajectory.time

    cm = dna1.atoms.center_of_mass()

    print (cm)

[Oliver Beckstein]

Hello Onyema Augustine Chimezie,

On Oct 12, 2021, at 8:24 AM, Onyema Augustine Chimezie <onyemaaugus...@gmail.com> wrote:

I ran a simulation on a protein-DNA system and I want to calculate the center of mass of the DNA fragment alone. So i used the command as shown below. It had no error but I can print or list the center of mass of the selected residue.

In the code below you’re calculating the center of mass over multiple residues, i.e., [x_com, y_com, z_com], i.e., a 1x3 array. I assume that this is really what you want.

If instead you want the centers of mass of all 147 residues as a 147x3 array then you can use dna1.center_of_mass(compound=“residues”) as described in the docs https://docs.mdanalysis.org/stable/documentation_pages/core/groups.html#MDAnalysis.core.groups.AtomGroup.center_of_mass . Note that if your molecule is broken across periodic boundaries from simulations run under periodic boundary conditions (PBC) then you will likely want to add the unwrap=True keyword to correctly calculate the COMs: dna1.center_of_mass(compound=“residues, unwrap=True). Whenever you have to deal with PBC you should carefully check your output!!!

In the following I just keep your original COM calculation but you MIGHT have to also add unwrap=True. Check.

I want to get the center of mass from first frame to last frame. 

You can store the center of mass in a list or numpy array by appending to a list and then converting the list to an array at the end:

import numpy as np

time = []

cm = []

dna1 = u.select_atoms("resid 1 to 147")

for ts in u.trajectory:

    time.append(u.trajectory.time)

    cm.append(dna1.atoms.center_of_mass())

time = np.array(time)

cm = np.array(cm)

There are fancier ways to solve this problem. For a start, you don’t need the intermediate list but you can instead allocate the array in advance:

import numpy as np

time = np.zeros(u.trajectory.n_frames, dtype=np.float32)

cm = np.zeros((u.trajectory.n_frames, 3), dtype=np.float32)

dna1 = u.select_atoms("resid 1 to 147")

for i, ts in enumerate(u.trajectory):

    time[i] = u.trajectory.time

    cm[i] = dna1.atoms.center_of_mass()

# now work with time and cm array

You can also create your own analysis class (similar to the built-in analysis tools in MDAnalysis.analysis) using “AnalysisFromFunction” as explained in more detail in the User Guide https://userguide.mdanalysis.org/stable/examples/analysis/custom_trajectory_analysis.html ):

from MDAnalysis.analysis.base import AnalysisFromFunction

def center_of_mass(ag):

    return ag.center_of_mass()

    

com = AnalysisFromFunction(center_of_mass, u.trajectory, dna1).run(verbose=True)

print(com.results.times)       # like time array before

print(com.resulys.timeseries)  # like cm array before 

Hope this helps.

Oliver

--

Oliver Beckstein (he/his/him)

email: orbe...@mdanalysis.org
twitter: @orbeckst

MDAnalysis – a NumFOCUS fiscally sponsored project

https://www.mdanalysis.org/

[Onyema Augustine Chimezie]

Thanks Oliver,

It really helped me. I am using your guide on other systems.

Thanks once again.