# DPD example
# counter-flowing Poiseuille flow
# [1] Backer, J.A., Lowe, C.P., Hoefsloot, H.C.J. & Iedema,
# P.D. J. Chem. Phys. 122, 154503 (2005).
# doi:10.1063/1.1883163
# [2] Espanol, P. & Warren, P. Europhys. Lett. 30, 191-196 (1995).
# doi:10.1209/0295-5075/30/4/001
# 
# 2D version of simulation in [1]
# the simulation produces spatially averages 
# vy.av (velocity), rho.av (density), simtemp.av (temperature), sxy.av (shear stress)
units		si
variable       ndim  equal 2
variable number_density equal 6
# domain size 
variable        xsize equal 12
variable        ysize equal 8
# temperature 
variable        kb equal 1.3806488e-23
variable        T equal 0.5/${kb}
# interaction parameters
variable        cutoff equal 1.0
variable        sigma equal 4.5
# gamma is not given in [1] 
# we get it using eq. (10) in [2]
variable        gamma equal ${sigma}^2/(2*${T}*${kb})
# body force (units of acceleration)
variable gy equal 0.055
# number of timesteps
variable ntime equal 40000
timestep	0.01

dimension       ${ndim}
atom_style	atomic
communicate	single vel yes

# lattice parameter
variable lsp    equal 1.0/${number_density}^(1.0/${ndim})
lattice sq      ${lsp} origin 0.5 0.5 0.0
variable        Npart equal ${xsize}*${ysize}*${number_density}

region		box block 0 ${xsize} 0 ${ysize} -0.1 0.1 units box
create_box	1 box
create_atoms	1 box

# TODO: find the way to generate the right number of atoms
# delete extra atoms
group extra_atoms id > ${Npart}
delete_atoms group extra_atoms

mass		1 1.0
velocity	all create ${T} 87287 loop geom

pair_style	dpd ${T} ${cutoff} 928948
pair_coeff	1 1 ${gamma} ${sigma}

neighbor	0.5 bin
neigh_modify    delay 0 every 1

fix		1 all nve

variable        Nfreq   equal  ${ntime}
variable        Nrepeat equal  round(0.9*${ntime})
fix av_vy all ave/spatial 1 ${Nrepeat} ${Nfreq} x center 0.2 vy file "vy.av"
fix av_rho all ave/spatial 1 ${Nrepeat} ${Nfreq} x center 0.5 density/number file "rho.av"

# simulation temperature (do not use vy, flow goes in this direction)
variable simtemp atom mass*vx^2
fix av_temp all ave/spatial 1 ${Nrepeat} ${Nfreq} x lower 0.5 v_simtemp file "simtemp.av"

# stresses are in units of pressure*volume must be divided by per-atom
# volume to have units of stress (pressue); components of the stress
# are in the following order xx(1), yy(2), zz(3), xy(4), xz(5), yz(6)
compute stress all stress/atom
variable stress_pressure atom c_stress[4]*${number_density}
fix av_xz_stress all ave/spatial 1 ${Nrepeat} ${Nfreq} x lower 0.5 v_stress_pressure file "sxy.av"

dump mdump all custom 10000 dump.dpd id type xs ys vx vy
variable fy atom mass*${gy}*((x>${xsize}/2.0)-(x<${xsize}/2.0))
fix reverce_periodic all addforce 0.0 v_fy 0.0 

run		${ntime}