NAMD Wiki: ReadingDCDinFortran

  You are encouraged to improve NamdWiki by adding content, correcting errors, or removing spam.

From time to time someone asks me how to read DCD files using fortran because of some 2006 mailing list posts. So, here is how I do it. I warn everyone that probably this is not the correct way to read the header, since I could never write a DCD file that VMD could read. But using the approach below you can read the coordinates with no problem and make your own analysis codes:

You must declare the variables, in fortran, like this:

double precision d
real t, x(maxframes,maxatoms), y(maxframes,maxatoms),
     z(maxframes,maxatoms), dummyr
integer nset, natom, dummyi
character*4 dummyc

(the dummy variables are there because I don't know what they mean in the header).

The header of the DCD file is then read with:

open(10,file=dcdfile,status='old',form='unformatted')
read(10) dummyc, nframes, (dummyi,i=1,8), dummyr, (dummyi,i=1,9)
read(10) dummyi, dummyr
read(10) natom

Then the coordinates are read with:

 do i = 1, nframes
   if(DCDUnitCell) read(10) (d, j=1, 6)
   read(10) (x(i,j),j=1,natom)
   read(10) (y(i,j),j=1,natom)
   read(10) (z(i,j),j=1,natom)
end do

Note that the first line of the four reading lines must be read something if the DCDUnitCell parameter was set to "on", which is the default value since namd 2.6, i think. Otherwise that line can be removed. That line contains the periodic cell size (read d as a vector to get the values).

If using fortran77 and 64-bit machines the program must be compiled with the -m32 option. Fortran90+ does not require that because (I think) the precision of the variables are not machine-dependent.

Some examples of fortran77 and fortran90 simple analysis programs that use this are available at http://lm-mdanalysis.googlecode.com

----------------------

Here is a complete f90 code that is based on the code above, which needed some tweaks to get running:

program BinRead

implicit none
character*80      :: filename
double precision  :: d
real              :: t,dummyr
real,allocatable  :: x(:,:),y(:,:),z(:,:)
integer           :: nset,natom,dummyi,i,j,nframes
character*4       :: dummyc

! getarg retrieves command line arguments, and sets first to "filename"
call getarg(1,filename)

open(10,file=trim(filename),status='old',form='unformatted')
read(10) dummyc, nframes, (dummyi,i=1,8), dummyr, (dummyi,i=1,9)
read(10) dummyi, dummyr
read(10) natom

allocate(x(nframes,natom))
allocate(y(nframes,natom))
allocate(z(nframes,natom))


do i = 1, nframes
   read(10) (d, j=1, 6)
   read(10) (x(i,j),j=1,natom)
   read(10) (y(i,j),j=1,natom)
   read(10) (z(i,j),j=1,natom)
end do

15 format(a17,i6,a12,i6,a6)

16 format(6f13.8)

print*
print 15, 'The dcd contains ',nframes,' frames and ',natom,' atoms'
print*


! print 16,x(1,21362),y(1,21362),z(1,21362),x(125,1),y(125,1),z(125,1)

deallocate(x)
deallocate(y)
deallocate(z)


close(10)

end program BinRead