forked from lijiext/lammps
133 lines
4.0 KiB
Plaintext
133 lines
4.0 KiB
Plaintext
|
###############################################################################
|
||
|
|
#
|
||
|
|
#
|
||
|
|
# This input script is a modified version of the example script lj_ehex.lmp
|
||
|
|
# which is part of the supplementary (open access) material of the paper
|
||
|
|
#
|
||
|
|
# P. Wirnsberger, D. Frenkel and C. Dellago,
|
||
|
|
# "An enhanced version of the heat exchange algorithm with excellent energy
|
||
|
|
# conservation properties", J. Chem. Phys. 143, 124104 (2015).
|
||
|
|
#
|
||
|
|
# The full article is available on arXiv: http://arxiv.org/pdf/1507.07081.
|
||
|
|
#
|
||
|
|
#
|
||
|
|
# Description:
|
||
|
|
# ------------
|
||
|
|
#
|
||
|
|
# This file is a LAMMPS input script for carrying out a NEMD simulation of
|
||
|
|
# Lennard-Jones using the eHEX/a algorithm. The run produces two files:
|
||
|
|
# "out.Tlj_ehex" contains the temperature profile and "out.Elj_ehex" the time
|
||
|
|
# evolution of the total energy.
|
||
|
|
#
|
||
|
|
###############################################################################
|
||
|
|
|
||
|
|
log log.lj_ehex
|
||
|
|
|
||
|
|
# heat flux
|
||
|
|
variable J equal 0.15
|
||
|
|
|
||
|
|
# timestep
|
||
|
|
variable dt equal 0.007
|
||
|
|
|
||
|
|
# cutoff radius for shifted LJ-potential
|
||
|
|
variable rc equal 3.0
|
||
|
|
|
||
|
|
# simulation time for the production run
|
||
|
|
variable tprod equal 5000
|
||
|
|
|
||
|
|
# total number of timesteps
|
||
|
|
variable Nprod equal floor(${tprod}/${dt})
|
||
|
|
|
||
|
|
# equilibrated steady state configuration
|
||
|
|
read_data "data.lj"
|
||
|
|
|
||
|
|
# use LJ shifted force pair style
|
||
|
|
pair_style lj/sf ${rc}
|
||
|
|
|
||
|
|
# with coefficients eps = 1, sigma = 1, and rc = 3.0
|
||
|
|
pair_coeff 1 1 1.0 1.0 ${rc}
|
||
|
|
|
||
|
|
# increase neigbor skin because of the large timestep
|
||
|
|
neighbor 0.8 bin
|
||
|
|
|
||
|
|
# options used for fix ave/time; sample the quantities every 10 steps
|
||
|
|
variable Nsamp equal 10
|
||
|
|
variable Nrepeat equal floor(${Nprod}/${Nsamp})
|
||
|
|
variable Nevery equal ${Nsamp}*${Nrepeat}
|
||
|
|
|
||
|
|
# box dimensions
|
||
|
|
variable Lz equal zhi-zlo
|
||
|
|
variable Lx equal xhi-xlo
|
||
|
|
variable Ly equal yhi-ylo
|
||
|
|
|
||
|
|
# reservoir width in z-direction
|
||
|
|
variable delta equal 2.
|
||
|
|
|
||
|
|
# specify z-extents of both reservoirs
|
||
|
|
variable zlo_Thi equal -${Lz}/4.-${delta}/2.
|
||
|
|
variable zhi_Thi equal ${zlo_Thi}+${delta}
|
||
|
|
variable zlo_Tlo equal ${zlo_Thi}+${Lz}/2.
|
||
|
|
variable zhi_Tlo equal ${zlo_Tlo}+${delta}
|
||
|
|
|
||
|
|
# resolution for fix ave/spatial
|
||
|
|
variable dz equal ${Lz}/60
|
||
|
|
|
||
|
|
# compute per-atom kinetic energy and temperature, respectively
|
||
|
|
# NOTE: In this example we ignored the centre of mass (com) velocities
|
||
|
|
# of the individual bins for simplicity. However, we took that
|
||
|
|
# into account for the publication.
|
||
|
|
compute ke all ke/atom
|
||
|
|
variable T atom c_ke/1.5
|
||
|
|
|
||
|
|
# specify the reservoirs
|
||
|
|
region Thi_region block INF INF INF INF ${zlo_Thi} ${zhi_Thi}
|
||
|
|
region Tlo_region block INF INF INF INF ${zlo_Tlo} ${zhi_Tlo}
|
||
|
|
|
||
|
|
# compute the temperature of the individual region
|
||
|
|
compute cTlo all temp/region Tlo_region
|
||
|
|
compute cThi all temp/region Thi_region
|
||
|
|
|
||
|
|
# calculate the energy flux from the specified heat flux
|
||
|
|
variable F equal ${J}*${Lx}*${Ly}*2.
|
||
|
|
|
||
|
|
# use fix ehex to create the gradient
|
||
|
|
# hot reservoir
|
||
|
|
fix fHi all ehex 1 +${F} region Thi_region
|
||
|
|
|
||
|
|
# cold reservoir
|
||
|
|
fix fLo all ehex 1 -${F} region Tlo_region
|
||
|
|
|
||
|
|
# use velocity Verlet for integration
|
||
|
|
fix fNVEGrad all nve
|
||
|
|
|
||
|
|
# calculate the centre of mass velocity of the entire box (vcmx, vcmy, vcmz)
|
||
|
|
variable vcmx equal "vcm(all,x)"
|
||
|
|
variable vcmy equal "vcm(all,y)"
|
||
|
|
variable vcmz equal "vcm(all,z)"
|
||
|
|
variable vcm2 equal v_vcmx*v_vcmx+v_vcmy*v_vcmy+v_vcmz*v_vcmz
|
||
|
|
|
||
|
|
# specify the timestep
|
||
|
|
timestep ${dt}
|
||
|
|
|
||
|
|
# frequency for console output
|
||
|
|
thermo 10000
|
||
|
|
|
||
|
|
# print timestep, temperature, total energy and v_com^2 to console
|
||
|
|
thermo_style custom step temp etotal v_vcm2
|
||
|
|
|
||
|
|
# calculate spatial average of temperature
|
||
|
|
compute cchT all chunk/atom bin/1d z lower ${dz}
|
||
|
|
fix fchT all ave/chunk ${Nsamp} ${Nrepeat} ${Nevery} cchT v_T file out.Tlj_ehex
|
||
|
|
|
||
|
|
# compute the total energy
|
||
|
|
compute cKe all ke
|
||
|
|
compute cPe all pe
|
||
|
|
variable E equal c_cKe+c_cPe
|
||
|
|
|
||
|
|
# track the time evolution of the total energy
|
||
|
|
fix fE all ave/time ${Nsamp} 1000 10000 v_E file out.Elj_ehex
|
||
|
|
|
||
|
|
# production run
|
||
|
|
run ${Nprod}
|
||
|
|
|