RMPCDMD
setup_thermal_fluid.f90
Go to the documentation of this file.
1 ! This file is part of RMPCDMD
2 ! Copyright (c) 2015-2016 Pierre de Buyl and contributors
3 ! License: BSD 3-clause (see file LICENSE)
4 
5 program setup_fluid
6  use common
7  use cell_system
8  use particle_system
9  use hilbert
10  use interaction
11  use hdf5
12  use h5md_module
13  use mpcd
14  use md
15  use threefry_module
16  use iso_c_binding
17  use omp_lib
18  implicit none
19 
20  type(cell_system_t) :: solvent_cells
21  type(particle_system_t) :: solvent
22 
23  type(profile_t) :: tz
24  type(histogram_t) :: rhoz
25 
26  type(h5md_file_t) :: datafile
27  integer(HID_T) :: fields_group
28  type(h5md_element_t) :: elem
29  type(h5md_element_t) :: e_solvent, e_solvent_v
30  type(h5md_element_t) :: elem_tz, elem_tz_count
31  type(h5md_element_t) :: elem_rhoz
32  type(h5md_element_t) :: elem_T
33  integer(HID_T) :: box_group, solvent_group
34 
35  integer, parameter :: N = 10800
36 
37  integer :: i, L(3), clock, error
38 
39  double precision :: v_com(3), wall_v(3,2)
40  double precision :: T
41 
42  double precision, parameter :: tau = 0.1d0
43  type(threefry_rng_t), allocatable :: state(:)
44  integer :: n_threads
45 
46  n_threads = omp_get_max_threads()
47  allocate(state(n_threads))
48  call threefry_rng_init(state,-8053552729396308626_c_int64_t)
49 
50  call h5open_f(error)
51 
52  l = [6, 6, 30]
53 
54  call solvent% init(n)
55 
56  do i=1, solvent% Nmax
57  solvent%vel(1,i) = threefry_normal(state(1))
58  solvent%vel(2,i) = threefry_normal(state(1))
59  solvent%vel(3,i) = threefry_normal(state(1))
60  end do
61  v_com = sum(solvent% vel, dim=2) / size(solvent% vel, dim=2)
62  solvent% vel = solvent% vel - spread(v_com, dim=2, ncopies=size(solvent% vel, dim=2))
63 
64  solvent% force = 0
65  solvent% species = 1
66  call solvent% random_placement(l*1.d0, state=state(1))
67 
68  call solvent_cells%init(l, 1.d0, has_walls=.true.)
69  solvent_cells% origin(3) = -0.5d0
70 
71  call solvent_cells%count_particles(solvent% pos)
72 
73  call datafile% create('data_setup_thermal_fluid.h5', 'RMPCDMD:setup_thermal_fluid', '0.0 dev', 'Pierre de Buyl')
74 
75  call tz% init(0.d0, solvent_cells% edges(3), l(3))
76  call rhoz% init(0.d0, solvent_cells% edges(3), l(3))
77 
78  call h5gcreate_f(datafile% particles, 'solvent', solvent_group, error)
79  call h5gcreate_f(solvent_group, 'box', box_group, error)
80  call h5md_write_attribute(box_group, 'dimension', 3)
81  call h5md_write_attribute(box_group, 'boundary', ['periodic', 'periodic', 'periodic'])
82  call elem% create_fixed(box_group, 'edges', l*1.d0)
83  call h5gclose_f(box_group, error)
84 
85  call e_solvent% create_time(solvent_group, 'position', solvent% pos, ior(h5md_time, h5md_store_time))
86  call e_solvent_v% create_time(solvent_group, 'velocity', solvent% vel, ior(h5md_time, h5md_store_time))
87 
88  call h5gcreate_f(datafile%id, 'fields', fields_group, error)
89  call elem_tz% create_time(fields_group, 'tz', tz% data, &
90  ior(h5md_time, h5md_store_time))
91  call elem_tz_count% create_time(fields_group, 'tz_count', &
92  tz% count, ior(h5md_time, h5md_store_time))
93  call elem_rhoz% create_time(fields_group, 'rhoz', rhoz% data, &
94  ior(h5md_time, h5md_store_time))
95  call h5gclose_f(fields_group, error)
96 
97  call elem_t% create_time(datafile% observables, 'temperature', t, &
98  ior(h5md_time, h5md_store_time))
99 
100  call solvent% sort(solvent_cells)
101  call solvent_cells%count_particles(solvent% pos)
102 
103  wall_v = 0
104 
105  do i = 1, 1000
106  call mpcd_stream_nogravity_zwall(solvent, solvent_cells, tau)
107  call md_vel(solvent, tau)
108  call random_number(solvent_cells% origin)
109  solvent_cells% origin = solvent_cells% origin - 1
110  call apply_pbc(solvent, solvent_cells%edges)
111  call solvent% sort(solvent_cells)
112  call wall_mpcd_step(solvent, solvent_cells, state, &
113  wall_temperature=[0.9d0, 1.1d0], wall_v=wall_v, wall_n=[10, 10])
114  end do
115 
116  do i = 1, 1000
117  call mpcd_stream_nogravity_zwall(solvent, solvent_cells, tau)
118  call md_vel(solvent, tau)
119 
120  call random_number(solvent_cells% origin)
121  solvent_cells% origin = solvent_cells% origin - 1
122  call apply_pbc(solvent, solvent_cells%edges)
123  call solvent% sort(solvent_cells)
124  call wall_mpcd_step(solvent, solvent_cells, state, &
125  wall_temperature=[0.9d0, 1.1d0], wall_v=wall_v, wall_n=[10, 10])
126  v_com = sum(solvent% vel, dim=2) / size(solvent% vel, dim=2)
127 
128  t = compute_temperature(solvent, solvent_cells, tz)
129  write(13,*) t, sum(solvent% vel**2)/(3*solvent% Nmax), v_com
130  call elem_t% append(t, i, i*tau)
131 
132  call compute_rho(solvent, rhoz)
133 
134  if (modulo(i, 50) == 0) then
135  call tz% norm()
136  call elem_tz% append(tz% data, i, i*tau)
137  call elem_tz_count% append(tz% count, i, i*tau)
138  call tz% reset()
139  rhoz% data = rhoz% data / (50.d0 * rhoz% dx)
140  call elem_rhoz% append(rhoz% data, i, i*tau)
141  rhoz% data = 0
142  end if
143 
144  end do
145 
146  call e_solvent% append(solvent% pos, i, i*tau)
147  call e_solvent_v% append(solvent% vel, i, i*tau)
148 
149  clock = 0
150  do i = 1 , solvent_cells% N
151  if ( solvent_cells% cell_count(i) > 0 ) clock = clock + 1
152  end do
153  write(*,*) clock, 'filled cells'
154  write(*,*) l(1)*l(2)*(l(3)+1), 'actual cells'
155 
156  call e_solvent% close()
157  call e_solvent_v% close()
158  call elem_tz% close()
159  call elem_rhoz% close()
160  call elem_t% close()
161 
162  call h5gclose_f(solvent_group, error)
163 
164  call h5gclose_f(datafile% observables, error)
165 
166  call datafile% close()
167 
168  call h5close_f(error)
169 
170 end program setup_fluid
md::md_vel
subroutine, public md_vel(particles, dt)
Definition: md.f90:220
mpcd::wall_mpcd_step
subroutine, public wall_mpcd_step(particles, cells, state, wall_temperature, wall_v, wall_n, thermostat, bulk_temperature, alpha, keep_cell_v)
Collisions in partially filled cells at the walls use the rule of Lamura et al (2001) ...
Definition: mpcd.f90:195
mpcd::compute_rho
subroutine, public compute_rho(particles, rhoz)
Compute density profile along z.
Definition: mpcd.f90:411
mpcd::compute_temperature
double precision function, public compute_temperature(particles, cells, tz)
Compute the temperature of a mpcd solvent.
Definition: mpcd.f90:343
particle_system
Data for particles.
Definition: particle_system.f90:10
setup_fluid
program setup_fluid
Definition: setup_simple_fluid.f90:5
md::apply_pbc
subroutine, public apply_pbc(particles, edges)
Definition: md.f90:203
particle_system::particle_system_t
Definition: particle_system.f90:22
hilbert
Compute compact Hilbert indices.
Definition: hilbert.f90:9
common
Utility routines.
Definition: common.f90:12
common::profile_t
Container for a profile, e.g. v(x)
Definition: common.f90:66
mpcd::mpcd_stream_nogravity_zwall
subroutine, public mpcd_stream_nogravity_zwall(particles, cells, dt)
Definition: mpcd.f90:703
hilbert::dim
integer, parameter dim
Definition: hilbert.f90:15
mpcd
Routines to perform MPCD dynamics.
Definition: mpcd.f90:18
md
Routines to perform Molecular Dynamics integration.
Definition: md.f90:13
common::histogram_t
Container for a histogram, e.g. p(x)
Definition: common.f90:86
cell_system
Spatial cells.
Definition: cell_system.f90:11
cell_system::cell_system_t
Definition: cell_system.f90:20
interaction
Lennard-Jones potential definition.
Definition: interaction.f90:10
Generated by   doxygen 1.8.13