to do that you need the bond fluctuations along your trajectory, then you can calculate this manually with respect to the equilibrium bond length. in VMD software, choose the two atoms you are interested in (press 2 and then select two atoms), go to graphics--> select labels--> select bonds, highlight the two selected atoms, in the same window, select save. the saved file will contain the bond length at every time step. then RMSD=sqrt (.......), you are done. 

I don't know why you do this. Many program can do this very well. I prefer AMBERTOOLS. It could handle most trajectory format. During the RMSD calculation, the program will adjust the molecular orientation to get the lowest RMSD.

to do that you need the bond fluctuations along your trajectory, then you can calculate this manually with respect to the equilibrium bond length. in VMD software, choose the two atoms you are interested in (press 2 and then select two atoms), go to graphics--> select labels--> select bonds, highlight the two selected atoms, in the same window, select save. the saved file will contain the bond length at every time step. then RMSD=sqrt (.......), you are done. 

Sufian, thanks for your answer, I will give that a go.

Hi James, 

As pointed out, you may do this in available software, but incase you want to do your own coding I will provide a bare-bones answer.

You can easily rig up a python, Mathematica, or even BASH script for this. Start by writing code to find each species you're interested in from the input file (e.g. find C1 and H1). This is done by either searching the coordinate file for either the label (C1 or H1), or by specifying its place in the ensemble. Calculate the displacement (bond length), d, between these atom coordinates in each frame (or equivalently for each time t), as

d(t)= ((C1x(t)-H1x(t))2+( C1y(t)-H1y(t))2+( C1z(t)-H1z(t))2)0.5 ,

then upon defining either the starting displacement (initial configuration) or the equilibrium displacement (final configuration), you may perform the root mean square displacement as a function of time t,

RMSD(t, 0 or eq) =(N-1Σi=1Ν(di(t) - d(0 or eq))2)0.5,

where N is the total number of frames (or total time).

Hope this helps!

Cheers, 

Quinn

Thanks very much Quinn that is a great answer.

I do want to code it, as it will become part of a bigger software package, I became confused when search online as I was often referred to algorithms involving rotating the molecules to a optimal superposition, which is not relevant in this case. Your answer is a great clarification thank you. 

To evaluate bond length variation along a MD trajectory, remember that your integration step should be of 1 fs or lower, and that you should not use bond constraints (this is a obvious statement, but just to be sure). Additionally, if your topology is correct, you should observe just a small variation along simulation, without significance for any practical purpose. For example, bond lengths are very frequently kept constrained in macromolecule simulations, without relevant negative impacts. But to validate a topology, geometric determinants are only the start, and should be succeeded by the calculation of physical chemical properties in vacuum or condensed phase, depending on our interests.

Hi Hugo, thanks for your comment, to answer some of your concerns:

We are only interested in the geometrical sampling space covered in the simulations currently, hence, this information is important to us and of significance; it will play a practical role in our research.

We are using steps of less than a femto second and do not apply bond constraints in the simulations. As you say, we see very small deviations in the bond length over the simulation.

Thank you for your comment. 

You can do it in VMD ( you get simple script from the net to do so) and in my understanding you can also use GROMACS tool for this!!