RMPCDMD
particle_system_io.f90
Go to the documentation of this file.
1 ! This file is part of RMPCDMD
2 ! Copyright (c) 2016-2017 Pierre de Buyl and contributors
3 ! License: BSD 3-clause (see file LICENSE)
4 
12 
14  use hdf5
15  use h5md_module
16  use particle_system
17  implicit none
18 
19  private
20 
21  public :: particle_system_io_t
23  public :: thermo_t
24 
25  type thermo_t
26  integer :: n_buffer
27  integer :: idx
28  double precision, allocatable, dimension(:) :: temperature
29  double precision, allocatable, dimension(:) :: potential_energy
30  double precision, allocatable, dimension(:) :: kinetic_energy
31  double precision, allocatable, dimension(:) :: internal_energy
32  double precision, allocatable, dimension(:,:) :: center_of_mass_velocity
33  contains
34  procedure :: init => thermo_init
35  procedure :: append => thermo_append
36  end type thermo_t
37 
39  integer :: mode
40  integer :: step
41  integer :: step_offset = -1
42  double precision :: time
43  double precision :: time_offset = -1
44  logical :: store
46 
48  integer :: nmax
49  integer :: error
50  integer(HID_T) :: group
51  type(h5md_element_t) :: box
52  type(h5md_element_t) :: position
53  type(h5md_element_t) :: image
54  type(h5md_element_t) :: velocity
55  type(h5md_element_t) :: force
56  type(h5md_element_t) :: id
57  type(h5md_element_t) :: species
58  type(particle_system_io_info_t) :: position_info
59  type(particle_system_io_info_t) :: image_info
60  type(particle_system_io_info_t) :: velocity_info
61  type(particle_system_io_info_t) :: force_info
62  type(particle_system_io_info_t) :: id_info
63  type(particle_system_io_info_t) :: species_info
64  contains
65  procedure :: init => ps_init
66  procedure :: close => ps_close
67  end type particle_system_io_t
68 
69 contains
70 
71  subroutine ps_init(this, h5md_file, name, ps)
72  class(particle_system_io_t), intent(inout) :: this
73  type(h5md_file_t), intent(inout) :: h5md_file
74  character(len=*), intent(in) :: name
75  type(particle_system_t), intent(in) :: ps
76  type(particle_system_io_info_t) :: info
77 
78  this% Nmax = ps% Nmax
79 
80  call h5gcreate_f(h5md_file% particles, name, this% group, this% error)
81 
82  if (this%position_info%store) then
83  info = this%position_info
84  if ( (info%step_offset < 0) .and. (info%time_offset > 0) .or. &
85  (info%step_offset > 0) .and. (info%time_offset < 0) ) then
86  stop 'use both or neither of step_offset and time_offset in ps_init position'
87  end if
88  if (iand(info%mode, h5md_time) == h5md_time) then
89  call this%position%create_time(this% group, 'position', ps% pos, info%mode)
90  else if (iand(info%mode, h5md_linear) == h5md_linear) then
91  if (info%step_offset > 0) then
92  call this%position%create_time(this% group, 'position', ps% pos, info%mode, info%step, &
93  info%time, step_offset=info%step_offset, time_offset=info%time_offset)
94  else
95  call this%position%create_time(this% group, 'position', ps% pos, info%mode, info%step, &
96  info%time)
97  end if
98  else if (iand(info%mode, h5md_fixed) == h5md_fixed) then
99  call this%position%create_fixed(this% group, 'position', ps% pos)
100  else
101  stop 'unknown storage for position in particle_system_io init'
102  end if
103  end if
104 
105  if (this%velocity_info%store) then
106  info = this%velocity_info
107  if ( (info%step_offset < 0) .and. (info%time_offset > 0) .or. &
108  (info%step_offset > 0) .and. (info%time_offset < 0) ) then
109  stop 'use both or neither of step_offset and time_offset in ps_init velocity'
110  end if
111  if (iand(info%mode, h5md_time) == h5md_time) then
112  call this%velocity%create_time(this% group, 'velocity', ps% vel, info%mode)
113  else if (iand(info%mode, h5md_linear) == h5md_linear) then
114  if (info%step_offset > 0) then
115  call this%velocity%create_time(this% group, 'velocity', ps% vel, info%mode, info%step, &
116  info%time, step_offset=info%step_offset, time_offset=info%time_offset)
117  else
118  call this%velocity%create_time(this% group, 'velocity', ps% vel, info%mode, info%step, &
119  info%time)
120  end if
121  else if (iand(info%mode, h5md_fixed) == h5md_fixed) then
122  call this%velocity%create_fixed(this% group, 'velocity', ps% vel)
123  else
124  stop 'unknown storage for velocity in particle_system_io init'
125  end if
126  end if
127 
128  if (this%force_info%store) then
129  info = this%force_info
130  if ( (info%step_offset < 0) .and. (info%time_offset > 0) .or. &
131  (info%step_offset > 0) .and. (info%time_offset < 0) ) then
132  stop 'use both or neither of step_offset and time_offset in ps_init velocity'
133  end if
134 
135  if (iand(info%mode, h5md_time) == h5md_time) then
136  call this%force%create_time(this% group, 'force', ps% vel, info%mode)
137  else if (iand(info%mode, h5md_linear) == h5md_linear) then
138  if (info%step_offset > 0) then
139  call this%force%create_time(this% group, 'force', ps% pos, info%mode, info%step, &
140  info%time, step_offset=info%step_offset, time_offset=info%time_offset)
141  else
142  call this%force%create_time(this% group, 'force', ps% pos, info%mode, info%step, &
143  info%time)
144  end if
145  else if (iand(info%mode, h5md_fixed) == h5md_fixed) then
146  call this%force%create_fixed(this% group, 'force', ps% vel)
147  else
148  stop 'unknown storage for force in particle_system_io init'
149  end if
150  end if
151 
152  if (this%image_info%store) then
153  info = this%image_info
154  if ( (info%step_offset < 0) .and. (info%time_offset > 0) .or. &
155  (info%step_offset > 0) .and. (info%time_offset < 0) ) then
156  stop 'use both or neither of step_offset and time_offset in ps_init image'
157  end if
158 
159  if (iand(info%mode, h5md_time) == h5md_time) then
160  call this%image%create_time(this% group, 'image', ps% image, info%mode)
161  else if (iand(info%mode, h5md_linear) == h5md_linear) then
162  if (info%step_offset > 0) then
163  call this%image%create_time(this% group, 'image', ps% image, info%mode, info%step, &
164  info%time, step_offset=info%step_offset, time_offset=info%time_offset)
165  else
166  call this%image%create_time(this% group, 'image', ps% image, info%mode, info%step, &
167  info%time)
168  end if
169  else if (iand(info%mode, h5md_fixed) == h5md_fixed) then
170  call this%image%create_fixed(this% group, 'image', ps% image)
171  else
172  stop 'unknown storage for image in particle_system_io init'
173  end if
174  end if
175 
176  if (this%id_info%store) then
177  info = this%id_info
178  if ( (info%step_offset < 0) .and. (info%time_offset > 0) .or. &
179  (info%step_offset > 0) .and. (info%time_offset < 0) ) then
180  stop 'use both or neither of step_offset and time_offset in ps_init id'
181  end if
182 
183  if (iand(info%mode, h5md_time) == h5md_time) then
184  call this%id%create_time(this% group, 'id', ps% id, info%mode)
185  else if (iand(info%mode, h5md_linear) == h5md_linear) then
186  if (info%step_offset > 0) then
187  call this%id%create_time(this% group, 'id', ps% id, info%mode, info%step, &
188  info%time, step_offset=info%step_offset, time_offset=info%time_offset)
189  else
190  call this%id%create_time(this% group, 'id', ps% id, info%mode, info%step, &
191  info%time)
192  end if
193  else if (iand(info%mode, h5md_fixed) == h5md_fixed) then
194  call this%id%create_fixed(this% group, 'id', ps% id)
195  else
196  stop 'unknown storage for id in particle_system_io init'
197  end if
198  end if
199 
200  if (this%species_info%store) then
201  info = this%species_info
202  if ( (info%step_offset < 0) .and. (info%time_offset > 0) .or. &
203  (info%step_offset > 0) .and. (info%time_offset < 0) ) then
204  stop 'use both or neither of step_offset and time_offset in ps_init species'
205  end if
206 
207  if (iand(info%mode, h5md_time) == h5md_time) then
208  call this%species%create_time(this% group, 'species', ps% species, info%mode)
209  else if (iand(info%mode, h5md_linear) == h5md_linear) then
210  if (info%step_offset > 0) then
211  call this%species%create_time(this% group, 'species', ps% species, info%mode, info%step, &
212  info%time, step_offset=info%step_offset, time_offset=info%time_offset)
213  else
214  call this%species%create_time(this% group, 'species', ps% species, info%mode, info%step, &
215  info%time)
216  end if
217  else if (iand(info%mode, h5md_fixed) == h5md_fixed) then
218  call this%species%create_fixed(this% group, 'species', ps% species)
219  else
220  stop 'unknown storage for species in particle_system_io init'
221  end if
222  end if
223 
224  end subroutine ps_init
225 
226  subroutine ps_close(this)
227  class(particle_system_io_t), intent(inout) :: this
228 
229  call this% position% close()
230  call h5gclose_f(this% group, this% error)
231 
232  end subroutine ps_close
233 
234  subroutine thermo_init(this, datafile, n_buffer, step, time, step_offset, time_offset)
235  class(thermo_t), intent(out) :: this
236  type(h5md_file_t), intent(inout) :: datafile
237  integer, intent(in) :: n_buffer
238  integer, intent(in) :: step
239  double precision, intent(in) :: time
240  integer, intent(in), optional :: step_offset
241  double precision, intent(in), optional :: time_offset
242 
243  type(h5md_element_t) :: e
244  double precision :: dummy, dummy_vec(3)
245  integer :: mode
246 
247  if (n_buffer <= 0) error stop 'n_buffer non-positive in thermo_init'
248 
249  if ( (present(step_offset) .and. .not. present(time_offset)) .or. &
250  (.not. present(step_offset) .and. present(time_offset)) ) &
251  then
252  stop 'in thermo_init, use both or neither of step_offset and time_offset'
253  end if
254 
255  this% n_buffer = n_buffer
256  this% idx = 0
257 
258  mode = ior(h5md_linear, h5md_store_time)
259 
260  allocate(this% temperature(n_buffer))
261  allocate(this% potential_energy(n_buffer))
262  allocate(this% kinetic_energy(n_buffer))
263  allocate(this% internal_energy(n_buffer))
264  allocate(this% center_of_mass_velocity(3,n_buffer))
265 
266  if (present(step_offset) .and. present(time_offset)) then
267  call e%create_time(datafile%observables, 'temperature', dummy, mode, step, time, &
268  step_offset=step_offset, time_offset=time_offset)
269  call e%close()
270  call e%create_time(datafile%observables, 'potential_energy', dummy, mode, step, time, &
271  step_offset=step_offset, time_offset=time_offset)
272  call e%close()
273  call e%create_time(datafile%observables, 'kinetic_energy', dummy, mode, step, time, &
274  step_offset=step_offset, time_offset=time_offset)
275  call e%close()
276  call e%create_time(datafile%observables, 'internal_energy', dummy, mode, step, time, &
277  step_offset=step_offset, time_offset=time_offset)
278  call e%close()
279  call e%create_time(datafile%observables, 'center_of_mass_velocity', dummy_vec, mode, step, time, &
280  step_offset=step_offset, time_offset=time_offset)
281  call e%close()
282  else
283  call e%create_time(datafile%observables, 'temperature', dummy, mode, step, time)
284  call e%close()
285  call e%create_time(datafile%observables, 'potential_energy', dummy, mode, step, time)
286  call e%close()
287  call e%create_time(datafile%observables, 'kinetic_energy', dummy, mode, step, time)
288  call e%close()
289  call e%create_time(datafile%observables, 'internal_energy', dummy, mode, step, time)
290  call e%close()
291  call e%create_time(datafile%observables, 'center_of_mass_velocity', dummy_vec, mode, step, time)
292  call e%close()
293  end if
294 
295  end subroutine thermo_init
296 
297  subroutine thermo_append(this, datafile, temperature, potential_energy, kinetic_energy, internal_energy, &
298  center_of_mass_velocity, add, force)
299  class(thermo_t), intent(inout) :: this
300  type(h5md_file_t), intent(inout) :: datafile
301  double precision, intent(in) :: temperature, potential_energy, kinetic_energy, internal_energy
302  double precision, intent(in) :: center_of_mass_velocity(3)
303  logical, intent(in), optional :: add, force
304 
305  integer :: i
306  type(h5md_element_t) :: e
307  logical :: do_add, do_append
308 
309  if (present(add)) then
310  do_add = add
311  else
312  do_add = .true.
313  end if
314 
315  if (do_add) then
316  i = this%idx + 1
317  this%idx = i
318 
319  this%temperature(i) = temperature
320  this%potential_energy(i) = potential_energy
321  this%kinetic_energy(i) = kinetic_energy
322  this%internal_energy(i) = internal_energy
323  this%center_of_mass_velocity(:,i) = center_of_mass_velocity
324  end if
325 
326  do_append = (this%idx == this%n_buffer)
327 
328  if (present(force)) then
329  if ((force) .and. (this%idx>0)) then
330  do_append = .true.
331  end if
332  end if
333 
334  if (do_append) then
335  call e%open_time(datafile%observables, 'temperature')
336  call e%append_buffer(this%temperature, force_size=this%idx)
337  call e%close()
338  call e%open_time(datafile%observables, 'potential_energy')
339  call e%append_buffer(this%potential_energy, force_size=this%idx)
340  call e%close()
341  call e%open_time(datafile%observables, 'kinetic_energy')
342  call e%append_buffer(this%kinetic_energy, force_size=this%idx)
343  call e%close()
344  call e%open_time(datafile%observables, 'internal_energy')
345  call e%append_buffer(this%internal_energy, force_size=this%idx)
346  call e%close()
347  call e%open_time(datafile%observables, 'center_of_mass_velocity')
348  call e%append_buffer(this%center_of_mass_velocity, force_size=this%idx)
349  call e%close()
350  this%idx = 0
351  end if
352 
353  end subroutine thermo_append
354 
355 end module particle_system_io
subroutine thermo_append(this, datafile, temperature, potential_energy, kinetic_energy, internal_energy, center_of_mass_velocity, add, force)
Data for particles.
Facilities for particle data I/O.
subroutine ps_init(this, h5md_file, name, ps)
subroutine ps_close(this)
subroutine thermo_init(this, datafile, n_buffer, step, time, step_offset, time_offset)