Homec4science
Diffusion lammps (lammm-master) Browse

lammps/examples/gaussFlow/in.GD
95755ba8b2e9lammm-master

Last Change
Raw File

in.GD
View Options

#LAMMPS input script
#in.GD
#see README for details
###############################################################################
#initialize variables
clear
#frequency for outputting info (timesteps)
variable dump_rate equal 50
variable thermo_rate equal 10
#equilibration time (timesteps)
variable equil equal 1000
#stabilization time (timesteps to reach steady-state)
variable stabil equal 1000
#data collection time (timesteps)
variable run equal 2000
#length of pipe
variable L equal 30
#width of pipe
variable d equal 20
#flux (mass/sigma*tau)
variable J equal 0.1
#simulation box dimensions
variable Lx equal 100
variable Ly equal 40
#bulk fluid density
variable dens equal 0.8
#lattice spacing for wall atoms
variable aWall equal 1.0 #1.7472
#timestep
variable ts equal 0.001
#temperature
variable T equal 2.0
#thermostat damping constant
variable tdamp equal ${ts}*100
units lj
dimension 2
atom_style atomic
###############################################################################
#create box
#create lattice with the spacing aWall
variable rhoWall equal ${aWall}^(-2)
lattice sq ${rhoWall}
#modify input dimensions to be multiples of aWall
variable L1 equal round($L/${aWall})*${aWall}
variable d1 equal round($d/${aWall})*${aWall}
variable Ly1 equal round(${Ly}/${aWall})*${aWall}
variable Lx1 equal round(${Lx}/${aWall})*${aWall}
#create simulation box
variable lx2 equal ${Lx1}/2
variable ly2 equal ${Ly1}/2
region simbox block -${lx2} ${lx2} -${ly2} ${ly2} 0 0.1 units box
create_box 2 simbox
#####################################################################
#set up potential
mass 1 1.0 #fluid atoms
mass 2 1.0 #wall atoms
pair_style lj/cut 2.5
pair_modify shift yes
pair_coeff 1 1 1.0 1.0 2.5
pair_coeff 1 2 1.0 1.0 1.12246
pair_coeff 2 2 0.0 0.0 0.0
timestep ${ts}
#####################################################################
#create atoms
#create wall atoms everywhere
create_atoms 2 box
#define region which is "walled off"
variable dhalf equal ${d1}/2
variable Lhalf equal ${L1}/2
region walltop block -${Lhalf} ${Lhalf} ${dhalf} EDGE -0.1 0.1 &
units box
region wallbot block -${Lhalf} ${Lhalf} EDGE -${dhalf} -0.1 0.1 &
units box
region outsidewall union 2 walltop wallbot side out
#remove wall atoms outside wall region
group outside region outsidewall
delete_atoms group outside
#remove wall atoms that aren't on edge of wall region
variable x1 equal ${Lhalf}-${aWall}
variable y1 equal ${dhalf}+${aWall}
region insideTop block -${x1} ${x1} ${y1} EDGE -0.1 0.1 units box
region insideBot block -${x1} ${x1} EDGE -${y1} -0.1 0.1 units box
region insideWall union 2 insideTop insideBot
group insideWall region insideWall
delete_atoms group insideWall
#define new lattice, to give correct fluid density
#y lattice const must be a multiple of aWall
variable atrue equal ${dens}^(-1/2)
variable ay equal round(${atrue}/${aWall})*${aWall}
#choose x lattice const to give correct density
variable ax equal (${ay}*${dens})^(-1)
#change Lx to be multiple of ax
variable Lx1 equal round(${Lx}/${ax})*${ax}
variable lx2 equal ${Lx1}/2
change_box all x final -${lx2} ${lx2} units box
#define new lattice
lattice custom ${dens} &
a1 ${ax} 0.0 0.0 a2 0.0 ${ay} 0.0 a3 0.0 0.0 1.0 &
basis 0.0 0.0 0.0
#fill in rest of box with bulk particles
variable delta equal 0.001
variable Ldelt equal ${Lhalf}+${delta}
variable dDelt equal ${dhalf}-${delta}
region left block EDGE -${Ldelt} EDGE EDGE -0.1 0.1 units box
region right block ${Ldelt} EDGE EDGE EDGE -0.1 0.1 units box
region pipe block -${Ldelt} ${Ldelt} -${dDelt} ${dDelt} -0.1 0.1 &
units box
region bulk union 3 left pipe right
create_atoms 1 region bulk
group bulk type 1
group wall type 2
#remove atoms that are too close to wall
delete_atoms overlap 0.9 bulk wall
neighbor 0.3 bin
neigh_modify delay 0 every 1 check yes
neigh_modify exclude group wall wall
velocity bulk create $T 78915 dist gaussian rot yes mom yes loop geom
#####################################################################
#set up PUT
#see Evans and Morriss, Phys. Rev. Lett. 56(20) 1986, p. 2172
#average number of particles per box, Evans and Morriss used 2.0
variable NperBox equal 8.0
#calculate box sizes
variable boxSide equal sqrt(${NperBox}/${dens})
variable nX equal round(lx/${boxSide})
variable nY equal round(ly/${boxSide})
variable dX equal lx/${nX}
variable dY equal ly/${nY}
#temperature of fluid (excluding wall)
compute myT bulk temp
#profile-unbiased temperature of fluid
compute myTp bulk temp/profile 1 1 0 xy ${nX} ${nY}
#thermo setup
thermo ${thermo_rate}
thermo_style custom step c_myT c_myTp etotal press
#dump initial configuration
dump 55 all custom 1 all.init.lammpstrj id type x y z vx vy vz
dump 56 wall custom 1 wall.init.lammpstrj id type x y z
dump_modify 55 sort id
dump_modify 56 sort id
run 0
undump 55
undump 56
#####################################################################
#equilibrate without GD
fix nvt bulk nvt temp $T $T ${tdamp}
fix_modify nvt temp myTp
fix 2 bulk enforce2d
run ${equil}
#####################################################################
#initialize the COM velocity and run to achieve steady-state
#calculate velocity to add: V=J/rho_total
variable Vadd equal $J*lx*ly/count(bulk)
#first remove any COM velocity, then add back the streaming velocity
velocity bulk zero linear
velocity bulk set ${Vadd} 0.0 0.0 units box sum yes mom no
fix GD bulk gaussFlow 1 0 0 #energy yes
#fix_modify GD energy yes
run ${stabil}
#####################################################################
#collect data
#print the applied force and total flux to ensure conservation of Jx
variable Fapp equal f_GD[1]
compute vxBulk bulk reduce sum vx
compute vyBulk bulk reduce sum vy
variable invVol equal 1.0/(lx*ly)
variable jx equal c_vxBulk*${invVol}
variable jy equal c_vyBulk*${invVol}
variable curr_step equal step
fix print_vCOM all print ${dump_rate} &
"${curr_step} ${Fapp} ${jx} ${jy}" file GD.out screen no &
title "timestep Fapp Jx Jy"
#compute IK1 pressure profile
#see Todd, Evans, and Davis, Phys. Rev. E 52(2) 1995, p. 1627
#use profile-unbiased temperature to remove the streaming velocity
#from the kinetic part of the pressure
compute spa bulk stress/atom myTp
#for the pressure profile, use the same grid as the PUT
compute chunkX bulk chunk/atom bin/1d x lower ${dX} units box
#output pressure profile and other profiles
#the pressure profile is (-1/2V)*(c_spa[1] + c_spa[2]), where
#V is the volume of a slice
fix profiles bulk ave/chunk 1 1 ${dump_rate} chunkX &
vx density/mass c_spa[1] c_spa[2] &
file x_profiles ave running overwrite
#compute velocity profile across the pipe with a finer grid
variable dYnew equal ${dY}/10
compute chunkY bulk chunk/atom bin/1d y center ${dYnew} units box &
region pipe
fix velYprof bulk ave/chunk 1 1 ${dump_rate} chunkY &
vx file Vy_profile ave running overwrite
#full trajectory
dump 7 bulk custom ${dump_rate} bulk.lammpstrj &
id type x y z
dump_modify 7 sort id
run ${run}
c4science · Help