dcdfort_trajectory Module Reference

Module that contains Trajectory type. More...

Data Types
type	trajectory
	Trajectory class. More...

Functions/Subroutines
subroutine	trajectory_open (this, filename, ndxfile)
	Trajectory class method which opens DCD file and optionally index file. More...
pure integer(kind=int32) function	trajectory_get_natoms (this, group)
	Trajectory class method which gets the number of atoms in the system or an index group. More...
integer(kind=int32) function	trajectory_skip_next (this, F)
	Trajectory class method which skips a specified number of frames. More...
integer(kind=int32) function	trajectory_read_next (this, F, ndxgrp, every)
	Trajectory class method which reads a specified number of frames into memory after using the open() method. More...
pure real(kind=real32) function, dimension(3)	trajectory_get_xyz (this, frame, atom, group)
	Trajectory class method which returns the coordinates of a particle. More...
subroutine	trajectory_read (this, dcdfile, ndxfile, ndxgrp, every, skip, last)
	Trajectory class method which opens, reads, and closes a trajectory file. More...
pure real(kind=real64) function, dimension(6)	trajectory_get_box (this, frame)
	Trajectory class method which returns the box from a simulation frame. More...
subroutine	trajectory_close (this)
	Trajectory class method which closes a DCD file which was opened with open() More...
pure real(kind=real64) function	trajectory_vol (this, frame)
	Trajectory class method which returns the volume of the box. More...

Detailed Description

Module that contains Trajectory type.

Function/Subroutine Documentation

trajectory_close()

subroutine dcdfort_trajectory::trajectory_close ( class(trajectory), intent(inout)  this )

private

Trajectory class method which closes a DCD file which was opened with open()

Parameters

[in,out]

this

Trajectory class

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

        call this%dcd%close()

trajectory_get_box()

pure real(kind=real64) function, dimension(6) dcdfort_trajectory::trajectory_get_box ( class(trajectory), intent(in)  this, integer(kind=int32), intent(in)  frame  )

private

Trajectory class method which returns the box from a simulation frame.

Parameters

[in,out]	this	Trajectory class
[in]	frame	Which frame to return the box dimensions

Returns

box dimensions of specified 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(kind=real64) :: trajectory_get_box(6)
        class(Trajectory), intent(in) :: this
        integer(kind=int32), intent(in) :: frame
        integer(kind=int32) :: i

        call trajectory_check_frame(this, frame)
        do i = 1, 6
            trajectory_get_box(i) = this%frameArray(frame)%box(i)
        end do

trajectory_get_natoms()

pure integer(kind=int32) function dcdfort_trajectory::trajectory_get_natoms ( class(trajectory), intent(in)  this, character(len=*), intent(in), optional  group  )

private

Trajectory class method which gets the number of atoms in the system or an index group.

Parameters

[in,out]	this	the Trajectory object
[in]	group	name of index group

Returns

number of atoms in index group (if specified) or in system

        implicit none
        integer(kind=int32) :: trajectory_get_natoms
        class(Trajectory), intent(in) :: this
        character(len=*), intent(in), optional :: group

        if (this%read_only_index_group .and. present(group)) then
            call error_stop_program("Do not specify an index group in natoms() when already specifying an & 
                &index group with read() or read_next().")
        end if

        trajectory_get_natoms = merge(this%ndx%get_natoms(group), this%NUMATOMS, present(group))

trajectory_get_xyz()

pure real(kind=real32) function, dimension(3) dcdfort_trajectory::trajectory_get_xyz ( class(trajectory), intent(in)  this, integer(kind=int32), intent(in)  frame, integer(kind=int32), intent(in)  atom, character(len=*), intent(in), optional  group  )

private

Trajectory class method which returns the coordinates of a particle.

Gets the coordinates of a particle from the read in trajectory. Returns a real array with 3 elements (x, y, and z). An optional index group can be specified; if so, the atomid is in relationship to the group.

Parameters

[in,out]	this	Trajectory class
[in]	frame
[in]	atom	atomid of particle to get
[in]	group	optional index group

Returns

coordinate of the particle

        implicit none
        real(kind=real32) :: trajectory_get_xyz(3)
        integer(kind=int32), intent(in) :: frame, atom
        integer(kind=int32) :: atom_tmp, natoms
        class(Trajectory), intent(in) :: this
        character(len=1024) :: message
        character(len=*), intent(in), optional :: group

        call trajectory_check_frame(this, frame)

        if (this%read_only_index_group .and. present(group)) then
            call error_stop_program("Do not specify an index group in x() when already specifying an & 
                &index group with read() or read_next().")
        end if

        atom_tmp = merge(this%ndx%get(group, atom), atom, present(group))
        natoms = merge(this%natoms(group), this%natoms(), present(group))

        if (atom > natoms .or. atom < 1) then
            write(message, "(a,i0,a,i0,a)") "Tried to access atom number ", atom_tmp, " when there are ", &
                natoms, ". Note that Fortran uses one-based indexing."
            call error_stop_program(trim(message))
        end if

        ! Faster to list individually
        trajectory_get_xyz(1) = this%frameArray(frame)%xyz(atom_tmp,1)
        trajectory_get_xyz(2) = this%frameArray(frame)%xyz(atom_tmp,2)
        trajectory_get_xyz(3) = this%frameArray(frame)%xyz(atom_tmp,3)

trajectory_open()

subroutine dcdfort_trajectory::trajectory_open	(	class(trajectory), intent(inout)	this,
		character(len=*), intent(in)	filename,
		character(len=*), intent(in), optional	ndxfile
	)

Trajectory class method which opens DCD file and optionally index file.

Parameters

[in,out]	this	the Trajectory object
[in]	filename	name of DCD file
[in]	ndxfile	name of GROMACS style index file

        implicit none
        class(Trajectory), intent(inout) :: this
        character(len=*), intent(in) :: filename
        character(len=*), intent(in), optional :: ndxfile

        call this%dcd%open(trim(filename))

        write(error_unit,'(a)') prompt//"Opened "//trim(filename)//" for reading."
        call this%dcd%read_header(this%nframes, this%istart, this%nevery, this%iend, this%timestep, this%NUMATOMS)

        write(error_unit,'(a,i0,a)') prompt, this%NUMATOMS, " atoms present in system."
        write(error_unit,'(a,i0,a)') prompt, this%NFRAMES, " frames present in trajectory file."
        write(error_unit,'(a,i0)') prompt//"First timestep in trajectory file: ", this%ISTART
        write(error_unit,'(a,i0)') prompt//"Last timestep in strajectory file: ", this%IEND
        write(error_unit,'(a,f12.6)') prompt//"Simulation timestep: ", this%timestep
        write(error_unit,'(a,i0,a)') prompt//"Trajectory written every ", this%nevery, " timesteps."

        this%FRAMES_REMAINING = this%NFRAMES

        this%N = this%NUMATOMS ! Save for use when user selects just one group
        if (present(ndxfile)) call this%ndx%read(ndxfile, this%NUMATOMS)

trajectory_read()

subroutine dcdfort_trajectory::trajectory_read ( class(trajectory), intent(inout)  this, character(len=*)  dcdfile, character(len=*), optional  ndxfile, character(len=*), optional  ndxgrp, integer(kind=int32), intent(in), optional  every, integer(kind=int32), intent(in), optional  skip, integer(kind=int32), intent(in), optional  last  )

private

Trajectory class method which opens, reads, and closes a trajectory file.

Parameters

[in,out]	this	Trajectory class
[in]	dcdfile	Name of DCD trajectory file
[in]	ndxfile	Name of optional index file
[in]	ndxgrp	Name of optional group. If specified, only that group will be read into memory; otherwise, all particles read into memory.
[in]	every	Read in every this many frames; default is 1
[in]	skip	Skip this many frames at the beginning of the trajectory file; default is 0
[in]	last	Last frame to read in

        implicit none
        class(Trajectory), intent(inout) :: this
        character(len=*), optional :: ndxfile, ndxgrp
        character(len=*) :: dcdfile
        integer(kind=int32) :: nread
        integer(kind=int32), intent(in), optional :: every, skip, last

        call this%open(dcdfile, ndxfile)

        if (present(every)) then
            write(error_unit,'(a,i0,a)') prompt//"Saving every ", every, " snapshots into memory."
        end if
        if (present(skip)) then
            write(error_unit,'(a,i0,a)') prompt//"Skipping first ", skip, " snapshots."
        end if
        if (present(last)) then
            write(error_unit,'(a,i0,a)') prompt//"Limiting to ", last, " snapshots."
        end if

        if (present(skip)) then
            nread = this%skip_next(skip)
            if (present(last)) then
                nread = this%read_next(last, ndxgrp, every)
            else
                nread = this%read_next(this%NFRAMES-nread, ndxgrp, every)
            end if
        else
            if (present(last)) then
                nread = this%read_next(last, ndxgrp, every)
            else
                nread = this%read_next(this%NFRAMES, ndxgrp, every)
            end if
        end if

        call this%close()

trajectory_read_next()

integer(kind=int32) function dcdfort_trajectory::trajectory_read_next ( class(trajectory), intent(inout)  this, integer(kind=int32), intent(in), optional  F, character(len=*), optional  ndxgrp, integer(kind=int32), intent(in), optional  every  )

private

Trajectory class method which reads a specified number of frames into memory after using the open() method.

Parameters

[in,out]	this	Trajectory class
[in]	F	frame number to read up to; if not specified, 1 frame is read
[in]	ndxgrp	read only this index group into memory
[in]	every	Read in every this many frames; default is to read in every frame

Returns

number of frames read in

        implicit none
        integer(kind=int32) :: trajectory_read_next
        class(Trajectory), intent(inout) :: this
        integer(kind=int32), intent(in), optional :: F, every
        character(len=*), optional :: ndxgrp
        integer(kind=int32) :: I, J, N, S, K
        real(kind=real32), allocatable :: xyz(:,:)

        ! If the user specified how many frames to read and it is greater than one, use it
        n = merge(f, 1, present(f))
        s = merge(every, 1, present(every)) - 1

        ! Are we near the end of the file?
        n = min(this%FRAMES_REMAINING, n)
        this%FRAMES_REMAINING = this%FRAMES_REMAINING - n

        if (present(every)) n = n / every

        if (allocated(this%frameArray)) deallocate(this%frameArray)
        allocate(this%frameArray(n))

        write(error_unit,*)

        this%read_only_index_group = .false.

        if (present(ndxgrp)) then

            allocate(xyz(this%N,3))
            this%NUMATOMS = this%natoms(trim(ndxgrp))
            do i = 1, n

                if (modulo(i*(s+1), 1000) .eq. 0) call print_frames_saved(i, ndxgrp)

                allocate(this%frameArray(i)%xyz(this%NUMATOMS,3))
                call this%dcd%skip_next(s)
                call this%dcd%read_next(xyz, this%frameArray(i)%box)

                do j = 1, size(this%ndx%group)
                    if (trim(this%ndx%group(j)%title) .eq. trim(ndxgrp)) then
                        ! Tends to be faster to list this out individually
                        do k = 1, this%NUMATOMS
                            this%frameArray(i)%xyz(k,1) = xyz(this%ndx%group(j)%LOC(k),1)
                            this%frameArray(i)%xyz(k,2) = xyz(this%ndx%group(j)%LOC(k),2)
                            this%frameArray(i)%xyz(k,3) = xyz(this%ndx%group(j)%LOC(k),3)
                        end do
                        exit
                    end if
                end do

            end do
            deallocate(xyz)

            this%read_only_index_group = .true.

        else

            do i = 1, n

                if (modulo(i*(s+1), 1000) .eq. 0) call print_frames_saved(i)

                allocate(this%frameArray(i)%xyz(this%NUMATOMS,3))
                call this%dcd%skip_next(s)
                call this%dcd%read_next(this%frameArray(i)%xyz, this%frameArray(i)%box)

            end do

        end if

        this%frames_read = n
        trajectory_read_next = n
        call print_frames_saved(n, ndxgrp)

trajectory_skip_next()

integer(kind=int32) function dcdfort_trajectory::trajectory_skip_next ( class(trajectory), intent(inout)  this, integer(kind=int32), intent(in), optional  F  )

private

Trajectory class method which skips a specified number of frames.

Parameters

[in,out]	this	Trajectory class
[in]	F	number of frames to skip; if not specified, 1 frame is skipped

Returns

number of frames skipped

        implicit none
        class(Trajectory), intent(inout) :: this
        integer(kind=int32), intent(in), optional :: F
        integer(kind=int32) :: trajectory_skip_next, N

        ! If the user specified how many frames to read and it is greater than one, use it
        n = merge(f, 1, present(f))

        ! Are we near the end of the file?
        n = min(this%FRAMES_REMAINING, n)
        this%FRAMES_REMAINING = this%FRAMES_REMAINING - n

        call this%dcd%skip_next(n)
        write(error_unit,'(a,i0,a)') prompt//"Skipped ", n, " frames."
        trajectory_skip_next = n

trajectory_vol()

pure real(kind=real64) function dcdfort_trajectory::trajectory_vol ( class(trajectory), intent(in)  this, integer(kind=int32), intent(in)  frame  )

private

Trajectory class method which returns the volume of the box.

Parameters

[in,out]	this	Trajectory class
[in]	frame	snapshot to get box volume of

Returns

the volume of the box of the frame specified

        implicit none
        class(Trajectory), intent(in) :: this
        integer(kind=int32), intent(in) :: frame
        real(kind=real64) :: trajectory_vol

        trajectory_vol = vol(this%box(frame))