if(!-var msf) getvar msf msg="No MSF file specified"

if(!-var mds) getvar mds msg="No MDS file specified"

if(!-var trj) getvar trj msg="No Trj file specified"

if(!-var pdb) getvar pdb msg="No ref PDB file specified"

#load force field

#include "loadff.inc"

#load md system

loadmsf file="$msf"

#load a reference for RMSD calculation

loadpdb file="$pdb" copy2ref1=true

#load MD core and setup

loadsetup file=$mds

#feed in trj

trjin file="$trj"

#remove files if already exists

if(-f rmsd.txt) exec rm -rf rmsd.txt

if(-f bbrmsd.txt) exec rm -rf bbrmsd.txt

exec rm -rf rmsd_h*.txt

exec rm -rf bbrmsd_h*.txt

#define mainpro

defcls main select="segn PRO* DNA* RNA* && not hydrogen"

defcls backbone select="(segn DNA* RNA* && name C3' C4' C5' P) \

                     || (segn PRO* && name C CA N O)"

defcls helix1 select="segn PRO* && seq 43:51"

defcls helix2 select="segn PRO* && seq 54:61"

defcls helix3 select="segn PRO* && seq 62:73"

defcls helix12 select="segn PRO* && seq 43:61"

defcls helix23 select="segn PRO* && seq 54:73"

ifrm=0

while($ifrm<$ntotframes){

  trjreadnextframe

  mdtime=$curmdtime

  if($mdtime<0) goto done #double check whether trjreadnextframe successfully or not

  #print progress for every 10ps

  tmod=$mdtime%10

  if($tmod==0) echo -hs "current MD time: $mdtime ps"

  #calculate RMSD for all heavy atoms comparing to initial structure

  calc rmsd fitgroup="main && helix1" calcgroup="main && helix1"

  echo "$mdtime $value" >> rmsd_h1.txt

  #calculate RMSD for main chain atoms

  calc rmsd fitgroup="backbone && helix1" calcgroup="backbone && helix1"

  echo "$mdtime $value" >> bbrmsd_h1.txt

  #calculate RMSD for all heavy atoms comparing to initial structure

  calc rmsd fitgroup="main && helix2" calcgroup="main && helix2"

  echo "$mdtime $value" >> rmsd_h2.txt

  #calculate RMSD for main chain atoms

  calc rmsd fitgroup="backbone && helix2" calcgroup="backbone && helix2"

  echo "$mdtime $value" >> bbrmsd_h2.txt

  #calculate RMSD for all heavy atoms comparing to initial structure

  calc rmsd fitgroup="main && helix3" calcgroup="main && helix3"

  echo "$mdtime $value" >> rmsd_h3.txt

  #calculate RMSD for main chain atoms

  calc rmsd fitgroup="backbone && helix3" calcgroup="backbone && helix3"

  echo "$mdtime $value" >> bbrmsd_h3.txt

  #calculate RMSD for all heavy atoms comparing to initial structure

  calc rmsd fitgroup="main && helix12" calcgroup="main && helix12"

  echo "$mdtime $value" >> rmsd_h12.txt

  #calculate RMSD for main chain atoms

  calc rmsd fitgroup="backbone && helix12" calcgroup="backbone && helix12"

  echo "$mdtime $value" >> bbrmsd_h12.txt

  #calculate RMSD for all heavy atoms comparing to initial structure

  calc rmsd fitgroup="main && helix23" calcgroup="main && helix23"

  echo "$mdtime $value" >> rmsd_h23.txt

  #calculate RMSD for main chain atoms

  calc rmsd fitgroup="backbone && helix23" calcgroup="backbone && helix23"

  echo "$mdtime $value" >> bbrmsd_h23.txt

  ifrm++

}

stop

:done

echo -hs "something went wrong when reading trjactory data"

abort