CHARMM Development Project Forums User Discussion & Questions Other User Discussion and Questions How to calculate RMSD?

Dear all :
I finished my charmm job, and I wanted to analyze the RMSD of the tracjectory file (*.dcd). So, my input file is as follows:
* Script file produced by INSIGHT
*
! Startup script for CHARMm
!
UPPER ! case for file to write
BOMBLEVEL -6
WRNLEV 0
PRNLEV 5
! Script to read parameter, psf and crd file
!
OPEN READ UNIT 21 CARD NAME "top_all27_prot_na.rtf"
READ RTF UNIT 21 CARD
CLOSE UNIT 21
OPEN READ UNIT 21 CARD NAME "charmm27_support.rtf"
READ RTF UNIT 21 CARD APPEND
CLOSE UNIT 21

OPEN READ UNIT 20 CARD NAME "par_all27_prot_na.prm"
READ PARA UNIT 20 CARD
CLOSE UNIT 20

OPEN READ UNIT 20 CARD NAME "hsa_5h2ond9_1.psf"
READ PSF UNIT 20 CARD
CLOSE UNIT 20

!OPEN READ UNIT 20 CARD NAME "hsa_5h2ond9_1.crd"
!READ COOR UNIT 20 CARD
!CLOSE UNIT 20

open unit 30 write form name "nonionnd9_rmsd.dat"

write title unit 30
*# Time RMSD
*


open unit 22 read unform name "hsa_5h2ond9_1.dcd"

traj iread 22 nread 1
set count = 1

label read

traj read

if count ne 1 goto skipfirst
define hsa select (segid HSA_ .and. .not. resname TIP3) end
coor copy comp select hsa end
label skipfirst

coor rms comp select hsa end
set rmsd = ?rms
write title unit 30
* @count @rmsd
incr count by 1
if count le 10000 goto read
stop

I used charmm on linux.
But it terminated in several seconds, and did't calculate all frames RMSD.
What's wrong with my input file?

thanks all !!

I cannot judge from your description what the problem is. There is however an easier way compute RMSD(t), see Script Archive, or correl.doc. BOMBLEVEL -6 also seems to be a rather dangerous setting.

I believe that you want the command;

COOR ORIE RMS ...

not

COOR RMS COMP ....

This will do the rms bestfit. If you really want an absolute RMSD then your command is fine. Assuming that you selected the correct number of atoms (did you?), then your script looks OK, but not optimal.

There is also an RMS time series in CORREL which is a better choice.