Section_howto.html
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Sat, Nov 30, 04:25

Section_howto.html
View Options

	<HTML>
	<CENTER><A HREF = "Section_accelerate.html">Previous Section</A> - <A HREF = "http://lammps.sandia.gov">LAMMPS WWW Site</A> - <A HREF = "Manual.html">LAMMPS Documentation</A> - <A HREF = "Section_commands.html#comm">LAMMPS Commands</A> - <A HREF = "Section_example.html">Next Section</A>
	</CENTER>






	<HR>

	<H3>6. How-to discussions
	</H3>
	<P>This section describes how to perform common tasks using LAMMPS.
	</P>
	6.1 <A HREF = "#howto_1">Restarting a simulation</A><BR>
	6.2 <A HREF = "#howto_2">2d simulations</A><BR>
	6.3 <A HREF = "#howto_3">CHARMM, AMBER, and DREIDING force fields</A><BR>
	6.4 <A HREF = "#howto_4">Running multiple simulations from one input script</A><BR>
	6.5 <A HREF = "#howto_5">Multi-replica simulations</A><BR>
	6.6 <A HREF = "#howto_6">Granular models</A><BR>
	6.7 <A HREF = "#howto_7">TIP3P water model</A><BR>
	6.8 <A HREF = "#howto_8">TIP4P water model</A><BR>
	6.9 <A HREF = "#howto_9">SPC water model</A><BR>
	6.10 <A HREF = "#howto_10">Coupling LAMMPS to other codes</A><BR>
	6.11 <A HREF = "#howto_11">Visualizing LAMMPS snapshots</A><BR>
	6.12 <A HREF = "#howto_12">Triclinic (non-orthogonal) simulation boxes</A><BR>
	6.13 <A HREF = "#howto_13">NEMD simulations</A><BR>
	6.14 <A HREF = "#howto_14">Finite-size spherical and aspherical particles</A><BR>
	6.15 <A HREF = "#howto_15">Output from LAMMPS (thermo, dumps, computes, fixes, variables)</A><BR>
	6.16 <A HREF = "#howto_16">Thermostatting, barostatting and computing temperature</A><BR>
	6.17 <A HREF = "#howto_17">Walls</A><BR>
	6.18 <A HREF = "#howto_18">Elastic constants</A><BR>
	6.19 <A HREF = "#howto_19">Library interface to LAMMPS</A><BR>
	6.20 <A HREF = "#howto_20">Calculating thermal conductivity</A><BR>
	6.21 <A HREF = "#howto_21">Calculating viscosity</A><BR>
	6.22 <A HREF = "#howto_22">Calculating a diffusion coefficient</A> <BR>

	<P>The example input scripts included in the LAMMPS distribution and
	highlighted in <A HREF = "Section_example.html">Section_example</A> also show how to
	setup and run various kinds of simulations.
	</P>
	<HR>

	<HR>

	<A NAME = "howto_1"></A><H4>6.1 Restarting a simulation
	</H4>
	<P>There are 3 ways to continue a long LAMMPS simulation. Multiple
	<A HREF = "run.html">run</A> commands can be used in the same input script. Each
	run will continue from where the previous run left off. Or binary
	restart files can be saved to disk using the <A HREF = "restart.html">restart</A>
	command. At a later time, these binary files can be read via a
	<A HREF = "read_restart.html">read_restart</A> command in a new script. Or they can
	be converted to text data files using the <A HREF = "Section_start.html#start_7">-r command-line
	switch</A> and read by a
	<A HREF = "read_data.html">read_data</A> command in a new script.
	</P>
	<P>Here we give examples of 2 scripts that read either a binary restart
	file or a converted data file and then issue a new run command to
	continue where the previous run left off. They illustrate what
	settings must be made in the new script. Details are discussed in the
	documentation for the <A HREF = "read_restart.html">read_restart</A> and
	<A HREF = "read_data.html">read_data</A> commands.
	</P>
	<P>Look at the <I>in.chain</I> input script provided in the <I>bench</I> directory
	of the LAMMPS distribution to see the original script that these 2
	scripts are based on. If that script had the line
	</P>
	<PRE>restart 50 tmp.restart
	</PRE>
	<P>added to it, it would produce 2 binary restart files (tmp.restart.50
	and tmp.restart.100) as it ran.
	</P>
	<P>This script could be used to read the 1st restart file and re-run the
	last 50 timesteps:
	</P>
	<PRE>read_restart tmp.restart.50
	</PRE>
	<PRE>neighbor 0.4 bin
	neigh_modify every 1 delay 1
	</PRE>
	<PRE>fix 1 all nve
	fix 2 all langevin 1.0 1.0 10.0 904297
	</PRE>
	<PRE>timestep 0.012
	</PRE>
	<PRE>run 50
	</PRE>
	<P>Note that the following commands do not need to be repeated because
	their settings are included in the restart file: <I>units, atom_style,
	special_bonds, pair_style, bond_style</I>. However these commands do
	need to be used, since their settings are not in the restart file:
	<I>neighbor, fix, timestep</I>.
	</P>
	<P>If you actually use this script to perform a restarted run, you will
	notice that the thermodynamic data match at step 50 (if you also put a
	"thermo 50" command in the original script), but do not match at step
	100. This is because the <A HREF = "fix_langevin.html">fix langevin</A> command
	uses random numbers in a way that does not allow for perfect restarts.
	</P>
	<P>As an alternate approach, the restart file could be converted to a data
	file as follows:
	</P>
	<PRE>lmp_g++ -r tmp.restart.50 tmp.restart.data
	</PRE>
	<P>Then, this script could be used to re-run the last 50 steps:
	</P>
	<PRE>units lj
	atom_style bond
	pair_style lj/cut 1.12
	pair_modify shift yes
	bond_style fene
	special_bonds 0.0 1.0 1.0
	</PRE>
	<PRE>read_data tmp.restart.data
	</PRE>
	<PRE>neighbor 0.4 bin
	neigh_modify every 1 delay 1
	</PRE>
	<PRE>fix 1 all nve
	fix 2 all langevin 1.0 1.0 10.0 904297
	</PRE>
	<PRE>timestep 0.012
	</PRE>
	<PRE>reset_timestep 50
	run 50
	</PRE>
	<P>Note that nearly all the settings specified in the original <I>in.chain</I>
	script must be repeated, except the <I>pair_coeff</I> and <I>bond_coeff</I>
	commands since the new data file lists the force field coefficients.
	Also, the <A HREF = "reset_timestep.html">reset_timestep</A> command is used to tell
	LAMMPS the current timestep. This value is stored in restart files,
	but not in data files.
	</P>
	<HR>

	<A NAME = "howto_2"></A><H4>6.2 2d simulations
	</H4>
	<P>Use the <A HREF = "dimension.html">dimension</A> command to specify a 2d simulation.
	</P>
	<P>Make the simulation box periodic in z via the <A HREF = "boundary.html">boundary</A>
	command. This is the default.
	</P>
	<P>If using the <A HREF = "create_box.html">create box</A> command to define a
	simulation box, set the z dimensions narrow, but finite, so that the
	create_atoms command will tile the 3d simulation box with a single z
	plane of atoms - e.g.
	</P>
	<PRE><A HREF = "create_box.html">create box</A> 1 -10 10 -10 10 -0.25 0.25
	</PRE>
	<P>If using the <A HREF = "read_data.html">read data</A> command to read in a file of
	atom coordinates, set the "zlo zhi" values to be finite but narrow,
	similar to the create_box command settings just described. For each
	atom in the file, assign a z coordinate so it falls inside the
	z-boundaries of the box - e.g. 0.0.
	</P>
	<P>Use the <A HREF = "fix_enforce2d.html">fix enforce2d</A> command as the last
	defined fix to insure that the z-components of velocities and forces
	are zeroed out every timestep. The reason to make it the last fix is
	so that any forces induced by other fixes will be zeroed out.
	</P>
	<P>Many of the example input scripts included in the LAMMPS distribution
	are for 2d models.
	</P>
	<P>IMPORTANT NOTE: Some models in LAMMPS treat particles as finite-size
	spheres, as opposed to point particles. In 2d, the particles will
	still be spheres, not disks, meaning their moment of inertia will be
	the same as in 3d.
	</P>
	<HR>

	<A NAME = "howto_3"></A><H4>6.3 CHARMM, AMBER, and DREIDING force fields
	</H4>
	<P>A force field has 2 parts: the formulas that define it and the
	coefficients used for a particular system. Here we only discuss
	formulas implemented in LAMMPS that correspond to formulas commonly
	used in the CHARMM, AMBER, and DREIDING force fields. Setting
	coefficients is done in the input data file via the
	<A HREF = "read_data.html">read_data</A> command or in the input script with
	commands like <A HREF = "pair_coeff.html">pair_coeff</A> or
	<A HREF = "bond_coeff.html">bond_coeff</A>. See <A HREF = "Section_tools.html">Section_tools</A>
	for additional tools that can use CHARMM or AMBER to assign force
	field coefficients and convert their output into LAMMPS input.
	</P>
	<P>See <A HREF = "#MacKerell">(MacKerell)</A> for a description of the CHARMM force
	field. See <A HREF = "#Cornell">(Cornell)</A> for a description of the AMBER force
	field.
	</P>




	<P>These style choices compute force field formulas that are consistent
	with common options in CHARMM or AMBER. See each command's
	documentation for the formula it computes.
	</P>
	<UL><LI><A HREF = "bond_harmonic.html">bond_style</A> harmonic
	<LI><A HREF = "angle_charmm.html">angle_style</A> charmm
	<LI><A HREF = "dihedral_charmm.html">dihedral_style</A> charmm
	<LI><A HREF = "pair_charmm.html">pair_style</A> lj/charmm/coul/charmm
	<LI><A HREF = "pair_charmm.html">pair_style</A> lj/charmm/coul/charmm/implicit
	<LI><A HREF = "pair_charmm.html">pair_style</A> lj/charmm/coul/long
	</UL>
	<UL><LI><A HREF = "special_bonds.html">special_bonds</A> charmm
	<LI><A HREF = "special_bonds.html">special_bonds</A> amber
	</UL>
	<P>DREIDING is a generic force field developed by the <A HREF = "http://www.wag.caltech.edu">Goddard
	group</A> at Caltech and is useful for
	predicting structures and dynamics of organic, biological and
	main-group inorganic molecules. The philosophy in DREIDING is to use
	general force constants and geometry parameters based on simple
	hybridization considerations, rather than individual force constants
	and geometric parameters that depend on the particular combinations of
	atoms involved in the bond, angle, or torsion terms. DREIDING has an
	<A HREF = "pair_hbond_dreiding.html">explicit hydrogen bond term</A> to describe
	interactions involving a hydrogen atom on very electronegative atoms
	(N, O, F).
	</P>
	<P>See <A HREF = "#Mayo">(Mayo)</A> for a description of the DREIDING force field
	</P>
	<P>These style choices compute force field formulas that are consistent
	with the DREIDING force field. See each command's
	documentation for the formula it computes.
	</P>
	<UL><LI><A HREF = "bond_harmonic.html">bond_style</A> harmonic
	<LI><A HREF = "bond_morse.html">bond_style</A> morse
	</UL>
	<UL><LI><A HREF = "angle_harmonic.html">angle_style</A> harmonic
	<LI><A HREF = "angle_cosine.html">angle_style</A> cosine
	<LI><A HREF = "angle_cosine_periodic.html">angle_style</A> cosine/periodic
	</UL>
	<UL><LI><A HREF = "dihedral_charmm.html">dihedral_style</A> charmm
	<LI><A HREF = "improper_umbrella.html">improper_style</A> umbrella
	</UL>
	<UL><LI><A HREF = "pair_buck.html">pair_style</A> buck
	<LI><A HREF = "pair_buck.html">pair_style</A> buck/coul/cut
	<LI><A HREF = "pair_buck.html">pair_style</A> buck/coul/long
	<LI><A HREF = "pair_lj.html">pair_style</A> lj/cut
	<LI><A HREF = "pair_lj.html">pair_style</A> lj/cut/coul/cut
	<LI><A HREF = "pair_lj.html">pair_style</A> lj/cut/coul/long
	</UL>
	<UL><LI><A HREF = "pair_hbond_dreiding.html">pair_style</A> hbond/dreiding/lj
	<LI><A HREF = "pair_hbond_dreiding.html">pair_style</A> hbond/dreiding/morse
	</UL>
	<UL><LI><A HREF = "special_bonds.html">special_bonds</A> dreiding
	</UL>
	<HR>

	<A NAME = "howto_4"></A><H4>6.4 Running multiple simulations from one input script
	</H4>
	<P>This can be done in several ways. See the documentation for
	individual commands for more details on how these examples work.
	</P>
	<P>If "multiple simulations" means continue a previous simulation for
	more timesteps, then you simply use the <A HREF = "run.html">run</A> command
	multiple times. For example, this script
	</P>
	<PRE>units lj
	atom_style atomic
	read_data data.lj
	run 10000
	run 10000
	run 10000
	run 10000
	run 10000
	</PRE>
	<P>would run 5 successive simulations of the same system for a total of
	50,000 timesteps.
	</P>
	<P>If you wish to run totally different simulations, one after the other,
	the <A HREF = "clear.html">clear</A> command can be used in between them to
	re-initialize LAMMPS. For example, this script
	</P>
	<PRE>units lj
	atom_style atomic
	read_data data.lj
	run 10000
	clear
	units lj
	atom_style atomic
	read_data data.lj.new
	run 10000
	</PRE>
	<P>would run 2 independent simulations, one after the other.
	</P>
	<P>For large numbers of independent simulations, you can use
	<A HREF = "variable.html">variables</A> and the <A HREF = "next.html">next</A> and
	<A HREF = "jump.html">jump</A> commands to loop over the same input script
	multiple times with different settings. For example, this
	script, named in.polymer
	</P>
	<PRE>variable d index run1 run2 run3 run4 run5 run6 run7 run8
	shell cd $d
	read_data data.polymer
	run 10000
	shell cd ..
	clear
	next d
	jump in.polymer
	</PRE>
	<P>would run 8 simulations in different directories, using a data.polymer
	file in each directory. The same concept could be used to run the
	same system at 8 different temperatures, using a temperature variable
	and storing the output in different log and dump files, for example
	</P>
	<PRE>variable a loop 8
	variable t index 0.8 0.85 0.9 0.95 1.0 1.05 1.1 1.15
	log log.$a
	read data.polymer
	velocity all create $t 352839
	fix 1 all nvt $t $t 100.0
	dump 1 all atom 1000 dump.$a
	run 100000
	clear
	next t
	next a
	jump in.polymer
	</PRE>
	<P>All of the above examples work whether you are running on 1 or
	multiple processors, but assumed you are running LAMMPS on a single
	partition of processors. LAMMPS can be run on multiple partitions via
	the "-partition" command-line switch as described in <A HREF = "Section_start.html#start_7">this
	section</A> of the manual.
	</P>
	<P>In the last 2 examples, if LAMMPS were run on 3 partitions, the same
	scripts could be used if the "index" and "loop" variables were
	replaced with <I>universe</I>-style variables, as described in the
	<A HREF = "variable.html">variable</A> command. Also, the "next t" and "next a"
	commands would need to be replaced with a single "next a t" command.
	With these modifications, the 8 simulations of each script would run
	on the 3 partitions one after the other until all were finished.
	Initially, 3 simulations would be started simultaneously, one on each
	partition. When one finished, that partition would then start
	the 4th simulation, and so forth, until all 8 were completed.
	</P>
	<HR>

	<A NAME = "howto_5"></A><H4>6.5 Multi-replica simulations
	</H4>
	<P>Several commands in LAMMPS run mutli-replica simulations, meaning
	that multiple instances (replicas) of your simulation are run
	simultaneously, with small amounts of data exchanged between replicas
	periodically.
	</P>
	<P>These are the relevant commands:
	</P>
	<UL><LI><A HREF = "neb.html">neb</A> for nudged elastic band calculations
	<LI><A HREF = "prd.html">prd</A> for parallel replica dynamics
	<LI><A HREF = "tad.html">tad</A> for temperature accelerated dynamics
	<LI><A HREF = "temper.html">temper</A> for parallel tempering
	<LI><A HREF = "fix_pimd.html">fix pimd</A> for path-integral molecular dynamics (PIMD)
	</UL>
	<P>NEB is a method for finding transition states and barrier energies.
	PRD and TAD are methods for performing accelerated dynamics to find
	and perform infrequent events. Parallel tempering or replica exchange
	runs different replicas at a series of temperature to facilitate
	rare-event sampling.
	</P>
	<P>These commands can only be used if LAMMPS was built with the REPLICA
	package. See the <A HREF = "Section_start.html#start_3">Making LAMMPS</A> section
	for more info on packages.
	</P>
	<P>PIMD runs different replicas whose individual particles are coupled
	together by springs to model a system or ring-polymers.
	</P>
	<P>This commands can only be used if LAMMPS was built with the USER-MISC
	package. See the <A HREF = "Section_start.html#start_3">Making LAMMPS</A> section
	for more info on packages.
	</P>
	<P>In all these cases, you must run with one or more processors per
	replica. The processors assigned to each replica are determined at
	run-time by using the <A HREF = "Section_start.html#start_7">-partition command-line
	switch</A> to launch LAMMPS on multiple
	partitions, which in this context are the same as replicas. E.g.
	these commands:
	</P>
	<PRE>mpirun -np 16 lmp_linux -partition 8x2 -in in.temper
	mpirun -np 8 lmp_linux -partition 8x1 -in in.neb
	</PRE>
	<P>would each run 8 replicas, on either 16 or 8 processors. Note the use
	of the <A HREF = "Section_start.html#start_7">-in command-line switch</A> to specify
	the input script which is required when running in multi-replica mode.
	</P>
	<P>Also note that with MPI installed on a machine (e.g. your desktop),
	you can run on more (virtual) processors than you have physical
	processors. Thus the above commands could be run on a
	single-processor (or few-processor) desktop so that you can run
	a multi-replica simulation on more replicas than you have
	physical processors.
	</P>
	<HR>

	<A NAME = "howto_6"></A><H4>6.6 Granular models
	</H4>
	<P>Granular system are composed of spherical particles with a diameter,
	as opposed to point particles. This means they have an angular
	velocity and torque can be imparted to them to cause them to rotate.
	</P>
	<P>To run a simulation of a granular model, you will want to use
	the following commands:
	</P>
	<UL><LI><A HREF = "atom_style.html">atom_style sphere</A>
	<LI><A HREF = "fix_nve_sphere.html">fix nve/sphere</A>
	<LI><A HREF = "fix_gravity.html">fix gravity</A>
	</UL>
	<P>This compute
	</P>
	<UL><LI><A HREF = "compute_erotate_sphere.html">compute erotate/sphere</A>
	</UL>
	<P>calculates rotational kinetic energy which can be <A HREF = "Section_howto.html#howto_15">output with
	thermodynamic info</A>.
	</P>
	<P>Use one of these 3 pair potentials, which compute forces and torques
	between interacting pairs of particles:
	</P>
	<UL><LI><A HREF = "pair_style.html">pair_style</A> gran/history
	<LI><A HREF = "pair_style.html">pair_style</A> gran/no_history
	<LI><A HREF = "pair_style.html">pair_style</A> gran/hertzian
	</UL>
	<P>These commands implement fix options specific to granular systems:
	</P>
	<UL><LI><A HREF = "fix_freeze.html">fix freeze</A>
	<LI><A HREF = "fix_pour.html">fix pour</A>
	<LI><A HREF = "fix_viscous.html">fix viscous</A>
	<LI><A HREF = "fix_wall_gran.html">fix wall/gran</A>
	</UL>
	<P>The fix style <I>freeze</I> zeroes both the force and torque of frozen
	atoms, and should be used for granular system instead of the fix style
	<I>setforce</I>.
	</P>
	<P>For computational efficiency, you can eliminate needless pairwise
	computations between frozen atoms by using this command:
	</P>
	<UL><LI><A HREF = "neigh_modify.html">neigh_modify</A> exclude
	</UL>
	<HR>

	<A NAME = "howto_7"></A><H4>6.7 TIP3P water model
	</H4>
	<P>The TIP3P water model as implemented in CHARMM
	<A HREF = "#MacKerell">(MacKerell)</A> specifies a 3-site rigid water molecule with
	charges and Lennard-Jones parameters assigned to each of the 3 atoms.
	In LAMMPS the <A HREF = "fix_shake.html">fix shake</A> command can be used to hold
	the two O-H bonds and the H-O-H angle rigid. A bond style of
	<I>harmonic</I> and an angle style of <I>harmonic</I> or <I>charmm</I> should also be
	used.
	</P>
	<P>These are the additional parameters (in real units) to set for O and H
	atoms and the water molecule to run a rigid TIP3P-CHARMM model with a
	cutoff. The K values can be used if a flexible TIP3P model (without
	fix shake) is desired. If the LJ epsilon and sigma for HH and OH are
	set to 0.0, it corresponds to the original 1983 TIP3P model
	<A HREF = "#Jorgensen">(Jorgensen)</A>.
	</P>
	<P>O mass = 15.9994<BR>
	H mass = 1.008<BR>
	O charge = -0.834<BR>
	H charge = 0.417<BR>
	LJ epsilon of OO = 0.1521<BR>
	LJ sigma of OO = 3.1507<BR>
	LJ epsilon of HH = 0.0460<BR>
	LJ sigma of HH = 0.4000<BR>
	LJ epsilon of OH = 0.0836<BR>
	LJ sigma of OH = 1.7753<BR>
	K of OH bond = 450<BR>
	r0 of OH bond = 0.9572<BR>
	K of HOH angle = 55<BR>
	theta of HOH angle = 104.52 <BR>
	</P>
	<P>These are the parameters to use for TIP3P with a long-range Coulombic
	solver (e.g. Ewald or PPPM in LAMMPS), see <A HREF = "#Price">(Price)</A> for
	details:
	</P>
	<P>O mass = 15.9994<BR>
	H mass = 1.008<BR>
	O charge = -0.830<BR>
	H charge = 0.415<BR>
	LJ epsilon of OO = 0.102<BR>
	LJ sigma of OO = 3.188<BR>
	LJ epsilon, sigma of OH, HH = 0.0<BR>
	K of OH bond = 450<BR>
	r0 of OH bond = 0.9572<BR>
	K of HOH angle = 55<BR>
	theta of HOH angle = 104.52 <BR>
	</P>
	<P>Wikipedia also has a nice article on <A HREF = "http://en.wikipedia.org/wiki/Water_model">water
	models</A>.
	</P>
	<HR>

	<A NAME = "howto_8"></A><H4>6.8 TIP4P water model
	</H4>
	<P>The four-point TIP4P rigid water model extends the traditional
	three-point TIP3P model by adding an additional site, usually
	massless, where the charge associated with the oxygen atom is placed.
	This site M is located at a fixed distance away from the oxygen along
	the bisector of the HOH bond angle. A bond style of <I>harmonic</I> and an
	angle style of <I>harmonic</I> or <I>charmm</I> should also be used.
	</P>
	<P>A TIP4P model is run with LAMMPS using either this command
	for a cutoff model:
	</P>
	<P><A HREF = "pair_lj.html">pair_style lj/cut/tip4p/cut</A>
	</P>
	<P>or these two commands for a long-range model:
	</P>
	<UL><LI><A HREF = "pair_lj.html">pair_style lj/cut/tip4p/long</A>
	<LI><A HREF = "kspace_style.html">kspace_style pppm/tip4p</A>
	</UL>
	<P>For both models, the bond lengths and bond angles should be held fixed
	using the <A HREF = "fix_shake.html">fix shake</A> command.
	</P>
	<P>These are the additional parameters (in real units) to set for O and H
	atoms and the water molecule to run a rigid TIP4P model with a cutoff
	<A HREF = "#Jorgensen">(Jorgensen)</A>. Note that the OM distance is specified in
	the <A HREF = "pair_style.html">pair_style</A> command, not as part of the pair
	coefficients.
	</P>
	<P>O mass = 15.9994<BR>
	H mass = 1.008<BR>
	O charge = -1.040<BR>
	H charge = 0.520<BR>
	r0 of OH bond = 0.9572<BR>
	theta of HOH angle = 104.52 <BR>
	OM distance = 0.15<BR>
	LJ epsilon of O-O = 0.1550<BR>
	LJ sigma of O-O = 3.1536<BR>
	LJ epsilon, sigma of OH, HH = 0.0<BR>
	Coulombic cutoff = 8.5 <BR>
	</P>
	<P>For the TIP4/Ice model (J Chem Phys, 122, 234511 (2005);
	http://dx.doi.org/10.1063/1.1931662) these values can be used:
	</P>
	<P>O mass = 15.9994<BR>
	H mass = 1.008<BR>
	O charge = -1.1794<BR>
	H charge = 0.5897<BR>
	r0 of OH bond = 0.9572<BR>
	theta of HOH angle = 104.52<BR>
	OM distance = 0.1577<BR>
	LJ epsilon of O-O = 0.21084<BR>
	LJ sigma of O-O = 3.1668<BR>
	LJ epsilon, sigma of OH, HH = 0.0<BR>
	Coulombic cutoff = 8.5 <BR>
	</P>
	<P>For the TIP4P/2005 model (J Chem Phys, 123, 234505 (2005);
	http://dx.doi.org/10.1063/1.2121687), these values can be used:
	</P>
	<P>O mass = 15.9994<BR>
	H mass = 1.008<BR>
	O charge = -1.1128<BR>
	H charge = 0.5564<BR>
	r0 of OH bond = 0.9572<BR>
	theta of HOH angle = 104.52<BR>
	OM distance = 0.1546<BR>
	LJ epsilon of O-O = 0.1852<BR>
	LJ sigma of O-O = 3.1589<BR>
	LJ epsilon, sigma of OH, HH = 0.0<BR>
	Coulombic cutoff = 8.5 <BR>
	</P>
	<P>These are the parameters to use for TIP4P with a long-range Coulombic
	solver (e.g. Ewald or PPPM in LAMMPS):
	</P>
	<P>O mass = 15.9994<BR>
	H mass = 1.008<BR>
	O charge = -1.0484<BR>
	H charge = 0.5242<BR>
	r0 of OH bond = 0.9572<BR>
	theta of HOH angle = 104.52<BR>
	OM distance = 0.1250<BR>
	LJ epsilon of O-O = 0.16275<BR>
	LJ sigma of O-O = 3.16435<BR>
	LJ epsilon, sigma of OH, HH = 0.0 <BR>
	</P>
	<P>Note that the when using the TIP4P pair style, the neighobr list
	cutoff for Coulomb interactions is effectively extended by a distance
	2 * (OM distance), to account for the offset distance of the
	fictitious charges on O atoms in water molecules. Thus it is
	typically best in an efficiency sense to use a LJ cutoff >= Coulomb
	cutoff + 2*(OM distance), to shrink the size of the neighbor list.
	This leads to slightly larger cost for the long-range calculation, so
	you can test the trade-off for your model. The OM distance and the LJ
	and Coulombic cutoffs are set in the <A HREF = "pair_lj.html">pair_style
	lj/cut/tip4p/long</A> command.
	</P>
	<P>Wikipedia also has a nice article on <A HREF = "http://en.wikipedia.org/wiki/Water_model">water
	models</A>.
	</P>
	<HR>

	<A NAME = "howto_9"></A><H4>6.9 SPC water model
	</H4>
	<P>The SPC water model specifies a 3-site rigid water molecule with
	charges and Lennard-Jones parameters assigned to each of the 3 atoms.
	In LAMMPS the <A HREF = "fix_shake.html">fix shake</A> command can be used to hold
	the two O-H bonds and the H-O-H angle rigid. A bond style of
	<I>harmonic</I> and an angle style of <I>harmonic</I> or <I>charmm</I> should also be
	used.
	</P>
	<P>These are the additional parameters (in real units) to set for O and H
	atoms and the water molecule to run a rigid SPC model.
	</P>
	<P>O mass = 15.9994<BR>
	H mass = 1.008<BR>
	O charge = -0.820<BR>
	H charge = 0.410<BR>
	LJ epsilon of OO = 0.1553<BR>
	LJ sigma of OO = 3.166<BR>
	LJ epsilon, sigma of OH, HH = 0.0<BR>
	r0 of OH bond = 1.0<BR>
	theta of HOH angle = 109.47 <BR>
	</P>
	<P>Note that as originally proposed, the SPC model was run with a 9
	Angstrom cutoff for both LJ and Coulommbic terms. It can also be used
	with long-range Coulombics (Ewald or PPPM in LAMMPS), without changing
	any of the parameters above, though it becomes a different model in
	that mode of usage.
	</P>
	<P>The SPC/E (extended) water model is the same, except
	the partial charge assignemnts change:
	</P>
	<P>O charge = -0.8476<BR>
	H charge = 0.4238 <BR>
	</P>
	<P>See the <A HREF = "#Berendsen">(Berendsen)</A> reference for more details on both
	the SPC and SPC/E models.
	</P>
	<P>Wikipedia also has a nice article on <A HREF = "http://en.wikipedia.org/wiki/Water_model">water
	models</A>.
	</P>
	<HR>

	<A NAME = "howto_10"></A><H4>6.10 Coupling LAMMPS to other codes
	</H4>
	<P>LAMMPS is designed to allow it to be coupled to other codes. For
	example, a quantum mechanics code might compute forces on a subset of
	atoms and pass those forces to LAMMPS. Or a continuum finite element
	(FE) simulation might use atom positions as boundary conditions on FE
	nodal points, compute a FE solution, and return interpolated forces on
	MD atoms.
	</P>
	<P>LAMMPS can be coupled to other codes in at least 3 ways. Each has
	advantages and disadvantages, which you'll have to think about in the
	context of your application.
	</P>
	<P>(1) Define a new <A HREF = "fix.html">fix</A> command that calls the other code. In
	this scenario, LAMMPS is the driver code. During its timestepping,
	the fix is invoked, and can make library calls to the other code,
	which has been linked to LAMMPS as a library. This is the way the
	<A HREF = "http://www.rpi.edu/~anderk5/lab">POEMS</A> package that performs constrained rigid-body motion on
	groups of atoms is hooked to LAMMPS. See the
	<A HREF = "fix_poems.html">fix_poems</A> command for more details. See <A HREF = "Section_modify.html">this
	section</A> of the documentation for info on how to add
	a new fix to LAMMPS.
	</P>


	<P>(2) Define a new LAMMPS command that calls the other code. This is
	conceptually similar to method (1), but in this case LAMMPS and the
	other code are on a more equal footing. Note that now the other code
	is not called during the timestepping of a LAMMPS run, but between
	runs. The LAMMPS input script can be used to alternate LAMMPS runs
	with calls to the other code, invoked via the new command. The
	<A HREF = "run.html">run</A> command facilitates this with its <I>every</I> option, which
	makes it easy to run a few steps, invoke the command, run a few steps,
	invoke the command, etc.
	</P>
	<P>In this scenario, the other code can be called as a library, as in
	(1), or it could be a stand-alone code, invoked by a system() call
	made by the command (assuming your parallel machine allows one or more
	processors to start up another program). In the latter case the
	stand-alone code could communicate with LAMMPS thru files that the
	command writes and reads.
	</P>
	<P>See <A HREF = "Section_modify.html">Section_modify</A> of the documentation for how
	to add a new command to LAMMPS.
	</P>
	<P>(3) Use LAMMPS as a library called by another code. In this case the
	other code is the driver and calls LAMMPS as needed. Or a wrapper
	code could link and call both LAMMPS and another code as libraries.
	Again, the <A HREF = "run.html">run</A> command has options that allow it to be
	invoked with minimal overhead (no setup or clean-up) if you wish to do
	multiple short runs, driven by another program.
	</P>
	<P>Examples of driver codes that call LAMMPS as a library are included in
	the examples/COUPLE directory of the LAMMPS distribution; see
	examples/COUPLE/README for more details:
	</P>
	<UL><LI>simple: simple driver programs in C++ and C which invoke LAMMPS as a
	library

	<LI>lammps_quest: coupling of LAMMPS and <A HREF = "http://dft.sandia.gov/Quest">Quest</A>, to run classical
	MD with quantum forces calculated by a density functional code

	<LI>lammps_spparks: coupling of LAMMPS and <A HREF = "http://www.sandia.gov/~sjplimp/spparks.html">SPPARKS</A>, to couple
	a kinetic Monte Carlo model for grain growth using MD to calculate
	strain induced across grain boundaries
	</UL>




	<P><A HREF = "Section_start.html#start_5">This section</A> of the documentation
	describes how to build LAMMPS as a library. Once this is done, you
	can interface with LAMMPS either via C++, C, Fortran, or Python (or
	any other language that supports a vanilla C-like interface). For
	example, from C++ you could create one (or more) "instances" of
	LAMMPS, pass it an input script to process, or execute individual
	commands, all by invoking the correct class methods in LAMMPS. From C
	or Fortran you can make function calls to do the same things. See
	<A HREF = "Section_python.html">Section_python</A> of the manual for a description
	of the Python wrapper provided with LAMMPS that operates through the
	LAMMPS library interface.
	</P>
	<P>The files src/library.cpp and library.h contain the C-style interface
	to LAMMPS. See <A HREF = "Section_howto.html#howto_19">Section_howto 19</A> of the
	manual for a description of the interface and how to extend it for
	your needs.
	</P>
	<P>Note that the lammps_open() function that creates an instance of
	LAMMPS takes an MPI communicator as an argument. This means that
	instance of LAMMPS will run on the set of processors in the
	communicator. Thus the calling code can run LAMMPS on all or a subset
	of processors. For example, a wrapper script might decide to
	alternate between LAMMPS and another code, allowing them both to run
	on all the processors. Or it might allocate half the processors to
	LAMMPS and half to the other code and run both codes simultaneously
	before syncing them up periodically. Or it might instantiate multiple
	instances of LAMMPS to perform different calculations.
	</P>
	<HR>

	<A NAME = "howto_11"></A><H4>6.11 Visualizing LAMMPS snapshots
	</H4>
	<P>LAMMPS itself does not do visualization, but snapshots from LAMMPS
	simulations can be visualized (and analyzed) in a variety of ways.
	</P>
	<P>LAMMPS snapshots are created by the <A HREF = "dump.html">dump</A> command which can
	create files in several formats. The native LAMMPS dump format is a
	text file (see "dump atom" or "dump custom") which can be visualized
	by the <A HREF = "Section_tools.html#xmovie">xmovie</A> program, included with the
	LAMMPS package. This produces simple, fast 2d projections of 3d
	systems, and can be useful for rapid debugging of simulation geometry
	and atom trajectories.
	</P>
	<P>Several programs included with LAMMPS as auxiliary tools can convert
	native LAMMPS dump files to other formats. See the
	<A HREF = "Section_tools.html">Section_tools</A> doc page for details. The first is
	the <A HREF = "Section_tools.html#charmm">ch2lmp tool</A>, which contains a
	lammps2pdb Perl script which converts LAMMPS dump files into PDB
	files. The second is the <A HREF = "Section_tools.html#arc">lmp2arc tool</A> which
	converts LAMMPS dump files into Accelrys' Insight MD program files.
	The third is the <A HREF = "Section_tools.html#cfg">lmp2cfg tool</A> which converts
	LAMMPS dump files into CFG files which can be read into the
	<A HREF = "http://mt.seas.upenn.edu/Archive/Graphics/A">AtomEye</A> visualizer.
	</P>
	<P>A Python-based toolkit distributed by our group can read native LAMMPS
	dump files, including custom dump files with additional columns of
	user-specified atom information, and convert them to various formats
	or pipe them into visualization software directly. See the <A HREF = "http://www.sandia.gov/~sjplimp/pizza.html">Pizza.py
	WWW site</A> for details. Specifically, Pizza.py can convert
	LAMMPS dump files into PDB, XYZ, <A HREF = "http://www.ensight.com">Ensight</A>, and VTK formats.
	Pizza.py can pipe LAMMPS dump files directly into the Raster3d and
	RasMol visualization programs. Pizza.py has tools that do interactive
	3d OpenGL visualization and one that creates SVG images of dump file
	snapshots.
	</P>
	<P>LAMMPS can create XYZ files directly (via "dump xyz") which is a
	simple text-based file format used by many visualization programs
	including <A HREF = "http://www.ks.uiuc.edu/Research/vmd">VMD</A>.
	</P>
	<P>LAMMPS can create DCD files directly (via "dump dcd") which can be
	read by <A HREF = "http://www.ks.uiuc.edu/Research/vmd">VMD</A> in conjunction with a CHARMM PSF file. Using this
	form of output avoids the need to convert LAMMPS snapshots to PDB
	files. See the <A HREF = "dump.html">dump</A> command for more information on DCD
	files.
	</P>
	<P>LAMMPS can create XTC files directly (via "dump xtc") which is GROMACS
	file format which can also be read by <A HREF = "http://www.ks.uiuc.edu/Research/vmd">VMD</A> for visualization.
	See the <A HREF = "dump.html">dump</A> command for more information on XTC files.
	</P>








	<HR>

	<A NAME = "howto_12"></A><H4>6.12 Triclinic (non-orthogonal) simulation boxes
	</H4>
	<P>By default, LAMMPS uses an orthogonal simulation box to encompass the
	particles. The <A HREF = "boundary.html">boundary</A> command sets the boundary
	conditions of the box (periodic, non-periodic, etc). The orthogonal
	box has its "origin" at (xlo,ylo,zlo) and is defined by 3 edge vectors
	starting from the origin given by <B>a</B> = (xhi-xlo,0,0); <B>b</B> =
	(0,yhi-ylo,0); <B>c</B> = (0,0,zhi-zlo). The 6 parameters
	(xlo,xhi,ylo,yhi,zlo,zhi) are defined at the time the simulation box
	is created, e.g. by the <A HREF = "create_box.html">create_box</A> or
	<A HREF = "read_data.html">read_data</A> or <A HREF = "read_restart.html">read_restart</A>
	commands. Additionally, LAMMPS defines box size parameters lx,ly,lz
	where lx = xhi-xlo, and similarly in the y and z dimensions. The 6
	parameters, as well as lx,ly,lz, can be output via the <A HREF = "thermo_style.html">thermo_style
	custom</A> command.
	</P>
	<P>LAMMPS also allows simulations to be performed in triclinic
	(non-orthogonal) simulation boxes shaped as a parallelepiped with
	triclinic symmetry. The parallelepiped has its "origin" at
	(xlo,ylo,zlo) and is defined by 3 edge vectors starting from the
	origin given by <B>a</B> = (xhi-xlo,0,0); <B>b</B> = (xy,yhi-ylo,0); <B>c</B> =
	(xz,yz,zhi-zlo). <I>xy,xz,yz</I> can be 0.0 or positive or negative values
	and are called "tilt factors" because they are the amount of
	displacement applied to faces of an originally orthogonal box to
	transform it into the parallelepiped. In LAMMPS the triclinic
	simulation box edge vectors <B>a</B>, <B>b</B>, and <B>c</B> cannot be arbitrary
	vectors. As indicated, <B>a</B> must lie on the positive x axis. <B>b</B> must
	lie in the xy plane, with strictly positive y component. <B>c</B> may have
	any orientation with strictly positive z component. The requirement
	that <B>a</B>, <B>b</B>, and <B>c</B> have strictly positive x, y, and z components,
	respectively, ensures that <B>a</B>, <B>b</B>, and <B>c</B> form a complete
	right-handed basis. These restrictions impose no loss of generality,
	since it is possible to rotate/invert any set of 3 crystal basis
	vectors so that they conform to the restrictions.
	</P>
	<P>For example, assume that the 3 vectors <B>A</B>,<B>B</B>,<B>C</B> are the edge
	vectors of a general parallelepiped, where there is no restriction on
	<B>A</B>,<B>B</B>,<B>C</B> other than they form a complete right-handed basis i.e.
	<B>A</B> x <B>B</B> . <B>C</B> > 0. The equivalent LAMMPS <B>a</B>,<B>b</B>,<B>c</B> are a linear
	rotation of <B>A</B>, <B>B</B>, and <B>C</B> and can be computed as follows:
	</P>
	<CENTER><IMG SRC = "Eqs/transform.jpg">
	</CENTER>
	<P>where A = \|<B>A</B>\| indicates the scalar length of <B>A</B>. The ^ hat symbol
	indicates the corresponding unit vector. <I>beta</I> and <I>gamma</I> are angles
	between the vectors described below. Note that by construction,
	<B>a</B>, <B>b</B>, and <B>c</B> have strictly positive x, y, and z components, respectively.
	If it should happen that
	<B>A</B>, <B>B</B>, and <B>C</B> form a left-handed basis, then the above equations
	are not valid for <B>c</B>. In this case, it is necessary
	to first apply an inversion. This can be achieved
	by interchanging two basis vectors or by changing the sign of one of them.
	</P>
	<P>For consistency, the same rotation/inversion applied to the basis vectors
	must also be applied to atom positions, velocities,
	and any other vector quantities.
	This can be conveniently achieved by first converting to
	fractional coordinates in the
	old basis and then converting to distance coordinates in the new basis.
	The transformation is given by the following equation:
	</P>
	<CENTER><IMG SRC = "Eqs/rotate.jpg">
	</CENTER>
	<P>where <I>V</I> is the volume of the box, <B>X</B> is the original vector quantity and
	<B>x</B> is the vector in the LAMMPS basis.
	</P>
	<P>There is no requirement that a triclinic box be periodic in any
	dimension, though it typically should be in at least the 2nd dimension
	of the tilt (y in xy) if you want to enforce a shift in periodic
	boundary conditions across that boundary. Some commands that work
	with triclinic boxes, e.g. the <A HREF = "fix_deform.html">fix deform</A> and <A HREF = "fix_nh.html">fix
	npt</A> commands, require periodicity or non-shrink-wrap
	boundary conditions in specific dimensions. See the command doc pages
	for details.
	</P>
	<P>The 9 parameters (xlo,xhi,ylo,yhi,zlo,zhi,xy,xz,yz) are defined at the
	time the simluation box is created. This happens in one of 3 ways.
	If the <A HREF = "create_box.html">create_box</A> command is used with a region of
	style <I>prism</I>, then a triclinic box is setup. See the
	<A HREF = "region.html">region</A> command for details. If the
	<A HREF = "read_data.html">read_data</A> command is used to define the simulation
	box, and the header of the data file contains a line with the "xy xz
	yz" keyword, then a triclinic box is setup. See the
	<A HREF = "read_data.html">read_data</A> command for details. Finally, if the
	<A HREF = "read_restart.html">read_restart</A> command reads a restart file which
	was written from a simulation using a triclinic box, then a triclinic
	box will be setup for the restarted simulation.
	</P>
	<P>Note that you can define a triclinic box with all 3 tilt factors =
	0.0, so that it is initially orthogonal. This is necessary if the box
	will become non-orthogonal, e.g. due to the <A HREF = "fix_nh.html">fix npt</A> or
	<A HREF = "fix_deform.html">fix deform</A> commands. Alternatively, you can use the
	<A HREF = "change_box.html">change_box</A> command to convert a simulation box from
	orthogonal to triclinic and vice versa.
	</P>
	<P>As with orthogonal boxes, LAMMPS defines triclinic box size parameters
	lx,ly,lz where lx = xhi-xlo, and similarly in the y and z dimensions.
	The 9 parameters, as well as lx,ly,lz, can be output via the
	<A HREF = "thermo_style.html">thermo_style custom</A> command.
	</P>
	<P>To avoid extremely tilted boxes (which would be computationally
	inefficient), LAMMPS normally requires that no tilt factor can skew
	the box more than half the distance of the parallel box length, which
	is the 1st dimension in the tilt factor (x for xz). This is required
	both when the simulation box is created, e.g. via the
	<A HREF = "create_box.html">create_box</A> or <A HREF = "read_data.html">read_data</A> commands,
	as well as when the box shape changes dynamically during a simulation,
	e.g. via the <A HREF = "fix_deform.html">fix deform</A> or <A HREF = "fix_nh.html">fix npt</A>
	commands.
	</P>
	<P>For example, if xlo = 2 and xhi = 12, then the x box length is 10 and
	the xy tilt factor must be between -5 and 5. Similarly, both xz and
	yz must be between -(xhi-xlo)/2 and +(yhi-ylo)/2. Note that this is
	not a limitation, since if the maximum tilt factor is 5 (as in this
	example), then configurations with tilt = ..., -15, -5, 5, 15, 25,
	... are geometrically all equivalent. If the box tilt exceeds this
	limit during a dynamics run (e.g. via the <A HREF = "fix_deform.html">fix deform</A>
	command), then the box is "flipped" to an equivalent shape with a tilt
	factor within the bounds, so the run can continue. See the <A HREF = "fix_deform.html">fix
	deform</A> doc page for further details.
	</P>
	<P>One exception to this rule is if the 1st dimension in the tilt
	factor (x for xy) is non-periodic. In that case, the limits on the
	tilt factor are not enforced, since flipping the box in that dimension
	does not change the atom positions due to non-periodicity. In this
	mode, if you tilt the system to extreme angles, the simulation will
	simply become inefficient, due to the highly skewed simulation box.
	</P>
	<P>The limitation on not creating a simulation box with a tilt factor
	skewing the box more than half the distance of the parallel box length
	can be overridden via the <A HREF = "box.html">box</A> command. Setting the <I>tilt</I>
	keyword to <I>large</I> allows any tilt factors to be specified.
	</P>
	<P>Box flips that may occur using the <A HREF = "fix_deform.html">fix deform</A> or
	<A HREF = "fix_nh.html">fix npt</A> commands can be turned off using the <I>flip no</I>
	option with either of the commands.
	</P>
	<P>Note that if a simulation box has a large tilt factor, LAMMPS will run
	less efficiently, due to the large volume of communication needed to
	acquire ghost atoms around a processor's irregular-shaped sub-domain.
	For extreme values of tilt, LAMMPS may also lose atoms and generate an
	error.
	</P>
	<P>Triclinic crystal structures are often defined using three lattice
	constants <I>a</I>, <I>b</I>, and <I>c</I>, and three angles <I>alpha</I>, <I>beta</I> and
	<I>gamma</I>. Note that in this nomenclature, the a, b, and c lattice
	constants are the scalar lengths of the edge vectors <B>a</B>, <B>b</B>, and <B>c</B>
	defined above. The relationship between these 6 quantities
	(a,b,c,alpha,beta,gamma) and the LAMMPS box sizes (lx,ly,lz) =
	(xhi-xlo,yhi-ylo,zhi-zlo) and tilt factors (xy,xz,yz) is as follows:
	</P>
	<CENTER><IMG SRC = "Eqs/box.jpg">
	</CENTER>
	<P>The inverse relationship can be written as follows:
	</P>
	<CENTER><IMG SRC = "Eqs/box_inverse.jpg">
	</CENTER>
	<P>The values of <I>a</I>, <I>b</I>, <I>c</I> , <I>alpha</I>, <I>beta</I> , and <I>gamma</I> can be printed
	out or accessed by computes using the
	<A HREF = "thermo_style.html">thermo_style custom</A> keywords
	<I>cella</I>, <I>cellb</I>, <I>cellc</I>, <I>cellalpha</I>, <I>cellbeta</I>, <I>cellgamma</I>,
	respectively.
	</P>
	<P>As discussed on the <A HREF = "dump.html">dump</A> command doc page, when the BOX
	BOUNDS for a snapshot is written to a dump file for a triclinic box,
	an orthogonal bounding box which encloses the triclinic simulation box
	is output, along with the 3 tilt factors (xy, xz, yz) of the triclinic
	box, formatted as follows:
	</P>
	<PRE>ITEM: BOX BOUNDS xy xz yz
	xlo_bound xhi_bound xy
	ylo_bound yhi_bound xz
	zlo_bound zhi_bound yz
	</PRE>
	<P>This bounding box is convenient for many visualization programs and is
	calculated from the 9 triclinic box parameters
	(xlo,xhi,ylo,yhi,zlo,zhi,xy,xz,yz) as follows:
	</P>
	<PRE>xlo_bound = xlo + MIN(0.0,xy,xz,xy+xz)
	xhi_bound = xhi + MAX(0.0,xy,xz,xy+xz)
	ylo_bound = ylo + MIN(0.0,yz)
	yhi_bound = yhi + MAX(0.0,yz)
	zlo_bound = zlo
	zhi_bound = zhi
	</PRE>
	<P>These formulas can be inverted if you need to convert the bounding box
	back into the triclinic box parameters, e.g. xlo = xlo_bound -
	MIN(0.0,xy,xz,xy+xz).
	</P>
	<P>One use of triclinic simulation boxes is to model solid-state crystals
	with triclinic symmetry. The <A HREF = "lattice.html">lattice</A> command can be
	used with non-orthogonal basis vectors to define a lattice that will
	tile a triclinic simulation box via the
	<A HREF = "create_atoms.html">create_atoms</A> command.
	</P>
	<P>A second use is to run Parinello-Rahman dyanamics via the <A HREF = "fix_nh.html">fix
	npt</A> command, which will adjust the xy, xz, yz tilt
	factors to compensate for off-diagonal components of the pressure
	tensor. The analalog for an <A HREF = "minimize.html">energy minimization</A> is
	the <A HREF = "fix_box_relax.html">fix box/relax</A> command.
	</P>
	<P>A third use is to shear a bulk solid to study the response of the
	material. The <A HREF = "fix_deform.html">fix deform</A> command can be used for
	this purpose. It allows dynamic control of the xy, xz, yz tilt
	factors as a simulation runs. This is discussed in the next section
	on non-equilibrium MD (NEMD) simulations.
	</P>
	<HR>

	<A NAME = "howto_13"></A><H4>6.13 NEMD simulations
	</H4>
	<P>Non-equilibrium molecular dynamics or NEMD simulations are typically
	used to measure a fluid's rheological properties such as viscosity.
	In LAMMPS, such simulations can be performed by first setting up a
	non-orthogonal simulation box (see the preceding Howto section).
	</P>
	<P>A shear strain can be applied to the simulation box at a desired
	strain rate by using the <A HREF = "fix_deform.html">fix deform</A> command. The
	<A HREF = "fix_nvt_sllod.html">fix nvt/sllod</A> command can be used to thermostat
	the sheared fluid and integrate the SLLOD equations of motion for the
	system. Fix nvt/sllod uses <A HREF = "compute_temp_deform.html">compute
	temp/deform</A> to compute a thermal temperature
	by subtracting out the streaming velocity of the shearing atoms. The
	velocity profile or other properties of the fluid can be monitored via
	the <A HREF = "fix_ave_spatial.html">fix ave/spatial</A> command.
	</P>
	<P>As discussed in the previous section on non-orthogonal simulation
	boxes, the amount of tilt or skew that can be applied is limited by
	LAMMPS for computational efficiency to be 1/2 of the parallel box
	length. However, <A HREF = "fix_deform.html">fix deform</A> can continuously strain
	a box by an arbitrary amount. As discussed in the <A HREF = "fix_deform.html">fix
	deform</A> command, when the tilt value reaches a limit,
	the box is flipped to the opposite limit which is an equivalent tiling
	of periodic space. The strain rate can then continue to change as
	before. In a long NEMD simulation these box re-shaping events may
	occur many times.
	</P>
	<P>In a NEMD simulation, the "remap" option of <A HREF = "fix_deform.html">fix
	deform</A> should be set to "remap v", since that is what
	<A HREF = "fix_nvt_sllod.html">fix nvt/sllod</A> assumes to generate a velocity
	profile consistent with the applied shear strain rate.
	</P>
	<P>An alternative method for calculating viscosities is provided via the
	<A HREF = "fix_viscosity.html">fix viscosity</A> command.
	</P>
	<HR>

	<A NAME = "howto_14"></A><H4>6.14 Finite-size spherical and aspherical particles
	</H4>
	<P>Typical MD models treat atoms or particles as point masses. Sometimes
	it is desirable to have a model with finite-size particles such as
	spheroids or ellipsoids or generalized aspherical bodies. The
	difference is that such particles have a moment of inertia, rotational
	energy, and angular momentum. Rotation is induced by torque coming
	from interactions with other particles.
	</P>
	<P>LAMMPS has several options for running simulations with these kinds of
	particles. The following aspects are discussed in turn:
	</P>
	<UL><LI>atom styles
	<LI>pair potentials
	<LI>time integration
	<LI>computes, thermodynamics, and dump output
	<LI>rigid bodies composed of finite-size particles
	</UL>
	<P>Example input scripts for these kinds of models are in the body,
	colloid, dipole, ellipse, line, peri, pour, and tri directories of the
	<A HREF = "Section_example.html">examples directory</A> in the LAMMPS distribution.
	</P>
	<H5>Atom styles
	</H5>
	<P>There are several <A HREF = "atom_style.html">atom styles</A> that allow for
	definition of finite-size particles: sphere, dipole, ellipsoid, line,
	tri, peri, and body.
	</P>
	<P>The sphere style defines particles that are spheriods and each
	particle can have a unique diameter and mass (or density). These
	particles store an angular velocity (omega) and can be acted upon by
	torque. The "set" command can be used to modify the diameter and mass
	of individual particles, after then are created.
	</P>
	<P>The dipole style does not actually define finite-size particles, but
	is often used in conjunction with spherical particles, via a command
	like
	</P>
	<PRE>atom_style hybrid sphere dipole
	</PRE>
	<P>This is because when dipoles interact with each other, they induce
	torques, and a particle must be finite-size (i.e. have a moment of
	inertia) in order to respond and rotate. See the <A HREF = "atom_style.html">atom_style
	dipole</A> command for details. The "set" command can be
	used to modify the orientation and length of the dipole moment of
	individual particles, after then are created.
	</P>
	<P>The ellipsoid style defines particles that are ellipsoids and thus can
	be aspherical. Each particle has a shape, specified by 3 diameters,
	and mass (or density). These particles store an angular momentum and
	their orientation (quaternion), and can be acted upon by torque. They
	do not store an angular velocity (omega), which can be in a different
	direction than angular momentum, rather they compute it as needed.
	The "set" command can be used to modify the diameter, orientation, and
	mass of individual particles, after then are created. It also has a
	brief explanation of what quaternions are.
	</P>
	<P>The line style defines line segment particles with two end points and
	a mass (or density). They can be used in 2d simulations, and they can
	be joined together to form rigid bodies which represent arbitrary
	polygons.
	</P>
	<P>The tri style defines triangular particles with three corner points
	and a mass (or density). They can be used in 3d simulations, and they
	can be joined together to form rigid bodies which represent arbitrary
	particles with a triangulated surface.
	</P>
	<P>The peri style is used with <A HREF = "pair_peri.html">Peridynamic models</A> and
	defines particles as having a volume, that is used internally in the
	<A HREF = "pair_peri.html">pair_style peri</A> potentials.
	</P>
	<P>The body style allows for definition of particles which can represent
	complex entities, such as surface meshes of discrete points,
	collections of sub-particles, deformable objects, etc. The body style
	is discussed in more detail on the <A HREF = "body.html">body</A> doc page.
	</P>
	<P>Note that if one of these atom styles is used (or multiple styles via
	the <A HREF = "atom_style.html">atom_style hybrid</A> command), not all particles in
	the system are required to be finite-size or aspherical.
	</P>
	<P>For example, in the ellipsoid style, if the 3 shape parameters are set
	to the same value, the particle will be a sphere rather than an
	ellipsoid. If the 3 shape parameters are all set to 0.0 or if the
	diameter is set to 0.0, it will be a point particle. In the line or
	tri style, if the lineflag or triflag is specified as 0, then it
	will be a point particle.
	</P>
	<P>Some of the pair styles used to compute pairwise interactions between
	finite-size particles also compute the correct interaction with point
	particles as well, e.g. the interaction between a point particle and a
	finite-size particle or between two point particles. If necessary,
	<A HREF = "pair_hybrid.html">pair_style hybrid</A> can be used to insure the correct
	interactions are computed for the appropriate style of interactions.
	Likewise, using groups to partition particles (ellipsoids versus
	spheres versus point particles) will allow you to use the appropriate
	time integrators and temperature computations for each class of
	particles. See the doc pages for various commands for details.
	</P>
	<P>Also note that for <A HREF = "dimension.html">2d simulations</A>, atom styles sphere
	and ellipsoid still use 3d particles, rather than as circular disks or
	ellipses. This means they have the same moment of inertia as the 3d
	object. When temperature is computed, the correct degrees of freedom
	are used for rotation in a 2d versus 3d system.
	</P>
	<H5>Pair potentials
	</H5>
	<P>When a system with finite-size particles is defined, the particles
	will only rotate and experience torque if the force field computes
	such interactions. These are the various <A HREF = "pair_style.html">pair
	styles</A> that generate torque:
	</P>
	<UL><LI><A HREF = "pair_gran.html">pair_style gran/history</A>
	<LI><A HREF = "pair_gran.html">pair_style gran/hertzian</A>
	<LI><A HREF = "pair_gran.html">pair_style gran/no_history</A>
	<LI><A HREF = "pair_dipole.html">pair_style dipole/cut</A>
	<LI><A HREF = "pair_gayberne.html">pair_style gayberne</A>
	<LI><A HREF = "pair_resquared.html">pair_style resquared</A>
	<LI><A HREF = "pair_brownian.html">pair_style brownian</A>
	<LI><A HREF = "pair_lubricate.html">pair_style lubricate</A>
	<LI><A HREF = "pair_line_lj.html">pair_style line/lj</A>
	<LI><A HREF = "pair_tri_lj.html">pair_style tri/lj</A>
	<LI><A HREF = "pair_body.html">pair_style body</A>
	</UL>
	<P>The granular pair styles are used with spherical particles. The
	dipole pair style is used with the dipole atom style, which could be
	applied to spherical or ellipsoidal particles. The GayBerne and
	REsquared potentials require ellipsoidal particles, though they will
	also work if the 3 shape parameters are the same (a sphere). The
	Brownian and lubrication potentials are used with spherical particles.
	The line, tri, and body potentials are used with line segment,
	triangular, and body particles respectively.
	</P>
	<H5>Time integration
	</H5>
	<P>There are several fixes that perform time integration on finite-size
	spherical particles, meaning the integrators update the rotational
	orientation and angular velocity or angular momentum of the particles:
	</P>
	<UL><LI><A HREF = "fix_nve_sphere.html">fix nve/sphere</A>
	<LI><A HREF = "fix_nvt_sphere.html">fix nvt/sphere</A>
	<LI><A HREF = "fix_npt_sphere.html">fix npt/sphere</A>
	</UL>
	<P>Likewise, there are 3 fixes that perform time integration on
	ellipsoidal particles:
	</P>
	<UL><LI><A HREF = "fix_nve_asphere.html">fix nve/asphere</A>
	<LI><A HREF = "fix_nvt_asphere.html">fix nvt/asphere</A>
	<LI><A HREF = "fix_npt_asphere.html">fix npt/asphere</A>
	</UL>
	<P>The advantage of these fixes is that those which thermostat the
	particles include the rotational degrees of freedom in the temperature
	calculation and thermostatting. The <A HREF = "fix_langevin">fix langevin</A>
	command can also be used with its <I>omgea</I> or <I>angmom</I> options to
	thermostat the rotational degrees of freedom for spherical or
	ellipsoidal particles. Other thermostatting fixes only operate on the
	translational kinetic energy of finite-size particles.
	</P>
	<P>These fixes perform constant NVE time integration on line segment,
	triangular, and body particles:
	</P>
	<UL><LI><A HREF = "fix_nve_line.html">fix nve/line</A>
	<LI><A HREF = "fix_nve_tri.html">fix nve/tri</A>
	<LI><A HREF = "fix_nve_body.html">fix nve/body</A>
	</UL>
	<P>Note that for mixtures of point and finite-size particles, these
	integration fixes can only be used with <A HREF = "group.html">groups</A> which
	contain finite-size particles.
	</P>
	<H5>Computes, thermodynamics, and dump output
	</H5>
	<P>There are several computes that calculate the temperature or
	rotational energy of spherical or ellipsoidal particles:
	</P>
	<UL><LI><A HREF = "compute_temp_sphere.html">compute temp/sphere</A>
	<LI><A HREF = "compute_temp_asphere.html">compute temp/asphere</A>
	<LI><A HREF = "compute_erotate_sphere.html">compute erotate/sphere</A>
	<LI><A HREF = "compute_erotate_asphere.html">compute erotate/asphere</A>
	</UL>
	<P>These include rotational degrees of freedom in their computation. If
	you wish the thermodynamic output of temperature or pressure to use
	one of these computes (e.g. for a system entirely composed of
	finite-size particles), then the compute can be defined and the
	<A HREF = "thermo_modify.html">thermo_modify</A> command used. Note that by default
	thermodynamic quantities will be calculated with a temperature that
	only includes translational degrees of freedom. See the
	<A HREF = "thermo_style.html">thermo_style</A> command for details.
	</P>
	<P>These commands can be used to output various attributes of finite-size
	particles:
	</P>
	<UL><LI><A HREF = "dump.html">dump custom</A>
	<LI><A HREF = "compute_property_atom.html">compute property/atom</A>
	<LI><A HREF = "dump.html">dump local</A>
	<LI><A HREF = "compute_body_local.html">compute body/local</A>
	</UL>
	<P>Attributes include the dipole moment, the angular velocity, the
	angular momentum, the quaternion, the torque, the end-point and
	corner-point coordinates (for line and tri particles), and
	sub-particle attributes of body particles.
	</P>
	<H5>Rigid bodies composed of finite-size particles
	</H5>
	<P>The <A HREF = "fix_rigid.html">fix rigid</A> command treats a collection of
	particles as a rigid body, computes its inertia tensor, sums the total
	force and torque on the rigid body each timestep due to forces on its
	constituent particles, and integrates the motion of the rigid body.
	</P>
	<P>If any of the constituent particles of a rigid body are finite-size
	particles (spheres or ellipsoids or line segments or triangles), then
	their contribution to the inertia tensor of the body is different than
	if they were point particles. This means the rotational dynamics of
	the rigid body will be different. Thus a model of a dimer is
	different if the dimer consists of two point masses versus two
	spheroids, even if the two particles have the same mass. Finite-size
	particles that experience torque due to their interaction with other
	particles will also impart that torque to a rigid body they are part
	of.
	</P>
	<P>See the "fix rigid" command for example of complex rigid-body models
	it is possible to define in LAMMPS.
	</P>
	<P>Note that the <A HREF = "fix_shake.html">fix shake</A> command can also be used to
	treat 2, 3, or 4 particles as a rigid body, but it always assumes the
	particles are point masses.
	</P>
	<P>Also note that body particles cannot be modeled with the <A HREF = "fix_rigid.html">fix
	rigid</A> command. Body particles are treated by LAMMPS
	as single particles, though they can store internal state, such as a
	list of sub-particles. Individual body partices are typically treated
	as rigid bodies, and their motion integrated with a command like <A HREF = "fix_nve_body.html">fix
	nve/body</A>. Interactions between pairs of body
	particles are computed via a command like <A HREF = "pair_body.html">pair_style
	body</A>.
	</P>
	<HR>

	<A NAME = "howto_15"></A><H4>6.15 Output from LAMMPS (thermo, dumps, computes, fixes, variables)
	</H4>
	<P>There are four basic kinds of LAMMPS output:
	</P>
	<UL><LI><A HREF = "thermo_style.html">Thermodynamic output</A>, which is a list
	of quantities printed every few timesteps to the screen and logfile.

	<LI><A HREF = "dump.html">Dump files</A>, which contain snapshots of atoms and various
	per-atom values and are written at a specified frequency.

	<LI>Certain fixes can output user-specified quantities to files: <A HREF = "fix_ave_time.html">fix
	ave/time</A> for time averaging, <A HREF = "fix_ave_spatial.html">fix
	ave/spatial</A> for spatial averaging, and <A HREF = "fix_print.html">fix
	print</A> for single-line output of
	<A HREF = "variable.html">variables</A>. Fix print can also output to the
	screen.

	<LI><A HREF = "restart.html">Restart files</A>.
	</UL>
	<P>A simulation prints one set of thermodynamic output and (optionally)
	restart files. It can generate any number of dump files and fix
	output files, depending on what <A HREF = "dump.html">dump</A> and <A HREF = "fix.html">fix</A>
	commands you specify.
	</P>
	<P>As discussed below, LAMMPS gives you a variety of ways to determine
	what quantities are computed and printed when the thermodynamics,
	dump, or fix commands listed above perform output. Throughout this
	discussion, note that users can also <A HREF = "Section_modify.html">add their own computes and fixes
	to LAMMPS</A> which can then generate values that can
	then be output with these commands.
	</P>
	<P>The following sub-sections discuss different LAMMPS command related
	to output and the kind of data they operate on and produce:
	</P>
	<UL><LI><A HREF = "#global">Global/per-atom/local data</A>
	<LI><A HREF = "#scalar">Scalar/vector/array data</A>
	<LI><A HREF = "#thermo">Thermodynamic output</A>
	<LI><A HREF = "#dump">Dump file output</A>
	<LI><A HREF = "#fixoutput">Fixes that write output files</A>
	<LI><A HREF = "#computeoutput">Computes that process output quantities</A>
	<LI><A HREF = "#fixoutput">Fixes that process output quantities</A>
	<LI><A HREF = "#compute">Computes that generate values to output</A>
	<LI><A HREF = "#fix">Fixes that generate values to output</A>
	<LI><A HREF = "#variable">Variables that generate values to output</A>
	<LI><A HREF = "#table">Summary table of output options and data flow between commands</A>
	</UL>
	<H5><A NAME = "global"></A>Global/per-atom/local data
	</H5>
	<P>Various output-related commands work with three different styles of
	data: global, per-atom, or local. A global datum is one or more
	system-wide values, e.g. the temperature of the system. A per-atom
	datum is one or more values per atom, e.g. the kinetic energy of each
	atom. Local datums are calculated by each processor based on the
	atoms it owns, but there may be zero or more per atom, e.g. a list of
	bond distances.
	</P>
	<H5><A NAME = "scalar"></A>Scalar/vector/array data
	</H5>
	<P>Global, per-atom, and local datums can each come in three kinds: a
	single scalar value, a vector of values, or a 2d array of values. The
	doc page for a "compute" or "fix" or "variable" that generates data
	will specify both the style and kind of data it produces, e.g. a
	per-atom vector.
	</P>
	<P>When a quantity is accessed, as in many of the output commands
	discussed below, it can be referenced via the following bracket
	notation, where ID in this case is the ID of a compute. The leading
	"c_" would be replaced by "f_" for a fix, or "v_" for a variable:
	</P>
	<DIV ALIGN=center><TABLE BORDER=1 >
	<TR><TD >c_ID </TD><TD > entire scalar, vector, or array</TD></TR>
	<TR><TD >c_ID[I] </TD><TD > one element of vector, one column of array</TD></TR>
	<TR><TD >c_ID[I][J] </TD><TD > one element of array
	</TD></TR></TABLE></DIV>

	<P>In other words, using one bracket reduces the dimension of the data
	once (vector -> scalar, array -> vector). Using two brackets reduces
	the dimension twice (array -> scalar). Thus a command that uses
	scalar values as input can typically also process elements of a vector
	or array.
	</P>
	<H5><A NAME = "thermo"></A>Thermodynamic output
	</H5>
	<P>The frequency and format of thermodynamic output is set by the
	<A HREF = "thermo.html">thermo</A>, <A HREF = "thermo_style.html">thermo_style</A>, and
	<A HREF = "thermo_modify.html">thermo_modify</A> commands. The
	<A HREF = "thermo_style.html">thermo_style</A> command also specifies what values
	are calculated and written out. Pre-defined keywords can be specified
	(e.g. press, etotal, etc). Three additional kinds of keywords can
	also be specified (c_ID, f_ID, v_name), where a <A HREF = "compute.html">compute</A>
	or <A HREF = "fix.html">fix</A> or <A HREF = "variable.html">variable</A> provides the value to be
	output. In each case, the compute, fix, or variable must generate
	global values for input to the <A HREF = "dump.html">thermo_style custom</A>
	command.
	</P>
	<P>Note that thermodynamic output values can be "extensive" or
	"intensive". The former scale with the number of atoms in the system
	(e.g. total energy), the latter do not (e.g. temperature). The
	setting for <A HREF = "thermo_modify.html">thermo_modify norm</A> determines whether
	extensive quantities are normalized or not. Computes and fixes
	produce either extensive or intensive values; see their individual doc
	pages for details. <A HREF = "variable.html">Equal-style variables</A> produce only
	intensive values; you can include a division by "natoms" in the
	formula if desired, to make an extensive calculation produce an
	intensive result.
	</P>
	<H5><A NAME = "dump"></A>Dump file output
	</H5>
	<P>Dump file output is specified by the <A HREF = "dump.html">dump</A> and
	<A HREF = "dump_modify.html">dump_modify</A> commands. There are several
	pre-defined formats (dump atom, dump xtc, etc).
	</P>
	<P>There is also a <A HREF = "dump.html">dump custom</A> format where the user
	specifies what values are output with each atom. Pre-defined atom
	attributes can be specified (id, x, fx, etc). Three additional kinds
	of keywords can also be specified (c_ID, f_ID, v_name), where a
	<A HREF = "compute.html">compute</A> or <A HREF = "fix.html">fix</A> or <A HREF = "variable.html">variable</A>
	provides the values to be output. In each case, the compute, fix, or
	variable must generate per-atom values for input to the <A HREF = "dump.html">dump
	custom</A> command.
	</P>
	<P>There is also a <A HREF = "dump.html">dump local</A> format where the user specifies
	what local values to output. A pre-defined index keyword can be
	specified to enumuerate the local values. Two additional kinds of
	keywords can also be specified (c_ID, f_ID), where a
	<A HREF = "compute.html">compute</A> or <A HREF = "fix.html">fix</A> or <A HREF = "variable.html">variable</A>
	provides the values to be output. In each case, the compute or fix
	must generate local values for input to the <A HREF = "dump.html">dump local</A>
	command.
	</P>
	<H5><A NAME = "fixoutput"></A>Fixes that write output files
	</H5>
	<P>Several fixes take various quantities as input and can write output
	files: <A HREF = "fix_ave_time.html">fix ave/time</A>, <A HREF = "fix_ave_spatial.html">fix
	ave/spatial</A>, <A HREF = "fix_ave_histo.html">fix ave/histo</A>,
	<A HREF = "fix_ave_correlate.html">fix ave/correlate</A>, and <A HREF = "fix_print.html">fix
	print</A>.
	</P>
	<P>The <A HREF = "fix_ave_time.html">fix ave/time</A> command enables direct output to
	a file and/or time-averaging of global scalars or vectors. The user
	specifies one or more quantities as input. These can be global
	<A HREF = "compute.html">compute</A> values, global <A HREF = "fix.html">fix</A> values, or
	<A HREF = "variable.html">variables</A> of any style except the atom style which
	produces per-atom values. Since a variable can refer to keywords used
	by the <A HREF = "thermo_style.html">thermo_style custom</A> command (like temp or
	press) and individual per-atom values, a wide variety of quantities
	can be time averaged and/or output in this way. If the inputs are one
	or more scalar values, then the fix generate a global scalar or vector
	of output. If the inputs are one or more vector values, then the fix
	generates a global vector or array of output. The time-averaged
	output of this fix can also be used as input to other output commands.
	</P>
	<P>The <A HREF = "fix_ave_spatial.html">fix ave/spatial</A> command enables direct
	output to a file of spatial-averaged per-atom quantities like those
	output in dump files, within 1d layers of the simulation box. The
	per-atom quantities can be atom density (mass or number) or atom
	attributes such as position, velocity, force. They can also be
	per-atom quantities calculated by a <A HREF = "compute.html">compute</A>, by a
	<A HREF = "fix.html">fix</A>, or by an atom-style <A HREF = "variable.html">variable</A>. The
	spatial-averaged output of this fix can also be used as input to other
	output commands.
	</P>
	<P>The <A HREF = "fix_ave_histo.html">fix ave/histo</A> command enables direct output
	to a file of histogrammed quantities, which can be global or per-atom
	or local quantities. The histogram output of this fix can also be
	used as input to other output commands.
	</P>
	<P>The <A HREF = "fix_ave_correlate.html">fix ave/correlate</A> command enables direct
	output to a file of time-correlated quantities, which can be global
	scalars. The correlation matrix output of this fix can also be used
	as input to other output commands.
	</P>
	<P>The <A HREF = "fix_print.html">fix print</A> command can generate a line of output
	written to the screen and log file or to a separate file, periodically
	during a running simulation. The line can contain one or more
	<A HREF = "variable.html">variable</A> values for any style variable except the atom
	style). As explained above, variables themselves can contain
	references to global values generated by <A HREF = "thermo_style.html">thermodynamic
	keywords</A>, <A HREF = "compute.html">computes</A>,
	<A HREF = "fix.html">fixes</A>, or other <A HREF = "variable.html">variables</A>, or to per-atom
	values for a specific atom. Thus the <A HREF = "fix_print.html">fix print</A>
	command is a means to output a wide variety of quantities separate
	from normal thermodynamic or dump file output.
	</P>
	<H5><A NAME = "computeoutput"></A>Computes that process output quantities
	</H5>
	<P>The <A HREF = "compute_reduce.html">compute reduce</A> and <A HREF = "compute_reduce.html">compute
	reduce/region</A> commands take one or more per-atom
	or local vector quantities as inputs and "reduce" them (sum, min, max,
	ave) to scalar quantities. These are produced as output values which
	can be used as input to other output commands.
	</P>
	<P>The <A HREF = "compute_slice.html">compute slice</A> command take one or more global
	vector or array quantities as inputs and extracts a subset of their
	values to create a new vector or array. These are produced as output
	values which can be used as input to other output commands.
	</P>
	<P>The <A HREF = "compute_property_atom.html">compute property/atom</A> command takes a
	list of one or more pre-defined atom attributes (id, x, fx, etc) and
	stores the values in a per-atom vector or array. These are produced
	as output values which can be used as input to other output commands.
	The list of atom attributes is the same as for the <A HREF = "dump.html">dump
	custom</A> command.
	</P>
	<P>The <A HREF = "compute_property_local.html">compute property/local</A> command takes
	a list of one or more pre-defined local attributes (bond info, angle
	info, etc) and stores the values in a local vector or array. These
	are produced as output values which can be used as input to other
	output commands.
	</P>
	<P>The <A HREF = "compute_atom_molecule.html">compute atom/molecule</A> command takes a
	list of one or more per-atom quantities (from a compute, fix, per-atom
	variable) and sums the quantities on a per-molecule basis. It
	produces a global vector or array as output values which can be used
	as input to other output commands.
	</P>
	<H5><A NAME = "fixoutput"></A>Fixes that process output quantities
	</H5>
	<P>The <A HREF = "fix_vector.html">fix vector</A> command can create global vectors as
	output from global scalars as input, accumulating them one element at
	a time.
	</P>
	<P>The <A HREF = "fix_ave_atom.html">fix ave/atom</A> command performs time-averaging
	of per-atom vectors. The per-atom quantities can be atom attributes
	such as position, velocity, force. They can also be per-atom
	quantities calculated by a <A HREF = "compute.html">compute</A>, by a
	<A HREF = "fix.html">fix</A>, or by an atom-style <A HREF = "variable.html">variable</A>. The
	time-averaged per-atom output of this fix can be used as input to
	other output commands.
	</P>
	<P>The <A HREF = "fix_store_state.html">fix store/state</A> command can archive one or
	more per-atom attributes at a particular time, so that the old values
	can be used in a future calculation or output. The list of atom
	attributes is the same as for the <A HREF = "dump.html">dump custom</A> command,
	including per-atom quantities calculated by a <A HREF = "compute.html">compute</A>,
	by a <A HREF = "fix.html">fix</A>, or by an atom-style <A HREF = "variable.html">variable</A>.
	The output of this fix can be used as input to other output commands.
	</P>
	<H5><A NAME = "compute"></A>Computes that generate values to output
	</H5>
	<P>Every <A HREF = "compute.html">compute</A> in LAMMPS produces either global or
	per-atom or local values. The values can be scalars or vectors or
	arrays of data. These values can be output using the other commands
	described in this section. The doc page for each compute command
	describes what it produces. Computes that produce per-atom or local
	values have the word "atom" or "local" in their style name. Computes
	without the word "atom" or "local" produce global values.
	</P>
	<H5><A NAME = "fix"></A>Fixes that generate values to output
	</H5>
	<P>Some <A HREF = "fix.html">fixes</A> in LAMMPS produces either global or per-atom or
	local values which can be accessed by other commands. The values can
	be scalars or vectors or arrays of data. These values can be output
	using the other commands described in this section. The doc page for
	each fix command tells whether it produces any output quantities and
	describes them.
	</P>
	<H5><A NAME = "variable"></A>Variables that generate values to output
	</H5>
	<P>Every <A HREF = "variable.html">variables</A> defined in an input script generates
	either a global scalar value or a per-atom vector (only atom-style
	variables) when it is accessed. The formulas used to define equal-
	and atom-style variables can contain references to the thermodynamic
	keywords and to global and per-atom data generated by computes, fixes,
	and other variables. The values generated by variables can be output
	using the other commands described in this section.
	</P>
	<H5><A NAME = "table"></A>Summary table of output options and data flow between commands
	</H5>
	<P>This table summarizes the various commands that can be used for
	generating output from LAMMPS. Each command produces output data of
	some kind and/or writes data to a file. Most of the commands can take
	data from other commands as input. Thus you can link many of these
	commands together in pipeline form, where data produced by one command
	is used as input to another command and eventually written to the
	screen or to a file. Note that to hook two commands together the
	output and input data types must match, e.g. global/per-atom/local
	data and scalar/vector/array data.
	</P>
	<P>Also note that, as described above, when a command takes a scalar as
	input, that could be an element of a vector or array. Likewise a
	vector input could be a column of an array.
	</P>
	<DIV ALIGN=center><TABLE BORDER=1 >
	<TR><TD >Command</TD><TD > Input</TD><TD > Output</TD><TD ></TD></TR>
	<TR><TD ><A HREF = "thermo_style.html">thermo_style custom</A></TD><TD > global scalars</TD><TD > screen, log file</TD><TD ></TD></TR>
	<TR><TD ><A HREF = "dump.html">dump custom</A></TD><TD > per-atom vectors</TD><TD > dump file</TD><TD ></TD></TR>
	<TR><TD ><A HREF = "dump.html">dump local</A></TD><TD > local vectors</TD><TD > dump file</TD><TD ></TD></TR>
	<TR><TD ><A HREF = "fix_print.html">fix print</A></TD><TD > global scalar from variable</TD><TD > screen, file</TD><TD ></TD></TR>
	<TR><TD ><A HREF = "print.html">print</A></TD><TD > global scalar from variable</TD><TD > screen</TD><TD ></TD></TR>
	<TR><TD ><A HREF = "compute.html">computes</A></TD><TD > N/A</TD><TD > global/per-atom/local scalar/vector/array</TD><TD ></TD></TR>
	<TR><TD ><A HREF = "fix.html">fixes</A></TD><TD > N/A</TD><TD > global/per-atom/local scalar/vector/array</TD><TD ></TD></TR>
	<TR><TD ><A HREF = "variable.html">variables</A></TD><TD > global scalars, per-atom vectors</TD><TD > global scalar, per-atom vector</TD><TD ></TD></TR>
	<TR><TD ><A HREF = "compute_reduce.html">compute reduce</A></TD><TD > per-atom/local vectors</TD><TD > global scalar/vector</TD><TD ></TD></TR>
	<TR><TD ><A HREF = "compute_slice.html">compute slice</A></TD><TD > global vectors/arrays</TD><TD > global vector/array</TD><TD ></TD></TR>
	<TR><TD ><A HREF = "compute_property_atom.html">compute property/atom</A></TD><TD > per-atom vectors</TD><TD > per-atom vector/array</TD><TD ></TD></TR>
	<TR><TD ><A HREF = "compute_property_local.html">compute property/local</A></TD><TD > local vectors</TD><TD > local vector/array</TD><TD ></TD></TR>
	<TR><TD ><A HREF = "compute_atom_molecule.html">compute atom/molecule</A></TD><TD > per-atom vectors</TD><TD > global vector/array</TD><TD ></TD></TR>
	<TR><TD ><A HREF = "fix_vector.html">fix vector</A></TD><TD > global scalars</TD><TD > global vector</TD><TD ></TD></TR>
	<TR><TD ><A HREF = "fix_ave_atom.html">fix ave/atom</A></TD><TD > per-atom vectors</TD><TD > per-atom vector/array</TD><TD ></TD></TR>
	<TR><TD ><A HREF = "fix_ave_time.html">fix ave/time</A></TD><TD > global scalars/vectors</TD><TD > global scalar/vector/array, file</TD><TD ></TD></TR>
	<TR><TD ><A HREF = "fix_ave_spatial.html">fix ave/spatial</A></TD><TD > per-atom vectors</TD><TD > global array, file</TD><TD ></TD></TR>
	<TR><TD ><A HREF = "fix_ave_histo.html">fix ave/histo</A></TD><TD > global/per-atom/local scalars and vectors</TD><TD > global array, file</TD><TD ></TD></TR>
	<TR><TD ><A HREF = "fix_ave_correlate.html">fix ave/correlate</A></TD><TD > global scalars</TD><TD > global array, file</TD><TD ></TD></TR>
	<TR><TD ><A HREF = "fix_store_state.html">fix store/state</A></TD><TD > per-atom vectors</TD><TD > per-atom vector/array</TD><TD ></TD></TR>
	<TR><TD >
	</TD></TR></TABLE></DIV>

	<HR>

	<A NAME = "howto_16"></A><H4>6.16 Thermostatting, barostatting, and computing temperature
	</H4>
	<P>Thermostatting means controlling the temperature of particles in an MD
	simulation. Barostatting means controlling the pressure. Since the
	pressure includes a kinetic component due to particle velocities, both
	these operations require calculation of the temperature. Typically a
	target temperature (T) and/or pressure (P) is specified by the user,
	and the thermostat or barostat attempts to equilibrate the system to
	the requested T and/or P.
	</P>
	<P>Temperature is computed as kinetic energy divided by some number of
	degrees of freedom (and the Boltzmann constant). Since kinetic energy
	is a function of particle velocity, there is often a need to
	distinguish between a particle's advection velocity (due to some
	aggregate motiion of particles) and its thermal velocity. The sum of
	the two is the particle's total velocity, but the latter is often what
	is wanted to compute a temperature.
	</P>
	<P>LAMMPS has several options for computing temperatures, any of which
	can be used in thermostatting and barostatting. These <A HREF = "compute.html">compute
	commands</A> calculate temperature, and the <A HREF = "compute_pressure.html">compute
	pressure</A> command calculates pressure.
	</P>
	<UL><LI><A HREF = "compute_temp.html">compute temp</A>
	<LI><A HREF = "compute_temp_sphere.html">compute temp/sphere</A>
	<LI><A HREF = "compute_temp_asphere.html">compute temp/asphere</A>
	<LI><A HREF = "compute_temp_com.html">compute temp/com</A>
	<LI><A HREF = "compute_temp_deform.html">compute temp/deform</A>
	<LI><A HREF = "compute_temp_partial.html">compute temp/partial</A>
	<LI><A HREF = "compute_temp_profile.html">compute temp/profile</A>
	<LI><A HREF = "compute_temp_ramp.html">compute temp/ramp</A>
	<LI><A HREF = "compute_temp_region.html">compute temp/region</A>
	</UL>
	<P>All but the first 3 calculate velocity biases (i.e. advection
	velocities) that are removed when computing the thermal temperature.
	<A HREF = "compute_temp_sphere.html">Compute temp/sphere</A> and <A HREF = "compute_temp_asphere.html">compute
	temp/asphere</A> compute kinetic energy for
	finite-size particles that includes rotational degrees of freedom.
	They both allow, as an extra argument, which is another temperature
	compute that subtracts a velocity bias. This allows the translational
	velocity of spherical or aspherical particles to be adjusted in
	prescribed ways.
	</P>
	<P>Thermostatting in LAMMPS is performed by <A HREF = "fix.html">fixes</A>, or in one
	case by a pair style. Several thermostatting fixes are available:
	Nose-Hoover (nvt), Berendsen, CSVR, Langevin, and direct rescaling
	(temp/rescale). Dissipative particle dynamics (DPD) thermostatting
	can be invoked via the <I>dpd/tstat</I> pair style:
	</P>
	<UL><LI><A HREF = "fix_nh.html">fix nvt</A>
	<LI><A HREF = "fix_nvt_sphere.html">fix nvt/sphere</A>
	<LI><A HREF = "fix_nvt_asphere.html">fix nvt/asphere</A>
	<LI><A HREF = "fix_nvt_sllod.html">fix nvt/sllod</A>
	<LI><A HREF = "fix_temp_berendsen.html">fix temp/berendsen</A>
	<LI><A HREF = "fix_temp_csvr.html">fix temp/csvr</A>
	<LI><A HREF = "fix_langevin.html">fix langevin</A>
	<LI><A HREF = "fix_temp_rescale.html">fix temp/rescale</A>
	<LI><A HREF = "pair_dpd.html">pair_style dpd/tstat</A>
	</UL>
	<P><A HREF = "fix_nh.html">Fix nvt</A> only thermostats the translational velocity of
	particles. <A HREF = "fix_nvt_sllod.html">Fix nvt/sllod</A> also does this, except
	that it subtracts out a velocity bias due to a deforming box and
	integrates the SLLOD equations of motion. See the <A HREF = "#howto_13">NEMD
	simulations</A> section of this page for further details. <A HREF = "fix_nvt_sphere.html">Fix
	nvt/sphere</A> and <A HREF = "fix_nvt_asphere.html">fix
	nvt/asphere</A> thermostat not only translation
	velocities but also rotational velocities for spherical and aspherical
	particles.
	</P>
	<P>DPD thermostatting alters pairwise interactions in a manner analagous
	to the per-particle thermostatting of <A HREF = "fix_langevin.html">fix
	langevin</A>.
	</P>
	<P>Any of the thermostatting fixes can use temperature computes that
	remove bias for two purposes: (a) computing the current temperature to
	compare to the requested target temperature, and (b) adjusting only
	the thermal temperature component of the particle's velocities. See
	the doc pages for the individual fixes and for the
	<A HREF = "fix_modify.html">fix_modify</A> command for instructions on how to assign
	a temperature compute to a thermostatting fix. For example, you can
	apply a thermostat to only the x and z components of velocity by using
	it in conjunction with <A HREF = "compute_temp_partial.html">compute
	temp/partial</A>.
	</P>
	<P>IMPORTANT NOTE: Only the nvt fixes perform time integration, meaning
	they update the velocities and positions of particles due to forces
	and velocities respectively. The other thermostat fixes only adjust
	velocities; they do NOT perform time integration updates. Thus they
	should be used in conjunction with a constant NVE integration fix such
	as these:
	</P>
	<UL><LI><A HREF = "fix_nve.html">fix nve</A>
	<LI><A HREF = "fix_nve_sphere.html">fix nve/sphere</A>
	<LI><A HREF = "fix_nve_asphere.html">fix nve/asphere</A>
	</UL>
	<P>Barostatting in LAMMPS is also performed by <A HREF = "fix.html">fixes</A>. Two
	barosttating methods are currently available: Nose-Hoover (npt and
	nph) and Berendsen:
	</P>
	<UL><LI><A HREF = "fix_nh.html">fix npt</A>
	<LI><A HREF = "fix_npt_sphere.html">fix npt/sphere</A>
	<LI><A HREF = "fix_npt_asphere.html">fix npt/asphere</A>
	<LI><A HREF = "fix_nh.html">fix nph</A>
	<LI><A HREF = "fix_press_berendsen.html">fix press/berendsen</A>
	</UL>
	<P>The <A HREF = "fix_nh.html">fix npt</A> commands include a Nose-Hoover thermostat
	and barostat. <A HREF = "fix_nh.html">Fix nph</A> is just a Nose/Hoover barostat;
	it does no thermostatting. Both <A HREF = "fix_nh.html">fix nph</A> and <A HREF = "fix_press_berendsen.html">fix
	press/bernendsen</A> can be used in conjunction
	with any of the thermostatting fixes.
	</P>
	<P>As with the thermostats, <A HREF = "fix_nh.html">fix npt</A> and <A HREF = "fix_nh.html">fix
	nph</A> only use translational motion of the particles in
	computing T and P and performing thermo/barostatting. <A HREF = "fix_npt_sphere.html">Fix
	npt/sphere</A> and <A HREF = "fix_npt_asphere.html">fix
	npt/asphere</A> thermo/barostat using not only
	translation velocities but also rotational velocities for spherical
	and aspherical particles.
	</P>
	<P>All of the barostatting fixes use the <A HREF = "compute_pressure.html">compute
	pressure</A> compute to calculate a current
	pressure. By default, this compute is created with a simple <A HREF = "compute_temp.html">compute
	temp</A> (see the last argument of the <A HREF = "compute_pressure.html">compute
	pressure</A> command), which is used to calculated
	the kinetic componenet of the pressure. The barostatting fixes can
	also use temperature computes that remove bias for the purpose of
	computing the kinetic componenet which contributes to the current
	pressure. See the doc pages for the individual fixes and for the
	<A HREF = "fix_modify.html">fix_modify</A> command for instructions on how to assign
	a temperature or pressure compute to a barostatting fix.
	</P>
	<P>IMPORTANT NOTE: As with the thermostats, the Nose/Hoover methods (<A HREF = "fix_nh.html">fix
	npt</A> and <A HREF = "fix_nh.html">fix nph</A>) perform time
	integration. <A HREF = "fix_press_berendsen.html">Fix press/berendsen</A> does NOT,
	so it should be used with one of the constant NVE fixes or with one of
	the NVT fixes.
	</P>
	<P>Finally, thermodynamic output, which can be setup via the
	<A HREF = "thermo_style.html">thermo_style</A> command, often includes temperature
	and pressure values. As explained on the doc page for the
	<A HREF = "thermo_style.html">thermo_style</A> command, the default T and P are
	setup by the thermo command itself. They are NOT the ones associated
	with any thermostatting or barostatting fix you have defined or with
	any compute that calculates a temperature or pressure. Thus if you
	want to view these values of T and P, you need to specify them
	explicitly via a <A HREF = "thermo_style.html">thermo_style custom</A> command. Or
	you can use the <A HREF = "thermo_modify.html">thermo_modify</A> command to
	re-define what temperature or pressure compute is used for default
	thermodynamic output.
	</P>
	<HR>

	<A NAME = "howto_17"></A><H4>6.17 Walls
	</H4>
	<P>Walls in an MD simulation are typically used to bound particle motion,
	i.e. to serve as a boundary condition.
	</P>
	<P>Walls in LAMMPS can be of rough (made of particles) or idealized
	surfaces. Ideal walls can be smooth, generating forces only in the
	normal direction, or frictional, generating forces also in the
	tangential direction.
	</P>
	<P>Rough walls, built of particles, can be created in various ways. The
	particles themselves can be generated like any other particle, via the
	<A HREF = "lattice.html">lattice</A> and <A HREF = "create_atoms.html">create_atoms</A> commands,
	or read in via the <A HREF = "read_data.html">read_data</A> command.
	</P>
	<P>Their motion can be constrained by many different commands, so that
	they do not move at all, move together as a group at constant velocity
	or in response to a net force acting on them, move in a prescribed
	fashion (e.g. rotate around a point), etc. Note that if a time
	integration fix like <A HREF = "fix_nve.html">fix nve</A> or <A HREF = "fix_nh.html">fix nvt</A>
	is not used with the group that contains wall particles, their
	positions and velocities will not be updated.
	</P>
	<UL><LI><A HREF = "fix_aveforce.html">fix aveforce</A> - set force on particles to average value, so they move together
	<LI><A HREF = "fix_setforce.html">fix setforce</A> - set force on particles to a value, e.g. 0.0
	<LI><A HREF = "fix_freeze.html">fix freeze</A> - freeze particles for use as granular walls
	<LI><A HREF = "fix_nve_noforce.html">fix nve/noforce</A> - advect particles by their velocity, but without force
	<LI><A HREF = "fix_move.html">fix move</A> - prescribe motion of particles by a linear velocity, oscillation, rotation, variable
	</UL>
	<P>The <A HREF = "fix_move.html">fix move</A> command offers the most generality, since
	the motion of individual particles can be specified with
	<A HREF = "variable.html">variable</A> formula which depends on time and/or the
	particle position.
	</P>
	<P>For rough walls, it may be useful to turn off pairwise interactions
	between wall particles via the <A HREF = "neigh_modify.html">neigh_modify
	exclude</A> command.
	</P>
	<P>Rough walls can also be created by specifying frozen particles that do
	not move and do not interact with mobile particles, and then tethering
	other particles to the fixed particles, via a <A HREF = "bond_style.html">bond</A>.
	The bonded particles do interact with other mobile particles.
	</P>
	<P>Idealized walls can be specified via several fix commands. <A HREF = "fix_wall_gran.html">Fix
	wall/gran</A> creates frictional walls for use with
	granular particles; all the other commands create smooth walls.
	</P>
	<UL><LI><A HREF = "fix_wall_reflect.html">fix wall/reflect</A> - reflective flat walls
	<LI><A HREF = "fix_wall.html">fix wall/lj93</A> - flat walls, with Lennard-Jones 9/3 potential
	<LI><A HREF = "fix_wall.html">fix wall/lj126</A> - flat walls, with Lennard-Jones 12/6 potential
	<LI><A HREF = "fix_wall.html">fix wall/colloid</A> - flat walls, with <A HREF = "pair_colloid.html">pair_style colloid</A> potential
	<LI><A HREF = "fix_wall.html">fix wall/harmonic</A> - flat walls, with repulsive harmonic spring potential
	<LI><A HREF = "fix_wall_region.html">fix wall/region</A> - use region surface as wall
	<LI><A HREF = "fix_wall_gran.html">fix wall/gran</A> - flat or curved walls with <A HREF = "pair_gran.html">pair_style granular</A> potential
	</UL>
	<P>The <I>lj93</I>, <I>lj126</I>, <I>colloid</I>, and <I>harmonic</I> styles all allow the
	flat walls to move with a constant velocity, or oscillate in time.
	The <A HREF = "fix_wall_region.html">fix wall/region</A> command offers the most
	generality, since the region surface is treated as a wall, and the
	geometry of the region can be a simple primitive volume (e.g. a
	sphere, or cube, or plane), or a complex volume made from the union
	and intersection of primitive volumes. <A HREF = "region.html">Regions</A> can also
	specify a volume "interior" or "exterior" to the specified primitive
	shape or <I>union</I> or <I>intersection</I>. <A HREF = "region.html">Regions</A> can also be
	"dynamic" meaning they move with constant velocity, oscillate, or
	rotate.
	</P>
	<P>The only frictional idealized walls currently in LAMMPS are flat or
	curved surfaces specified by the <A HREF = "fix_wall_gran.html">fix wall/gran</A>
	command. At some point we plan to allow regoin surfaces to be used as
	frictional walls, as well as triangulated surfaces.
	</P>
	<HR>

	<A NAME = "howto_18"></A><H4>6.18 Elastic constants
	</H4>
	<P>Elastic constants characterize the stiffness of a material. The formal
	definition is provided by the linear relation that holds between the
	stress and strain tensors in the limit of infinitesimal deformation.
	In tensor notation, this is expressed as s_ij = C_ijkl * e_kl, where
	the repeated indices imply summation. s_ij are the elements of the
	symmetric stress tensor. e_kl are the elements of the symmetric strain
	tensor. C_ijkl are the elements of the fourth rank tensor of elastic
	constants. In three dimensions, this tensor has 3^4=81 elements. Using
	Voigt notation, the tensor can be written as a 6x6 matrix, where C_ij
	is now the derivative of s_i w.r.t. e_j. Because s_i is itself a
	derivative w.r.t. e_i, it follows that C_ij is also symmetric, with at
	most 7*6/2 = 21 distinct elements.
	</P>
	<P>At zero temperature, it is easy to estimate these derivatives by
	deforming the simulation box in one of the six directions using the
	<A HREF = "change_box.html">change_box</A> command and measuring the change in the
	stress tensor. A general-purpose script that does this is given in the
	examples/elastic directory described in <A HREF = "Section_example.html">this
	section</A>.
	</P>
	<P>Calculating elastic constants at finite temperature is more
	challenging, because it is necessary to run a simulation that perfoms
	time averages of differential properties. One way to do this is to
	measure the change in average stress tensor in an NVT simulations when
	the cell volume undergoes a finite deformation. In order to balance
	the systematic and statistical errors in this method, the magnitude of
	the deformation must be chosen judiciously, and care must be taken to
	fully equilibrate the deformed cell before sampling the stress
	tensor. Another approach is to sample the triclinic cell fluctuations
	that occur in an NPT simulation. This method can also be slow to
	converge and requires careful post-processing <A HREF = "#Shinoda">(Shinoda)</A>
	</P>
	<HR>

	<A NAME = "howto_19"></A><H4>6.19 Library interface to LAMMPS
	</H4>
	<P>As described in <A HREF = "Section_start.html#start_5">Section_start 5</A>, LAMMPS
	can be built as a library, so that it can be called by another code,
	used in a <A HREF = "Section_howto.html#howto_10">coupled manner</A> with other
	codes, or driven through a <A HREF = "Section_python.html">Python interface</A>.
	</P>
	<P>All of these methodologies use a C-style interface to LAMMPS that is
	provided in the files src/library.cpp and src/library.h. The
	functions therein have a C-style argument list, but contain C++ code
	you could write yourself in a C++ application that was invoking LAMMPS
	directly. The C++ code in the functions illustrates how to invoke
	internal LAMMPS operations. Note that LAMMPS classes are defined
	within a LAMMPS namespace (LAMMPS_NS) if you use them from another C++
	application.
	</P>
	<P>Library.cpp contains these 4 functions:
	</P>
	<PRE>void lammps_open(int, char , MPI_Comm, void )
	void lammps_close(void *)
	void lammps_file(void , char )
	char lammps_command(void , char *)
	</PRE>
	<P>The lammps_open() function is used to initialize LAMMPS, passing in a
	list of strings as if they were <A HREF = "Section_start.html#start_7">command-line
	arguments</A> when LAMMPS is run in
	stand-alone mode from the command line, and a MPI communicator for
	LAMMPS to run under. It returns a ptr to the LAMMPS object that is
	created, and which is used in subsequent library calls. The
	lammps_open() function can be called multiple times, to create
	multiple instances of LAMMPS.
	</P>
	<P>LAMMPS will run on the set of processors in the communicator. This
	means the calling code can run LAMMPS on all or a subset of
	processors. For example, a wrapper script might decide to alternate
	between LAMMPS and another code, allowing them both to run on all the
	processors. Or it might allocate half the processors to LAMMPS and
	half to the other code and run both codes simultaneously before
	syncing them up periodically. Or it might instantiate multiple
	instances of LAMMPS to perform different calculations.
	</P>
	<P>The lammps_close() function is used to shut down an instance of LAMMPS
	and free all its memory.
	</P>
	<P>The lammps_file() and lammps_command() functions are used to pass a
	file or string to LAMMPS as if it were an input script or single
	command in an input script. Thus the calling code can read or
	generate a series of LAMMPS commands one line at a time and pass it
	thru the library interface to setup a problem and then run it,
	interleaving the lammps_command() calls with other calls to extract
	information from LAMMPS, perform its own operations, or call another
	code's library.
	</P>
	<P>Other useful functions are also included in library.cpp. For example:
	</P>
	<PRE>void lammps_extract_global(void , char *)
	void lammps_extract_atom(void , char *)
	void lammps_extract_compute(void , char *, int, int)
	void lammps_extract_fix(void , char *, int, int, int, int)
	void lammps_extract_variable(void , char , char )
	int lammps_get_natoms(void *)
	void lammps_get_coords(void , double )
	void lammps_put_coords(void , double )
	</PRE>
	<P>These can extract various global or per-atom quantities from LAMMPS as
	well as values calculated by a compute, fix, or variable. The "get"
	and "put" operations can retrieve and reset atom coordinates.
	See the library.cpp file and its associated header file library.h for
	details.
	</P>
	<P>The key idea of the library interface is that you can write any
	functions you wish to define how your code talks to LAMMPS and add
	them to src/library.cpp and src/library.h, as well as to the <A HREF = "Section_python.html">Python
	interface</A>. The routines you add can access or
	change any LAMMPS data you wish. The examples/COUPLE and python
	directories have example C++ and C and Python codes which show how a
	driver code can link to LAMMPS as a library, run LAMMPS on a subset of
	processors, grab data from LAMMPS, change it, and put it back into
	LAMMPS.
	</P>
	<HR>

	<A NAME = "howto_20"></A><H4>6.20 Calculating thermal conductivity
	</H4>
	<P>The thermal conductivity kappa of a material can be measured in at
	least 4 ways using various options in LAMMPS. See the examples/KAPPA
	directory for scripts that implement the 4 methods discussed here for
	a simple Lennard-Jones fluid model. Also, see <A HREF = "Section_howto.html#howto_21">this
	section</A> of the manual for an analogous
	discussion for viscosity.
	</P>
	<P>The thermal conducitivity tensor kappa is a measure of the propensity
	of a material to transmit heat energy in a diffusive manner as given
	by Fourier's law
	</P>
	<P>J = -kappa grad(T)
	</P>
	<P>where J is the heat flux in units of energy per area per time and
	grad(T) is the spatial gradient of temperature. The thermal
	conductivity thus has units of energy per distance per time per degree
	K and is often approximated as an isotropic quantity, i.e. as a
	scalar.
	</P>
	<P>The first method is to setup two thermostatted regions at opposite
	ends of a simulation box, or one in the middle and one at the end of a
	periodic box. By holding the two regions at different temperatures
	with a <A HREF = "Section_howto.html#howto_13">thermostatting fix</A>, the energy
	added to the hot region should equal the energy subtracted from the
	cold region and be proportional to the heat flux moving between the
	regions. See the paper by <A HREF = "#Ikeshoji">Ikeshoji and Hafskjold</A> for
	details of this idea. Note that thermostatting fixes such as <A HREF = "fix_nh.html">fix
	nvt</A>, <A HREF = "fix_langevin.html">fix langevin</A>, and <A HREF = "fix_temp_rescale.html">fix
	temp/rescale</A> store the cumulative energy they
	add/subtract.
	</P>
	<P>Alternatively, as a second method, the <A HREF = "fix_heat.html">fix heat</A>
	command can used in place of thermostats on each of two regions to
	add/subtract specified amounts of energy to both regions. In both
	cases, the resulting temperatures of the two regions can be monitored
	with the "compute temp/region" command and the temperature profile of
	the intermediate region can be monitored with the <A HREF = "fix_ave_spatial.html">fix
	ave/spatial</A> and <A HREF = "compute_ke_atom.html">compute
	ke/atom</A> commands.
	</P>
	<P>The third method is to perform a reverse non-equilibrium MD simulation
	using the <A HREF = "fix_thermal_conductivity.html">fix thermal/conductivity</A>
	command which implements the rNEMD algorithm of Muller-Plathe.
	Kinetic energy is swapped between atoms in two different layers of the
	simulation box. This induces a temperature gradient between the two
	layers which can be monitored with the <A HREF = "fix_ave_spatial.html">fix
	ave/spatial</A> and <A HREF = "compute_ke_atom.html">compute
	ke/atom</A> commands. The fix tallies the
	cumulative energy transfer that it performs. See the <A HREF = "fix_thermal_conductivity.html">fix
	thermal/conductivity</A> command for
	details.
	</P>
	<P>The fourth method is based on the Green-Kubo (GK) formula which
	relates the ensemble average of the auto-correlation of the heat flux
	to kappa. The heat flux can be calculated from the fluctuations of
	per-atom potential and kinetic energies and per-atom stress tensor in
	a steady-state equilibrated simulation. This is in contrast to the
	two preceding non-equilibrium methods, where energy flows continuously
	between hot and cold regions of the simulation box.
	</P>
	<P>The <A HREF = "compute_heat_flux.html">compute heat/flux</A> command can calculate
	the needed heat flux and describes how to implement the Green_Kubo
	formalism using additional LAMMPS commands, such as the <A HREF = "fix_ave_correlate.html">fix
	ave/correlate</A> command to calculate the needed
	auto-correlation. See the doc page for the <A HREF = "compute_heat_flux.html">compute
	heat/flux</A> command for an example input script
	that calculates the thermal conductivity of solid Ar via the GK
	formalism.
	</P>
	<HR>

	<A NAME = "howto_21"></A><H4>6.21 Calculating viscosity
	</H4>
	<P>The shear viscosity eta of a fluid can be measured in at least 4 ways
	using various options in LAMMPS. See the examples/VISCOSITY directory
	for scripts that implement the 4 methods discussed here for a simple
	Lennard-Jones fluid model. Also, see <A HREF = "Section_howto.html#howto_20">this
	section</A> of the manual for an analogous
	discussion for thermal conductivity.
	</P>
	<P>Eta is a measure of the propensity of a fluid to transmit momentum in
	a direction perpendicular to the direction of velocity or momentum
	flow. Alternatively it is the resistance the fluid has to being
	sheared. It is given by
	</P>
	<P>J = -eta grad(Vstream)
	</P>
	<P>where J is the momentum flux in units of momentum per area per time.
	and grad(Vstream) is the spatial gradient of the velocity of the fluid
	moving in another direction, normal to the area through which the
	momentum flows. Viscosity thus has units of pressure-time.
	</P>
	<P>The first method is to perform a non-equlibrium MD (NEMD) simulation
	by shearing the simulation box via the <A HREF = "fix_deform.html">fix deform</A>
	command, and using the <A HREF = "fix_nvt_sllod.html">fix nvt/sllod</A> command to
	thermostat the fluid via the SLLOD equations of motion.
	Alternatively, as a second method, one or more moving walls can be
	used to shear the fluid in between them, again with some kind of
	thermostat that modifies only the thermal (non-shearing) components of
	velocity to prevent the fluid from heating up.
	</P>
	<P>In both cases, the velocity profile setup in the fluid by this
	procedure can be monitored by the <A HREF = "fix_ave_spatial.html">fix
	ave/spatial</A> command, which determines
	grad(Vstream) in the equation above. E.g. the derivative in the
	y-direction of the Vx component of fluid motion or grad(Vstream) =
	dVx/dy. The Pxy off-diagonal component of the pressure or stress
	tensor, as calculated by the <A HREF = "compute_pressure.html">compute pressure</A>
	command, can also be monitored, which is the J term in the equation
	above. See <A HREF = "Section_howto.html#howto_13">this section</A> of the manual
	for details on NEMD simulations.
	</P>
	<P>The third method is to perform a reverse non-equilibrium MD simulation
	using the <A HREF = "fix_viscosity.html">fix viscosity</A> command which implements
	the rNEMD algorithm of Muller-Plathe. Momentum in one dimension is
	swapped between atoms in two different layers of the simulation box in
	a different dimension. This induces a velocity gradient which can be
	monitored with the <A HREF = "fix_ave_spatial.html">fix ave/spatial</A> command.
	The fix tallies the cummulative momentum transfer that it performs.
	See the <A HREF = "fix_viscosity.html">fix viscosity</A> command for details.
	</P>
	<P>The fourth method is based on the Green-Kubo (GK) formula which
	relates the ensemble average of the auto-correlation of the
	stress/pressure tensor to eta. This can be done in a steady-state
	equilibrated simulation which is in contrast to the two preceding
	non-equilibrium methods, where momentum flows continuously through the
	simulation box.
	</P>
	<P>Here is an example input script that calculates the viscosity of
	liquid Ar via the GK formalism:
	</P>
	<PRE># Sample LAMMPS input script for viscosity of liquid Ar
	</PRE>
	<PRE>units real
	variable T equal 86.4956
	variable V equal vol
	variable dt equal 4.0
	variable p equal 400 # correlation length
	variable s equal 5 # sample interval
	variable d equal $p*$s # dump interval
	</PRE>
	<PRE># convert from LAMMPS real units to SI
	</PRE>
	<PRE>variable kB equal 1.3806504e-23 # [J/K/</B> Boltzmann
	variable atm2Pa equal 101325.0
	variable A2m equal 1.0e-10
	variable fs2s equal 1.0e-15
	variable convert equal ${atm2Pa}${atm2Pa}${fs2s}${A2m}${A2m}*${A2m}
	</PRE>
	<PRE># setup problem
	</PRE>
	<PRE>dimension 3
	boundary p p p
	lattice fcc 5.376 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
	region box block 0 4 0 4 0 4
	create_box 1 box
	create_atoms 1 box
	mass 1 39.948
	pair_style lj/cut 13.0
	pair_coeff * * 0.2381 3.405
	timestep ${dt}
	thermo $d
	</PRE>
	<PRE># equilibration and thermalization
	</PRE>
	<PRE>velocity all create $T 102486 mom yes rot yes dist gaussian
	fix NVT all nvt temp $T $T 10 drag 0.2
	run 8000
	</PRE>
	<PRE># viscosity calculation, switch to NVE if desired
	</PRE>
	<PRE>#unfix NVT
	#fix NVE all nve
	</PRE>
	<PRE>reset_timestep 0
	variable pxy equal pxy
	variable pxz equal pxz
	variable pyz equal pyz
	fix SS all ave/correlate $s $p $d &
	v_pxy v_pxz v_pyz type auto file S0St.dat ave running
	variable scale equal ${convert}/(${kB}$T)$V$s${dt}
	variable v11 equal trap(f_SS[3])*${scale}
	variable v22 equal trap(f_SS[4])*${scale}
	variable v33 equal trap(f_SS[5])*${scale}
	thermo_style custom step temp press v_pxy v_pxz v_pyz v_v11 v_v22 v_v33
	run 100000
	variable v equal (v_v11+v_v22+v_v33)/3.0
	variable ndens equal count(all)/vol
	print "average viscosity: $v [Pa.s/</B> @ $T K, ${ndens} /A^3"
	</PRE>
	<HR>

	<A NAME = "howto_22"></A><H4>6.22 Calculating a diffusion coefficient
	</H4>
	<P>The diffusion coefficient D of a material can be measured in at least
	2 ways using various options in LAMMPS. See the examples/DIFFUSE
	directory for scripts that implement the 2 methods discussed here for
	a simple Lennard-Jones fluid model.
	</P>
	<P>The first method is to measure the mean-squared displacement (MSD) of
	the system, via the <A HREF = "compute_msd.html">compute msd</A> command. The slope
	of the MSD versus time is proportional to the diffusion coefficient.
	The instantaneous MSD values can be accumulated in a vector via the
	<A HREF = "fix_vector.html">fix vector</A> command, and a line fit to the vector to
	compute its slope via the <A HREF = "variable.html">variable slope</A> function, and
	thus extract D.
	</P>
	<P>The second method is to measure the velocity auto-correlation function
	(VACF) of the system, via the <A HREF = "compute_vacf.html">compute vacf</A>
	command. The time-integral of the VACF is proportional to the
	diffusion coefficient. The instantaneous VACF values can be
	accumulated in a vector via the <A HREF = "fix_vector.html">fix vector</A> command,
	and time integrated via the <A HREF = "variable.html">variable trap</A> function,
	and thus extract D.
	</P>
	<HR>

	<HR>

	<A NAME = "Berendsen"></A>

	<P><B>(Berendsen)</B> Berendsen, Grigera, Straatsma, J Phys Chem, 91,
	6269-6271 (1987).
	</P>
	<A NAME = "Cornell"></A>

	<P><B>(Cornell)</B> Cornell, Cieplak, Bayly, Gould, Merz, Ferguson,
	Spellmeyer, Fox, Caldwell, Kollman, JACS 117, 5179-5197 (1995).
	</P>
	<A NAME = "Horn"></A>

	<P><B>(Horn)</B> Horn, Swope, Pitera, Madura, Dick, Hura, and Head-Gordon,
	J Chem Phys, 120, 9665 (2004).
	</P>
	<A NAME = "Ikeshoji"></A>

	<P><B>(Ikeshoji)</B> Ikeshoji and Hafskjold, Molecular Physics, 81, 251-261
	(1994).
	</P>
	<A NAME = "MacKerell"></A>

	<P><B>(MacKerell)</B> MacKerell, Bashford, Bellott, Dunbrack, Evanseck, Field,
	Fischer, Gao, Guo, Ha, et al, J Phys Chem, 102, 3586 (1998).
	</P>
	<A NAME = "Mayo"></A>

	<P><B>(Mayo)</B> Mayo, Olfason, Goddard III, J Phys Chem, 94, 8897-8909
	(1990).
	</P>
	<A NAME = "Jorgensen"></A>

	<P><B>(Jorgensen)</B> Jorgensen, Chandrasekhar, Madura, Impey, Klein, J Chem
	Phys, 79, 926 (1983).
	</P>
	<A NAME = "Price"></A>

	<P><B>(Price)</B> Price and Brooks, J Chem Phys, 121, 10096 (2004).
	</P>
	<A NAME = "Shinoda"></A>

	<P><B>(Shinoda)</B> Shinoda, Shiga, and Mikami, Phys Rev B, 69, 134103 (2004).
	</P>
	</HTML>

Section_howto.htmlNo OneTemporaryActions

File Metadata

Section_howto.htmlView Options

Event Timeline

Section_howto.html
No OneTemporary
Actions

Section_howto.html
View Options