CHARMM Development Project Forums CHARMM Community Script Archive rmsd-rgyr.inp

*FILENAME: rmsd-rgyr.inp
*PURPOSE: compute rmsd vs initial structure and radius of gyration from trajectory
*AUTHOR: Lennart Nilsson, Karolinska Institutet, October 2003
*
!unix environment variable CHM_HOME points to CHARMM installation directory
read rtf card name $CHM_HOME/toppar/top_all22_prot.inp
read para card name $CHM_HOME/toppar/par_all22_prot.inp
read psf card name my_psf.psf
read coor card name my_reference.crd
! put reference coordinates into comparison set as well
coor copy comp
open unit 11 write form name rmsd.dat
write title unit 11
* time rmsd
*
open unit 12 write form name rgyr.dat
write title unit 12
* time rgyr
*
! assume this is a trajectory where overall translation/rotation has been removed and
! without water or other unnecessary things
open unit 51 read unform name traj1.cor
open unit 52 read unform name traj2.cor
correl maxtime 10000
enter v1 rms
enter v2 gyra

traj firstu 51 nunit 2 begin 10500 skip 500 stop 55000

write v1 dumb time unit 11
*hi
*

write v2 dumb time unit 12
*hi
*

end

dear Lennart,

I try to plot the autocorrelation of the radius of gyration of a pdb file vs. the radius of gyration of the production trajectory.

The purpose of this is to find out the simulation time i need to explore the conformational space of my protein.

Because I'm new in this field, I was glad to see that input files are distributed (that's great!!), but I couldn't manage/adapt your script for my purpose: I was only able to get the rms output and no rgyr output..

----> here i add your rmsd-rgyr.inp modified for my machine

Quote:

---

*FILENAME: rmsd-rgyr.inp
*PURPOSE: compute rmsd vs initial structure and radius of gyration from trajectory
*AUTHOR: Lennart Nilsson, Karolinska Institutet, October 2003
*
*Adapted by rhz, December 2004
*


!----------Alias

SET toppar /usr/local/c31b1/toppar ! Location of top and parameters

!---------------

read rtf card name @toppar/top_all22_prot.inp
read para card name @toppar/par_all22_prot.inp

read psf card name 1l2y.psf
read coor card name 1l2y.crd


! put reference coordinates into comparison set as well

coor copy comp tcage_minimized.pdb
open unit 11 write form name rmsd.dat
write title unit 11
* time rmsd
*

open unit 12 write form name rgyr.dat
write title unit 12
* time rgyr
*

! assume this is a trajectory where overall translation/rotation has been removed and
! without water or other unnecessary things

open unit 51 read unform name 1l2y.dcd


correl maxtimesteps 100000 maxseries 4 maxatoms 1000

enter v1 rms
enter v2 gyra

traj firstu 51 nunit 2 begin 5200 skip 500 stop 100000 sele all end


! write down rgyr data

open unit 51 write form name rgyr_1l2y.dat


write v1 dumb time unit 11
*hi
*
write v2 dumb time unit 12
*hi
*

end

---

I've read the correl.doc of CHARMM, but I can't find out how to extract the information about the radius of gyration..

If you could help me out it would be very very nice

Does the script work for you without your modifications? What are the error messages&symptoms?

I don't know if this helps
1/ You seem to have a trajectory that is contained in two files ("nunit 2"), but you have only opened one of them; you need to open the second one on unit 52

2/ The "open unit 51 write form name rgyr_1l2y.dat" should be removed - units 11 and 12 are already opened and will be used to write the rmsd and rgyr data.

yes, the script works well (no skull) , but the output is not as desired..

anyway:

1/ you are right: I have only 1 trajectory, so I corrected it like you wrote = nunit 1

2/ the command line: open unit 51 write form name rgyr_1l2y.dat is removed


what I've done now is to add the CORFUN function like:

Quote:

---

correl maxtimesteps 100000 maxseries 4 maxatoms 1000

enter v1 rms
enter v2 gyra

traj firstu 51 nunit 1 begin 5200 skip 500 stop 100000 sele all end

corfun v2 v2 fft p1 nltc
write correl plot unit 21
* corfun v2 v2 fft p1 nltc
*

---

if I understand it right, it means that it will do an autocorrelation between v2(0) and v2(t). What I don't understand is: how may I get the right output? since I need to plot the [rgyr(0).rgyr(t)] in function of the the production time and get -instead- only the fort.21 & fort.22 output..

In any case, may I ask you what method would you use to calculate the necessary simulation time for a given protein?

best regards, r_hz

If you are only interested in rgyr, remove the "enter v1 rms".
You can compute the correlation function for a specific part of your
trajectory using BEGIN and STOP keywords to the TRAJ command.
I do not know what kind of output you want, but unless you use the PLT2 plotting program your "write corr ..." will not be very useful. Try this instead:

write corr dumb time unit 21
* this title is not used at all (to get a header do a write title before
* invoking the correl module)
*

From correl.doc:
"The PLOT option will create a BINARY file for plotting by PLT2.
The first line of the title is used as the plot title, but this may be
reset in PLT2.
The DUMB options will simply write out the values with no title
or header to a card file, one value to a line. If the TIME option is
specified, the time value will preceed the time series values"

Note also that with some versions of the Portland fortran compiler there is a bug in the optimizer, which makes the fft option compute an erroneous correlation function; the workaround is to use direct instead of fft.

The quick answer about the necessary length of a simulation is that you should continue until a steady state is reached in a plot of the quantity you are interested in vs time, and then the production part has to be long enough that averages can be obtained with small enough errors to make your results significant. There is a chapter on this in Allen&Tildesley "Computer Simulations of Liquids".

Thank you very very much for your enlightments !!

r_hz

Dear all,
I used the script to caculate rmsd . But the error has happened . Could you tell me how to solve it . THX
yeng-tseng

my script
___________________________________________________
*FILENAME: cpt-equil.inp
*PURPOSE: setup and run Constant Pressure&Temperature (rectangular box) MD simulation
*AUTHOR: Wang Yeng-tseng
*
!unix environment variable CHM_HOME points to CHARMM installation directory
read rtf card name top_all22_prot.inp
read para card name par_all22_prot.inp
read psf card name 4pti-solv.psf
read coor card name 4pti-solv.crd
! put reference coordinates into comparison set as well
coor copy comp
open unit 11 write form name rmsd.dat
write title unit 11
* time rmsd
*

! assume this is a trajectory where overall translation/rotation has been removed and
! without water or other unnecessary things
open unit 51 read unform name 4pti-heat.cor

correl maxtime 10000
enter v1 rms


traj firstu 51 nunit 1 begin 1 end 5500000

write v1 dumb time unit 11
*hi
*


end
stop
________________________________________________

error
_________________________________________________
CORREL> traj firstu 51 nunit 1 begin 1 end 5500000

***** LEVEL 0 WARNING FROM *****
***** SELEction keyword missing
******************************************
BOMLEV ( 0) IS REACHED - TERMINATING. WRNLEV IS 5

If you had not changed the STOP on the traj command to an END it should have worked, if also the psf matches your trajectory; if the trajectory contains solvent and protein the resulting rmsd(t) will rather meaningless.

Dear Lennart ,
I have already removed the water form *.cor or *.dcd
file. But the error message still happened , could you tell me why ?
THX

__________________________________________________
*FILENAME: rmsd-rgyr.inp
*PURPOSE: compute rmsd vs initial structure and radius of gyration from trajectory
*AUTHOR: Lennart Nilsson, Karolinska Institutet, October 2003
*
!unix environment variable CHM_HOME points to CHARMM installation directory
read rtf card name top_all22_prot.inp
read para card name par_all22_prot.inp
read psf card name 4pti_no_wat.psf
read coor card name 4pti_no_wat.crd
! put reference coordinates into comparison set as well
coor copy comp
open unit 11 write form name rmsd.dat
write title unit 11
* time rmsd
*

! assume this is a trajectory where overall translation/rotation has been removed and
! without water or other unnecessary things
open unit 51 read unform name 4pti-equil_nwat.cor
correl maxtime 10000
enter v1 rms

traj firstu 51 nunit 1 begin 1 skip 1 end 500

write v1 dumb time unit 11
*hi
*
end


___________________________________________________
General atom nonbond list generation found:
185128 ATOM PAIRS WERE FOUND FOR ATOM LIST
7480 GROUP PAIRS REQUIRED ATOM SEARCHES


CORREL> enter v1 rms

CORREL>

CORREL> traj firstu 51 nunit 1 begin 1 skip 1 end 500

***** LEVEL 0 WARNING FROM *****
***** SELEction keyword missing
******************************************
BOMLEV ( 0) IS REACHED - TERMINATING. WRNLEV IS 5

I can only repeat what I said in the previous post: "If you had not changed the STOP on the TRAJ command to an END it should have worked".
Note also that BEGIN, SKIP and STOP should be followed by numbers specifying integration steps (counted from the STARt of the simulation), not the sequential number of the coordinate set in the trajectory; you probably do not want 1 to 500.