RMPCDMD
single_sphere_thermo_trap.f90
Go to the documentation of this file.
1 ! This file is part of RMPCDMD
2 ! Copyright (c) 2016 Pierre de Buyl and contributors
3 ! License: BSD 3-clause (see file LICENSE)
4 
30 
32  use rmpcdmd_module
33  use hdf5
34  use h5md_module
35  use threefry_module
36  use parsetext
37  use iso_c_binding
38  use omp_lib
39  implicit none
40 
41  integer, parameter :: N_species = 1
42  integer, parameter :: N_species_colloids = 1
43 
44  type(threefry_rng_t), allocatable :: state(:)
45 
46  type(cell_system_t) :: solvent_cells
47  type(particle_system_t) :: solvent
48  type(particle_system_t) :: colloids
49  type(neighbor_list_t) :: neigh
50  type(lj_params_t) :: solvent_colloid_lj
51 
52  type(profile_t) :: tz
53  type(histogram_t) :: rhoz
54  type(profile_t) :: vx
55  type(h5md_element_t) :: elem_tz, elem_tz_count
56  type(h5md_element_t) :: elem_rhoz
57  double precision, allocatable :: v_xz(:,:,:), v_xyz(:,:,:,:)
58  integer, allocatable :: v_xz_count(:,:), v_xyz_count(:,:,:)
59  type(h5md_element_t) :: v_xz_el
60 
61  integer :: rho
62  integer :: N
63  integer :: error
64 
65  double precision :: epsilon(n_species,n_species_colloids)
66  double precision :: sigma(n_species,n_species_colloids), sigma_cut(n_species,n_species_colloids)
67  double precision :: mass(n_species_colloids)
68  double precision :: max_cut
69 
70  double precision :: v_com(3), wall_v(3,2), wall_t(2)
71 
72  double precision :: e1, e2
73  double precision :: tau, dt , T, alpha
74  double precision :: skin, co_max, so_max
75  integer :: N_MD_steps, N_loop
76  integer :: vxz_interval
77  integer :: N_therm
78  integer :: n_extra_sorting
79  double precision :: kin_e, temperature
80 
81  type(h5md_file_t) :: hfile
82  type(h5md_element_t) :: dummy_element
83  integer(HID_T) :: fields_group
84  type(h5md_element_t) :: vx_el
85  type(thermo_t) :: thermo_data
86  type(particle_system_io_t) :: sphere_io
87  type(particle_system_io_t) :: solvent_io
88  integer(HID_T) :: box_group
89 
90  type(pto) :: config
91  integer :: L(3), n_threads
92  integer :: i, j, ijk
93 
94  type(timer_t) :: flag_timer, change_timer, varia
95  integer(HID_T) :: timers_group
96 
97  ! trap parameters
98  double precision :: k, trap_center(3)
99  logical :: move
100  type(args_t) :: args
101 
102  args = get_input_args()
103  call ptparse(config, args%input_file, 11)
104 
105  call flag_timer%init('flag')
106  call change_timer%init('change')
107  call varia%init('varia')
108 
109  n_threads = omp_get_max_threads()
110  allocate(state(n_threads))
111  call threefry_rng_init(state, args%seed)
112 
113  call h5open_f(error)
114 
115  l = ptread_ivec(config, 'L', 3)
116  if (modulo(l(2),2) /= 0) error stop 'non-even Ly is not supported'
117  rho = ptread_i(config, 'rho')
118  n = rho *l(1)*l(2)*l(3)
119  tau =ptread_d(config, 'tau')
120  alpha = ptread_d(config,'alpha')
121 
122  n_md_steps = ptread_i(config, 'N_MD')
123  dt = tau / n_md_steps
124  n_loop = ptread_i(config, 'N_loop')
125  n_therm = ptread_i(config, 'N_therm')
126  vxz_interval = ptread_i(config, 'vxz_interval')
127 
128  wall_t = ptread_dvec(config, 'wall_T', 2)
129  t = ptread_d(config, 'T')
130 
131  k = ptread_d(config, 'k')
132 
133  sigma = ptread_d(config, 'sigma')
134  epsilon = ptread_d(config, 'epsilon')
135 
136  sigma_cut = sigma*2**(1.d0/6.d0)
137  max_cut = maxval(sigma_cut)
138 
139  call solvent_colloid_lj% init(epsilon, sigma, sigma_cut)
140 
141  mass(1) = rho * sigma(1,1)**3 * 4 * 3.14159265/3
142 
143  call solvent% init(n,n_species)
144 
145  move = ptread_l(config, 'move')
146  call colloids% init(1, n_species_colloids, mass)
147  colloids% species(1) = 1
148  colloids% vel = 0
149 
150  call hfile%create(args%output_file, 'RMPCDMD::single_sphere_thermo_trap', &
151  rmpcdmd_revision, 'Pierre de Buyl')
152  call thermo_data%init(hfile, n_buffer=50, step=n_md_steps, time=n_md_steps*dt)
153 
154  call ptkill(config)
155 
156  sphere_io%id_info%store = .false.
157  sphere_io%force_info%store = .true.
158  sphere_io%force_info%mode = ior(h5md_linear,h5md_store_time)
159  sphere_io%force_info%step = n_md_steps
160  sphere_io%force_info%time = n_md_steps*dt
161  sphere_io%position_info%store = .true.
162  sphere_io%position_info%mode = ior(h5md_linear,h5md_store_time)
163  sphere_io%position_info%step = n_md_steps
164  sphere_io%position_info%time = n_md_steps*dt
165  sphere_io%image_info%store = .true.
166  sphere_io%image_info%mode = ior(h5md_linear,h5md_store_time)
167  sphere_io%image_info%step = n_md_steps
168  sphere_io%image_info%time = n_md_steps*dt
169  sphere_io%velocity_info%store = .true.
170  sphere_io%velocity_info%mode = ior(h5md_linear,h5md_store_time)
171  sphere_io%velocity_info%step = n_md_steps
172  sphere_io%velocity_info%time = n_md_steps*dt
173  sphere_io%species_info%store = .true.
174  sphere_io%species_info%mode = h5md_fixed
175  call sphere_io%init(hfile, 'sphere', colloids)
176 
177  solvent_io%force_info%store = .false.
178  solvent_io%id_info%store = .false.
179  solvent_io%position_info%store = .true.
180  solvent_io%position_info%mode = ior(h5md_linear,h5md_store_time)
181  solvent_io%position_info%step = n_loop*n_md_steps
182  solvent_io%position_info%time = n_loop*n_md_steps*dt
183  solvent_io%image_info%store = .true.
184  solvent_io%image_info%mode = ior(h5md_linear,h5md_store_time)
185  solvent_io%image_info%step = n_loop*n_md_steps
186  solvent_io%image_info%time = n_loop*n_md_steps*dt
187  solvent_io%velocity_info%store = .true.
188  solvent_io%velocity_info%mode = ior(h5md_linear,h5md_store_time)
189  solvent_io%velocity_info%step = n_loop*n_md_steps
190  solvent_io%velocity_info%time = n_loop*n_md_steps*dt
191  solvent_io%species_info%store = .true.
192  solvent_io%species_info%mode = ior(h5md_linear,h5md_store_time)
193  solvent_io%species_info%step = n_loop*n_md_steps
194  solvent_io%species_info%time = n_loop*n_md_steps*dt
195  call solvent_io%init(hfile, 'solvent', solvent)
196 
197  solvent% force = 0
198  solvent% species = 1
199  call solvent_cells%init(l, 1.d0, has_walls=.true.)
200 
201  trap_center = solvent_cells%edges/2
202  colloids%pos(:,1) = trap_center
203  colloids%vel = 0
204  colloids%force = 0
205 
206  call vx% init(0.d0, solvent_cells% edges(3), l(3))
207  call tz% init(0.d0, solvent_cells% edges(3), l(3))
208  call rhoz% init(0.d0, solvent_cells% edges(3), l(3))
209 
210  allocate(v_xz_count(l(3), l(1)))
211  allocate(v_xyz_count(l(3), l(2), l(1)))
212  allocate(v_xz(2, l(3), l(1)))
213  allocate(v_xyz(3, l(3), l(2), l(1)))
214 
215  call h5gcreate_f(hfile%id, 'fields', fields_group, error)
216  call vx_el%create_time(fields_group, 'vx', vx%data, ior(h5md_linear,h5md_store_time), &
217  step=n_md_steps, time=n_md_steps*dt)
218  call elem_tz% create_time(fields_group, 'tz', tz% data, ior(h5md_time, h5md_store_time))
219  call elem_tz_count% create_time(fields_group, 'tz_count', tz% count, ior(h5md_time, h5md_store_time))
220  call elem_rhoz% create_time(fields_group, 'rhoz', rhoz% data, ior(h5md_time, h5md_store_time))
221  call v_xz_el%create_time(fields_group, 'v_xz', v_xz, ior(h5md_time, h5md_store_time))
222 
223  call h5gcreate_f(sphere_io%group, 'box', box_group, error)
224  call h5md_write_attribute(box_group, 'dimension', 3)
225  call dummy_element%create_fixed(box_group, 'edges', solvent_cells%edges)
226  call h5gclose_f(box_group, error)
227 
228  call h5gcreate_f(solvent_io%group, 'box', box_group, error)
229  call h5md_write_attribute(box_group, 'dimension', 3)
230  call dummy_element%create_fixed(box_group, 'edges', solvent_cells%edges)
231  call h5gclose_f(box_group, error)
232 
233  call solvent% random_placement(solvent_cells% edges, colloids, solvent_colloid_lj, state(1))
234  solvent% pos_old = solvent% pos
235 
236  do i=1, solvent% Nmax
237  solvent% vel(1,i) = threefry_normal(state(1))*sqrt(t)
238  solvent% vel(2,i) = threefry_normal(state(1))*sqrt(t)
239  solvent% vel(3,i) = threefry_normal(state(1))*sqrt(t)
240  end do
241 
242  call solvent% sort(solvent_cells)
243 
244  skin = 1
245  call neigh% init(colloids% Nmax, 25*floor((max_cut+skin)**2)*rho)
246 
247  n_extra_sorting = 0
248 
249  call neigh% make_stencil(solvent_cells, max_cut+skin+0.1)
250  call neigh% update_list(colloids, solvent, max_cut+skin, solvent_cells)
251 
252  solvent% force = 0
253  colloids% force = 0
254  e1 = compute_force(colloids, solvent, neigh, solvent_cells% edges, solvent_colloid_lj)
255  e2 = compute_force_harmonic_trap(colloids, k, trap_center)
256  solvent% force_old = solvent% force
257  colloids% force_old = colloids% force
258 
259  i = 0
260  wall_v = 0
261  v_xz_count = 0
262  v_xz = 0
263  v_xyz_count = 0
264  v_xyz = 0
265 
266  write(*,*) 'Running for', n_loop+n_therm, 'loops'
267  !start RMPCDMD
268  setup: do i = 1, n_loop+n_therm
269  if (modulo(i,50) == 0) write(*,'(i09)',advance='no') i
270  md_loop: do j = 1, n_md_steps
271  call mpcd_stream_nogravity_zwall(solvent, solvent_cells, dt)
272 
273 
274  if (move) colloids% pos = colloids% pos + dt * colloids% vel + &
275  dt**2 * colloids% force / (2*colloids% mass(1))
276 
277  so_max = solvent% maximum_displacement()
278  co_max = colloids% maximum_displacement()
279 
280  if ( (co_max >= skin*0.1) .or. (so_max >= skin*0.9) ) then
281  call apply_pbc(colloids, solvent_cells% edges)
282  call apply_pbc(solvent, solvent_cells% edges)
283  call solvent% sort(solvent_cells)
284  call neigh% update_list(colloids, solvent, max_cut+skin, solvent_cells)
285  call varia%tic()
286  solvent% pos_old = solvent% pos
287  colloids% pos_old = colloids% pos
288  call varia%tac()
289  n_extra_sorting = n_extra_sorting + 1
290  end if
291 
292  call switch(solvent% force, solvent% force_old)
293  call switch(colloids% force, colloids% force_old)
294 
295  !$omp parallel do
296  do ijk = 1, solvent%Nmax
297  solvent% force(:,ijk) = 0
298  end do
299  colloids% force = 0
300  e1 = compute_force(colloids, solvent, neigh, solvent_cells% edges, solvent_colloid_lj)
301  e2 = compute_force_harmonic_trap(colloids, k, trap_center)
302 
303  call md_vel(solvent, dt)
304  if (move) colloids% vel = colloids% vel + &
305  dt * ( colloids% force + colloids% force_old ) / (2 * colloids% mass(1))
306 
307  end do md_loop
308 
309  call solvent_cells%random_shift(state(1))
310 
311  call apply_pbc(colloids, solvent_cells% edges)
312  call apply_pbc(solvent, solvent_cells% edges)
313  call solvent% sort(solvent_cells)
314  call neigh% update_list(colloids, solvent, max_cut+skin, solvent_cells)
315  solvent%pos_old = solvent% pos
316  colloids%pos_old = colloids% pos
317 
318  call rescale_at_walls
319  call wall_mpcd_step(solvent, solvent_cells, state, &
320  wall_temperature=wall_t, wall_v=wall_v, wall_n=[rho, rho], alpha=alpha)
321 
322  if (i>n_therm) then
323  temperature = compute_temperature(solvent, solvent_cells, tz)
324  kin_e = colloids% mass(1)*sum(colloids% vel(:,1)**2)/2 + &
325  sum(solvent% vel**2)/2
326  v_com = (sum(solvent% vel, dim=2) + mass(1)*colloids%vel(:,1)) / (solvent%Nmax + mass(1))
327  call thermo_data%append(hfile, temperature, e1, kin_e, e1+kin_e, v_com)
328  call compute_rho(solvent, rhoz)
329 
330  call varia%tic()
331  call compute_vx(solvent, vx)
332  call vx% norm()
333  call vx_el%append(vx%data)
334  call vx% reset()
335  call tz% norm()
336  call elem_tz% append(tz% data, i, i*tau)
337  call elem_tz_count% append(tz% count, i, i*tau)
338  call tz% reset()
339  rhoz% data = rhoz% data / rhoz% dx
340  call elem_rhoz% append(rhoz% data, i, i*tau)
341  rhoz% data = 0
343  if (modulo(i-n_therm, vxz_interval)==0) then
344  call div_vxz()
345  call v_xz_el%append(v_xz, i, i*tau)
346  v_xz_count = 0
347  v_xz = 0
348  end if
349  call varia%tac()
350 
351  call sphere_io%position%append(colloids%pos)
352  call sphere_io%velocity%append(colloids%vel)
353  call sphere_io%force%append(colloids%force)
354  call sphere_io%image%append(colloids%image)
355  end if
356 
357  end do setup
358 
359  call thermo_data%append(hfile, temperature, e1+e2, kin_e, e1+e2+kin_e, v_com, add=.false., force=.true.)
360 
361  write(*,*) ''
362  write(*,*) 'n extra sorting', n_extra_sorting
363 
364  call solvent_io%position%append(solvent%pos)
365  call solvent_io%velocity%append(solvent%vel)
366  call solvent_io%image%append(solvent%image)
367  call solvent_io%species%append(solvent%species)
368 
369  call div_vxyz()
370  call dummy_element%create_fixed(fields_group, 'v_xyz', v_xyz)
371  call h5gclose_f(fields_group, error)
372 
373  call h5gcreate_f(hfile%id, 'timers', timers_group, error)
374  call h5md_write_dataset(timers_group, solvent%time_stream%name, solvent%time_stream%total)
375  call h5md_write_dataset(timers_group, solvent%time_md_vel%name, solvent%time_md_vel%total)
376  call h5md_write_dataset(timers_group, solvent%time_step%name, solvent%time_step%total)
377  call h5md_write_dataset(timers_group, solvent%time_count%name, solvent%time_count%total)
378  call h5md_write_dataset(timers_group, solvent%time_sort%name, solvent%time_sort%total)
379  call h5md_write_dataset(timers_group, solvent%time_ct%name, solvent%time_ct%total)
380  call h5md_write_dataset(timers_group, solvent%time_max_disp%name, solvent%time_max_disp%total)
381  call h5md_write_dataset(timers_group, solvent%time_apply_pbc%name, solvent%time_apply_pbc%total)
382  call h5md_write_dataset(timers_group, neigh%time_update%name, neigh%time_update%total)
383  call h5md_write_dataset(timers_group, varia%name, varia%total)
384  call h5md_write_dataset(timers_group, neigh%time_force%name, neigh%time_force%total)
385 
386  call h5md_write_dataset(timers_group, 'total', solvent%time_stream%total + &
387  solvent%time_step%total + solvent%time_count%total + solvent%time_sort%total + &
388  solvent%time_ct%total + solvent%time_md_vel%total + solvent%time_max_disp%total + &
389  flag_timer%total + change_timer%total + solvent%time_apply_pbc%total+ &
390  neigh%time_update%total + varia%total + neigh%time_force%total)
391 
392  call h5gclose_f(timers_group, error)
393 
394  call elem_tz%close()
395  call elem_tz_count%close()
396  call elem_rhoz%close()
397  call sphere_io%close()
398  call hfile%close()
399  call h5close_f(error)
400 
401 contains
402 
403  function compute_force_harmonic_trap(p, k, x0) result(e)
404  type(particle_system_t), intent(inout) :: p
405  double precision, intent(in) :: k, x0(3)
406  double precision :: e
407 
408  integer :: i
409 
410  e = 0
411  do i = 1, p%Nmax
412  p%force(:,i) = p%force(:,i) - k * (p%pos(:,i)-x0)
413  e = e + k * sum( (p%pos(:,i)-x0)**2 ) / 2
414  end do
415 
416  end function compute_force_harmonic_trap
417 
418  subroutine rescale_at_walls
420  integer :: i
421  integer :: cell_idx, n, start, wall_idx, cell(3)
422  double precision :: local_v(3), local_k, local_T, factor
423 
424  do cell_idx = 1, solvent_cells% N
425  if (solvent_cells% cell_count(cell_idx) <= 1) cycle
426 
427  start = solvent_cells% cell_start(cell_idx)
428  n = solvent_cells% cell_count(cell_idx)
429 
430  ! Find whether we are in a wall cell
431  cell = compact_h_to_p(cell_idx - 1, solvent_cells% M) + 1
432  if (cell(3) == 1) then
433  wall_idx = 1
434  else if (cell(3) == solvent_cells% L(3)) then
435  wall_idx = 2
436  local_t = 1.1
437  else
438  wall_idx = -1
439  end if
440  if (wall_idx==-1) cycle
441  local_t = wall_t(wall_idx)
442 
443  local_v = 0
444  do i = start, start + n - 1
445  local_v = local_v + solvent% vel(:, i)
446  end do
447  local_v = local_v / n
448 
449  local_k = 0
450  do i = start, start + n - 1
451  local_k = local_k + sum((solvent% vel(:, i)-local_v)**2)/2
452  end do
453  local_k = local_k/(3*(n-1))
454  factor = sqrt(local_t/(2*local_k))
455 
456  do i = start, start + n - 1
457  solvent% vel(:, i) = local_v + factor*(solvent% vel(:, i)-local_v)
458  end do
459 
460  end do
461 
462  end subroutine rescale_at_walls
463 
464  subroutine compute_vxz_and_vxyz
465  integer :: i, s, ix, iy, iz
466 
467  do i = 1, solvent%Nmax
468  s = solvent%species(i)
469  if (s <= 0) cycle
470  ix = modulo(floor(solvent%pos(1,i)/solvent_cells%a), l(1)) + 1
471  iy = modulo(floor(solvent%pos(2,i)/solvent_cells%a), l(2)) + 1
472  iz = modulo(floor(solvent%pos(3,i)/solvent_cells%a), l(3)) + 1
473  v_xyz_count(iz, iy, ix) = v_xyz_count(iz, iy, ix) + 1
474  v_xyz(:, iz, iy, ix) = v_xyz(:, iz, iy, ix) + solvent%vel(:, i)
475  if ( (iy == l(2)/2) .or. (iy == 1+l(2)/2) ) then
476  v_xz_count(iz, ix) = v_xz_count(iz, ix) + 1
477  v_xz(1, iz, ix) = v_xz(1, iz, ix) + solvent%vel(1, i)
478  v_xz(2, iz, ix) = v_xz(2, iz, ix) + solvent%vel(3, i)
479  end if
480  end do
481 
482  end subroutine compute_vxz_and_vxyz
483 
484  subroutine div_vxz()
486  integer :: i, j
487  do i = 1, l(1)
488  do j = 1, l(3)
489  if (v_xz_count(j, i) > 0) then
490  v_xz(:, j, i) = v_xz(:, j, i) / v_xz_count(j, i)
491  end if
492  end do
493  end do
494 
495  end subroutine div_vxz
496 
497  subroutine div_vxyz()
499  integer :: i, j, k
500  do i = 1, l(1)
501  do j = 1, l(2)
502  do k = 1, l(3)
503  if (v_xyz_count(k, j, i) > 0) then
504  v_xyz(:, k, j, i) = v_xyz(:, k, j, i) / v_xyz_count(k, j, i)
505  end if
506  end do
507  end do
508  end do
509 
510  end subroutine div_vxyz
511 
512 end program single_sphere_thermo_trap
double precision function compute_force_harmonic_trap(p, k, x0)
subroutine rescale_at_walls
subroutine compute_vxz_and_vxyz
subroutine div_vxyz()
subroutine div_vxz()
program single_sphere_thermo_trap
Simulate a thermal gradient with an embedded colloidal sphere.