RMPCDMD
setup_bulk_decay.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 
6  use common
7  use cell_system
9  use hilbert
10  use interaction
11  use hdf5
12  use h5md_module
13  use mpcd
14  use threefry_module
15  use iso_c_binding
16  use omp_lib
17  implicit none
18 
19  type(threefry_rng_t), target, allocatable :: state(:)
20 
21  type(cell_system_t) :: solvent_cells
22  type(particle_system_t) :: solvent
23 
24  type(profile_t) :: tz
25 
26  type(h5md_file_t) :: datafile
27  type(h5md_element_t) :: elem
28  type(h5md_element_t) :: e_solvent, e_solvent_v, e_solvent_spec, e_solvent_id
29  type(h5md_element_t) :: e_solvent_image
30  integer(HID_T) :: box_group, solvent_group
31 
32  integer :: rho, N, L(3)
33  integer :: i, j, k
34  integer :: error, clock, n_threads
35 
36  double precision, parameter :: tau=1.d0
37  double precision :: bulk_rate
38  double precision :: v_com(3)
39 
40  integer, parameter :: n_species = 2
41  integer, dimension(n_species) :: n_solvent
42  type(h5md_element_t) :: n_solvent_el
43 
44  n_threads = omp_get_max_threads()
45  allocate(state(n_threads))
46  call threefry_rng_init(state, -4881001983349061041_c_int64_t)
47 
48  call h5open_f(error)
49 
50  l = 30
51  rho = 9
52  bulk_rate = 0.01
53  n = rho*l(1)*l(2)*l(3)
54 
55  call solvent% init(n, n_species)
56 
57  do i=1, solvent% Nmax
58  solvent% vel(1,i) = threefry_normal(state(1))
59  solvent% vel(2,i) = threefry_normal(state(1))
60  solvent% vel(3,i) = threefry_normal(state(1))
61  end do
62  v_com = sum(solvent% vel, dim=2) / size(solvent% vel, dim=2)
63  solvent% vel = solvent% vel - spread(v_com, dim=2, ncopies=size(solvent% vel, dim=2))
64 
65  solvent% force = 0
66  solvent% species = 2
67  call solvent% random_placement(l*1.d0, state=state(1))
68 
69  call solvent_cells%init(l, 1.d0)
70  solvent_cells%is_reac = .true.
71 
72  call solvent_cells%count_particles(solvent% pos)
73 
74  call datafile% create('bulk_decay.h5', 'RMPCDMD:setup_bulk_decay', 'N/A', 'Pierre de Buyl')
75 
76  call tz% 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  call e_solvent_spec% create_time(solvent_group, 'species', solvent% species, ior(h5md_time, h5md_store_time))
88  call e_solvent_id% create_time(solvent_group, 'id', solvent% id, ior(h5md_time, h5md_store_time))
89  call e_solvent_image% create_time(solvent_group, 'image', solvent% image, ior(h5md_time, h5md_store_time))
90 
91  call n_solvent_el%create_time(datafile%observables, 'n_solvent', &
92  n_solvent, h5md_linear, step=1, time=tau)
93 
94  call solvent% sort(solvent_cells)
95  solvent_cells% origin(1) = threefry_double(state(1)) - 1
96  solvent_cells% origin(2) = threefry_double(state(1)) - 1
97  solvent_cells% origin(3) = threefry_double(state(1)) - 1
98 
99  do i = 1, 10
100  call simple_mpcd_step(solvent, solvent_cells, state)
101  call mpcd_stream_periodic(solvent, solvent_cells, tau)
102  solvent_cells% origin(1) = threefry_double(state(1)) - 1
103  solvent_cells% origin(2) = threefry_double(state(1)) - 1
104  solvent_cells% origin(3) = threefry_double(state(1)) - 1
105  call solvent% sort(solvent_cells)
106  end do
107 
108  solvent% image = 0
109 
110  call update_n_solvent
111 
112  do i = 1, 200
113  call simple_mpcd_step(solvent, solvent_cells, state)
114  v_com = sum(solvent% vel, dim=2) / size(solvent% vel, dim=2)
115 
116  call mpcd_stream_periodic(solvent, solvent_cells, tau)
117 
118  solvent_cells% origin(1) = threefry_double(state(1)) - 1
119  solvent_cells% origin(2) = threefry_double(state(1)) - 1
120  solvent_cells% origin(3) = threefry_double(state(1)) - 1
121  call solvent% sort(solvent_cells)
122 
123  call bulk_reaction(solvent, solvent_cells, 2, 1, bulk_rate, tau, state)
124 
125  call update_n_solvent
126 
127  call e_solvent% append(solvent% pos, i, i*tau)
128  call e_solvent_image% append(solvent% image, i, i*tau)
129  call e_solvent_v% append(solvent% vel, i, i*tau)
130  call e_solvent_spec% append(solvent% species, i, i*tau)
131  call e_solvent_id% append(solvent% id, i, i*tau)
132 
133  write(*,*) i, compute_temperature(solvent, solvent_cells)
134  write(13,*) compute_temperature(solvent, solvent_cells, tz), sum(solvent% vel**2)/(3*solvent% Nmax), v_com
135 
136  end do
137 
138  clock = 0
139  do i = 1 , solvent_cells% N
140  if ( solvent_cells% cell_count(i) > 0 ) clock = clock + 1
141  end do
142  write(*,*) clock, 'filled cells'
143  write(*,*) l(1)*l(2)*l(3), 'actual cells'
144 
145  call e_solvent% close()
146  call e_solvent_v% close()
147  call n_solvent_el%close()
148 
149  call h5gclose_f(solvent_group, error)
150 
151  call tz% norm()
152  call elem% create_fixed(datafile% observables, 'tz', tz% data)
153  call elem% create_fixed(datafile% observables, 'tz_count', tz% count)
154 
155  call h5gclose_f(datafile% observables, error)
156 
157  call datafile% close()
158 
159  call h5close_f(error)
160 
161 contains
162 
163  subroutine update_n_solvent
164  n_solvent = 0
165  do k = 1, solvent%Nmax
166  j = solvent%species(k)
167  if (j <= 0) continue
168  n_solvent(j) = n_solvent(j) + 1
169  end do
170  call n_solvent_el%append(n_solvent)
171  end subroutine update_n_solvent
172 
173 end program setup_bulk_decay
double precision function, public compute_temperature(particles, cells, tz)
Compute the temperature of a mpcd solvent.
Definition: mpcd.f90:343
subroutine, public mpcd_stream_periodic(particles, cells, dt)
Stream MPCD particles in a periodic system.
Definition: mpcd.f90:441
Data for particles.
subroutine, public simple_mpcd_step(particles, cells, state, alpha, thermostat, T, hydro)
Perform a collision.
Definition: mpcd.f90:72
subroutine, public bulk_reaction(p, c, from, to, rate, tau, state)
Apply a bulk unimolecular RMPCD reaction.
Definition: mpcd.f90:606
subroutine update_n_solvent
Compute compact Hilbert indices.
Definition: hilbert.f90:9
Utility routines.
Definition: common.f90:12
Container for a profile, e.g. v(x)
Definition: common.f90:66
integer, parameter dim
Definition: hilbert.f90:15
Routines to perform MPCD dynamics.
Definition: mpcd.f90:18
program setup_bulk_decay
Spatial cells.
Definition: cell_system.f90:11
Lennard-Jones potential definition.
Definition: interaction.f90:10