Hi,
On Fri, Oct 4, 2013 at 11:33 PM, Uday Midya <umidya.iitkgp.gmail.com> wrote:
> open(10,file='file.dcd',status='old',form='unformatted')
> read(10)dummyc,nframes,(dummyi,i=1,8),dummyr,(dummyi,i=1,9)
Admittedly I'm not a DCD expert, but it seems to me that this code
shouldn't work, at least for Charmm DCD files. As far as I know the
beginning of the DCD header format is this:
Header Block Size Begin: integer (32/64 bit) which equals 84
4 bytes ('C' 'O' 'R' 'D')
80 bytes (20 * 4 byte integers)
Header Block Size End: integer (32/64 bit) which equals 84
Since you are compiling with 'm32' everything is forced to be a 4 byte
integer, which is fine. However, you're only reading 1 integer before
you read in the # of frames which is wrong; you need to read 2 (the
'84' and 'CORD'). I'm not sure if NAMD does things differently but
this is the kind of DCD that ptraj and cpptraj both write. You will
need to put one more dummy variable before nframes.
> read(10)dummyi,dummyr
> read(10)natom
> do 1000 ii=1,nframes
> read(inunit)bx,bx1,by,by1,bz1,bz
> read(inunit)(x(j),j=1,natom)
> read(inunit)(y(j),j=1,natom)
> read(inunit)(z(j),j=1,natom)
> 1000 continue
Continued from above, the rest of the DCD header/coords format is:
Title Block Size Begin: integer (32/64 bit)
Ntitle: 4 byte integer
Titles: Ntitle * 80 bytes
Title Block Size End: integer (32/64 bit)
Atom Block Size Begin: integer (32/64 bit) which equals 4
Natom: 4 byte integer
Atom Block Size End: integer (32/64 bit) which equals 4
< Stuff regarding free/fixed atoms if there are fixed atoms >
< Box for frame 0 if box coords >
< X coords for frame 0>
< Y coords for frame 0>
< Z coords for frame 0>
...
It looks to me like you're not reading enough variables prior to
getting 'natom'; the titles in particular are completely skipped.
Also, the box coordinates and X/Y/Z coordinate blocks should be
bounded by integers. I'm not sure if NAMD writes a different format
DCD file, but the code you pasted would definitely not work with the
DCD format output from ptraj/cpptraj. Check the code in cpptraj files
Traj_CharmmDcd.cpp writeDcdHeader() and writeFrame() to see how
cpptraj writes out DCD trajectories (it's fairly well-commented). Also
there are issues with endianness, but you won't have to worry about
this as long as the machine you are reading with has the same
endianness as the machine that wrote the trajectory.
> Is there such 10-12 line fortran code to read binary mdcrd file (NetCDF).
More recent versions of Amber have NetCDF reading code in
$AMBERHOME/src/sander/trajene.F90 (look inside the #ifdef BINTRAJ
blocks). You can get some more info regarding the NetCDF format in
Amber from
http://ambermd.org/netcdf/nctraj.html. Also, code for
writing NetCDF trajectories in Fortran are in:
$AMBERHOME/src/pmemd/src/bintraj.F90
$AMBERHOME/src/sander/bintraj.F90
Hope this helps,
-Dan
--
-------------------------
Daniel R. Roe, PhD
Department of Medicinal Chemistry
University of Utah
30 South 2000 East, Room 201
Salt Lake City, UT 84112-5820
http://home.chpc.utah.edu/~cheatham/
(801) 587-9652
(801) 585-6208 (Fax)
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Mon Oct 07 2013 - 13:30:04 PDT
