single_dimer_pbc.f90

Go to the documentation of this file.

! This file is part of RMPCDMD
! Copyright (c) 2015-2017 Pierre de Buyl and contributors
! License: BSD 3-clause (see file LICENSE)


program single_dimer_pbc
  use rmpcdmd_module
  use hdf5
  use h5md_module
  use threefry_module
  use parsetext
  use iso_c_binding
  use omp_lib
  implicit none

  type(cell_system_t) :: solvent_cells
  type(particle_system_t) :: solvent
  type(particle_system_t) :: colloids
  type(neighbor_list_t) :: neigh
  type(lj_params_t) :: solvent_colloid_lj
  type(lj_params_t) :: colloid_lj

  integer, parameter :: N_species = 2

  integer :: rho
  integer :: N
  integer :: error

  double precision :: sigma_N, sigma_C, max_cut
  double precision :: epsilon(2,2)
  double precision :: sigma(2,2), sigma_cut(2,2)
  double precision :: mass(2)

  double precision :: e1, e2
  double precision :: tau, dt , T
  double precision :: d,prob
  double precision :: skin, co_max, so_max
  integer :: N_MD_steps, N_loop, collide_every
  integer :: n_extra_sorting

  type(pto) :: config

  integer :: i, L(3)
  integer :: j, k
  type(timer_t), target :: varia, main, time_flag, time_refuel, time_change
  double precision :: total_time
  type(timer_list_t) :: timer_list
  integer(HID_T) :: timers_group

  type(threefry_rng_t), allocatable :: state(:)
  integer :: n_threads
  type(h5md_file_t) :: hfile
  type(h5md_element_t) :: dummy_element
  integer(HID_T) :: fields_group, params_group
  integer(HID_T) :: box_group
  type(thermo_t) :: thermo_data
  double precision :: temperature, kin_e
  double precision :: v_com(3)
  double precision :: com_pos(3)
  double precision :: unit_r(3)
  type(particle_system_io_t) :: dimer_io
  type(particle_system_io_t) :: solvent_io
  double precision :: bulk_rate
  logical :: bulk_rmpcd

  integer, dimension(N_species) :: n_solvent
  type(h5md_element_t) :: n_solvent_el

  integer, parameter :: block_length = 8
  type(axial_correlator_t) :: axial_cf
  integer(HID_T) :: correlator_group

  integer :: equilibration_loops
  integer :: colloid_sampling
  logical :: sampling
  logical :: collision_step
  logical :: do_hydro, do_thermostat

  type(histogram_t) :: z_hist
  type(h5md_element_t) :: z_hist_el
  double precision :: cyl_shell_rmin, cyl_shell_rmax
  type(args_t) :: args

  args = get_input_args()
  call ptparse(config, args%input_file, 11)

  n_threads = omp_get_max_threads()
  allocate(state(n_threads))
  call threefry_rng_init(state, args%seed)

  call main%init('main')
  call varia%init('varia')
  call time_flag%init('flag')
  call time_refuel%init('refuel')
  call time_change%init('change')

  call timer_list%init(16)
  call timer_list%append(varia)
  call timer_list%append(time_flag)
  call timer_list%append(time_refuel)
  call timer_list%append(time_change)

  call h5open_f(error)
  call hfile%create(args%output_file, 'RMPCDMD::single_dimer_pbc', &
       rmpcdmd_revision, 'Pierre de Buyl')
  call h5gcreate_f(hfile%id, 'parameters', params_group, error)
  call hdf5_util_write_dataset(params_group, 'seed', args%seed)
  prob = ptread_d(config,'probability', loc=params_group)
  bulk_rmpcd = ptread_l(config, 'bulk_rmpcd', loc=params_group)
  bulk_rate = ptread_d(config, 'bulk_rate', loc=params_group)

  l = ptread_ivec(config, 'L', 3, loc=params_group)
  rho = ptread_i(config, 'rho', loc=params_group)
  n = rho *l(1)*l(2)*l(3)

  t = ptread_d(config, 'T', loc=params_group)
  d = ptread_d(config, 'd', loc=params_group)
  
  dt = ptread_d(config, 'dt', loc=params_group)
  n_md_steps = ptread_i(config, 'N_MD', loc=params_group)

  collide_every = ptread_i(config, 'collide_every', loc=params_group)
  tau = dt*n_md_steps*collide_every

  n_loop = ptread_i(config, 'N_loop', loc=params_group)
  equilibration_loops = ptread_i(config, 'equilibration_loops', loc=params_group)

  colloid_sampling = ptread_i(config, 'colloid_sampling', loc=params_group)
  if (modulo(n_md_steps, colloid_sampling) /= 0) then
     error stop 'colloid_sampling must divide N_MD with no remainder'
  end if

  sigma_c = ptread_d(config, 'sigma_C', loc=params_group)
  sigma_n = ptread_d(config, 'sigma_N', loc=params_group)

  ! solvent index first, colloid index second, in solvent_colloid_lj
  epsilon(:,1) = ptread_dvec(config, 'epsilon_C', 2, loc=params_group)
  epsilon(:,2) = ptread_dvec(config, 'epsilon_N', 2, loc=params_group)

  sigma(:,1) = sigma_c
  sigma(:,2) = sigma_n
  sigma_cut = sigma*2**(1.d0/6.d0)
  max_cut = maxval(sigma_cut)

  call solvent_colloid_lj% init(epsilon, sigma, sigma_cut)

  epsilon(1,1) = ptread_d(config, 'epsilon_C_C', loc=params_group)
  epsilon(1,2) = ptread_d(config, 'epsilon_N_C', loc=params_group)
  epsilon(2,1) = epsilon(1,2)
  epsilon(2,2) = ptread_d(config, 'epsilon_N_N', loc=params_group)

  sigma(1,1) = 2*sigma_c
  sigma(1,2) = sigma_c + sigma_n
  sigma(2,1) = sigma_c + sigma_n
  sigma(2,2) = 2*sigma_n
  sigma_cut = sigma*2**(1.d0/6.d0)

  call colloid_lj% init(epsilon, sigma, sigma_cut)

  mass(1) = rho * sigma_c**3 * 4 * 3.14159265/3
  mass(2) = rho * sigma_n**3 * 4 * 3.14159265/3

  call solvent% init(n,2, system_name='solvent') !there will be 2 species of solvent particles

  call colloids% init(2,2, mass, system_name='colloids') !there will be 2 species of colloids

  call thermo_data%init(hfile, n_buffer=50, step=n_md_steps, time=n_md_steps*dt)

  do_hydro = ptread_l(config, 'do_hydro', loc=params_group)
  do_thermostat = ptread_l(config, 'do_thermostat', loc=params_group)

  call h5gclose_f(params_group, error)
  call ptkill(config)

  call axial_cf%init(block_length, n_loop, n_loop*n_md_steps)
  
  colloids% species(1) = 1
  colloids% species(2) = 2
  colloids% vel = 0
  colloids% force = 0

  dimer_io%force_info%store = .false.
  dimer_io%id_info%store = .false.
  dimer_io%position_info%store = .true.
  dimer_io%position_info%mode = ior(h5md_linear,h5md_store_time)
  dimer_io%position_info%step = colloid_sampling
  dimer_io%position_info%time = colloid_sampling*dt
  dimer_io%image_info%store = .true.
  dimer_io%image_info%mode = ior(h5md_linear,h5md_store_time)
  dimer_io%image_info%step = colloid_sampling
  dimer_io%image_info%time = colloid_sampling*dt
  dimer_io%velocity_info%store = .true.
  dimer_io%velocity_info%mode = ior(h5md_linear,h5md_store_time)
  dimer_io%velocity_info%step = colloid_sampling
  dimer_io%velocity_info%time = colloid_sampling*dt
  dimer_io%species_info%store = .true.
  dimer_io%species_info%mode = h5md_fixed
  call dimer_io%init(hfile, 'dimer', colloids)

  solvent_io%force_info%store = .false.
  solvent_io%id_info%store = .false.
  solvent_io%position_info%store = .true.
  solvent_io%position_info%mode = ior(h5md_linear,h5md_store_time)
  solvent_io%position_info%step = n_loop*n_md_steps
  solvent_io%position_info%time = n_loop*n_md_steps*dt
  solvent_io%image_info%store = .true.
  solvent_io%image_info%mode = ior(h5md_linear,h5md_store_time)
  solvent_io%image_info%step = n_loop*n_md_steps
  solvent_io%image_info%time = n_loop*n_md_steps*dt
  solvent_io%velocity_info%store = .true.
  solvent_io%velocity_info%mode = ior(h5md_linear,h5md_store_time)
  solvent_io%velocity_info%step = n_loop*n_md_steps
  solvent_io%velocity_info%time = n_loop*n_md_steps*dt
  solvent_io%species_info%store = .true.
  solvent_io%species_info%mode = ior(h5md_linear,h5md_store_time)
  solvent_io%species_info%step = n_loop*n_md_steps
  solvent_io%species_info%time = n_loop*n_md_steps*dt
  call solvent_io%init(hfile, 'solvent', solvent)

  do k = 1, solvent%Nmax
     solvent% vel(1,k) = threefry_normal(state(1))*sqrt(t)
     solvent% vel(2,k) = threefry_normal(state(1))*sqrt(t)
     solvent% vel(3,k) = threefry_normal(state(1))*sqrt(t)
  end do
  solvent% vel = solvent% vel - spread(sum(solvent% vel, dim=2)/solvent% Nmax, 2, solvent% Nmax)
  solvent% force = 0
  solvent% species = 1

  call solvent_cells%init(l, 1.d0)
  colloids% pos(:,1) = solvent_cells% edges/2.0
  colloids% pos(:,2) = solvent_cells% edges/2.0 
  colloids% pos(1,2) = colloids% pos(1,2) + d

  call n_solvent_el%create_time(hfile%observables, 'n_solvent', &
       n_solvent, h5md_linear, step=n_md_steps, &
       time=n_md_steps*dt)

  call h5gcreate_f(dimer_io%group, 'box', box_group, error)
  call h5md_write_attribute(box_group, 'dimension', 3)
  call dummy_element%create_fixed(box_group, 'edges', solvent_cells%edges)
  call h5gclose_f(box_group, error)

  call h5gcreate_f(solvent_io%group, 'box', box_group, error)
  call h5md_write_attribute(box_group, 'dimension', 3)
  call dummy_element%create_fixed(box_group, 'edges', solvent_cells%edges)
  call h5gclose_f(box_group, error)

  call solvent% random_placement(solvent_cells% edges, colloids, solvent_colloid_lj, state(1))

  call solvent% sort(solvent_cells)

  call neigh% init(colloids% Nmax, int(300*max(sigma_c,sigma_n)**3))

  skin = 1
  n_extra_sorting = 0

  call neigh% make_stencil(solvent_cells, max_cut+skin)

  call neigh% update_list(colloids, solvent, max_cut+skin, solvent_cells, solvent_colloid_lj)

  cyl_shell_rmin = max(sigma_n, sigma_c)
  cyl_shell_rmax = cyl_shell_rmin + 2
  call h5gcreate_f(hfile%id, 'fields', fields_group, error)
  call z_hist%init(-5.d0, 5.d0+d, 100)
  call z_hist_el%create_time(fields_group, 'cylindrical_shell_histogram', z_hist%data, &
       h5md_linear, step=n_md_steps, time=n_md_steps*dt)
  call h5md_write_attribute(z_hist_el%id, 'r_min', cyl_shell_rmin)
  call h5md_write_attribute(z_hist_el%id, 'r_max', cyl_shell_rmax)
  call h5md_write_attribute(z_hist_el%id, 'xmin', z_hist%xmin)
  call h5md_write_attribute(z_hist_el%id, 'dx', z_hist%dx)
  call h5gclose_f(fields_group, error)

  e1 = compute_force(colloids, solvent, neigh, solvent_cells% edges, solvent_colloid_lj)
  e2 = compute_force_n2(colloids, solvent_cells% edges, colloid_lj)
  solvent% force_old = solvent% force
  colloids% force_old = colloids% force

  call h5fflush_f(hfile%id, h5f_scope_global_f, error)

  write(*,*) 'Box size', l
  write(*,*) 'Fluid particles', n
  if (do_hydro) then
     write(*,*) 'MPCD tau', tau
  else
     write(*,*) 'RS resampling time', tau
  end if
  write(*,*) 'Dimer mass', mass

  write(*,*) 'Running for', equilibration_loops, '+', n_loop, 'loops'

  call main%tic()
  sampling = .false.
  do i = 0, n_loop+equilibration_loops
     if (i==equilibration_loops) sampling = .true.
     if (modulo(i,20) == 0) write(*,'(i05)',advance='no') i
     md_loop: do j = 1, n_md_steps
        call md_pos(solvent, dt)

        call varia%tic()
        ! Extra copy for rattle
        colloids% pos_rattle = colloids% pos
        do k=1, colloids% Nmax
           colloids% pos(:,k) = colloids% pos(:,k) + dt * colloids% vel(:,k) + &
                dt**2 * colloids% force(:,k) / (2 * colloids% mass(k))
        end do

        call rattle_dimer_pos(colloids, d, dt, solvent_cells% edges)
        call varia%tac()

        so_max = solvent% maximum_displacement()
        co_max = colloids% maximum_displacement()

        if ( (co_max >= skin*0.1d0) .or. (so_max >= skin*0.89d0) ) then
           call varia%tic()
           call apply_pbc(solvent, solvent_cells% edges)
           call apply_pbc(colloids, solvent_cells% edges)
           call varia%tac()
           call solvent% sort(solvent_cells)
           call neigh% update_list(colloids, solvent, max_cut + skin, solvent_cells, solvent_colloid_lj)
           call varia%tic()
           !$omp parallel do
           do k = 1, solvent%Nmax
              solvent% pos_old(:,k) = solvent% pos(:,k)
           end do
           colloids% pos_old = colloids% pos
           n_extra_sorting = n_extra_sorting + 1
           call varia%tac()
        end if

        call varia%tic()
        call switch(solvent% force, solvent% force_old)
        call switch(colloids% force, colloids% force_old)

        !$omp parallel do
        do k = 1, solvent%Nmax
           solvent% force(:,k) = 0
        end do
        colloids% force = 0
        call varia%tac()
        e1 = compute_force(colloids, solvent, neigh, solvent_cells% edges, solvent_colloid_lj)
        e2 = compute_force_n2(colloids, solvent_cells% edges, colloid_lj)

        call md_vel(solvent, dt)

        call varia%tic()
        do k=1, colloids% Nmax
           colloids% vel(:,k) = colloids% vel(:,k) + &
             dt * ( colloids% force(:,k) + colloids% force_old(:,k) ) / (2 * colloids% mass(k))
        end do

        call rattle_dimer_vel(colloids, d, dt, solvent_cells% edges)
        call varia%tac()

        if (sampling) then
           v_com = sum(colloids%vel, dim=2)/2
           unit_r = rel_pos(colloids%pos(:,1), colloids%pos(:,2), solvent_cells%edges)
           unit_r = unit_r / norm2(unit_r)
           call axial_cf%add_fast((i-equilibration_loops)*n_md_steps+j-1, v_com, unit_r)
        end if

        if ((sampling) .and. (modulo(j, colloid_sampling)==0)) then
           call dimer_io%position%append(colloids%pos)
           call dimer_io%velocity%append(colloids%vel)
           call dimer_io%image%append(colloids%image)
        end if

        call time_flag%tic()
        call flag_particles_nl
        call time_flag%tac()
        call time_change%tic()
        call change_species
        call time_change%tac()

     end do md_loop

     collision_step = (modulo(i, collide_every) == 0)

     call varia%tic()

     if (sampling) then
        temperature = compute_temperature(solvent, solvent_cells)
        kin_e = (colloids% mass(1)*sum(colloids% vel(:,1)**2) + &
             colloids% mass(2)*sum(colloids% vel(:,2)**2))/2 + &
             sum(solvent% vel**2)/2
        v_com = (sum(solvent% vel, dim=2) + &
             mass(1)*colloids%vel(:,1) + mass(2)*colloids%vel(:,2)) / &
             (solvent%Nmax + mass(1) + mass(2))

        call thermo_data%append(hfile, temperature, e1+e2, kin_e, e1+e2+kin_e, v_com)

        com_pos = ( colloids%pos(:,1)+colloids%image(:,1)*solvent_cells%edges + &
             colloids%pos(:,2)+colloids%image(:,2)*solvent_cells%edges)/2
        call axial_cf%add(i-equilibration_loops, com_pos, unit_r)

     end if

     if (collision_step) call solvent_cells%random_shift(state(1))
     call varia%tac()

     call apply_pbc(solvent, solvent_cells% edges)
     call apply_pbc(colloids, solvent_cells% edges)
     call solvent% sort(solvent_cells)
     call neigh% update_list(colloids, solvent, max_cut+skin, solvent_cells, solvent_colloid_lj)
     call varia%tic()
     !$omp parallel do
     do k = 1, solvent%Nmax
        solvent% pos_old(:,k) = solvent% pos(:,k)
     end do
     colloids% pos_old = colloids% pos
     call varia%tac()

     if (collision_step) call simple_mpcd_step(solvent, solvent_cells, state, &
          thermostat=do_thermostat, t=t, hydro=do_hydro)
     
     call time_refuel%tic()
     if (bulk_rmpcd) then
        call bulk_reaction(solvent, solvent_cells, 2, 1, bulk_rate, tau, state)
     else
        call refuel
     end if

     call time_refuel%tac()

     n_solvent = 0
     do k = 1, solvent%Nmax
        j = solvent%species(k)
        if (j <= 0) continue
        n_solvent(j) = n_solvent(j) + 1
     end do
     call n_solvent_el%append(n_solvent)

     z_hist%data = 0
     call compute_cylindrical_shell_histogram(z_hist, colloids%pos(:,1), colloids%pos(:,2), &
          solvent_cells%edges, 2, cyl_shell_rmin, cyl_shell_rmax, solvent)
     call z_hist_el%append(z_hist%data)

  end do
  call main%tac()
  write(*,*) ''

  write(*,*) 'n extra sorting', n_extra_sorting

  ! create a group for block correlators and write the data

  call h5gcreate_f(hfile%id, 'block_correlators', correlator_group, error)
  call axial_cf%write(correlator_group, n_md_steps, n_md_steps*dt, 1, dt)
  call h5gclose_f(correlator_group, error)

  ! write solvent coordinates for last step

  call solvent_io%position%append(solvent%pos)
  call solvent_io%velocity%append(solvent%vel)
  call solvent_io%image%append(solvent%image)
  call solvent_io%species%append(solvent%species)

  ! Store timing data
  call timer_list%append(solvent%time_stream)
  call timer_list%append(solvent%time_step)
  call timer_list%append(solvent%time_count)
  call timer_list%append(solvent%time_sort)
  call timer_list%append(solvent%time_ct)
  call timer_list%append(solvent%time_md_pos)
  call timer_list%append(solvent%time_md_vel)
  call timer_list%append(solvent%time_max_disp)
  call timer_list%append(colloids%time_self_force)
  call timer_list%append(neigh%time_update)
  call timer_list%append(neigh%time_force)
  call timer_list%append(colloids%time_max_disp)

  call h5gcreate_f(hfile%id, 'timers', timers_group, error)
  call timer_list%write(timers_group, total_time)
  call h5md_write_dataset(timers_group, 'total', total_time)
  call h5md_write_dataset(timers_group, main%name, main%total)
  call h5gclose_f(timers_group, error)

  call dimer_io%close()
  call solvent_io%close()
  call z_hist_el%close()
  call n_solvent_el%close()
  call hfile%close()
  call h5close_f(error)

contains

  subroutine flag_particles_nl
    double precision :: dist_to_C_sq
    integer :: r, s
    double precision :: x(3)

    do s = 1,neigh% n(1)
       r = neigh%list(s,1)
       if (solvent% species(r) == 1) then
          x = rel_pos(colloids% pos(:,1),solvent% pos(:,r),solvent_cells% edges)
          dist_to_c_sq = dot_product(x, x)
          if (dist_to_c_sq < solvent_colloid_lj%cut_sq(1,1)) then
             solvent% flags(r) = ibset(solvent% flags(r), reac_bit)
          end if
       end if
    end do

  end subroutine flag_particles_nl
  
  subroutine change_species
    double precision :: dist_to_C_sq
    double precision :: dist_to_N_sq
    integer :: m
    double precision :: x(3)
    integer :: thread_id

    !$omp parallel private(thread_id)
    thread_id = omp_get_thread_num() + 1
    !$omp do private(x, dist_to_C_sq, dist_to_N_sq)
    do m = 1, solvent% Nmax
       if (btest(solvent% flags(m), reac_bit)) then
          x = rel_pos(colloids% pos(:,1), solvent% pos(:,m), solvent_cells% edges)
          dist_to_c_sq = dot_product(x, x)
          x = rel_pos(colloids% pos(:,2), solvent% pos(:,m), solvent_cells% edges)
          dist_to_n_sq = dot_product(x, x)
          if ( &
               (dist_to_c_sq > solvent_colloid_lj%cut_sq(1,1)) &
               .and. &
               (dist_to_n_sq > solvent_colloid_lj%cut_sq(1,2)) &
               ) &
               then
             if (threefry_double(state(thread_id)) <= prob) then
                solvent% species(m) = 2
             end if
             solvent%flags(m) = ibclr(solvent%flags(m), reac_bit)
          end if
       end if
    end do
    !$omp end do
    !$omp end parallel
  end subroutine change_species
  
  subroutine refuel
    double precision :: dist_to_C_sq
    double precision :: dist_to_N_sq
    double precision :: far
    double precision :: x(3)
    integer :: n

    far = (l(1)*0.45)**2

    !$omp parallel do private(x, dist_to_C_sq, dist_to_N_sq)
    do n = 1,solvent% Nmax
       if (solvent% species(n) == 2) then
          x = rel_pos(colloids% pos(:,1), solvent% pos(:,n), solvent_cells% edges)
          dist_to_c_sq = dot_product(x, x)
          x= rel_pos(colloids% pos(:,2), solvent% pos(:,n), solvent_cells% edges)
          dist_to_n_sq = dot_product(x, x)
          if ((dist_to_c_sq > far) .and. (dist_to_n_sq > far)) then
             solvent% species(n) = 1
          end if
       end if
    end do
  end subroutine refuel

end program single_dimer_pbc