CHARMM Development Project
Posted By: lennart rmsd-nucleotide.str - 10/23/03 10:43 AM
*FILENAME: rmsd-nucleotide.str
*PURPOSE: compute atomic RMS diff between main and comp coordinates
* averaged per nucleotide for nucleic acid non hydrogen atoms
*AUTHOR: Lennart Nilsson, Karolinska Institutet (October 8, 2003)
! Only the relevant segment is present (no solvent, ions,...)
! MAIN and COMP coordinate sets contain the structures to be compared
! Unit 21 is open for formatted write to the file where the result is written
! CHARMM variable R1 set to the (numerical!) resid of the first nucleotide
! NB: standard top_all22_na atom names are assumed for the backbone atoms
! This file calculates the RMS difference between main and comparison
! coordinate sets after least-squares superposition using all heavy atoms
! and averages the rms difference per residue.
! N- anc C-terminal atoms are included with the sidechains.
! The average rmsd for the whole residue, backbone and sidechain is written to unit 21
! Uses scalar storage arrays 1,2,3,4,5,6
!read rtf card name /applic/charmm/toppar/top_all27_na.rtf
!read para card name /applic/charmm/toppar/par_all27_na.prm
!read psf name my_structure.psf
!read coor card name my_first_set.crd
!read coor card name my_second_set.crd comp
!!delete atom sele .not. segid dna end ! we don't want rmsd for the waters
!set R1 14 ! this DNA fragment starts with nucleotide 14
!open unit 21 write form name rmsd.dat
!stream rmsd-nucleotide.str

! superpose
define heavy sele .not. hydrogen end
! the sugar and phosphate groups make up the backbone
define backbone sele heavy .and. (type *P .or. type *' .or. type *T) end

! save main coordinates
scalar sca3 = x
scalar sca4 = y
scalar sca5 = z
scalar sca6 = wmain

coor orie rms sele heavy end
! get difference
coor diff
! and square it
scalar x ipow 2
scalar x store 1
scalar y ipow 2
scalar y +store 1
scalar z ipow 2
scalar z +store 1
scalar wmain recall 1
! now take square root of sum to get RMS
scalar wmain sqrt
! average per nucleotide, all non-hydrogen atoms
scalar x = wmain
scalar x aver byres sele heavy end
! average per nucleotide, backbone atoms
scalar y = wmain
scalar y aver byres -
sele backbone end
! average per nucleotide, base atoms
scalar z = wmain
scalar z aver byres -
sele heavy .and. .not. backbone end

! now print it out
bomlev -1 ! in order not to crash when there are no sidechain atoms in GLY
set i 1
set n @r1
write title unit 21
label LOOP
coor stat sele ires @i .and. type C1' end
format (F12.4)
set rall ?xave
incr rall by 0.0 ! trick to get the format applied to the variable
trim rall to 12
set rback ?yave
incr rback by 0.0
trim rback to 12

coor stat sele ires @i .and. (type N1 .or. type N9) end
set rside ?zave
if ?NSEL eq 0 set rside 0.0
incr rside by 0.0
trim rside to 12
set r @n
format (I6)
incr r by 0.0
trim r to 6
write title unit 21
* @r @rall @rback @rside
incr n by 1
incr i by 1
if i le ?NRES goto loop

format ! best to reset to default ....
! restore main coordinates
scalar x = sca3
scalar y = sca4
scalar z = sca5
scalar wmain = sca6

© CHARMM forums