*FILENAME: rmsf-nucleotide.str
*PURPOSE: average atomic RMS fluctuations (around average structure from trajectory)
* 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,...)
! A "COOR DYNA" COMMAND has been executed prior to running this script
! in order to get the relevant fluctuations into the wmain array,
! eg (where only the DNA segment will be used, and we also orient to remove
! unwanted overall rotation/translation):
!! NB: Both selections are needed, and the orient has to be first!!
! coor dyna <traj-spec> orient sele segid dna end sele segid dna end
! dele atom sele .not. segid dna end ! this removes all but the DNA
! 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 nucleotide averaged RMS fluctuations
! The average rmsf for the whole nucleotide, backbone and base is written to unit 21
! Uses scalar storage arrays 3,4,5,6
! Defines selections HEAVY and BACKBONE
!read rtf card name /applic/charmm/toppar/top_all27_na.rtf
!read para card name /applic/charmm/toppar/par_all27_na.prm
!read psf card name my_structure.psf
!read coor card name reference.crd comp ! for the orientation operation
!open unit 51 read unform name traj1.cor
!open unit 52 read unform name traj2.cor
!coor dyna firstu 51 nunit 52 begin 15500 skip 500 end 350000 -
! orient select segid main end sele segid main end
!!delete atom sele .not. segid main end ! we don't want rmsd for the waters
!set R1 14 ! this fragment starts with nucleotide 14
!open unit 21 write form name rmsf.dat
!stream rmsf-nucleotide.str

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

! average per nucleotide, all atoms
scalar x = wmain
scalar x aver byres sele all 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

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


Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden