dcdfort
Modern Fortran library for analyzing DCD trajectories
Data Types | Functions/Subroutines
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)
 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]thisdcdreader object
169 
170  implicit none
171  class(dcdfile), intent(inout) :: this
172 
173  close(this%u)
174 

◆ dcdfile_open()

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

Opens file to read from.

Parameters
[in,out]thisdcdreader object
[in]filenamename of of DCD file to read from
58 
59  implicit none
60  character (len=*), intent(in) :: filename
61  class(dcdfile), intent(inout) :: this
62  integer :: filesize
63  logical :: ex
64 
65  inquire(file=trim(filename), exist=ex, size=this%filesize)
66 
67  if (ex .eqv. .false.) then
68  call error_stop_program(trim(filename)//" does not exist.")
69  end if
70 
71  open(newunit=this%u, file=trim(filename), form="unformatted", access="stream")
72 

◆ 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]thisdcdreader object
[in]nframesrnumber of frames (snapshots) in file
[out]istartfirst timestep of trajectory file
[out]neveryhow often (in timesteps) file was written to
[out]iendlast timestep of trajectory file
[out]timesteptimestep of simulation
[out]natomsnumber of atoms in each snapshot
85 
86  implicit none
87 
88  integer :: dummy
89  integer, intent(out) :: nframes, istart, nevery, iend, natoms
90  integer :: i, ntitle, n, framesize, nframes2, curr_pos
91  character (len=4) :: cord_string
92  character (len=80) :: title_string
93  real, intent(out) :: timestep
94  class(dcdfile), intent(inout) :: this
95 
96  read(this%u) dummy
97  if (dummy .ne. 84) then
98  call error_stop_program("This dcd file format is not supported, or the file header is corrupt.")
99  end if
100 
101  read(this%u) cord_string
102  if (cord_string .ne. "CORD") then
103  call error_stop_program("This dcd file format is not supported, or the file header is corrupt.")
104  end if
105 
106  ! Number of snapshots in file
107  read(this%u) nframes
108 
109  ! Timestep of first snapshot
110  read(this%u) istart
111 
112  ! Save snapshots every this many steps
113  read(this%u) nevery
114 
115  ! Timestep of last snapshot
116  read(this%u) iend
117 
118  do i = 1, 5
119  read(this%u) dummy
120  end do
121 
122  ! Simulation timestep
123  read(this%u) timestep
124 
125  do i = 1, 12
126  read(this%u) dummy
127  end do
128 
129  read(this%u) ntitle
130  do i = 1, ntitle
131  read(this%u) title_string
132 
133  n = 1
134  do while (n .le. 80 .and. title_string(n:n) .ne. c_null_char)
135  n = n + 1
136  end do
137 
138  write(error_unit,'(a)') trim(title_string(1:n))
139  end do
140 
141  read(this%u) dummy
142  read(this%u) dummy
143  if (dummy .ne. 4) then
144  call error_stop_program("This dcd file format is not supported, or the file header is corrupt.")
145  end if
146 
147  ! Number of atoms in each snapshot
148  read(this%u) natoms
149 
150  read(this%u) dummy
151 
152  inquire(unit=this%u, pos=curr_pos)
153 
154  ! Each frame has natoms*3 (4 bytes each), plus 6 box dimensions (8 bytes each)
155  framesize = natoms*12 + 48
156  ! Just an estimate
157  nframes2 = (this%filesize-curr_pos)/framesize
158  if ( nframes2 .ne. nframes) then
159  write(error_unit,'(a,i0,a,i0,a)') "WARNING: Header indicates ", nframes, &
160  &" frames, but file size indicates ", nframes2, "."
161  nframes = nframes2
162  end if
163 

◆ 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]thisdcdreader object
[in,out]xyzcoordinates of all atoms in snapshot
[in,out]boxdimensions of snapshot
182 
183  implicit none
184  real, allocatable, intent(inout) :: xyz(:,:)
185  real(8), intent(inout) :: box(6)
186  integer :: dummy
187  class(dcdfile), intent(inout) :: this
188 
189  ! Should be 48
190  read(this%u) dummy
191 
192  read(this%u) box(1) ! A
193  read(this%u) box(6) ! gamma
194  read(this%u) box(2) ! B
195  read(this%u) box(5) ! beta
196  read(this%u) box(4) ! alpha
197  read(this%u) box(3) ! C
198  if (box(4) >= -1.0 .and. box(4) <= 1.0 .and. box(5) >= -1.0 .and. box(5) <= 1.0 .and. &
199  box(6) >= -1.0 .and. box(6) <= 1.0) then
200  box(4) = 90.0 - asin(box(4)) * 180.0 / pi
201  box(5) = 90.0 - asin(box(5)) * 180.0 / pi
202  box(6) = 90.0 - asin(box(6)) * 180.0 / pi
203  end if
204 
205  ! 48 again
206  read(this%u) dummy
207  read(this%u) dummy
208 
209  read(this%u) xyz(1,:)
210 
211  ! 48 again
212  read(this%u) dummy
213  read(this%u) dummy
214 
215  read(this%u) xyz(2,:)
216 
217  ! 48 again
218  read(this%u) dummy
219  read(this%u) dummy
220 
221  read(this%u) xyz(3,:)
222 
223  read(this%u) dummy
224 

◆ dcdfile_skip_next()

subroutine dcdfort_reader::dcdfile_skip_next ( class(dcdfile), intent(inout)  this)

Skips reading this frame into memory.

Parameters
[in,out]thisdcdreader object
230 
231  implicit none
232  real :: real_dummy
233  integer :: dummy, i
234  real(8) :: box_dummy(6)
235  class(dcdfile), intent(inout) :: this
236 
237  ! Should be 48
238  read(this%u) dummy
239 
240  read(this%u) box_dummy
241 
242  ! 48 again
243  read(this%u) dummy
244  read(this%u) dummy
245 
246  do i = 1, dummy/4
247  read(this%u) real_dummy
248  end do
249 
250  ! 48 again
251  read(this%u) dummy
252  read(this%u) dummy
253 
254  do i = 1, dummy/4
255  read(this%u) real_dummy
256  end do
257 
258  ! 48 again
259  read(this%u) dummy
260  read(this%u) dummy
261 
262  do i = 1, dummy/4
263  read(this%u) real_dummy
264  end do
265 
266  read(this%u) dummy
267 
Generated on Fri Oct 13 2017 20:42:11 for dcdfort by   doxygen 1.8.13