dcdfort_index Module Reference

Module that contains IndexFile type. More...

Data Types
type	indexfile
	IndexFile class. More...
type	ndxgroups
	ndxgroups class More...

Functions/Subroutines
subroutine	indexfile_read (this, filename, N)
	Open, read in, and process the GROMACS-style index file, storing information in memory. More...
integer function	indexfile_get (this, group_name, I)
	Gets the number of atoms in a group. More...

Detailed Description

Module that contains IndexFile type.

This module contains the IndexFile class which can be used to read in and process GROMACS-style index files. It should not be necessary to create an object from this class independently from a Trajectory object, since it is included in the Trajectory class.

Function/Subroutine Documentation

indexfile_get()

integer function dcdfort_index::indexfile_get ( class(indexfile), intent(inout)  this, character (len=*), intent(in)  group_name, integer, intent(in), optional  I  )

private

Gets the number of atoms in a group.

Parameters

[in,out]	this	IndexFile class
[in]	group_name	index group title or name (in brackets in the index file)
[in]	I	the location in the group

Returns

If an atom is specified, integer returns the overall index for that atom; otherwise, returns number of atoms in group

        implicit none
        integer :: indexfile_get
        class(indexfile), intent(inout) :: this
        character (len=*), intent(in) :: group_name
        integer, intent(in), optional :: i
        integer :: j

        do j = 1, size(this%group)

            if (trim(this%group(j)%title) .eq. trim(group_name)) then

                indexfile_get = merge(this%group(j)%LOC(i), this%group(j)%NUMATOMS, present(i))
                return

            end if

        end do

        ! If user asked to get the number of atoms in an index group, and that index group is not 
        ! in the index file, just return 0. If the user specified an atom number, then throw an error,
        ! since the overall index cannot be returned in that case
        if (.not. present(i)) then
            if (this%group_warning) then
                write(error_unit, '(a)') prompt//"WARNING: No atoms found in index group '"//trim(group_name)//"'."
                write(error_unit, '(a)') prompt//"This warning will not appear again for any other index groups."
                this%group_warning = .false.
            end if
            indexfile_get = 0
        else
            indexfile_get = -1
            write(error_unit, '(a)') prompt//"ERROR: "//trim(group_name)//" is not in index file. The groups available are:"
            do j = 1, size(this%group)
                write(error_unit,'(a10,a,i0,a)') this%group(j)%title, " (", this%group(j)%NUMATOMS, ")"
            end do
            stop 1
        end if

indexfile_read()

subroutine dcdfort_index::indexfile_read	(	class(indexfile), intent(inout)	this,
		character (len=*), intent(in)	filename,
		integer, intent(in), optional	N
	)

Open, read in, and process the GROMACS-style index file, storing information in memory.

Parameters

[in,out]	this	IndexFile class
[in]	filename	name of GROMACS-style index file
[in]	N	number of atoms in the entire system; used as a check

        implicit none
        class(indexfile), intent(inout) :: this
        character (len=*), intent(in) :: filename
        character (len=2048) :: line, ncols_string, fmt_string
        integer :: index_file_unit, io_status, ngrps, i, j, ncols
        integer, allocatable :: indices_tmp(:), title_loc(:), num_array(:)
        logical :: ex
        integer, intent(in), optional :: n

        write(error_unit,'(a)') prompt//"Reading in "//trim(filename)//"."

        ! Does the file exist?
        inquire(file=trim(filename), exist=ex)
        if (ex .eqv. .false.) call error_stop_program(trim(filename)//" does not exist.")

        ! Is in index file?
        open(newunit=index_file_unit, file=trim(filename), status="old")
        read(index_file_unit, '(a)', iostat=io_status) line
        if (index(line, "[") .eq. 0) call error_stop_program(trim(filename)//" is not a valid index file.")

        ! How many groups are in it?
        rewind index_file_unit
        io_status = 0
        ngrps = 0
        do while (io_status .eq. 0)
            read(index_file_unit, '(a)', iostat=io_status) line
            if (io_status .ne. 0) goto 100
            if (index(line, "[") .ne. 0) ngrps = ngrps + 1
        end do

100     continue

        if (allocated(this%group)) deallocate(this%group)
        allocate(this%group(ngrps), title_loc(ngrps+1)) ! Add one to include end of file

        ! Now find the title locations and save their names
        rewind index_file_unit
        i = 1
        j = 1
        io_status = 0
        do while (io_status .eq. 0)

            read(index_file_unit, '(a)', iostat=io_status) line
            if (io_status .ne. 0) goto 200
            if (index(line, "[") .ne. 0) then
                this%group(i)%title = trim(line(index(line, "[")+2:index(line, "]")-2))
                title_loc(i) = j
                i = i + 1
            end if
            j = j + 1

        end do

200     continue

        if (this%group(1)%title .ne. "System") call error_stop_program("Index file does not have 'System' group as first group.")

        ! Index files can have a varying number of columns. This attempts to
        ! detect the correct number by reading in the second line of the file,
        ! which should be a list of indices for the "System" group.
        ncols = 100
        allocate(num_array(ncols))
        io_status = 5000
        do while (io_status .ne. 0)
            ncols = ncols - 1
            write(ncols_string, '(i0)') ncols
            write(fmt_string, '(a)') '('//trim(ncols_string)//'i0)'
            rewind index_file_unit
            read(index_file_unit, '(a)', iostat=io_status) line
            read(index_file_unit, '(a)', iostat=io_status) line
            read(line, *, iostat=io_status) num_array(1:ncols)
        end do

        title_loc(i) = j ! End of file location

        ! Now finally get all of the indices for each group
        ! Allocate for total number of atoms in system, since that is the maximum
        allocate(indices_tmp(n))
        do i = 1, ngrps

            ! Initial guess only how many items are in the group
            ! Add 1, bc loop subtracts 1 at the beginning
            this%group(i)%NUMATOMS = (title_loc(i+1)-title_loc(i)-1)*ncols + 1

            if (n < this%group(i)%NUMATOMS) this%group(i)%NUMATOMS = n + 1
            io_status = 5000

            do while (io_status .ne. 0)

                ! Our guess was too large if we made it back here, go to the beginning and reduce our guess by 1, try again
                rewind index_file_unit
                this%group(i)%NUMATOMS = this%group(i)%NUMATOMS - 1
                if (this%group(i)%NUMATOMS .le. 0) then 
                    this%group(i)%NUMATOMS = 0
                    goto 300
                end if

                ! Read all the way to the group
                do j = 1, title_loc(i); read(index_file_unit, '(a)', iostat=io_status) line; end do

                ! Attempt to read into array
                read(index_file_unit, *, iostat=io_status) indices_tmp(1:this%group(i)%NUMATOMS)

            end do

            ! Specifying array bounds for array to be allocated is not required for F2008 but is required for F2003
            allocate(this%group(i)%LOC(1:this%group(i)%NUMATOMS), source=indices_tmp(1:this%group(i)%NUMATOMS))

300         cycle

        end do
        deallocate(indices_tmp)

        do i = 1, ngrps-1
            do j = i+1, ngrps
                if (i .ne. j) then
                    if (this%group(i)%title .eq. this%group(j)%title) then
                        write(error_unit,*)
                        write(error_unit,'(a, a, a)') prompt//"WARNING: Index group ", this%group(i)%title, &
                            " was specified more than once in index file."
                        write(error_unit,*)
                    end if
                end if
            end do
        end do

        ! Error checking to see if index file goes with the trajectory file

        ! If the number of atoms is not the same in the System group (first group) and dcd file
        if (this%group(1)%numatoms .ne. n .or. this%group(1)%loc(this%group(1)%numatoms) .ne. n) then
            call error_stop_program("Index file does not match dcd file.")
        end if

        if (present(n)) then
            do i = 1, ngrps

                ! If number of atoms in index group is larger than number of atoms in dcd file
                if (this%group(i)%numatoms .gt. n) call error_stop_program("Index file does not match dcd file.")

                ! If a location number is greater than number of atoms in dcd file
                do j = 1, this%group(i)%numatoms
                    if (this%group(i)%loc(j) .gt. n) call error_stop_program("Index file does not match dcd file.")
                end do

            end do
        end if

        close(index_file_unit)

        write(error_unit,'(a,i0,a)') prompt//"Read in ", ngrps, " index groups."