dcdfort_reader Module Reference

Module that contains dcdreader class. More...

Data Types
type	dcdfile
	dcdwriter class More...

Functions/Subroutines
subroutine	dcdfile_open (this, filename)
	Opens file to read from. More...
subroutine	dcdfile_read_header (this, nframes, istart, nevery, iend, timestep, natoms)
	Reads header of open DCD file. More...
subroutine	dcdfile_close (this)
	Closes DCD file. More...
subroutine	dcdfile_read_next (this, xyz, box)
	Reads next frame into memory. More...
subroutine	dcdfile_skip_next (this, n)
	Skips reading this frame into memory. More...

Detailed Description

Module that contains dcdreader class.

Function/Subroutine Documentation

dcdfile_close()

subroutine dcdfort_reader::dcdfile_close

(

class(dcdfile), intent(inout)

this

)

Closes DCD file.

Parameters

[in,out]

this

dcdreader object

        implicit none
        class(dcdfile), intent(inout) :: this

        close(this%u)

dcdfile_open()

subroutine dcdfort_reader::dcdfile_open	(	class(dcdfile), intent(inout)	this,
		character (len=*), intent(in)	filename
	)

Opens file to read from.

Also detects if the DCD file is compatible with this library and swaps endinness if necessary.

Parameters

[in,out]	this	dcdreader object
[in]	filename	name of of DCD file to read from

        implicit none
        character (len=*), parameter :: magic_string = "CORD"
        integer, parameter :: magic_number = 84
        character (len=*), intent(in) :: filename
        class(dcdfile), intent(inout) :: this
        integer :: line1, charmm_version, has_extra_block, four_dimensions
        character (len=4) :: line2
        logical :: ex

        ! Does file exist?
        inquire(file=trim(filename), exist=ex, size=this%filesize)
        if (ex .eqv. .false.) then
            call error_stop_program(trim(filename)//" does not exist.")
        end if

        ! Open file in native endinness
        open(newunit=this%u, file=trim(filename), form="unformatted", access="stream")

        ! Read in magic number of magic string
        read(this%u,pos=1) line1
        read(this%u) line2

        ! Ensure the magic number and string are correct, if not we'll swap the endinness
        if (line1 .ne. magic_number .or. line2 .ne. magic_string) then

            ! Try converting to the reverse endianness
            close(this%u)
            open(newunit=this%u, file=trim(filename), form="unformatted", access="stream", convert="swap")

            read(this%u,pos=1) line2
            read(this%u) line2

            ! We tried both native and reverse endiness and didn't have magic number or string
            if (line1 .ne. magic_number .or. line2 .ne. magic_string) then
                call error_stop_program("This DCD file format is not supported, or the file header is corrupt.")
            end if

        end if

        ! Check if the file identifies as CHARMM (LAMMPS pretends to be CHARMM v. 24)
        read(this%u, pos=85) charmm_version
        if (charmm_version .eq. 0) then
            call error_stop_program("DCD file indicates it is not CHARMM. Only CHARMM-style DCD files are supported.")
        end if

        ! We only support files with the extra unitcell block
        read(this%u, pos=49) has_extra_block
        if (has_extra_block .ne. 1) then
            call error_stop_program("DCD file indicates it does not have unit cell information. Only DCD files with&
                & unit cell information are supported.")
        end if

        ! We don't support files with four dimensions
        read(this%u) four_dimensions
        if (four_dimensions .eq. 1) then
            call error_stop_program("DCD file indicates it has four dimensions. Only DCD files with three dimensions&
                & are supported.")
        end if

dcdfile_read_header()

subroutine dcdfort_reader::dcdfile_read_header	(	class(dcdfile), intent(inout)	this,
		integer, intent(out)	nframes,
		integer, intent(out)	istart,
		integer, intent(out)	nevery,
		integer, intent(out)	iend,
		real, intent(out)	timestep,
		integer, intent(out)	natoms
	)

Reads header of open DCD file.

Should be called after open()

Parameters

[in,out]	this	dcdreader object
[in]	nframes	rnumber of frames (snapshots) in file
[out]	istart	first timestep of trajectory file
[out]	nevery	how often (in timesteps) file was written to
[out]	iend	last timestep of trajectory file
[out]	timestep	timestep of simulation
[out]	natoms	number of atoms in each snapshot

        implicit none

        integer, intent(out) :: nframes, istart, nevery, iend, natoms
        integer :: i, ntitle, n, framesize, nframes2, dummy
        character (len=80) :: title_string
        real, intent(out) :: timestep
        class(dcdfile), intent(inout) :: this

        read(this%u, pos=9) nframes, istart, nevery, iend

        read(this%u, pos=45) timestep

        read(this%u, pos=97) ntitle
        if (ntitle > 0) then
            write(error_unit,'(a)') prompt//"The following titles were found:"
        end if
        do i = 1, ntitle
            read(this%u) title_string

            n = 1
            do while (n .le. 80 .and. title_string(n:n) .ne. c_null_char)
                n = n + 1
            end do

            write(error_unit,'(a)') prompt//"  "//trim(title_string(1:n))
        end do

        read(this%u) dummy, dummy
        if (dummy .ne. 4) then
            call error_stop_program("This DCD file format is not supported, or the file header is corrupt.")
        end if

        ! Number of atoms in each snapshot
        read(this%u) natoms, dummy
        if (dummy .ne. 4) then
            call error_stop_program("This DCD file format is not supported, or the file header is corrupt.")
        end if

        ! Each frame has natoms*3 (4 bytes each) = natoms*12
        ! plus 6 box dimensions (8 bytes each) = 48
        ! Additionally there are 32 bytes of file information in each frame
        this%framesize = natoms*12 + 80
        ! Header is 276 bytes
        nframes2 = (this%filesize-276)/this%framesize
        if ( nframes2 .ne. nframes) then
            write(error_unit,'(a,i0,a,i0,a)') prompt//"WARNING: Header indicates ", nframes, &
                &" frames, but file size indicates ", nframes2, "." 
            nframes = nframes2
        end if

dcdfile_read_next()

subroutine dcdfort_reader::dcdfile_read_next	(	class(dcdfile), intent(inout)	this,
		real, dimension(:,:), intent(inout), allocatable	xyz,
		real(8), dimension(6), intent(inout)	box
	)

Reads next frame into memory.

Parameters

[in,out]	this	dcdreader object
[in,out]	xyz	coordinates of all atoms in snapshot. First dimension is atom number, second dimension is x, y, or z !coordinate.
[in,out]	box	dimensions of snapshot. Array containing 6 elements, ordered as A, B, C, alpha, beta, gamma. A = length of unit cell vector along x-axis; B = length of unit cell vector in xy-plane; C = length of unit cell vector in yz-plane; alpha = cosine of angle between B and C; beta = cosine of angle between A and C; gamma = cosine of angle between A and B;

        implicit none
        real, allocatable, intent(inout) :: xyz(:,:)
        real(8), intent(inout) :: box(6)
        integer :: dummy(2)
        class(dcdfile), intent(inout) :: this
    
        ! Should be 48
        read(this%u) dummy(1)

        !            A       gamma   B       beta    alpha   C
        read(this%u) box(1), box(6), box(2), box(5), box(4), box(3)

        ! 48, then no. of bytes for x coordinates, x coordinates (repeat for y and z coordinates)
        read(this%u) dummy, xyz(:,1), dummy, xyz(:,2), dummy, xyz(:,3)

        read(this%u) dummy(1)

dcdfile_skip_next()

subroutine dcdfort_reader::dcdfile_skip_next	(	class(dcdfile), intent(inout)	this,
		integer, intent(in), optional	n
	)

Skips reading this frame into memory.

Parameters

[in,out]	this	dcdreader object
[in]	n	number of frames to skip

        implicit none
        real :: real_dummy
        integer :: dummy, i
        integer(8) :: pos, newpos
        integer, intent(in), optional :: n
        real(8) :: box_dummy(6)
        class(dcdfile), intent(inout) :: this
   
        ! Where are we?
        inquire(unit=this%u, pos=pos)

        ! We subtract 4 bytes so that the next read of the 4-byte integer will line things up properly for the next read
        if (.not. present(n)) then
            newpos = pos + this%framesize - 4
        else
            newpos = pos + this%framesize*n - 4
        end if
        
        read(this%u, pos=newpos) dummy