RMPCDMD
three_bead_enzyme.f90
Go to the documentation of this file.
1 ! This file is part of RMPCDMD
2 ! Copyright (c) 2018 Pierre de Buyl
3 ! License: BSD 3-clause (see file LICENSE)
4 
35 
37  use rmpcdmd_module
38  use hdf5
39  use h5md_module
40  use threefry_module
41  use parsetext
42  use iso_c_binding
43  use omp_lib
44  implicit none
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  type(lj_params_t) :: colloid_lj
52 
53  integer, parameter :: N_species = 3
54  integer, parameter :: N_species_colloids = 2
55  integer :: N_colloids, N_enzymes
56 
57  integer :: rho
58  integer :: N
59  integer :: error
60 
61  double precision :: sigma_N, sigma_E, max_cut
62  double precision :: epsilon(n_species,n_species_colloids)
63  double precision :: sigma(n_species,n_species_colloids), sigma_cut(n_species,n_species_colloids)
64  double precision, allocatable :: mass(:)
65 
66  double precision :: elastic_k, link_d(2)
67  double precision :: angle_k, link_angles(2)
68 
69  double precision :: e1, e2, e3
70  double precision :: tau, dt , T
71  double precision :: skin, co_max, so_max
72  integer :: N_MD_steps, N_loop
73  integer :: n_extra_sorting
74 
75  type(pto) :: config
76 
77  integer :: i, L(3)
78  double precision :: current_time
79  integer :: j, k
80  type(timer_t), target :: varia, main, time_flag, time_change
81  type(timer_t), target :: time_select, time_release, time_reset
82  double precision :: total_time
83  type(timer_list_t) :: timer_list
84  integer(HID_T) :: timers_group
85 
86  type(threefry_rng_t), allocatable :: state(:)
87  integer :: n_threads
88  type(h5md_file_t) :: hfile
89  type(h5md_element_t) :: dummy_element
90  integer(HID_T) :: fields_group, params_group
91  integer(HID_T) :: kinetics_group
92  integer(HID_T) :: box_group
93  integer(HID_T) :: dummy_id
94  type(thermo_t) :: thermo_data
95  double precision :: temperature, kin_e
96  double precision :: v_com(3)
97  double precision :: tmp_vec(3), normal_vec(3), tmp_q(4)
98  type(particle_system_io_t) :: enzyme_io
99  double precision :: bulk_rate(2)
100  logical :: bulk_rmpcd
101  logical :: immediate_chemistry
102 
103  integer :: enzyme_i
104  integer :: enz_1, enz_2, enz_3
105  double precision :: dist
106  double precision :: enzyme_capture_radius
107  double precision :: proba_p, proba_s
108  double precision :: rate_release_p, rate_release_s, total_rate
109  double precision :: substrate_fraction, product_relative_fraction
110  logical :: placement_too_close
111 
112  !! State variables for the enzyme kinetics
113  integer, allocatable :: bound_molecule_id(:)
114  double precision, allocatable :: next_reaction_time(:)
115  logical, allocatable :: enzyme_bound(:)
116  double precision, allocatable :: link_angle(:)
117 
118  integer, dimension(N_species) :: n_solvent
119  type(h5md_element_t) :: n_solvent_el
120  type(h5md_element_t) :: link_angle_el
121 
122  integer, dimension(N_species) :: n_plus, n_minus
123  type(h5md_element_t) :: n_plus_el, n_minus_el
124 
125  !! Histogram variable
126  type(histogram_t), allocatable :: radial_hist(:)
127  double precision, allocatable :: dummy_hist(:,:,:)
128 
129  !! Kinetics data
130  type(enzyme_kinetics_t), allocatable :: kinetics_data(:)
131 
132  integer :: equilibration_loops
133  integer :: colloid_sampling
134  logical :: sampling
135  type(args_t) :: args
136 
137  args = get_input_args()
138  call ptparse(config, args%input_file, 11)
139 
140  n_threads = omp_get_max_threads()
141  allocate(state(n_threads))
142  call threefry_rng_init(state, args%seed)
143 
144  call main%init('main')
145  call varia%init('varia')
146  call time_flag%init('flag')
147  call time_change%init('change')
148  call time_select%init('select_substrate')
149  call time_release%init('release_substrate')
150  call time_reset%init('reset_bit')
151 
152  call timer_list%init(20)
153  call timer_list%append(varia)
154  call timer_list%append(time_flag)
155  call timer_list%append(time_change)
156  call timer_list%append(time_select)
157  call timer_list%append(time_release)
158  call timer_list%append(time_reset)
159 
160  call h5open_f(error)
161  call hfile%create(args%output_file, 'RMPCDMD::three_bead_enzyme', &
162  rmpcdmd_revision, 'Pierre de Buyl')
163  call h5gcreate_f(hfile%id, 'parameters', params_group, error)
164  call hdf5_util_write_dataset(params_group, 'seed', args%seed)
165 
166  bulk_rmpcd = ptread_l(config, 'bulk_rmpcd', loc=params_group)
167  immediate_chemistry = ptread_l(config, 'immediate_chemistry', loc=params_group)
168  bulk_rate = ptread_dvec(config, 'bulk_rate', 2, loc=params_group)
169 
170  proba_s = ptread_d(config, 'proba_s', loc=params_group)
171  proba_p = ptread_d(config, 'proba_p', loc=params_group)
172  rate_release_s = ptread_d(config, 'rate_release_s', loc=params_group)
173  rate_release_p = ptread_d(config, 'rate_release_p', loc=params_group)
174  total_rate = rate_release_s + rate_release_p
175  enzyme_capture_radius = ptread_d(config, 'enzyme_capture_radius', loc=params_group)
176  substrate_fraction = ptread_d(config, 'substrate_fraction', loc=params_group)
177  product_relative_fraction = ptread_d(config, 'product_relative_fraction', loc=params_group)
178 
179  l = ptread_ivec(config, 'L', 3, loc=params_group)
180  rho = ptread_i(config, 'rho', loc=params_group)
181 
182  t = ptread_d(config, 'T', loc=params_group)
183  link_d = ptread_dvec(config, 'link_d', 2, loc=params_group)
184  link_angles = ptread_dvec(config, 'link_angles', 2, loc=params_group)
185  elastic_k = ptread_d(config, 'elastic_k', loc=params_group)
186  angle_k = ptread_d(config, 'angle_k', loc=params_group)
187 
188  tau = ptread_d(config, 'tau', loc=params_group)
189  n_md_steps = ptread_i(config, 'N_MD', loc=params_group)
190  colloid_sampling = ptread_i(config, 'colloid_sampling', loc=params_group)
191  if (colloid_sampling <= 0) then
192  error stop 'colloid_sampling must be positive'
193  end if
194  dt = tau / n_md_steps
195  n_loop = ptread_i(config, 'N_loop', loc=params_group)
196  equilibration_loops = ptread_i(config, 'equilibration_loops', loc=params_group)
197 
198  n_enzymes = ptread_i(config, 'N_enzymes', loc=params_group)
199  n_colloids = 3*n_enzymes
200 
201  allocate(link_angle(n_enzymes))
202  allocate(next_reaction_time(n_enzymes))
203  allocate(bound_molecule_id(n_enzymes))
204  allocate(enzyme_bound(n_enzymes))
205 
206  link_angle = link_angles(1)
207 
208  sigma_e = ptread_d(config, 'sigma_E', loc=params_group)
209  sigma_n = ptread_d(config, 'sigma_N', loc=params_group)
210 
211  n = rho * ( l(1)*l(2)*l(3) - int(n_enzymes*4*pi/3*(2*sigma_n**3 + sigma_e**3)) )
212 
213  ! solvent index first, colloid index second, in solvent_colloid_lj
214  epsilon(:,1) = ptread_dvec(config, 'epsilon_E', 3, loc=params_group)
215  epsilon(:,2) = ptread_dvec(config, 'epsilon_N', 3, loc=params_group)
216 
217  sigma(:,1) = sigma_e
218  sigma(:,2) = sigma_n
219  sigma_cut = sigma*2**(1.d0/6.d0)
220  max_cut = maxval(sigma_cut)
221 
222  if (max_cut > enzyme_capture_radius) then
223  write(*,*) 'enzyme_capture_radius must be greater than max_cut'
224  stop
225  end if
226 
227  call solvent_colloid_lj% init(epsilon, sigma, sigma_cut)
228 
229  epsilon(:,:) = ptread_d(config, 'epsilon_colloid', loc=params_group)
230 
231  sigma(1,1) = 2*sigma_e
232  sigma(1,2) = sigma_e + sigma_n
233  sigma(2,1) = sigma_e + sigma_n
234  sigma(2,2) = 2*sigma_n
235  sigma = sigma*1.1d0
236  sigma_cut = sigma*2**(1.d0/6.d0)
237 
238  call colloid_lj% init(epsilon(1:n_species_colloids,1:n_species_colloids), &
239  sigma(1:n_species_colloids,1:n_species_colloids), &
240  sigma_cut(1:n_species_colloids,1:n_species_colloids))
241 
242  allocate(mass(3*n_enzymes))
243  do i = 1, n_enzymes
244  mass(3*(i-1)+1) = rho * sigma_n**3 * 4 * pi/3
245  mass(3*(i-1)+2) = rho * sigma_e**3 * 4 * pi/3
246  mass(3*(i-1)+3) = rho * sigma_n**3 * 4 * pi/3
247  end do
248 
249 
250  allocate(radial_hist(n_colloids))
251  do i = 1, n_colloids
252  call radial_hist(i)%init(min(sigma_e, sigma_n)/4, max(sigma_e, sigma_n)*3, 100, n_species=n_species)
253  end do
254 
255  call solvent% init(n, n_species, system_name='solvent')
256 
257  call colloids% init(n_colloids, n_species_colloids, mass, system_name='colloids')
258  deallocate(mass)
259 
260  call thermo_data%init(hfile, n_buffer=50, step=n_md_steps, time=n_md_steps*dt, &
261  step_offset=n_md_steps, time_offset=n_md_steps*dt)
262 
263  call h5gclose_f(params_group, error)
264  call ptkill(config)
265 
266  allocate(kinetics_data(n_enzymes))
267  do enzyme_i = 1, n_enzymes
268  call kinetics_data(enzyme_i)%bind_substrate%init(block_size=32)
269  call kinetics_data(enzyme_i)%release_substrate%init(block_size=32)
270  call kinetics_data(enzyme_i)%bind_product%init(block_size=32)
271  call kinetics_data(enzyme_i)%release_product%init(block_size=32)
272  end do
273 
274  do enzyme_i = 1, n_enzymes
275  colloids% species(3*(enzyme_i-1)+1) = 2
276  colloids% species(3*(enzyme_i-1)+2) = 1
277  colloids% species(3*(enzyme_i-1)+3) = 2
278  end do
279 
280  colloids% vel = 0
281  colloids% force = 0
282 
283  enzyme_io%force_info%store = .false.
284  enzyme_io%id_info%store = .false.
285  enzyme_io%position_info%store = .true.
286  enzyme_io%position_info%mode = ior(h5md_linear,h5md_store_time)
287  enzyme_io%position_info%step = colloid_sampling
288  enzyme_io%position_info%step_offset = colloid_sampling
289  enzyme_io%position_info%time = colloid_sampling*n_md_steps*dt
290  enzyme_io%position_info%time_offset = colloid_sampling*n_md_steps*dt
291  enzyme_io%image_info%store = .true.
292  enzyme_io%image_info%mode = ior(h5md_linear,h5md_store_time)
293  enzyme_io%image_info%step = enzyme_io%position_info%step
294  enzyme_io%image_info%step_offset = enzyme_io%position_info%step_offset
295  enzyme_io%image_info%time = enzyme_io%position_info%time
296  enzyme_io%image_info%time_offset = enzyme_io%position_info%time_offset
297  enzyme_io%velocity_info%store = .true.
298  enzyme_io%velocity_info%mode = ior(h5md_linear,h5md_store_time)
299  enzyme_io%velocity_info%step = enzyme_io%position_info%step
300  enzyme_io%velocity_info%step_offset = enzyme_io%position_info%step_offset
301  enzyme_io%velocity_info%time = enzyme_io%position_info%time
302  enzyme_io%velocity_info%time_offset = enzyme_io%position_info%time_offset
303  enzyme_io%species_info%store = .true.
304  enzyme_io%species_info%mode = h5md_fixed
305  call enzyme_io%init(hfile, 'enzyme', colloids)
306 
307  do k = 1, solvent%Nmax
308  solvent% vel(1,k) = threefry_normal(state(1))*sqrt(t)
309  solvent% vel(2,k) = threefry_normal(state(1))*sqrt(t)
310  solvent% vel(3,k) = threefry_normal(state(1))*sqrt(t)
311  if (dble(k)/dble(solvent%Nmax) <= substrate_fraction) then
312  if (dble(k)/dble(solvent%Nmax)/substrate_fraction <= product_relative_fraction) then
313  solvent%species(k) = 2
314  else
315  solvent%species(k) = 1
316  end if
317  else
318  ! neutral species
319  solvent%species(k) = 3
320  end if
321  end do
322  solvent% vel = solvent% vel - spread(sum(solvent% vel, dim=2)/solvent% Nmax, 2, solvent% Nmax)
323  solvent% force = 0
324 
325  call solvent_cells%init(l, 1.d0)
326 
327  ! place enzymes without overlap
328  do enzyme_i = 1, n_enzymes
329  enz_1 = 3*(enzyme_i-1)+1
330  enz_2 = enz_1+1
331  enz_3 = enz_1+2
332 
333  placement_too_close = .true.
334  do while (placement_too_close)
335 
336  colloids% pos(:,enz_2) = [ threefry_double(state(1))*solvent_cells%edges(1), &
337  threefry_double(state(1))*solvent_cells%edges(2), &
338  threefry_double(state(1))*solvent_cells%edges(3) ]
339 
340  ! Pick one vector to place bead 1 and another one normal to the first
341  tmp_vec = rand_sphere(state(1))
342  normal_vec = rand_sphere(state(1))
343  normal_vec = normal_vec - tmp_vec * dot_product(normal_vec, tmp_vec)
344  ! Place normal_vec at an angle link_angles(1) w.r.t. tmp_vec
345  tmp_q = qnew(s=cos(link_angles(1)/2), v=sin(link_angles(1)/2)*normal_vec/norm2(normal_vec))
346  normal_vec = qvector(qmul(tmp_q, qmul(qnew(v=tmp_vec), qconj(tmp_q))))
347 
348  colloids%pos(:,enz_1) = colloids%pos(:,enz_2) + tmp_vec*link_d(1)
349  colloids%pos(:,enz_3) = colloids%pos(:,enz_2) + normal_vec*link_d(2)
350 
351  placement_too_close = .false.
352 
353  do j = 1, 3*(enzyme_i-1)
354  tmp_vec = colloids%pos(:,j)
355  do i = 1, 3
356  dist = norm2(rel_pos(colloids%pos(:,(enzyme_i-1)*3+i), &
357  tmp_vec, solvent_cells%edges))
358  if ( dist < &
359  colloid_lj%sigma(colloids%species((enzyme_i-1)*3+i), colloids%species(j)) &
360  ) then
361  placement_too_close = .true.
362  end if
363  end do
364  end do
365  end do
366  end do
367 
368  call n_solvent_el%create_time(hfile%observables, 'n_solvent', &
369  n_solvent, ior(h5md_linear, h5md_store_time), step=n_md_steps, &
370  time=n_md_steps*dt, step_offset=n_md_steps, time_offset=n_md_steps*dt)
371 
372  call n_plus_el%create_time(hfile%observables, 'n_plus', &
373  n_plus, ior(h5md_linear, h5md_store_time), step=n_md_steps, &
374  time=n_md_steps*dt, step_offset=n_md_steps, time_offset=n_md_steps*dt)
375 
376  call n_minus_el%create_time(hfile%observables, 'n_minus', &
377  n_minus, ior(h5md_linear, h5md_store_time), step=n_md_steps, &
378  time=n_md_steps*dt, step_offset=n_md_steps, time_offset=n_md_steps*dt)
379 
380  call link_angle_el%create_time(hfile%observables, 'link_angle', &
381  link_angle, ior(h5md_linear, h5md_store_time), step=n_md_steps, &
382  time=n_md_steps*dt, step_offset=n_md_steps, time_offset=n_md_steps*dt)
383 
384  call h5gcreate_f(enzyme_io%group, 'box', box_group, error)
385  call h5md_write_attribute(box_group, 'dimension', 3)
386  call dummy_element%create_fixed(box_group, 'edges', solvent_cells%edges)
387  call h5gclose_f(box_group, error)
388 
389 
390 
391  call solvent% random_placement(solvent_cells% edges, colloids, solvent_colloid_lj, state(1))
392 
393  call solvent% sort(solvent_cells)
394 
395  call neigh% init(colloids% Nmax, int(rho*30*max(sigma_e,sigma_n)**3))
396 
397  skin = 1
398  n_extra_sorting = 0
399 
400  call neigh% make_stencil(solvent_cells, enzyme_capture_radius+skin)
401 
402  call neigh% update_list(colloids, solvent, enzyme_capture_radius+skin, solvent_cells, solvent_colloid_lj)
403 
404  e1 = compute_force(colloids, solvent, neigh, solvent_cells% edges, solvent_colloid_lj)
405  e2 = compute_force_n2(colloids, solvent_cells% edges, colloid_lj)
406  e3 = compute_bead_force()
407 
408  solvent% force_old = solvent% force
409  colloids% force_old = colloids% force
410 
411  call h5fflush_f(hfile%id, h5f_scope_global_f, error)
412 
413  write(*,*) 'Box size', l
414  write(*,*) n, 'fluid particles'
415 
416  write(*,*) 'Running for', equilibration_loops, '+', n_loop, 'loops'
417 
418  call main%tic()
419  sampling = .false.
420  enzyme_bound = .false.
421  next_reaction_time = huge(next_reaction_time(1))
422 
423  do i = 0, n_loop+equilibration_loops
424  if (i==equilibration_loops) sampling = .true.
425  if (modulo(i,256) == 0) write(*,'(i08)',advance='no') i
426 
427  ! reset substrate/product creation/destruction counters
428  n_plus = 0
429  n_minus = 0
430 
431  md_loop: do j = 1, n_md_steps
432 
433  current_time = (i-equilibration_loops)*tau + (j-1)*dt
434 
435  call md_pos(solvent, dt)
436  do k=1, colloids% Nmax
437  colloids% pos(:,k) = colloids% pos(:,k) + dt * colloids% vel(:,k) + &
438  dt**2 * colloids% force(:,k) / (2 * colloids%mass(k))
439  end do
440 
441  so_max = solvent% maximum_displacement()
442  co_max = colloids% maximum_displacement()
443 
444  ! select substrate for binding
445  ! requires neigbor list with enzyme_capture_radius + skin
446  if (sampling .or. immediate_chemistry) then
447  call select_substrate
448  call time_release%tic()
449  ! Check if unbinding should occur
450  do enzyme_i = 1, n_enzymes
451  if (current_time >= next_reaction_time(enzyme_i)) then
452  if (threefry_double(state(1)) < rate_release_s / total_rate) then
453  ! release substrate
454  call unbind_molecule(enzyme_i, 1)
455  call fix_neighbor_list(solvent%id_to_idx(bound_molecule_id(enzyme_i)))
456  n_plus(1) = n_plus(1) + 1
457  if (sampling) call kinetics_data(enzyme_i)%release_substrate%append(current_time)
458  else
459  ! release product
460  call unbind_molecule(enzyme_i, 2)
461  call fix_neighbor_list(solvent%id_to_idx(bound_molecule_id(enzyme_i)))
462  n_plus(2) = n_plus(2) + 1
463  if (sampling) call kinetics_data(enzyme_i)%release_product%append(current_time)
464  end if
465  next_reaction_time(enzyme_i) = huge(next_reaction_time(enzyme_i))
466  end if
467  end do
468  call time_release%tac()
469  end if
470 
471 
472  if ( co_max + so_max >= skin ) then
473  call apply_pbc(solvent, solvent_cells% edges)
474  call apply_pbc(colloids, solvent_cells% edges)
475  call solvent% sort(solvent_cells)
476  call neigh% update_list(colloids, solvent, enzyme_capture_radius + skin, solvent_cells, solvent_colloid_lj)
477  call varia%tic()
478  !$omp parallel do
479  do k = 1, solvent%Nmax
480  solvent% pos_old(:,k) = solvent% pos(:,k)
481  end do
482  colloids% pos_old = colloids% pos
483  n_extra_sorting = n_extra_sorting + 1
484  call varia%tac()
485  end if
486 
487  call varia%tic()
488  call switch(solvent% force, solvent% force_old)
489  call switch(colloids% force, colloids% force_old)
490 
491  !$omp parallel do
492  do k = 1, solvent%Nmax
493  solvent% force(:,k) = 0
494  end do
495  colloids% force = 0
496  call varia%tac()
497  e1 = compute_force(colloids, solvent, neigh, solvent_cells% edges, solvent_colloid_lj)
498  e2 = compute_force_n2(colloids, solvent_cells% edges, colloid_lj)
499  e3 = compute_bead_force()
500 
501  call md_vel(solvent, dt)
502 
503  call varia%tic()
504  do k=1, colloids% Nmax
505  colloids% vel(:,k) = colloids% vel(:,k) + &
506  dt * ( colloids% force(:,k) + colloids% force_old(:,k) ) / (2 * colloids%mass(k))
507  end do
508  call varia%tac()
509 
510 
511  end do md_loop
512 
513 
514 
515  call solvent_cells%random_shift(state(1))
516  call apply_pbc(solvent, solvent_cells% edges)
517  call apply_pbc(colloids, solvent_cells% edges)
518  call solvent% sort(solvent_cells)
519  call neigh% update_list(colloids, solvent, enzyme_capture_radius+skin, solvent_cells, solvent_colloid_lj)
520  call varia%tic()
521  !$omp parallel do
522  do k = 1, solvent%Nmax
523  solvent% pos_old(:,k) = solvent% pos(:,k)
524  end do
525  colloids% pos_old = colloids% pos
527  call varia%tac()
528 
529  !! TODO: add thermostat and hydro option
530  call simple_mpcd_step(solvent, solvent_cells, state)
531 
532  if (bulk_rmpcd) then
533  call bulk_reaction(solvent, solvent_cells, 1, 2, bulk_rate(1), tau, state)
534  call bulk_reaction(solvent, solvent_cells, 2, 1, bulk_rate(2), tau, state)
535  end if
536 
537  call varia%tic()
538  !$omp parallel do
539  do k = 1, solvent%Nmax
540  solvent% pos_old(:,k) = solvent% pos(:,k)
541  end do
542  colloids% pos_old = colloids% pos
543  call varia%tac()
544 
545  call varia%tic()
546 
547  if (sampling) then
548  temperature = compute_temperature(solvent, solvent_cells)
549  kin_e = 0
550  v_com = 0
551  !$omp parallel do private(k) reduction(+:kin_e) reduction(+:v_com)
552  do k = 1, solvent%Nmax
553  if (solvent%species(k)>0) then
554  kin_e = kin_e + sum(solvent%vel(:,k)**2) / 2
555  v_com = v_com + solvent%vel(:,k)
556  end if
557  end do
558  do k = 1, colloids%Nmax
559  v_com = v_com + colloids%mass(k) * colloids%vel(:,k)
560  kin_e = kin_e + colloids%mass(k) * sum(colloids%vel(:,k)**2) / 2
561  end do
562 
563  call thermo_data%append(hfile, temperature, e1+e2+e3, kin_e, e1+e2+e3+kin_e, v_com)
564 
565  ! update radial histogram
566  !$omp parallel do
567  do k = 1, colloids%Nmax
568  call compute_radial_histogram(radial_hist(k), colloids%pos(:,k), solvent_cells%edges, solvent, solvent_cells)
569  end do
570  n_solvent = 0
571  !$omp parallel do private(k,j) reduction(+:n_solvent)
572  do k = 1, solvent%Nmax
573  j = solvent%species(k)
574  if (j <= 0) cycle
575  n_solvent(j) = n_solvent(j) + 1
576  end do
577  call n_solvent_el%append(n_solvent)
578  call n_plus_el%append(n_plus)
579  call n_minus_el%append(n_minus)
580  call link_angle_el%append(link_angle)
581  if (modulo(i, colloid_sampling)==0) then
582  call enzyme_io%position%append(colloids%pos)
583  call enzyme_io%velocity%append(colloids%vel)
584  call enzyme_io%image%append(colloids%image)
585  end if
586 
587  end if
588 
589  call varia%tac()
590 
591  end do
592  call main%tac()
593  write(*,*) ''
594 
595  write(*,*) 'n extra sorting', n_extra_sorting
596 
597  call thermo_data%append(hfile, temperature, e1+e2+e3, kin_e, e1+e2+e3+kin_e, &
598  v_com, add=.false., force=.true.)
599 
600  ! write solvent coordinates for last step
601 
602  call h5gcreate_f(hfile%id, 'fields', fields_group, error)
603 
604  ! copy radial histogram data
605  allocate(dummy_hist(size(radial_hist(1)%data, 1), &
606  size(radial_hist(1)%data, 2), size(radial_hist)))
607  do i = 1, size(radial_hist)
608  call correct_radial_histogram(radial_hist(i))
609  dummy_hist(:,:,i) = radial_hist(i)%data / dble(n_loop)
610  end do
611 
612  call dummy_element%create_fixed(fields_group, 'radial_histogram', dummy_hist)
613 
614  call h5oopen_f(fields_group, 'radial_histogram', dummy_id, error)
615  call h5md_write_attribute(dummy_id, 'xmin', radial_hist(1)%xmin)
616  call h5md_write_attribute(dummy_id, 'dx', radial_hist(1)%dx)
617  call h5oclose_f(dummy_id, error)
618 
619  call h5gclose_f(fields_group, error)
620 
621 
622  ! Store timing data
623  call timer_list%append(solvent%time_stream)
624  call timer_list%append(solvent%time_step)
625  call timer_list%append(solvent%time_count)
626  call timer_list%append(solvent%time_sort)
627  call timer_list%append(solvent%time_ct)
628  call timer_list%append(solvent%time_md_pos)
629  call timer_list%append(solvent%time_md_vel)
630  call timer_list%append(solvent%time_max_disp)
631  call timer_list%append(solvent%time_apply_pbc)
632  call timer_list%append(colloids%time_self_force)
633  call timer_list%append(colloids%time_apply_pbc)
634  call timer_list%append(neigh%time_update)
635  call timer_list%append(neigh%time_force)
636  call timer_list%append(colloids%time_max_disp)
637 
638  call h5gcreate_f(hfile%id, 'timers', timers_group, error)
639  call timer_list%write(timers_group, total_time)
640  call h5md_write_dataset(timers_group, 'total', total_time)
641  call h5md_write_dataset(timers_group, main%name, main%total)
642  call h5gclose_f(timers_group, error)
643 
644  call h5gcreate_f(hfile%id, 'kinetics', kinetics_group, error)
645  do i = 1, n_enzymes
646  call h5gcreate_f(kinetics_group, numbered_string('enzyme_', i, 4), dummy_id, error)
647  k = kinetics_data(i)%bind_substrate%current_idx
648  call h5md_write_dataset(dummy_id, 'bind_substrate', kinetics_data(i)%bind_substrate%data(1:k))
649  k = kinetics_data(i)%release_substrate%current_idx
650  call h5md_write_dataset(dummy_id, 'release_substrate', kinetics_data(i)%release_substrate%data(1:k))
651  k = kinetics_data(i)%bind_product%current_idx
652  call h5md_write_dataset(dummy_id, 'bind_product', kinetics_data(i)%bind_product%data(1:k))
653  k = kinetics_data(i)%release_product%current_idx
654  call h5md_write_dataset(dummy_id, 'release_product', kinetics_data(i)%release_product%data(1:k))
655  call h5gclose_f(dummy_id, error)
656  end do
657  call h5gclose_f(kinetics_group, error)
658 
659  call enzyme_io%close()
660  call n_solvent_el%close()
661  call n_plus_el%close()
662  call n_minus_el%close()
663  call link_angle_el%close()
664  call hfile%close()
665  call h5close_f(error)
666 
667 contains
668 
669  function compute_bead_force() result(en)
670  double precision :: en
671 
672  integer :: i, i_enzyme
673  double precision :: f(3), r12(3), r
674  double precision :: r32(3), d12, d32, costheta
675 
676  en = 0
677 
678  do i_enzyme = 1, n_enzymes
679  do i = 1, 2
680  r12 = rel_pos(colloids%pos(:,(i_enzyme-1)*3+i), colloids%pos(:,(i_enzyme-1)*3+i+1), solvent_cells%edges)
681  r = norm2(r12)
682  en = en + elastic_k*(r-link_d(i))**2/2
683  f = -elastic_k*(r-link_d(i))*r12/r
684  colloids%force(:,(i_enzyme-1)*3+i) = colloids%force(:,(i_enzyme-1)*3+i) + f
685  colloids%force(:,(i_enzyme-1)*3+i+1) = colloids%force(:,(i_enzyme-1)*3+i+1) - f
686  end do
687 
688  r12 = rel_pos(colloids%pos(:,(i_enzyme-1)*3+1), colloids%pos(:,(i_enzyme-1)*3+2), solvent_cells%edges)
689  d12 = norm2(r12)
690  r32 = rel_pos(colloids%pos(:,(i_enzyme-1)*3+3), colloids%pos(:,(i_enzyme-1)*3+2), solvent_cells%edges)
691  d32 = norm2(r32)
692  costheta = dot_product(r12, r32)/(d12*d32)
693  f = -fprime(costheta, link_angle(i_enzyme)) * (r32/(d32*d12) - costheta * r12/d12**2)
694  colloids%force(:,(i_enzyme-1)*3+1) = colloids%force(:,(i_enzyme-1)*3+1) + f
695  colloids%force(:,(i_enzyme-1)*3+2) = colloids%force(:,(i_enzyme-1)*3+2) - f
696  f = -fprime(costheta, link_angle(i_enzyme)) * (r12/(d12*d32) - costheta * r32/d32**2)
697  colloids%force(:,(i_enzyme-1)*3+3) = colloids%force(:,(i_enzyme-1)*3+3) + f
698  colloids%force(:,(i_enzyme-1)*3+2) = colloids%force(:,(i_enzyme-1)*3+2) - f
699 
700  en = en + angle_k * (acos(costheta)-link_angle(i_enzyme))**2 / 2
701  end do
702 
703  end function compute_bead_force
704 
705  function compute_bead_energy(enzyme_idx, factor) result(en)
706  integer, intent(in) :: enzyme_idx
707  integer, intent(in) :: factor
708  double precision :: en
709 
710  double precision :: r12(3), d12, r32(3), d32, costheta, f(3)
711 
712  r12 = rel_pos(colloids%pos(:,(enzyme_idx-1)*3+1), colloids%pos(:,(enzyme_idx-1)*3+2), solvent_cells%edges)
713  d12 = norm2(r12)
714  r32 = rel_pos(colloids%pos(:,(enzyme_idx-1)*3+3), colloids%pos(:,(enzyme_idx-1)*3+2), solvent_cells%edges)
715  d32 = norm2(r32)
716  costheta = dot_product(r12, r32)/(d12*d32)
717 
718  f = -factor*fprime(costheta, link_angle(enzyme_idx)) * (r32/(d32*d12) - costheta * r12/d12**2)
719  colloids%force(:,(enzyme_idx-1)*3+1) = colloids%force(:,(enzyme_idx-1)*3+1) + f
720  colloids%force(:,(enzyme_idx-1)*3+2) = colloids%force(:,(enzyme_idx-1)*3+2) - f
721  f = -factor*fprime(costheta, link_angle(enzyme_idx)) * (r12/(d12*d32) - costheta * r32/d32**2)
722  colloids%force(:,(enzyme_idx-1)*3+3) = colloids%force(:,(enzyme_idx-1)*3+3) + f
723  colloids%force(:,(enzyme_idx-1)*3+2) = colloids%force(:,(enzyme_idx-1)*3+2) - f
724 
725  en = angle_k * (acos(costheta)-link_angle(enzyme_idx))**2 / 2
726 
727  end function compute_bead_energy
728 
730  function fprime(x, theta_0) result(r)
731  double precision, intent(in) :: x
732  double precision, intent(in) :: theta_0
733  double precision :: r
734 
735  r = - angle_k * (acos(x) - theta_0) / sqrt(1-x**2)
736 
737  end function fprime
738 
739 
744  subroutine select_substrate
746  integer :: i, m, s_sp
747  double precision :: dist, x_enzyme(3)
748  double precision :: total_p, xi
749  double precision :: proba_something
750 
751  integer, parameter :: list_size = 256
752  integer :: n_s, n_p
753  integer :: idx
754  integer :: list_s(list_size), list_p(list_size)
755  logical :: reaction_happens
756 
757  integer :: enzyme_i, enz_2
758  integer :: s_found, p_found
759 
760  call time_select%tic()
761 
762  select_substrate_loop: do enzyme_i = 1, n_enzymes
763 
764  if (enzyme_bound(enzyme_i)) cycle select_substrate_loop
765 
766  enz_2 = (enzyme_i-1)*3 + 2
767 
768  n_s = 0
769  n_p = 0
770  s_found = 0
771  p_found = 0
772 
773  x_enzyme = modulo(colloids%pos(:,enz_2), solvent_cells%edges)
774 
775  select_loop: do i = 1, neigh%n(enz_2)
776  m = neigh%list(i, enz_2)
777 
778  s_sp = solvent%species(m)
779 
780  ! skip neutral fluid particles
781  if ((s_sp == 0) .or. (s_sp == 3)) cycle select_loop
782 
783  ! only pick particles that are *not* doing MD but are in the range following the
784  ! position update
785 
786  dist = norm2(rel_pos(x_enzyme, solvent%pos(:,m), solvent_cells%edges))
787 
788  if (.not. btest(solvent%flags(m), md_bit) .and. &
789  dist <= solvent_colloid_lj%cut(s_sp, colloids%species(enz_2)) &
790  ) then
791  ! select for current round
792  solvent%flags(m) = ibset(solvent%flags(m), enzyme_region_bit)
793  if (s_sp==1) then
794  n_s = n_s + 1
795  if (n_s > list_size) then
796  stop 'exceed size of list_s in select_substrate'
797  end if
798  list_s(n_s) = m
799  else if (s_sp==2) then
800  n_p = n_p + 1
801  if (n_p > list_size) then
802  stop 'exceed size of list_p in select_substrate'
803  end if
804  list_p(n_p) = m
805  end if
806  end if
807  end do select_loop
808 
809  total_p = n_s*proba_s + n_p*proba_p
810  proba_something = 1 - ((1-proba_s)**n_s)*((1-proba_p)**n_p)
811  reaction_happens = (threefry_double(state(1)) < proba_something)
812 
813  ! Pick either a substrate or a product if a reaction happens
814 
815  if (reaction_happens) then
816  xi = threefry_double(state(1))*total_p
817 
818  if (xi < n_s*proba_s) then
819  ! pick in substrates
820  idx = floor(xi * n_s / total_p) + 1
821  m = list_s(idx)
822  n_minus(1) = n_minus(1) + 1
823  ! use list_s(idx)
824  else
825  ! pick in products
826  idx = floor(xi * n_p / total_p) + 1
827  m = list_p(idx)
828  n_minus(2) = n_minus(2) + 1
829  end if
830  call bind_molecule(enzyme_i, m)
831  end if
832 
833  end do select_substrate_loop
834 
835  call time_select%tac()
836 
837  end subroutine select_substrate
838 
839 
841  subroutine bind_molecule(enzyme_idx, idx)
842  integer, intent(in) :: enzyme_idx
843  integer, intent(in) :: idx
844  double precision :: com_v(3), w_ab(3), excess_kinetic_energy, mu_ab
845  double precision :: en1, en2
846  integer :: p(3), cell_idx
847  integer :: enz_2
848 
849  ! adjust velocity
850  ! compute excess kinetic energy
851 
852  enz_2 = 3*(enzyme_idx-1)+2
853 
854  com_v = (colloids%vel(:,enz_2)*colloids%mass(enz_2) + solvent%vel(:,idx)) &
855  / (colloids%mass(enz_2) + 1)
856  w_ab = colloids%vel(:,enz_2) - solvent%vel(:,idx)
857 
858  mu_ab = 1/(1+1/colloids%mass(enz_2))
859 
860  excess_kinetic_energy = mu_ab * dot_product(w_ab, w_ab) / 2
861 
862  colloids%vel(:,enz_2) = com_v
863 
864  p = solvent_cells%cartesian_indices(solvent%pos(:, idx))
865 
866  cell_idx = compact_p_to_h(p, solvent_cells%M) + 1
867 
868  ! transfer mass
869  colloids%mass(enz_2) = colloids%mass(enz_2) + 1
870 
871  bound_molecule_id(enzyme_idx) = solvent%id(idx)
872  if (sampling) then
873  if (solvent%species(idx) == 1) then
874  call kinetics_data(enzyme_idx)%bind_substrate%append(current_time)
875  else
876  call kinetics_data(enzyme_idx)%bind_product%append(current_time)
877  end if
878  end if
879 
880  solvent%species(idx) = 0
881 
882  ! compute conformational energy difference when changing the angle
883  en1 = compute_bead_energy(enzyme_idx, factor=-1)
884  link_angle(enzyme_idx) = link_angles(2)
885  en2 = compute_bead_energy(enzyme_idx, factor=1)
886 
887  call add_energy_to_cell(cell_idx, excess_kinetic_energy + en1 - en2)
888 
889  e3 = e3 + en2 - en1
890 
891  enzyme_bound(enzyme_idx) = .true.
892 
893  ! sample from Poisson process
894  next_reaction_time(enzyme_idx) = current_time - log(threefry_double(state(1)))/total_rate
895 
896  end subroutine bind_molecule
897 
899  subroutine unbind_molecule(enzyme_idx, to_species)
900  integer, intent(in) :: enzyme_idx
901  integer, intent(in) :: to_species
902 
903  integer :: i, p(3), cell_idx
904  double precision :: x_new(3), dist, en1, en2, min_r, delta_r, r
905  logical :: too_close
906  integer :: enz_1, enz_2, enz_3
907 
908  enz_1 = 3*(enzyme_idx-1)+1
909  enz_2 = enz_1+1
910  enz_3 = enz_1+2
911 
912  min_r = solvent_colloid_lj%cut(to_species, 1)
913  delta_r = enzyme_capture_radius - min_r
914 
915  ! transfer mass
916  colloids%mass(enz_2) = colloids%mass(enz_2) - 1
917 
918  ! Place the molecule outside of the interaction range of all colloids
919  too_close = .true.
920  placement_loop: do while (too_close)
921  r = min_r + 1d-3
922  x_new = colloids%pos(:,enz_2) + rand_sphere(state(1))*r
923  x_new = modulo(x_new, solvent_cells%edges)
924  p = solvent_cells%cartesian_indices(x_new)
925  if (solvent_cells%cell_count(compact_p_to_h(p, solvent_cells%M)+1) < 3) cycle placement_loop
926  too_close = .false.
927  do i = 1, 3*n_enzymes
928  dist = norm2(rel_pos(colloids%pos(:,i), x_new, solvent_cells%edges))
929  too_close = too_close .or. &
930  (dist < solvent_colloid_lj%cut(to_species, colloids%species(i)))
931  end do
932  end do placement_loop
933 
934  ! Retrieve location in particle array
935  i = solvent%id_to_idx(bound_molecule_id(enzyme_idx))
936 
937  solvent%species(i) = to_species
938  solvent%pos(:,i) = x_new
939  solvent%flags(i) = ibset(solvent%flags(i), enzyme_region_bit)
940 
941  ! Use c.o.m. velocity so that no kinetic energy exchange must take place
942  solvent%vel(:,i) = colloids%vel(:,enz_2)
943 
944  ! compute conformational energy difference when changing the angle
945  en1 = compute_bead_energy(enzyme_idx, factor=-1)
946  link_angle(enzyme_idx) = link_angles(1)
947  en2 = compute_bead_energy(enzyme_idx, factor=1)
948 
949  p = solvent_cells%cartesian_indices(solvent%pos(:, i))
950 
951  cell_idx = compact_p_to_h(p, solvent_cells%M) + 1
952  call add_energy_to_cell(cell_idx, en1 - en2)
953 
954  e3 = e3 + en2 - en1
955 
956  enzyme_bound(enzyme_idx) = .false.
957 
958  end subroutine unbind_molecule
959 
960  subroutine add_energy_to_cell(cell_idx, energy)
961  integer, intent(in) :: cell_idx
962  double precision, intent(in) :: energy
963 
964  integer :: i, start, n, n_effective
965  integer :: cell(3), actual_cell(3), actual_idx, cell_shift(3)
966  integer :: counter
967  double precision :: com_v(3), kin, factor
968  double precision :: remaining_energy, xi(3), change
969 
970  integer, parameter :: max_counter = 100
971 
972  start = solvent_cells% cell_start(cell_idx)
973  n = solvent_cells% cell_count(cell_idx)
974 
975  counter = 1
976  remaining_energy = energy
977  cell = compact_h_to_p(cell_idx - 1, solvent_cells%M) + 1
978 
979  ! If energy is positive, go for one cell.
980  ! todo: always go around several cells
981 
982  do counter = 1, max_counter
983  xi(1) = -1 + 3*threefry_double(state(1))
984  xi(2) = -1 + 3*threefry_double(state(1))
985  xi(3) = -1 + 3*threefry_double(state(1))
986 
987  cell_shift = floor(xi)
988  actual_cell = modulo(cell + cell_shift , solvent_cells% L)
989  actual_idx = compact_p_to_h(actual_cell-1, solvent_cells%M) + 1
990 
991  start = solvent_cells% cell_start(actual_idx)
992  n = solvent_cells% cell_count(actual_idx)
993 
994  n_effective = 0
995  com_v = 0
996  do i = start, start + n - 1
997  if (solvent%species(i) > 0) then
998  com_v = com_v + solvent% vel(:, i)
999  n_effective = n_effective + 1
1000  end if
1001  end do
1002 
1003  if (n_effective <= 1) cycle
1004 
1005  com_v = com_v / n_effective
1006 
1007  kin = 0
1008  do i = start, start + n - 1
1009  if (solvent%species(i) > 0) then
1010  kin = kin + sum((solvent%vel(:,i)-com_v)**2)
1011  end if
1012  end do
1013  kin = kin / 2
1014 
1015  if (energy > 0) then
1016  change = energy
1017  else
1018  change = max(-kin/2, remaining_energy)
1019  end if
1020 
1021  remaining_energy = remaining_energy - change
1022 
1023  factor = sqrt((kin + change)/kin)
1024  do i = start, start + n - 1
1025  solvent%vel(:,i) = com_v + factor*(solvent%vel(:,i) - com_v)
1026  end do
1027 
1028  if (abs(remaining_energy) < 1d-9) exit
1029  end do
1030 
1031  if (abs(remaining_energy) > 1d-4) then
1032  write(*,*) ''
1033  write(*,*) 'energy', energy
1034  write(*,*) 'counter', counter
1035  stop 'error in add_energy_to_cell'
1036  end if
1037 
1038 
1039  end subroutine add_energy_to_cell
1040 
1041  subroutine reset_enzyme_region_bit
1043  integer :: cell_idx
1044  integer :: i, start, n
1045 
1046  call time_reset%tic()
1047 
1048  !$omp parallel do private(cell_idx, i, start, n)
1049  do cell_idx = 1, solvent_cells%N
1050 
1051  start = solvent_cells%cell_start(cell_idx)
1052  n = solvent_cells%cell_count(cell_idx)
1053 
1054  if (.not. solvent_cells%is_md(cell_idx)) then
1055  do i = start, start + n - 1
1056  solvent%flags(i) = ibclr(solvent%flags(i), enzyme_region_bit)
1057  end do
1058  end if
1059 
1060  end do
1061 
1062  call time_reset%tac()
1063 
1064  end subroutine reset_enzyme_region_bit
1065 
1067  subroutine fix_neighbor_list(idx)
1068  integer, intent(in) :: idx
1069 
1070  integer :: s, colloid_i, n
1071  double precision :: x(3), d
1072 
1073  s = solvent%species(idx)
1074  x = solvent%pos(:,idx)
1075 
1076  do colloid_i = 1, 3*n_enzymes
1077  d = norm2(rel_pos(x, colloids%pos(:,colloid_i), solvent_cells%edges))
1078  if (d <= solvent_colloid_lj%cut(s, colloids%species(colloid_i)) + skin) then
1079  ! add to nlist
1080  n = neigh%n(colloid_i)
1081  if ( n < neigh%Nmax ) then
1082  n = n+1
1083  neigh%list(n, colloid_i) = idx
1084  neigh%n(colloid_i) = n
1085  end if
1086  end if
1087  end do
1088 
1089  end subroutine fix_neighbor_list
1090 
1091 end program three_bead_enzyme
subroutine select_substrate
Select substrate molecules for binding to enzyme.
double precision function compute_bead_energy(enzyme_idx, factor)
double precision function fprime(x, theta_0)
derivative of the angular harmonic arccos term
subroutine bind_molecule(enzyme_idx, idx)
Bind solvent molecule idx to enzyme.
subroutine fix_neighbor_list(idx)
Check all colloids&#39; neighbor lists to possibly add solvent particle idx.
subroutine add_energy_to_cell(cell_idx, energy)
subroutine reset_enzyme_region_bit
subroutine unbind_molecule(enzyme_idx, to_species)
Bind solvent molecule idx to enzyme.
program three_bead_enzyme
Simulate a three bead enzyme model.
double precision function compute_bead_force()