A question about pytraj.align()

[1heman]

Hi,
thanks to making-it-rain, I have been introduced to AmberTools and pytraj.
Prior to this, I have been a GROMACS user.

Unlike GROMACS, where we have to perform 2 operations on the trajectory namely centering and fitting, there seems to be only one here: pytraj.align.

I am particularly not sure how accurate this function is. I have been simulating protein-ligand systems. In many cases, in a visualizer, I have found the ligand to be all over the place but the RMSD of the ligand and the protein speak otherwise.

This is how I have been employing the function:
pt.align(traj_load, ref=0)
Note that I have not given any mask assuming that this would align the trajectory w.r.t. all atoms of the initial frame.
I would also be very pleased to know how to select only the atoms of the protein and ligand for this alignment process.

Thank you,
Hemant

[drroe]

I am particularly not sure how accurate this function is. I have been simulating protein-ligand systems. In many cases, in a visualizer, I have found the ligand to be all over the place but the RMSD of the ligand and the protein speak otherwise.

We will be better able to help if you could post your exact input (commands etc).

I would also be very pleased to know how to select only the atoms of the protein and ligand for this alignment process.

See section 3.2 on Atom mask selection syntax in the CPPTRAJ manual.

[1heman]

I am using a colab notebook created by making-it-rain

I am not sure how helpful these commands by themselves will be but here they are:

#@title **Concatenate and align the trajectory**
Skip = "1" #@param ["1", "2", "5", "10", "20", "50"]
stride_traj = Skip
Output_format = "dcd" #@param ["dcd", "pdb", "trr", "xtc"]

simulation_time_analysis = stride_time_ps*nstride
simulation_ns = float(Stride_Time)*int(Number_of_strides)
number_frames = int(simulation_time_analysis)/int(Write_the_trajectory)
number_frames_analysis = number_frames/int(stride_traj)


traj_end = os.path.join(workDir, str(Jobname) + "_all.dcd")
traj_end2 = os.path.join(workDir, str(Jobname) + "_all." + str(Output_format))
template =  os.path.join(workDir, str(Jobname) + '_%s.dcd')

flist = [template % str(i) for i in range(1, nstride + 1)]
#print(flist)

trajlist = pt.load(flist, pdb, stride=stride_traj)
traj_image = trajlist.iterframe(autoimage=True, rmsfit=0)
traj_write = pt.write_traj(traj_end, traj_image, overwrite=True)
traj_load = pt.load(traj_end, pdb)
traj_align = pt.align(traj_load, ref=0) 
traj_write = pt.write_traj(traj_end, traj_align, overwrite=True, options='dcd')
traj_write = pt.write_traj(traj_end2, traj_align, overwrite=True, options=Output_format)
traj_load = pt.load(traj_end, os.path.join(workDir, "SYS_gaff2.prmtop"))
print(traj_load)


traj_end_check = os.path.exists(traj_end2)

if traj_end_check == True:
  print("Trajectory concatenated successfully! :-)")
else:
  print("ERROR: Check your inputs! ")

There, in the Load, view and check the trajectory subsection, the above lines are contained.

The full notebook can be found here

[hainm]

hi @1heman: I am sorry for my delay reply (I normally check github notification via my email but github did not seem to send this thread to my email).

I am particularly not sure how accurate this function is. I have been simulating protein-ligand systems. In many cases, in a > visualizer, I have found the ligand to be all over the place but the RMSD of the ligand and the protein speak otherwise.

To make sure the code correct, please specify a protein segment, saying residue :1-12, write the aligned trajectory and visualize it.

pt.align(traj_load, ref=0, mask=':1-12')

Note that I have not given any mask assuming that this would align the trajectory w.r.t. all atoms of the initial frame.

Yes, all atoms.

[hainm]

@1heman: FYIpytraj does have center function (but like you said, pytraj.align (cpptraj under the hood) combines both center and fitting)

pytraj/pytraj/all_actions.py

Lines 2164 to 2206 in 30d43ce

	def center(traj=None,
	mask="",
	center='box',
	mass=False,
	top=None,
	frame_indices=None):
	"""Center coordinates in `mask` to specified point.
	Parameters
	----------
	traj : Trajectory-like or Frame iterator
	mask : str, mask
	center : str, {'box', 'origin', array-like}, default 'box'
	if 'origin', center on coordinate origin (0, 0, 0)
	if 'box', center on box center
	if array-like, center on that point
	mass : bool, default: False
	if True, use mass weighted
	top : Topology, optional, default: None
	Examples
	--------
	>>> import pytraj as pt
	>>> traj = pt.datafiles.load_tz2_ortho()
	>>> # load all frames to memory so we can 'mutate' them
	>>> traj = traj[:]
	>>> # all atoms, center to box center (x/2, y/2, z/2)
	>>> traj = pt.center(traj)
	>>> # center at origin, use @CA
	>>> traj = pt.center(traj, '@CA', center='origin')
	>>> # center to box center, use mass weighted
	>>> traj = pt.center(traj, mass=True)
	>>> traj = pt.center(traj, ':1', mass=True)
	Returns
	-------
	updated traj
	See also
	--------
	pytraj.translate