NAMD Wiki: ReadingDCDinFortran
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:
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
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