#360 - 10/23/03 06:47 AM
rmsd-rgyr.inp
|
Forum Member
Registered: 09/25/03
Posts: 4755
Loc: ~ 59N, 15E
|
*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
_________________________
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
|
Top
|
|
|
|
#361 - 12/01/04 07:34 AM
Re: rmsd-rgyr.inp
[Re: lennart]
|
Forum Member
Registered: 12/01/04
Posts: 7
Loc: lausanne
|
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 
|
Top
|
|
|
|
#362 - 12/01/04 08:15 AM
Re: rmsd-rgyr.inp
[Re: r_hz]
|
Forum Member
Registered: 09/25/03
Posts: 4755
Loc: ~ 59N, 15E
|
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
|
Top
|
|
|
|
#363 - 12/01/04 09:54 AM
Re: rmsd-rgyr.inp
[Re: lennart]
|
Forum Member
Registered: 12/01/04
Posts: 7
Loc: lausanne
|
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
Edited by r_hz (12/01/04 10:13 AM)
|
Top
|
|
|
|
#364 - 12/01/04 10:24 AM
Re: rmsd-rgyr.inp
[Re: r_hz]
|
Forum Member
Registered: 09/25/03
Posts: 4755
Loc: ~ 59N, 15E
|
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".
Edited by lennart (12/01/04 11:02 AM)
_________________________
Lennart Nilsson Karolinska Institutet Stockholm, Sweden
|
Top
|
|
|
|
#366 - 02/24/05 07:08 AM
Re: rmsd-rgyr.inp
[Re: lennart]
|
Forum Member
Registered: 12/30/03
Posts: 80
|
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
|
Top
|
|
|
|
#367 - 02/24/05 08:20 AM
Re: rmsd-rgyr.inp
[Re: label]
|
Forum Member
Registered: 09/25/03
Posts: 4755
Loc: ~ 59N, 15E
|
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
|
Top
|
|
|
|
#368 - 02/25/05 04:06 AM
Re: rmsd-rgyr.inp
[Re: lennart]
|
Forum Member
Registered: 12/30/03
Posts: 80
|
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
|
Top
|
|
|
|
#369 - 02/25/05 04:13 AM
Re: rmsd-rgyr.inp
[Re: label]
|
Forum Member
Registered: 09/25/03
Posts: 4755
Loc: ~ 59N, 15E
|
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
|
Top
|
|
|
|
|
|