RMPCDMD
particle_system.f90
Go to the documentation of this file.
1 ! This file is part of RMPCDMD
2 ! Copyright (c) 2015-2017 Pierre de Buyl and contributors
3 ! License: BSD 3-clause (see file LICENSE)
4 
9 
11  use common
12  use interaction
13  implicit none
14  private
15 
16  public :: particle_system_t
18  public :: compute_radial_histogram
19  public :: correct_radial_histogram
21 
23  integer :: nmax
24  integer :: n_species
25  double precision, allocatable :: mass(:)
26  double precision, contiguous, pointer :: pos1(:,:)
27  double precision, contiguous, pointer :: pos2(:,:)
28  double precision, contiguous, pointer :: pos(:,:)
29  double precision, contiguous, pointer :: pos_old(:,:)
30  double precision, contiguous, pointer :: pos_pointer(:,:)
31  double precision, allocatable :: pos_rattle(:,:)
32  double precision, contiguous, pointer :: vel1(:,:)
33  double precision, contiguous, pointer :: vel2(:,:)
34  double precision, contiguous, pointer :: vel(:,:)
35  double precision, contiguous, pointer :: vel_old(:,:)
36  double precision, contiguous, pointer :: vel_pointer(:,:)
37  double precision, contiguous, pointer :: force1(:,:)
38  double precision, contiguous, pointer :: force2(:,:)
39  double precision, contiguous, pointer :: force3(:,:)
40  double precision, contiguous, pointer :: force4(:,:)
41  double precision, contiguous, pointer :: force(:,:)
42  double precision, contiguous, pointer :: force_store(:,:)
43  double precision, contiguous, pointer :: force_old(:,:)
44  double precision, contiguous, pointer :: force_old_store(:,:)
45  double precision, contiguous, pointer :: force_pointer(:,:)
46  integer, contiguous, pointer :: id1(:)
47  integer, contiguous, pointer :: id2(:)
48  integer, contiguous, pointer :: id(:)
49  integer, contiguous, pointer :: id_old(:)
50  integer, contiguous, pointer :: id_pointer(:)
51  integer, contiguous, pointer :: species1(:)
52  integer, contiguous, pointer :: species2(:)
53  integer, contiguous, pointer :: species(:)
54  integer, contiguous, pointer :: species_old(:)
55  integer, contiguous, pointer :: species_pointer(:)
56  integer, contiguous, pointer :: image1(:,:)
57  integer, contiguous, pointer :: image2(:,:)
58  integer, contiguous, pointer :: image(:,:)
59  integer, contiguous, pointer :: image_old(:,:)
60  integer, contiguous, pointer :: image_pointer(:,:)
61  integer, contiguous, pointer :: flags1(:)
62  integer, contiguous, pointer :: flags2(:)
63  integer, contiguous, pointer :: flags(:)
64  integer, contiguous, pointer :: flags_old(:)
65  integer, allocatable :: id_to_idx(:)
66  type(timer_t), pointer :: time_stream, time_step, time_sort, time_count, time_ct
67  type(timer_t), pointer :: time_md_pos, time_md_vel, time_self_force, time_max_disp
68  type(timer_t), pointer :: time_apply_pbc
69  type(timer_t), pointer :: time_rattle_pos, time_rattle_vel
70  type(timer_t), pointer :: time_elastic
71  contains
72  procedure :: init
73  procedure :: init_from_file
74  procedure :: random_placement
75  procedure :: sort
76  procedure :: maximum_displacement
77  end type particle_system_t
78 
79 contains
80 
81  subroutine init(this, Nmax, n_species, mass, system_name)
82  class(particle_system_t), intent(out) :: this
83  integer, intent(in) :: Nmax
84  integer, intent(in), optional :: n_species
85  double precision, intent(in), optional :: mass(:)
86  character(len=*), intent(in), optional :: system_name
87 
88  integer :: i
89  character(len=24) :: system_name_var
90 
91  this% Nmax = nmax
92  if (present(n_species)) then
93  this% n_species = n_species
94  else
95  this% n_species = 1
96  end if
97  if (present(mass)) then
98  this% mass = mass
99  end if
100 
101  allocate(this% pos1(3, nmax))
102  allocate(this% pos2(3, nmax))
103  this% pos => this% pos1
104  this% pos_old => this% pos2
105  allocate(this% pos_rattle(3, nmax))
106 
107  allocate(this% vel1(3, nmax))
108  allocate(this% vel2(3, nmax))
109  this% vel => this% vel1
110  this% vel_old => this% vel2
111 
112  allocate(this% force1(3, nmax))
113  allocate(this% force2(3, nmax))
114  allocate(this% force3(3, nmax))
115  allocate(this% force4(3, nmax))
116  this% force => this% force1
117  this% force_store => this% force2
118  this% force_old => this% force3
119  this% force_old_store => this% force4
120 
121  allocate(this% id1(nmax))
122  allocate(this% id2(nmax))
123  this% id => this% id1
124  this% id_old => this% id2
125 
126  do i = 1, this% Nmax
127  this% id(i) = i
128  end do
129 
130  allocate(this% species1(nmax))
131  allocate(this% species2(nmax))
132  this% species => this% species1
133  this% species_old => this% species2
134 
135  allocate(this% image1(3, nmax))
136  allocate(this% image2(3, nmax))
137  this% image => this% image1
138  this% image_old => this% image2
139 
140  this% image = 0
141 
142  allocate(this% flags1(nmax))
143  allocate(this% flags2(nmax))
144  this% flags => this% flags1
145  this% flags_old => this% flags2
146 
147  this% flags = 0
148 
149  allocate(this%id_to_idx(nmax))
150 
151  if (present(system_name)) then
152  system_name_var = system_name
153  else
154  system_name_var = ''
155  end if
156 
157  allocate(this%time_stream)
158  call this%time_stream%init('stream', trim(system_name_var))
159  allocate(this%time_step)
160  call this%time_step%init('step', trim(system_name_var))
161  allocate(this%time_sort)
162  call this%time_sort%init('sort', trim(system_name_var))
163  allocate(this%time_count)
164  call this%time_count%init('count', trim(system_name_var))
165  allocate(this%time_ct)
166  call this%time_ct%init('compute_temperature', trim(system_name_var))
167  allocate(this%time_md_pos)
168  call this%time_md_pos%init('md_pos', trim(system_name_var))
169  allocate(this%time_md_vel)
170  call this%time_md_vel%init('md_vel', trim(system_name_var))
171  allocate(this%time_self_force)
172  call this%time_self_force%init('self force', trim(system_name_var))
173  allocate(this%time_max_disp)
174  call this%time_max_disp%init('maximum disp', trim(system_name_var))
175  allocate(this%time_rattle_pos)
176  call this%time_rattle_pos%init('rattle_pos', trim(system_name_var))
177  allocate(this%time_rattle_vel)
178  call this%time_rattle_vel%init('rattle_vel', trim(system_name_var))
179  allocate(this%time_elastic)
180  call this%time_elastic%init('elastic_network', trim(system_name_var))
181  allocate(this%time_apply_pbc)
182  call this%time_apply_pbc%init('apply_pbc', trim(system_name_var))
183 
184  end subroutine init
185 
186  subroutine init_from_file(this, filename, group_name, mode, idx)
187  use h5md_module
188  use hdf5
189  implicit none
190  class(particle_system_t), intent(out) :: this
191  character(len=*), intent(in) :: filename, group_name
192  integer, intent(in) :: mode
193  integer, intent(in), optional :: idx
194 
195  type(h5md_file_t) :: f
196  type(h5md_element_t) :: e
197  integer(HID_T) :: g, mem_space, file_space
198  integer(HSIZE_T) :: dims(3), maxdims(3), start(3)
199  integer :: rank
200  double precision, allocatable :: pos(:,:)
201  integer, allocatable :: species(:)
202  logical :: do_species
203 
204  call f% open(filename, h5f_acc_rdonly_f)
205 
206  call h5md_check_valid(f% particles, 'particles group not found')
207 
208  call h5md_check_exists(f% particles, group_name, 'group '//group_name//' not found')
209 
210  call h5gopen_f(f% particles, group_name, g, f% error)
211 
212  if (mode == h5md_fixed) then
213  call e% read_fixed(g, 'position', pos)
214  call h5oexists_by_name_f(g, 'species', do_species, f%error)
215  call this% init(size(pos, 2), 2)
216  if (do_species) call e% read_fixed(g, 'species', species)
217  this% pos = pos
218  deallocate(pos)
219  if (do_species) then
220  this%species = species
221  deallocate(species)
222  end if
223  else if ( (mode == h5md_time) .or. (mode == h5md_linear) ) then
224  if (.not. present(idx) ) error stop 'missing idx in init_from_file'
225  call e% open_time(g, 'position')
226  call this% init(e% Nmax)
227  dims = [3, e% Nmax, 1]
228  call h5screate_simple_f(3, dims, mem_space, e% error)
229 
230  call h5dget_space_f(e% v, file_space, e% error)
231  call h5sget_simple_extent_ndims_f(file_space, rank, e% error)
232  if (rank /= 3) then
233  error stop 'invalid dataset rank'
234  end if
235  call h5sget_simple_extent_dims_f(file_space, dims, maxdims, e% error)
236  start = [0, 0, idx]
237  dims(3) = 1
238  call h5sselect_hyperslab_f(file_space, h5s_select_set_f, start, dims, e% error)
239  call h5dread_f(e% v, h5t_native_double, this% pos, dims, e% error, mem_space, file_space)
240  call h5sclose_f(file_space, e% error)
241  call h5sclose_f(mem_space, e% error)
242  call h5dclose_f(e% v, e% error)
243  end if
244 
245  call h5gclose_f(g, f% error)
246  call f% close()
247 
248  end subroutine init_from_file
249 
250  subroutine random_placement(this, L, other, lj_params, state)
251  use threefry_module
252  implicit none
253  class(particle_system_t), intent(inout) :: this
254  double precision, intent(in) :: L(3)
255  type(particle_system_t), intent(inout), optional :: other
256  type(lj_params_t), intent(in), optional :: lj_params
257  type(threefry_rng_t), intent(inout) :: state
258 
259  integer :: i, j, N_obstacles
260  integer :: s1, s2
261  double precision :: x(3), rsq
262  logical :: tooclose
263 
264  if (present(other) .or. present(lj_params)) then
265  if ( .not. ( present(other) .and. present(lj_params) ) ) then
266  error stop 'other and lj_params must be present/absent together'
267  end if
268  n_obstacles = other% Nmax
269  else
270  n_obstacles = 0
271  end if
272 
273  if (n_obstacles .gt. 0) then
274  do i = 1, this% Nmax
275  s1 = this% species(i)
276  tooclose = .true.
277  do while ( tooclose )
278  x(1) = threefry_double(state) * l(1)
279  x(2) = threefry_double(state) * l(2)
280  x(3) = threefry_double(state) * l(3)
281  tooclose = .false.
282  do j = 1, n_obstacles
283  s2 = other% species(j)
284  rsq = sum(rel_pos(x, other% pos(:, j), l)**2)
285  ! request (s1, s2) as indices for lj are solvent first and colloids second
286  if ( rsq < lj_params% cut_sq(s1, s2) ) then
287  tooclose = .true.
288  exit
289  end if
290  end do
291  end do
292  this% pos(:, i) = x
293  end do
294  else
295  do i = 1, this%Nmax
296  do j = 1, 3
297  this% pos(j, i) = threefry_double(state)*l(j)
298  end do
299  end do
300  end if
301 
302  end subroutine random_placement
303 
304  subroutine sort(this, cells)
306  use hilbert
307  implicit none
308  class(particle_system_t), intent(inout) :: this
309  type(cell_system_t), intent(inout) :: cells
310 
311  integer :: i, idx, start, p(3)
312  integer :: L(3)
313  logical :: nowalls
314 
315  l = cells% L
316  nowalls = .not. cells% has_walls
317 
318  call this%time_count%tic()
319  call cells%count_particles(this% pos)
320  call this%time_count%tac()
321 
322  call this%time_sort%tic()
323  !$omp parallel do private(p, idx, start)
324  do i=1, this% Nmax
325  if (this% species(i) == 0) continue
326  p = cells%cartesian_indices(this%pos(:,i))
327  idx = compact_p_to_h(p, cells% M) + 1
328  !$omp atomic capture
329  start = cells% cell_start(idx)
330  cells% cell_start(idx) = cells% cell_start(idx) + 1
331  !$omp end atomic
332  this% pos_old(:, start) = this% pos(:, i)
333  this% image_old(:,start) = this% image(:,i)
334  this% vel_old(:, start) = this% vel(:, i)
335  this% force_store(:, start) = this% force(:, i)
336  this% force_old_store(:, start) = this% force_old(:, i)
337  this% id_old(start) = this% id(i)
338  this% id_to_idx(this%id(i)) = start
339  this% species_old(start) = this% species(i)
340  this% flags_old(start) = this% flags(i)
341  end do
342 
343  cells% cell_start = cells% cell_start - cells% cell_count
344  call this%time_sort%tac()
345 
346  call switch(this% pos, this% pos_old)
347  call switch(this% image, this% image_old)
348  call switch(this% vel, this% vel_old)
349  call switch(this% force, this% force_store)
350  call switch(this% force_old, this% force_old_store)
351  call switch(this% id, this% id_old)
352  call switch(this% species, this% species_old)
353  call switch(this% flags, this% flags_old)
354 
355  end subroutine sort
356 
357  function maximum_displacement(this) result(r)
358  implicit none
359  class(particle_system_t), intent(inout) :: this
360  double precision :: r
361 
362  integer :: i
363  double precision :: r_sq, rmax_sq
364 
365  rmax_sq = 0
366  call this%time_max_disp%tic()
367  !$omp parallel do private(i, r_sq) reduction(MAX:rmax_sq)
368  do i = 1, this% Nmax
369  if (this% species(i) >= 0) then
370  r_sq = sum((this%pos(:, i)-this%pos_old(:, i))**2)
371  rmax_sq = max(rmax_sq, r_sq)
372  end if
373  end do
374  call this%time_max_disp%tac()
375 
376  r = sqrt(rmax_sq)
377 
378  end function maximum_displacement
379 
380  function cell_maximum_displacement(cells, p, delta_t) result(r)
382  implicit none
383  class(cell_system_t), intent(in) :: cells
384  class(particle_system_t), intent(in) :: p
385  double precision, intent(in) :: delta_t
386  double precision :: r
387 
388  integer :: i, j
389  double precision :: r_sq, rmax_sq, local_rmax_sq
390 
391  rmax_sq = 0
392 
393  !$omp parallel do private(i, j, r_sq) reduction(MAX:rmax_sq)
394  do i = 1, cells%N
395  if (.not. cells%is_md(i)) cycle
396  local_rmax_sq = 0
397  do j = cells%cell_start(i), cells%cell_start(i) + cells%cell_count(i) - 1
398  if (p% species(j) >= 0) then
399  r_sq = sum((p%pos(:, j)-p%pos_old(:, j))**2)
400  if (r_sq > local_rmax_sq) then
401  local_rmax_sq = r_sq
402  end if
403  end if
404  end do
405  rmax_sq = max(rmax_sq, local_rmax_sq)
406  end do
407 
408  r = max(sqrt(rmax_sq), cells%max_v*delta_t)
409 
410  end function cell_maximum_displacement
411 
412  subroutine compute_cylindrical_shell_histogram(hist, x1, x2, L, bin_species, r_min, r_max, solvent)
413  type(histogram_t), intent(inout) :: hist
414  double precision, intent(in), dimension(3) :: x1, x2, L
415  integer, intent(in) :: bin_species
416  double precision, intent(in) :: r_min, r_max
417  type(particle_system_t), intent(in) :: solvent
418 
419  integer :: i, s
420 
421  double precision :: unit_z(3), x(3), delta_x(3), norm, elevation
422  double precision :: r_min_sq, r_max_sq
423 
424  r_min_sq = r_min**2
425  r_max_sq = r_max**2
426 
427  unit_z = rel_pos(x2, x1, l)
428  norm = sqrt(dot_product(unit_z, unit_z))
429  unit_z = unit_z/norm
430 
431  do i = 1, solvent%Nmax
432  s = solvent%species(i)
433  if (s/=bin_species) cycle
434  x = solvent%pos(:,i)
435  delta_x = rel_pos(x, x1, l)
436  elevation = dot_product(unit_z, delta_x)
437  delta_x = delta_x - elevation*unit_z
438  ! only count particles in cylindrical shell
439  norm = dot_product(delta_x,delta_x)
440  if ( (norm>r_min_sq) .and. (norm<r_max_sq)) then
441  call hist%bin(elevation)
442  end if
443  end do
444 
446 
447  subroutine compute_radial_histogram(hist, x1, L, solvent, cells)
449  use cell_system
450  implicit none
451  type(histogram_t), intent(inout) :: hist
452  double precision, intent(in), dimension(3) :: x1, L
453  type(particle_system_t), intent(in) :: solvent
454  type(cell_system_t), intent(in) :: cells
455 
456  integer :: max_i, i, j, k, cell_idx
457  integer :: i_particle, s, start, count
458  integer :: p(3), p_loop(3)
459 
460  double precision :: delta_x(3), norm
461  double precision :: r_max
462 
463  r_max = hist%xmin + hist%dx*hist%n
464 
465  p = cells%cartesian_indices(x1)
466 
467  max_i = floor(r_max/cells%a) + 1
468  p_loop = -max_i
469 
470  do i = -max_i, max_i
471  do j = -max_i, max_i
472  do k = -max_i, max_i
473  cell_idx = compact_p_to_h(modulo(p + [i, j, k], cells%L), cells%M) + 1
474  start = cells%cell_start(cell_idx)
475  count = cells%cell_count(cell_idx)
476 
477  do i_particle = start, start + count - 1
478 
479  s = solvent%species(i_particle)
480  if (s <= 0) cycle
481  delta_x = rel_pos(solvent%pos(:,i_particle), x1, l)
482  norm = norm2(delta_x)
483  if (norm < r_max) then
484  call hist%bin(norm, s)
485  end if
486 
487  end do
488 
489  end do
490  end do
491  end do
492 
493  end subroutine compute_radial_histogram
494 
495  subroutine correct_radial_histogram(hist)
496  type(histogram_t), intent(inout) :: hist
497 
498  integer :: i
499  double precision :: r
500 
501  do i = 1, hist%n
502  r = hist%xmin + (i-0.5d0) * hist%dx
503  hist%data(:,i) = hist%data(:,i) / (4*pi*r**2*hist%dx)
504  end do
505 
506  end subroutine correct_radial_histogram
507 
508 end module particle_system
subroutine init(this, Nmax, n_species, mass, system_name)
pure integer function, public compact_p_to_h(p, m)
Definition: hilbert.f90:242
subroutine, public compute_radial_histogram(hist, x1, L, solvent, cells)
subroutine sort(this, cells)
Data for particles.
double precision function maximum_displacement(this)
subroutine, public correct_radial_histogram(hist)
Compute compact Hilbert indices.
Definition: hilbert.f90:9
Utility routines.
Definition: common.f90:12
subroutine init_from_file(this, filename, group_name, mode, idx)
double precision function, public cell_maximum_displacement(cells, p, delta_t)
pure double precision function, dimension(3), public rel_pos(x, y, L)
Return x-y distance with minimum image convention.
Definition: common.f90:163
Container for a histogram, e.g. p(x)
Definition: common.f90:86
subroutine random_placement(this, L, other, lj_params, state)
Spatial cells.
Definition: cell_system.f90:11
subroutine, public compute_cylindrical_shell_histogram(hist, x1, x2, L, bin_species, r_min, r_max, solvent)
Lennard-Jones potential definition.
Definition: interaction.f90:10
double precision, parameter, public pi
Definition: common.f90:47