Page 1 of 2 1 2 >
Topic Options
#346 - 10/23/03 06:46 AM rmsf-residue.str
lennart Online   content

Forum Member

Registered: 09/25/03
Posts: 4741
Loc: ~ 59N, 15E
*FILENAME: rmsf-residue.str
*PURPOSE: average atomic RMS fluctuations (around average structure from trajectory)
* per residue for protein non hydrogen atoms
*AUTHOR: Lennart Nilsson, Karolinska Institutet (October 8, 2003)
*
!ASSUMPTIONS:
! 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 protein 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 prot end sele segid prot end
! dele atom sele .not. segid prot end ! this removes all but the protein
! 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 residue
! NB: standard top_all22_na atom names are assumed for the backbone atoms
!RESULT:
! This file calculates nucleotide averaged RMS fluctuations
! The average rmsf for the whole residue, backbone and sidechainis written to unit 21
! N- and C-terminal atoms are included with the sidchains
! Bomblev is set to -1 to avoid stopping on empty selections (Gly has no sidechain)
!LOCAL VARIABLES:
! R, RALL,RBACK,RSIDE,I,N
! Uses scalar storage arrays 3,4,5,6
! Defines selections HEAVY and BACKBONE
!
!USAGE EXAMPLE:
!read rtf card name /applic/charmm/toppar/top_all22_prot.inp
!read para card name /applic/charmm/toppar/par_all22_prot.inp
!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 residue 14
!open unit 21 write form name rmsf.dat
!stream rmsf-residue.str

define heavy sele .not. hydrogen end
define backbone sele heavy .and. (type N .or. type c .or. type ca .or. type o ) end
! save main coordinates
scalar sca3 = x
scalar sca4 = y
scalar sca5 = z
scalar sca6 = wmain

! average per residue, all atoms
scalar x = wmain
scalar x aver byres sele all end
! average per residue, backbone atoms
scalar y = wmain
scalar y aver byres sele backbone end
! average per residue, sidchain atoms
scalar z = wmain
scalar z aver byres -
sele heavy .and. .not. backbone end


! now print it out
set bomlev -1
set i 1
set n @r1
write title unit 21
* IRES RMS_ALL RMS_BACKB RMS_SIDE
*
label LOOP
coor stat sele ires @i .and. type CA 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 CB 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

return


Edited by lennart (04/14/10 05:40 AM)
Edit Reason: Corrected error in read psf in example
_________________________
Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden

Top
#347 - 02/01/05 06:45 AM Re: rmsf-residue.str [Re: lennart]
cathy Offline
Forum Member

Registered: 02/01/05
Posts: 31
thank you for the scipt you provided. i can run this script onto my system. but i am quite new in using charmm. i wonder how i can use the data i got. for example, how to plot the rmsf of the CA atoms as a function of residue number.


Edited by cathyhuang (02/01/05 08:46 AM)

Top
#348 - 02/02/05 01:41 AM Re: rmsf-residue.str [Re: cathy]
lennart Online   content

Forum Member

Registered: 09/25/03
Posts: 4741
Loc: ~ 59N, 15E
If you modify the definition of "backbone" to only include the CA atoms you will get the CA RMSF in the backbone column; the sidechain column will in this case contain data that is influenced also by non-sidechain atoms.
_________________________
Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden

Top
#349 - 05/10/06 02:35 PM Re: rmsf-residue.str [Re: lennart]
tong Offline
Forum Member

Registered: 01/16/06
Posts: 22
Thank you for a helpful script. I've used this script to calculate the rmsf of my protein. It looks fine except the fixed residues belonging to reservoir zone (I did stochastic boundary condition and fixed some amino acids in reservoir zone).

!========================================
open read formatted unit 27 name "../NEWTRA/traj_mergeall.psf"
read psf card unit 27

! Read reference structure
! here use starting coor

OPEN READ UNIT 20 CARD NAME "../../re-min-ab.crd"
READ COOR UNIT 20 CARD COMP
CLOSE UNIT 20

open unit 51 unfo read name "../NEWTRA/traj_mergeall.dcd"

coor dyna firstu 51 nunit 1 skip 500 -
orient select segid prot end sele segid prot end
delete atom sele segid SOLV end ! we don't want rmsd for the waters
set R1 1 ! this fragment starts with residue 1

open unit 21 write form name rmsf_t1.dat
stream rmsf-residue.str

stop
!==============================================

And this is the x.dat

IRES RMSF_ALL RMSF_BB RMSF_SC
1 0.8221 0.6295 0.6206
2 0.5446 0.4162 0.4724
3 0.3746 0.3382 0.3206
4 0.3856 0.3393 0.0000
5 0.6509 0.3170 0.4725



The rmsf of fixed residues:

113 15.8395 18.0510 14.9435
132 4.1948 5.2551 3.7975
134 8.0672 8.0798 0.0000
135 7.1833 6.9330 0.0000
136 9.1360 7.6044 9.6288
137 6.3501 6.0880 0.0000
138 5.2755 4.8729 0.0000
139 4.0817 4.6815 3.7308


The rmsf values of all other residues look fine only the fixed residues give the problem because I've got the high rmsf values of these residues. I just wonder because in fact the rmsf of fixed aminos should be 0 because they were fixed during the simulation.

Is that something wrong with my input file? How can I get the correct value of the rmsf values of the fixed residues?
Any help is greatly appreciated. Thanks in advance.

Top
#350 - 05/10/06 04:06 PM Re: rmsf-residue.str [Re: lennart]
lennart Online   content

Forum Member

Registered: 09/25/03
Posts: 4741
Loc: ~ 59N, 15E
Have not checked the code, but it may be related to the orienting of the coordinates - it is anyhow not importants. I would use CONS HARM rather than CONS FIX.
_________________________
Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden

Top
#351 - 05/10/06 04:28 PM Re: rmsf-residue.str [Re: lennart]
tong Offline
Forum Member

Registered: 01/16/06
Posts: 22
Thank you Prof. Lennart for your quickly reply.
Anyway, do you think is there the way to get the correct rmsf of the fixed residues? Or can I ignore the high rmsf of these fixed residues (use 0 instead of high values)? Thank you.

Top
#352 - 05/10/06 04:38 PM Re: rmsf-residue.str [Re: tong]
lennart Online   content

Forum Member

Registered: 09/25/03
Posts: 4741
Loc: ~ 59N, 15E
The correct value is known (hint: fixed residues do NOT move). If you absolutely must get those known valuse in the resulting "coordinate" file you can use scalar commands (scalar.doc) to set them as you wish.
_________________________
Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden

Top
#353 - 05/10/06 04:57 PM Re: rmsf-residue.str [Re: lennart]
tong Offline
Forum Member

Registered: 01/16/06
Posts: 22
Thank you so much Prof. Lennart.
After I applied the scalar commands, it works and I've got the correct values of the fixed residues.
Thank you.

Top
#354 - 01/11/08 03:12 AM Re: rmsf-residue.str [Re: lennart]
irina Offline
Forum Member

Registered: 03/14/06
Posts: 2
Dear Lennart,

I use your script to calculate RMSF for my protein, bellow
the input script. Every time I have an error message (see after input), Somehow the script can not be executed properly on selection of the third atom. Do you have any ideas why? I checked the reference, it has the third atom. I also did not do any changes in the script.
Thank you!
Irina

set direcCHARMM /data/qqq
open unit 1 read form name @direcCHARMM/charmm27.top
read rtf card unit 1
close unit 1
open unit 1 read form name @direcCHARMM/charmm27.par
read param card unit 1
close unit 1
open unit 1 read form name @direcCHARMM/qqq.psf
read psf card unit 1
close unit 1

open unit 1 read card name @direcCHARMM/qqq.pdb
read coor pdb unit 1 comp
close unit 1

open read unit 20 unform name @direcCHARMM/"dyn11.dcd"
coor dyna nunits 1 firstu 20 begin 100 skip 1000 stop 500000 -
orient select segid prot end sele segid prot end
set R1 1
open unit 21 write form name rmsf.dat
stream rmsf-residue.str

stop



CHARMM> if i le ?NRES goto loop
RDCMND substituted energy or value "?NRES" to "348"
Comparing " 3" and "348".
IF test evaluated as true. Performing command

CHARMM> coor stat sele ires @i .and. type CA end
Parameter: I -> " 3"
SELRPN> 1 atoms have been selected out of 5517
STATISTICS FOR 1 SELECTED ATOMS:
XMIN = 1.100687 XMAX = 1.100687 XAVE = 1.100687
YMIN = 1.053481 YMAX = 1.053481 YAVE = 1.053481
ZMIN = 1.032451 ZMAX = 1.032451 ZAVE = 1.032451
WMIN = 1.032451 WMAX = 1.032451 WAVE = 1.032451

CHARMM> format (F12.4)

CHARMM> set rall ?xave
RDCMND substituted energy or value "?XAVE" to " 1.1007"
Parameter: RALL <- "1.1007"

CHARMM> incr rall by 0.0 ! trick to get the format applied to the variable
Parameter: RALL <- " 1.1007"

CHARMM> trim rall to 12
Parameter: RALL <- "1.1007"

CHARMM> set rback ?yave
RDCMND substituted energy or value "?YAVE" to " 1.0535"
Parameter: RBACK <- "1.0535"

CHARMM> incr rback by 0.0
Parameter: RBACK <- " 1.0535"

CHARMM> trim rback to 12
Parameter: RBACK <- "1.0535"

CHARMM>

CHARMM> coor stat sele ires @i .and. type CB end
Parameter: I -> " 3"
SELRPN> 0 atoms have been selected out of 5517

***** LEVEL 0 WARNING FROM <CORMAN> *****
***** ZERO ATOMS SELECTED
******************************************
BOMLEV ( 0) IS REACHED - TERMINATING. WRNLEV IS 5



/---------\
/ \
/ \
/ \
! XXXX XXXX !
! XXXX XXXX !

Top
#355 - 01/11/08 03:18 AM Re: rmsf-residue.str [Re: irina]
lennart Online   content

Forum Member

Registered: 09/25/03
Posts: 4741
Loc: ~ 59N, 15E
This is probably a glycine, which does not have a CB atom. Are you absolutely sure that there is a CB in this residue? If this is really the case, then it is not possible to do anything more about your problem with the information you have provided.
The script is designed to handle this case; it has a BOMLEV -1 command, which you seem to have deleted (without changing the script?).
_________________________
Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden

Top
Page 1 of 2 1 2 >

Moderator:  chmgr, John Legato, petrella