Error when analysing trajectory object

[kellegra]

Hi there,

I am trying to analyze a trajectory object. Running a single trajectory passes and all plots can be generated. But when using the trajectory object by using tr.analysis_all_namd() it throws out the error shown in the bottom. Not sure how to deal with this, since other trajectories work. The trajectories that worked, were generated following the gromacs-ramd workflow. The trajectories that do not work were generated following the the initial RAMD paper from 2018. But not sure, why this would end up in this issue.

I hope, you have any ideas how I could fix this problem.

Best regards, Gerald

`ValueError                                Traceback (most recent call last)
Input In [15], in <cell line: 13>()
     11 step =3 
     12 # We start the analysis of the namd data and define ...
---> 13 tr.analysis_all_namd(WB_analysis = True,auxi_selection = auxi_selection,step_analysis=step, start_analysis=end)

File ~/tramd/MD-IFP-1.2/MD-IFP/Scripts/Trajectories.py:926, in trajectories.analysis_all_namd(self, WB_analysis, step_analysis, start_analysis, RE, Lipids, auxi_selection)
    924 print("\n\n>>>>>>>>>>>>>>>>>","Replica: ",repl, "\n")
    925 step = max(step_analysis, 1)
--> 926 length,start,rmsd_prot,rmsd_lig, rmsd_auxi,Rgr_prot, Rgr_lig,com_lig,df_prop,df_HB,df_WB= self.analysis_traj(repl,start_analysis,step,WB_analysis, RE,Lipids,auxi_selection,reference ="ref")
    927 df_prop_complete = table_combine(df_HB,df_WB,df_prop,sel_ligands,self.namd.contact_collection)   
    928 self.namd.length.append((self.timestep/1000)*length)

File ~/tramd/MD-IFP-1.2/MD-IFP/Scripts/Trajectories.py:882, in trajectories.analysis_traj(self, traj, start_analysis, step_analysis, WB_analysis, RE, Lipids, auxi_selection, reference)
    880 Rgr_prot.append(u_mem.select_atoms("protein").radius_of_gyration())
    881 Rgr_lig.append(u_mem.select_atoms("resname "+sel_ligands).radius_of_gyration())
--> 882 rmsd = superimpose_traj(ref,u_mem,selection_rmsd)  
    883 if(rmsd[0] > 10.0): print("for the frame %s protein RMSD is very large: %s" %(i,rmsd[0]))
    884 if(rmsd[1] > 10.0): print("for the frame %s ligand RMSD is very large: %s" %(i,rmsd[1]))

File ~/tramd/MD-IFP-1.2/MD-IFP/Scripts/Trajectories.py:1364, in superimpose_traj(ref, u, sel_list)
   1362 u.atoms.translate(-u_CA.center_of_mass())
   1363 u0 =  u_CA.positions - u_CA.center_of_mass()
-> 1364 R, rmsd = align.rotation_matrix(u0,ref0)  # compute rotation matrix
   1365 u.atoms.rotate(R)
   1366 u.atoms.translate(ref_CA.center_of_mass()) # translate back to the old center of mass position

File ~/.conda/envs/md/lib/python3.8/site-packages/MDAnalysis/analysis/align.py:273, in rotation_matrix(a, b, weights)
    270 b = np.asarray(b, dtype=np.float64)
    272 if a.shape != b.shape:
--> 273     raise ValueError("'a' and 'b' must have same shape")
    275 if np.allclose(a, b) and weights is None:
    276     return np.eye(3, dtype=np.float64), 0.0

ValueError: 'a' and 'b' must have same shape

[kellegra]

Okay, I checked what the align.py of the MDAnalysis package is actually doing. It seems that there is an issue with the selection I did for the auxi_selection used in tr.analysis_all_namd() function. It fails with the reported error when using

auxi_selection = ["backbone"]

I changed the selection to

auxi_selection = ["((resid 1:353) and (name C CA N))"]

what made the whole thing working. I think there is a issue with the selection "backbone", since my protein has two phosphorylated residues, but not sure if this is really causing the problem.

In the end, the different selection made it.