*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
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
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.
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
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?
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".
Last edited by lennart; 12/01/0404:02 PM.
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
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.
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
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.
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
from the example 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 end 55000
why there is two *.cor files , both these files have same values repeated twice or values of nsavc starting from 10500 to 55000 partly filled in each other.
The latter is the usual case; perhaps steps 500 thru 30000 are in the first file, and steps 30500 thru 60000 are in the second, so both are needed to analyze 10500 thru 55000.
I ran the script using psf, pdb, and traj made from your merge script and deleted water. The error I get running the rmsd-rgyr script is: At line 959 of file dynio.f Fortran runtime error: End of file
The output file is attached. This is the input file:
Code:
!*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
stream toppar.str
read psf card name "new_psf.psf"
open unit 2 card read name -
"test.pdb"
read coor pdb resid unit 2
! 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 -
"traj_all.cor"
correl maxtime 10000000
enter v1 rms
enter v2 gyra
traj firstu 51
write v1 dumb time unit 11
*hi
*
write v2 dumb time unit 12
*hi
*
end
***** LEVEL -2 WARNING FROM ***** ***** Apparent problem with endianness ****************************************** BOMLEV ( 0) IS REACHED - TERMINATING. WRNLEV IS 2
