LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
[lammps-users] Pressure calculated using temperature compute that removes bias results in lower average pressure than measured at equilibrium
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[lammps-users] Pressure calculated using temperature compute that removes bias results in lower average pressure than measured at equilibrium


From: Daniel Rankin <a1608247@...3703...>
Date: Wed, 22 Nov 2017 12:36:06 +1030

Dear LAMMPS users

I am simulating a pressure-driven flow of a 2D Lennard-Jones fluid through a pore in between two membranes using a similar script to the one in examples/USER/misc/flow_gauss in the lammps directory, but using the pump method described in Zhu, F., E. Tajkhorshid, and K. Schulten. 2004. Biophys. J. 86:50–57. I have measured the pressure profile in the bulk regions of the reservoirs using the "compute stress/atom" command before and after applying a pressure difference. When there is no pressure difference applied to the system I use the "compute temp" command to calculate the temperature of the fluid and then use this compute ID as a temp-ID in the "compute stress/atom" command, which I use to calculate the pressure profile. When I apply a pressure difference using the pump method (see the reference above) I use a temperature compute that removes a centre of mass velocity field from the particles before computing the temperature ("compute temp/profile"). Using this temperature compute in "compute stress/atom" results in a lower average pressure (averaging over both reservoirs) compared with the equilibrium pressure profile. The difference is more severe when I measure the pressure in the reservoirs of a 3D Lennard-Jones system for a different type of simulation but using the same "compute temp/profile" command.

My questions are:

1. Do I need to correct the fluid pressure calculated using "compute stress/atom" for the extra degrees of freedom resulting from using "compute temp/profile"?
2. How do I do this?

I am using the lammps version from 11Aug17.

Thanks, Daniel.
#LAMMPS input script
#Based on in.GD in LAMMPS example folder (USER/misc/flow_gauss)
#Daniel Rankin November 2017

###############################################################################
#initialize variables
clear

#frequency for outputting info (timesteps)
variable        dump_rate       equal 50
variable        thermo_rate     equal 10

#equilibration time (timesteps)
variable        equil           equal 1000000

#stabilization time (timesteps to reach steady-state)
variable        stabil          equal 1000000

#data collection time (timesteps)
variable        run             equal 1000000

#length of pipe
variable        L               equal 30

#width of pipe
variable        d               equal 20

#flux (mass/sigma*tau)
variable        dp               equal 0.05

#simulation box dimensions
variable        Lx              equal 400
variable        Ly              equal 60

#total width (sigma) of region where forces are applied to fluid
variable	fdist		equal 2

#distance from start of simulation cell to apply forces
#to fluid particles
variable	LL		equal -0.5*(${Lx}-${fdist})

#distance from end of simulation cell to apply forces
#to fluid particles
variable	LR		equal 0.5*(${Lx}-${fdist})

#coordinates for determining bulk regions of reservoirs
variable	resL1		equal ${LL}+10
variable	resL2		equal -0.5*$L-15
variable	resR1		equal 0.5*$L+15
variable	resR2		equal ${LR}-10

#bulk fluid density
variable        dens            equal 0.72

#lattice spacing for wall atoms
variable        aWall           equal 1.0 #1.7472

#timestep
variable        ts              equal 0.005

#temperature
variable        T               equal 0.5

#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 4.0
pair_coeff      1 1 1.0 1.0
pair_coeff      1 2 1.0 1.0
pair_coeff      2 2 1.0 1.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.6 bin
neigh_modify    delay 10

velocity        bulk create $T 8858125 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
thermo_modify 	temp myT

#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 pump method

fix             nvt bulk nvt temp $T $T ${tdamp}
fix_modify      nvt temp myTp
fix             2 bulk enforce2d

#for the pressure profile, use the same grid as the PUT
compute         chunkX bulk chunk/atom bin/1d x lower ${dX} units box

#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     spaeq bulk stress/atom myT

#set up bulk regions of reservoirs
region		resL block ${resL1} ${resL2} EDGE EDGE -0.1 0.1 units box
region		resR block ${resR1} ${resR2} EDGE EDGE -0.1 0.1 units box

#bulk reservoir area
variable 	res_bulk_area equal (${resL2}-${resL1})*ly

#group particles dynamically into reservoirs
group		gresL dynamic bulk region resL every 1
group		gresR dynamic bulk region resR every 1

#calculate pressures drop across reservoirs
compute		pLeq gresL reduce sum c_spaeq[1] c_spaeq[2]
compute		pReq gresR reduce sum c_spaeq[1] c_spaeq[2]
variable    pressLeq equal -(c_pLeq[1]+c_pLeq[2])/(2*v_res_bulk_area)
variable    pressReq equal -(c_pReq[1]+c_pReq[2])/(2*v_res_bulk_area)
variable	dpeq equal v_pressReq-v_pressLeq
fix			pressLeq bulk ave/time 10 1000 10000 v_pressLeq &
			file pressLeq.dat
fix			pressReq bulk ave/time 10 1000 10000 v_pressReq &
			file pressReq.dat
fix			dpeq bulk ave/time 10 1000 10000 v_dpeq file dpeq.dat

#variable for calculating pressure profile
variable	peq atom -(c_spaeq[1]+c_spaeq[2])/(2*${dX}*${Ly})

#output pressure profile
#the pressure profile is (-1/2V)*(c_spa[1] + c_spa[2]), where
#V is the volume of a slice

fix				peq bulk ave/chunk 1 1 ${dump_rate} chunkX &
				v_peq norm none file peq.dat ave running overwrite

run             ${equil}

unfix		pressLeq
unfix		pressReq
unfix		dpeq
unfix		peq

#####################################################################
#run to achieve stabilisation

#apply pressure gradient via pump method

#apply forces to fluid within regions of width 1 sigma at the ends of
#the simulation cell
region		rforce1 block EDGE ${LL} EDGE EDGE -0.1 0.1 units box
region		rforce2 block ${LR} EDGE EDGE EDGE -0.1 0.1 units box
region		rforce union 2 rforce1 rforce2
group		Fbulk dynamic bulk region rforce every 1 
run				1
variable 	fapp equal v_dp*ly/count(Fbulk)
fix			Fapp Fbulk addforce ${fapp} 0.0 0.0

#thermo setup
thermo          ${thermo_rate}
thermo_style    custom step c_myT c_myTp etotal press
thermo_modify 	temp myTp

run             ${stabil}

#####################################################################
#collect data

#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

#calculate pressures drop across reservoirs
compute		pL gresL reduce sum c_spa[1] c_spa[2]
compute		pR gresR reduce sum c_spa[1] c_spa[2]
variable    pressL equal -(c_pL[1]+c_pL[2])/(2*v_res_bulk_area)
variable    pressR equal -(c_pR[1]+c_pR[2])/(2*v_res_bulk_area)
variable	dp equal v_pressR-v_pressL
fix			pressL bulk ave/time 10 1000 10000 v_pressL file pressL.dat
fix			pressR bulk ave/time 10 1000 10000 v_pressR file pressR.dat
fix			dp bulk ave/time 10 1000 10000 v_dp file dp.dat

#variable for calculating pressure profile
variable	p atom -(c_spa[1]+c_spa[2])/(2*${dX}*${Ly})

#profiles

fix             dens bulk ave/chunk 1 1 ${dump_rate} chunkX &
				density/mass file dens.dat ave running overwrite
                
#output pressure profile
#the pressure profile is (-1/2V)*(c_spa[1] + c_spa[2]), where
#V is the volume of a slice

fix				p bulk ave/chunk 1 1 ${dump_rate} chunkX &
				v_p norm none file p.dat 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 -${dDelt} ${dYnew} &
				bound x -${Ldelt} ${Ldelt} bound y -${dDelt} ${dDelt} &
				units box region pipe
fix             velYprof bulk ave/chunk 1 1 ${dump_rate} chunkY &
                vx file Vy_profile.dat ave running overwrite              

#trajectory
dump            atomdump  all atom 10000 lj.lammpstrj

run             ${run}