# ====================================================================
# DPD thermostated liquid of LJ chains sheared between (111) fcc walls
# ====================================================================

# Couette flow of a DPD thermostated (1) liquid of LJ chains (2) confined
# between (111) fcc walls. The liquid viscosity and the liquid/solid slip
# length (3) can be computed from the recorded shear stress and velocity
# profile.
# Initial configuration (data.shear) is generated by the perl script chain.pl.
# 3 different shear velocities (U = +-0.1, +-0.2, +-0.3) are applied
# successively, to check the linear response of the system.

# This script can be used to test the influence of various parameters on the
# viscosity, in particular gamma coefficient of the DPD thermostat, chain
# length ($nchain in chain.pl), and on the slip length, in particular epsilon
# coefficient of the LJ potential between liquid and solid atoms.

# Detailed analysis procedure: Once the equilibrium confinement height has
# been established from rho.shear (denoting h the half-distance between the
# liquid/solid interfaces), viscosity can be computed as the ratio between the
# shear stress sigma_xy (from sigma.shear) and the shear rate gammadot (from
# the velocity profiles vx.U*). The slip length b can then be computed as:
# b = U/gammadot - h.

# (1) http://lammps.sandia.gov/doc/pair_dpd.html
# (2) http://link.aip.org/link/doi/10.1063/1.3469860
# (3) http://pubs.rsc.org/en/content/articlelanding/2007/sm/b616490k


# =======================================
# gamma coefficient of the DPD thermostat
# =======================================
variable gamma equal 1.0

# ======================================================================
# epsilon coefficient of the LJ potential between liquid and solid atoms
# ======================================================================
variable epsilon equal 0.5

# ========================
# temperature and pressure
# ========================
variable T equal 0.73
variable p equal 0.007

units lj
boundary p f p
communicate single vel yes

atom_style bond
pair_style hybrid/overlay lj/cut/opt 2.5 dpd/tstat $T $T 2.5 83628
bond_style harmonic
special_bonds fene

read_data data.shear

mass * 1.0
pair_coeff 1 1 lj/cut/opt 1.0 1.0
pair_coeff 1 2 lj/cut/opt ${epsilon} 1.0
pair_coeff 2 2 none
pair_coeff 1 1 dpd/tstat ${gamma}
bond_coeff 1 1500.0 1.0
neigh_modify delay 5

group liq type 1
group wall type 2
variable ymiddle equal 0.5*(yhi+ylo)
region hi block INF INF ${ymiddle} INF INF INF units box
group hi region hi
group wallhi intersect wall hi
group walllo subtract wall wallhi

# ============================================================
# upper wall is used as a piston to impose the liquid pressure
# ============================================================
fix setfhi wallhi setforce 0.0 NULL 0.0
variable Sy equal lx*lz
variable fy equal $p*${Sy}/count(wallhi)
fix avefhi wallhi aveforce NULL -${fy} NULL
fix setflo walllo setforce 0.0 0.0 0.0
fix nveall all nve

# =================================================
# output density, temperature and velocity profiles
# =================================================
variable Tatom atom 0.5*mass*(vy*vy+vz*vz)
fix 1 all ave/spatial 10 10000 200000 y lower 0.1 density/number v_Tatom units box file rho.eq
fix 2 all ave/spatial 10 80000 1200000 y lower 0.1 density/number v_Tatom units box file rho.U0.1
fix 3 all ave/spatial 10 80000 2200000 y lower 0.1 density/number v_Tatom units box file rho.U0.2
fix 4 all ave/spatial 10 80000 3200000 y lower 0.1 density/number v_Tatom units box file rho.U0.3
fix 5 all ave/spatial 10 80000 1200000 y lower 1.0 vx units box file vx.U0.1
fix 6 all ave/spatial 10 80000 2200000 y lower 1.0 vx units box file vx.U0.2
fix 7 all ave/spatial 10 80000 3200000 y lower 1.0 vx units box file vx.U0.3

# =============================================
# output pressure and shear stress on the walls
# =============================================
variable shi equal -f_setfhi[1]/${Sy}
variable phi equal f_avefhi[2]/${Sy}
variable slo equal f_setflo[1]/${Sy}
variable plo equal -f_setflo[2]/${Sy}
fix sigma all ave/time 10 1000 10000 v_phi v_plo v_shi v_slo file sigma.shear

# =====================================================
# dump atom positions for visualization, e.g. using VMD
# =====================================================
#dump 1 all atom 10000 dump.shear.lammpstrj
#dump_modify 1 flush yes
# =====================================
# dump images of the atom configuration
# =====================================
#dump 2 all image 4000 dump.*.ppm type type view 0.0 0.0 zoom 1.6
#dump_modify 2 flush yes

# ==========================================================
# output liquid temperature (computed only in the directions
# perpendicular to the flow) and position of the upper wall
# ==========================================================
compute Tliq liq temp/partial 0 1 1
variable ywall equal xcm(wallhi,y)
thermo_style custom step cpu c_Tliq v_ywall etotal
thermo 1000
thermo_modify flush yes

timestep 0.003  # cf. stiff bonds 
velocity liq create $T 87654 dist gaussian mom yes rot yes
run 200000
unfix 1

variable U equal 0.1
velocity wallhi set $U NULL NULL units box
velocity walllo set -$U NULL NULL units box
run 1000000
unfix 2
unfix 5

variable U equal 0.2
velocity wallhi set $U NULL NULL units box
velocity walllo set -$U NULL NULL units box
run 1000000
unfix 3
unfix 6

variable U equal 0.3
velocity wallhi set $U NULL NULL units box
velocity walllo set -$U NULL NULL units box
run 1000000
unfix 4
unfix 7