Merge pull request #1709 from lammps/units-example

KAPPA: compute thermal conductivity via several methods

MC: using LAMMPS in a Monte Carlo mode to relax the energy of a system

SPIN: examples for features of the SPIN package

UNITS: examples that run the same simulation in lj, real, metal units

USER: examples for USER packages and USER-contributed commands

VISCOSITY: compute viscosity via several methods :tb(s=:)

flow:     Couette and Poiseuille flow in a 2d channel

friction: frictional contact of spherical asperities between 2d surfaces

gcmc:     Grand Canonical Monte Carlo (GCMC) via the fix gcmc command

gjf:      use of fix langevin Gronbech-Jensen/Farago option

granregion: use of fix wall/region/gran as boundary on granular particles

hugoniostat: Hugoniostat shock dynamics

hyper:    global and local hyperdynamics of diffusion on Pt surface

energy-evaluation engine in a iterative Monte Carlo energy-relaxation

loop.

The UNITS directory contains examples of input scripts modeling the

same Lennard-Jones liquid model, written in 3 different unit systems:

lj, real, and metal.  So that you can see how to scale/unscale input

and output values read/written by LAMMPS to verify you are performing

the same simulation in different unit systems.

The USER directory contains subdirectories of user-provided example

scripts for ser packages.  See the README files in those directories

for more info.  See the doc/Section_start.html file for more info

This directory has 3 scripts which show how to run the same problem

using the 3 most common units system used in LAMMPS: lj, real, and

metal units.  As stated on the units command doc page:

"Any simulation you perform for one choice of units can be duplicated

with any other unit setting LAMMPS supports. ...  To perform the same

simulation in a different set of units you must change all the

unit-based input parameters in your input script and other input files

(data file, potential files, etc) correctly to the new units.  And you

must correctly convert all output from the new units to the old units

when comparing to the original results.  That is often not simple to

do."

These examples are meant to illustrate how to do this for a simple

Lennard-Jones liquid (argon).  All of the scripts have a set of

variables defined at the top which can be changed as command line

arguments (e.g. -v cutoff 3.0).  All 3 scripts give identical output,

modulo round-offs due to the finite precision of the conversion

factors used, either internally in LAMMPS or in the scripts.  If there

were run for a long time, the trajectories would diverge, but they

would still give statistically identical results.

The LJ script is the simplest; it is similar to the bench/in.lj

script.

The real and metal scripts each have a set of variables at the top

which define scale factors for converting quantities like distance,

energy, pressure from reduced LJ units to real or metal units.  Once

these are defined the rest of the input script is very similar to the

LJ script.  The approprate scale factor is applied to every input.

Output quantities are printed in both the native real/metal units and

unscaled back to LJ units.  So that you can see the outputs are the

same if you examine the log files.  Comments about this comparison

are at the bottom of the real and metal scripts.

These 3 scripts are provided, because converting from lj reduced units

to physical units (e.g. real or metal) or vice versa is the trickiest

case.  Converting input scripts between 2 sets of physical units

(e.g. reak <--> metal) is much easier.  But you can use the same ideas

as in these scripts; just define a set of scale/unscale factors.

See Allen & Tildesley's Computer Simulation of Liquids, Appendix B for

a nice discussion of reduced units.  It will explain the conversion

formulas used in the real and metal scripts.

Hopefully, if you study these scripts, you should be able to convert

an input script of your own, written in one set of units, to an

identical input script in an alternate set of units.  Where

"identical" means it runs the same simulation in a statistical sense.

You can find the full set of scale factors LAMMPS uses internally for

different unit systems it supports, at the top of the src/udpate.cpp

file.  A couple of those values are used in the real and metal

scripts.

# Ar in lj units

# simulation params in reduced units

# settable from command line

# epsilon = sigma = mass = 1.0

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

# script

units		lj

atom_style	atomic

lattice		fcc ${rhostar}

region		box block 0 $x 0 $y 0 $z

create_box	1 box

create_atoms	1 box

mass		1 1.0

velocity	all create ${tinitial} 12345

pair_style	lj/cut ${cutoff}

pair_coeff	1 1 1.0 1.0

neighbor	${skin} bin

neigh_modify	delay 0 every 20 check no

fix		1 all nve

timestep	${dt}

thermo		10

run		100

# Ar in metal 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 8.617343e-5          # kB in eV/K

variable        avogadro index 6.02214129e23  # Avogadro's number

# Ar properties in metal 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 eV

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 eV

# tmpscale = scale factor on temperature, converts reduced temp to degrees K

# tscale = scale factor on time, converts reduced time to ps

#   formula is t = t* / sqrt(epsilon/mass/sigma^2), but need t in fs

#   use epsilon (Joule), mass (kg/atom), sigma (meter) to get t in seconds

# pscale = scale factor on pressure, converts reduced pressure to bars

#   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 bars

variable        eVtoJoule index 1.602e-19     # convert eV to Joules

variable        NtMtoAtm equal 1.0e-5         # convert Nt/meter^2 to bars

variable        tmpscale equal ${epskb}

variable        epsilonJ equal ${epsilon}*${eVtoJoule}

variable        massKgAtom equal ${mass}/1000.0/${avogadro}

variable        sigmaM equal ${sigma}/1.0e10

variable        sigmaMsq equal ${sigmaM}*${sigmaM}

variable        tscale equal 1.0e12/sqrt(${epsilonJ}/${massKgAtom}/${sigmaMsq})

variable        sigmaM3 equal ${sigmaM}*${sigmaM}*${sigmaM}

variable        pscale equal ${NtMtoAtm}/(${sigmaM3}/(${epsilonJ}))

# 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		metal

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 metal units

# columns 5-9 = temp,energy.press in reduced units, compare to in.ar.lj

# need to include metal 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}