forked from lijiext/lammps
99 lines
3.4 KiB
Plaintext
99 lines
3.4 KiB
Plaintext
# Ar in real units
|
|
|
|
# simulation params in reduced units
|
|
# settable from command line
|
|
# epsilon, sigma, mass set below
|
|
|
|
variable x index 5
|
|
variable y index 5
|
|
variable z index 5
|
|
variable rhostar index 0.8842
|
|
variable dt index 0.005
|
|
variable cutoff index 2.5
|
|
variable skin index 0.3
|
|
variable tinitial index 1.0
|
|
variable nthermo index 10
|
|
variable nsteps index 100
|
|
|
|
# physical constants from update.cpp
|
|
|
|
variable kb index 0.0019872067 # kB in Kcal/mole/K
|
|
variable avogadro index 6.02214129e23 # Avogadro's number
|
|
|
|
# Ar properties in real units
|
|
|
|
variable epskb index 117.7 # LJ epsilon/kB in degrees K
|
|
variable sigma index 3.504 # LJ sigma in Angstroms
|
|
variable epsilon equal ${epskb}*${kb} # LJ epsilon in Kcal/mole
|
|
variable mass index 39.95 # mass in g/mole
|
|
|
|
# scale factors
|
|
|
|
# sigma = scale factor on distance, converts reduced distance to Angs
|
|
# epsilon = scale factor on energy, converts reduced energy to Kcal/mole
|
|
# tmpscale = scale factor on temperature, converts reduced temp to degrees K
|
|
# tscale = scale factor on time, converts reduced time to fs
|
|
# formula is t = t* / sqrt(epsilon/mass/sigma^2), but need t in fs
|
|
# use epsilon (Joule/mole), mass (kg/mole), sigma (meter) to get t in seconds
|
|
# pscale = scale factor on pressure, converts reduced pressure to atmospheres
|
|
# formula is P = P* / (sigma^3/epsilon), but need P in atmospheres
|
|
# use sigma (meter), epsilon (Joule) to get P in nt/meter^2, convert to atms
|
|
|
|
variable KcaltoJoule index 4.1868e3 # convert Kcals to Joules
|
|
variable NtMtoAtm equal 1.0/1.0135e5 # convert Nt/meter^2 to Atmospheres
|
|
|
|
variable tmpscale equal ${epskb}
|
|
variable epsJmole equal ${epsilon}*${KcaltoJoule}
|
|
variable massKgmole equal ${mass}/1000.0
|
|
variable sigmaM equal ${sigma}/1.0e10
|
|
variable sigmaMsq equal ${sigmaM}*${sigmaM}
|
|
variable tscale equal 1.0e15/sqrt(${epsJmole}/${massKgmole}/${sigmaMsq})
|
|
variable sigmaM3 equal ${sigmaM}*${sigmaM}*${sigmaM}
|
|
variable pscale equal ${NtMtoAtm}/(${sigmaM3}/(${epsJmole}/${avogadro}))
|
|
|
|
# variables
|
|
# alat = lattice constant in Angs (at reduced density rhostar)
|
|
# temp = reduced temperature for output
|
|
# epair,emol,etotal = reduced epair,emol,etotal energies for output
|
|
# press = reduced pressure for output
|
|
|
|
variable alat equal (4.0*${sigma}*${sigma}*${sigma}/${rhostar})^(1.0/3.0)
|
|
variable temp equal temp/${tmpscale}
|
|
variable epair equal epair/${epsilon}
|
|
variable emol equal emol/${epsilon}
|
|
variable etotal equal etotal/${epsilon}
|
|
variable press equal press/${pscale}
|
|
|
|
# same script as in.ar.lj
|
|
|
|
units real
|
|
atom_style atomic
|
|
|
|
lattice fcc ${alat}
|
|
region box block 0 $x 0 $y 0 $z
|
|
create_box 1 box
|
|
create_atoms 1 box
|
|
mass 1 ${mass}
|
|
|
|
velocity all create $(v_tinitial*v_epskb) 12345
|
|
|
|
pair_style lj/cut $(v_cutoff*v_sigma)
|
|
pair_coeff 1 1 ${epsilon} ${sigma}
|
|
|
|
neighbor $(v_skin*v_sigma) bin
|
|
neigh_modify delay 0 every 20 check no
|
|
|
|
fix 1 all nve
|
|
|
|
timestep $(v_dt*v_tscale)
|
|
|
|
# columns 2,3,4 = temp,pe,press in real units
|
|
# columns 5-9 = temp,energy.press in reduced units, compare to in.ar.lj
|
|
# need to include real unit output to enable use of reduced variables
|
|
|
|
thermo_style custom step temp pe press v_temp v_epair v_emol v_etotal v_press
|
|
thermo_modify norm yes
|
|
thermo ${nthermo}
|
|
|
|
run ${nsteps}
|