Section_howto.html
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Wed, Jul 10, 12:41

Section_howto.html
View Options

	<HTML>
	<CENTER><A HREF = "Section_commands.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>4. How-to discussions
	</H3>
	<P>The following sections describe what commands can be used to perform
	certain kinds of LAMMPS simulations.
	</P>
	4.1 <A HREF = "#4_1">Restarting a simulation</A><BR>
	4.2 <A HREF = "#4_2">2d simulations</A><BR>
	4.3 <A HREF = "#4_3">CHARMM and AMBER force fields</A><BR>
	4.4 <A HREF = "#4_4">Running multiple simulations from one input script</A><BR>
	4.5 <A HREF = "#4_5">Parallel tempering</A><BR>
	4.6 <A HREF = "#4_6">Granular models</A><BR>
	4.7 <A HREF = "#4_7">TIP3P water model</A><BR>
	4.8 <A HREF = "#4_8">TIP4P water model</A><BR>
	4.9 <A HREF = "#4_9">SPC water model</A><BR>
	4.10 <A HREF = "#4_10">Coupling LAMMPS to other codes</A><BR>
	4.11 <A HREF = "#4_11">Visualizing LAMMPS snapshots</A><BR>
	4.12 <A HREF = "#4_12">Non-orthogonal simulation boxes</A><BR>
	4.13 <A HREF = "#4_13">NEMD simulations</A><BR>
	4.14 <A HREF = "#4_14">Aspherical particles</A><BR>
	4.15 <A HREF = "#4_15">Output from LAMMPS</A> <BR>

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

	<A NAME = "4_1"></A><H4>4.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 and read by a
	<A HREF = "read_data.html">read_data</A> command in a new script. <A HREF = "Section_tools.html">This
	section</A> discusses the <I>restart2data</I> tool that is
	used to perform the conversion.
	</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 using this tool:
	</P>
	<PRE>restart2data 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 = "4_2"></A><H4>4.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>
	<HR>

	<A NAME = "4_3"></A><H4>4.3 CHARMM and AMBER force fields
	</H4>
	<P>There are many different ways to compute forces in the <A HREF = "http://www.scripps.edu/brooks">CHARMM</A>
	and <A HREF = "http://amber.scripps.edu">AMBER</A> molecular dynamics codes, only some of which are
	available as options in LAMMPS. 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. 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">this section</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_style.html">bond_style</A> harmonic
	<LI><A HREF = "angle_style.html">angle_style</A> charmm
	<LI><A HREF = "dihedral_style.html">dihedral_style</A> charmm
	<LI><A HREF = "pair_style.html">pair_style</A> lj/charmm/coul/charmm
	<LI><A HREF = "pair_style.html">pair_style</A> lj/charmm/coul/charmm/implicit
	<LI><A HREF = "pair_style.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>
	<HR>

	<A NAME = "4_4"></A><H4>4.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
	cd $d
	read_data data.polymer
	run 10000
	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
	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#2_6">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 = "4_5"></A><H4>4.5 Parallel tempering
	</H4>
	<P>The <A HREF = "temper.html">temper</A> command can be used to perform a parallel
	tempering or replica-exchange simulation where multiple copies of a
	simulation are run at different temperatures on different sets of
	processors, and Monte Carlo temperature swaps are performed between
	pairs of copies.
	</P>
	<P>Use the -procs and -in <A HREF = "Section_start.html#2_6">command-line switches</A>
	to launch LAMMPS on multiple partitions.
	</P>
	<P>In your input script, define a set of temperatures, one for each
	processor partition, using the <A HREF = "variable.html">variable</A> command:
	</P>
	<PRE>variable t proc 300.0 310.0 320.0 330.0
	</PRE>
	<P>Define a fix of style <A HREF = "fix_nvt.html">nvt</A> or <A HREF = "fix_langevin.html">langevin</A>
	to control the temperature of each simulation:
	</P>
	<PRE>fix myfix all nvt $t $t 100.0
	</PRE>
	<P>Use the <A HREF = "temper.html">temper</A> command in place of a <A HREF = "run.html">run</A>
	command to perform a simulation where tempering exchanges will take
	place:
	</P>
	<PRE>temper 100000 100 $t myfix 3847 58382
	</PRE>
	<HR>

	<A NAME = "4_6"></A><H4>4.6 Granular models
	</H4>
	<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</A> granular
	<LI><A HREF = "fix_nve_gran.html">fix nve/gran</A>
	<LI><A HREF = "fix_gravity.html">fix gravity</A>
	<LI><A HREF = "thermo_style.html">thermo_style</A> gran
	</UL>
	<P>Use one of these 3 pair potentials:
	</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_gran_diag.html">fix gran/diag</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 = "4_7"></A><H4>4.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>
	</P>
	<P>O charge = -0.834<BR>
	H charge = 0.417 <BR>
	</P>
	<P>LJ epsilon of OO = 0.1521<BR>
	LJ sigma of OO = 3.188<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>
	</P>
	<P>K of OH bond = 450<BR>
	r0 of OH bond = 0.9572 <BR>
	</P>
	<P>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 (Ewald or PPPM in LAMMPS):
	</P>
	<P>O mass = 15.9994<BR>
	H mass = 1.008 <BR>
	</P>
	<P>O charge = -0.830<BR>
	H charge = 0.415 <BR>
	</P>
	<P>LJ epsilon of OO = 0.102<BR>
	LJ sigma of OO = 3.1507<BR>
	LJ epsilon, sigma of OH, HH = 0.0 <BR>
	</P>
	<P>K of OH bond = 450<BR>
	r0 of OH bond = 0.9572 <BR>
	</P>
	<P>K of HOH angle = 55<BR>
	theta of HOH angle = 104.52 <BR>
	</P>
	<HR>

	<A NAME = "4_8"></A><H4>4.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>Currently, only a four-point model for long-range Coulombics is
	implemented via the LAMMPS <A HREF = "pair_lj.html">pair style
	lj/cut/coul/long/tip4p</A>. We plan to add a cutoff
	version in the future. 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>
	</P>
	<P>O charge = -1.040<BR>
	H charge = 0.520 <BR>
	</P>
	<P>r0 of OH bond = 0.9572<BR>
	theta of HOH angle = 104.52 <BR>
	</P>
	<P>OM distance = 0.15 <BR>
	</P>
	<P>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>
	</P>
	<P>These are the parameters to use for TIP4P with a long-range Coulombic
	solver (Ewald or PPPM in LAMMPS):
	</P>
	<P>O mass = 15.9994<BR>
	H mass = 1.008 <BR>
	</P>
	<P>O charge = -1.0484<BR>
	H charge = 0.5242 <BR>
	</P>
	<P>r0 of OH bond = 0.9572<BR>
	theta of HOH angle = 104.52 <BR>
	</P>
	<P>OM distance = 0.1250 <BR>
	</P>
	<P>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>
	<HR>

	<A NAME = "4_9"></A><H4>4.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 with long-range
	Coulombics (Ewald or PPPM in LAMMPS).
	</P>
	<P>O mass = 15.9994<BR>
	H mass = 1.008 <BR>
	</P>
	<P>O charge = -0.820<BR>
	H charge = 0.410 <BR>
	</P>
	<P>LJ epsilon of OO = 0.1553<BR>
	LJ sigma of OO = 3.166<BR>
	LJ epsilon, sigma of OH, HH = 0.0 <BR>
	</P>
	<P>r0 of OH bond = 1.0<BR>
	theta of HOH angle = 109.47 <BR>
	</P>
	<HR>

	<A NAME = "4_10"></A><H4>4.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">this section</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><A HREF = "Section_start.html#2_4">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, or Fortran (or any other
	language that supports a vanilla C-like interface, e.g. a scripting
	language). 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. Library.cpp and library.h contain such a C interface with the
	functions:
	</P>
	<PRE>void lammps_open(int, char , MPI_Comm, void );
	void lammps_close(void *);
	void lammps_file(void , char );
	char lammps_command(doivd , char *);
	</PRE>
	<P>The functions contain C++ code you could write in a C++ application
	that was invoking LAMMPS directly. Note that LAMMPS classes are
	defined wihin a LAMMPS namespace (LAMMPS_NS) if you use them
	from another C++ application.
	</P>
	<P>Two of the routines in library.cpp are of particular note. The
	lammps_open() function initiates LAMMPS and takes an MPI communicator
	as an argument. It returns a pointer to a LAMMPS "object". As with
	C++, the lammps_open() function can be called mutliple 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.
	</P>
	<P>Library.cpp contains a lammps_command() function to which the caller
	passes a single LAMMPS command (a string). Thus the calling code can
	read or generate a series of LAMMPS commands (e.g. an input script)
	one line at a time and pass it thru the library interface to setup a
	problem and then run it.
	</P>
	<P>A few other sample functions are included in library.cpp, but the key
	idea is that you can write any functions you wish to define an
	interface for how your code talks to LAMMPS and add them to
	library.cpp and library.h. The routines you add can access any LAMMPS
	data. The examples/couple directory has example C++ and C codes which
	show how a stand-alone code can link 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 = "4_11"></A><H4>4.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 geoemtry
	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's 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://164.107.79.177/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.cs.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 = "4_12"></A><H4>4.12 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). If the box size
	is xprd by yprd by zprd then the 3 mutually orthogonal edge vectors of
	an orthogonal simulation box are a = (xprd,0,0), b = (0,yprd,0), and c
	= (0,0,zprd).
	</P>
	<P>LAMMPS also allows non-orthogonal simulation boxes (triclinic
	symmetry) to be defined with 3 additional "tilt" parameters which
	change the edge vectors of the simulation box to be a = (xprd,0,0), b
	= (xy,yprd,0), and c = (xz,yz,zprd). The xy, xz, and yz parameters
	can be positive or negative. The simulation box must be periodic in
	both dimensions associated with a tilt factor. For example, if xz !=
	0.0, then the x and z dimensions must be periodic.
	</P>
	<P>To avoid extremely tilted boxes (which would be computationally
	inefficient), 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). 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 all equivalent.
	</P>
	<P>You tell LAMMPS to use a non-orthogonal box when the simulation box is
	defined. 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 non-orthogonal domain 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 non-orthogonal domain 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
	non-orthogonal box will be enabled for the restarted simulation.
	</P>
	<P>Note that you can define a non-orthogonal box with all 3 tilt factors
	= 0.0, so that it is initially orthogonal. This is necessary if the
	box will become non-orthogonal. Alternatively, you can use the
	<A HREF = "change_box.html">change_box</A> command to convert a simulation box from
	orthogonal to non-orthogonal and vice versa.
	</P>
	<P>One use of non-orthogonal 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
	non-orthogonal simulation box via the <A HREF = "create_atoms.html">create_atoms</A>
	command. Note that while the box edge vectors a,b,c cannot be
	arbitrary vectors (e.g. a must be aligned with the x axis), it is
	possible to rotate any crystal's basis vectors so that they meet these
	restrictions.
	</P>
	<P>A second use of non-orthogonal boxes 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, and yz tilt factors as a simulation runs.
	</P>
	<P>Another use of non-orthogonal boxes is to perform non-equilibrium MD
	(NEMD) simulations, as discussed in the next section.
	</P>
	<HR>

	<A NAME = "4_13"></A><H4>4.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 preceeding Howto section).
	</P>
	<P>A shear strain can be applied to the simuaation 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 re-shaped 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 = "4_14"></A><H4>4.14 Aspherical particles
	</H4>
	<P>LAMMPS supports ellipsoidal particles via the <A HREF = "atom_style.html">atom_style
	ellipsoid</A> and <A HREF = "shape.html">shape</A> commands. The
	latter command defines the 3 axes (diameters) of a general ellipsoid.
	The <A HREF = "pair_gayberne.html">pair_style gayberne</A> command can be used to
	define a Gay-Berne (GB) potential for how ellipsoidal particles
	interact with each other and with spherical particles. The GB
	potential is like a Lennard-Jones (LJ) potential, generalized for
	orientiation-dependent interactions.
	</P>
	<P>The orientation of ellipsoidal particles is stored as a quaternion.
	See the <A HREF = "set.html">set</A> command for a brief explanation of quaternions
	and how the orientation of such particles can be initialized. The
	data file read by the <A HREF = "read_data.html">read_data</A> command contains
	quaternions for each atom in the Atoms section if <A HREF = "atom_style.html">atom_style
	ellipsoid</A> is being used. The <A HREF = "compute_temp_asphere.html">compute
	temp/asphere</A> command can be used to
	calculate the temperature of a group of ellipsoidal particles, taking
	account of rotational degrees of freedom. The motion of the particles
	can be integrated via the <A HREF = "fix_nve_asphere.html">fix nve/asphere</A>, <A HREF = "fix_nvt_asphere.html">fix
	nvt/asphere</A>, or <A HREF = "fix_npt_asphere.html">fix
	npt/asphere</A> commands. All of these commands are
	part of the ASPHERE package in LAMMPS.
	</P>
	<P>Computationally, the cost for two ellipsoidal particles to interact is
	30 times (or more) expensive than for 2 spherical LJ particles. Thus
	if you are modeling a system with many spherical particles (e.g. as
	the solvent), then you should insure sphere-sphere interactions are
	computed with a cheaper potential than GB. This can be done by
	setting the particle's 3 shape parameters to all be equal (a sphere).
	Additionally, the corresponding GB potential coefficients can be set
	so the GB potential will treat the pair of particles as LJ spheres.
	Details are given in the doc page for the <A HREF = "pair_gayberne.html">pair_style
	gayberne</A>. Alternatively, the <A HREF = "pair_hybrid.html">pair_style
	hybrid</A> potential can be used, with the sphere-sphere
	interactions computed by another pair potential, such as <A HREF = "pair_lj.html">pair_style
	lj/cut</A>.
	</P>
	<HR>

	<A NAME = "4_15"></A><H4>4.15 Output from LAMMPS
	</H4>
	<P>Aside from <A HREF = "restart.html">restart files</A>, there are two basic kinds of
	LAMMPS output. The first is <A HREF = "thermo_style.html">thermodynamic output</A>,
	which is a list of quantities printed every few timesteps to the
	screen and logfile. The second is <A HREF = "dump.html">dump files</A>, which
	contain snapshots of atoms and various per-atom values and are written
	at a specified frequency. A simulation prints one set of
	thermodynamic output; it may generate zero, or one, or multiple dump
	files. LAMMPS gives you a variety of ways to determine what
	quantities are computed and printed when thermodynamic info or dump
	files are output. There are also three fixes which can do their own
	output of user-defined quantities: <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>. These are described below.
	</P>
	<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 = "themo_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) which include time-averaged versions of
	temperature, pressure, and a few other variables (tave, pave, etc).
	Three addtional 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(s) to be output. Each of
	these are described in turn.
	</P>
	<P>In LAMMPS, a <A HREF = "compute.html">compute</A> comes in two flavors: ones that
	compute one or more global values (e.g. temperature, kinetic energy
	tensor) and ones that compute one or more per-atom values. There is a
	<A HREF = "compute_sum.html">compute sum</A> command which sums per-atom quantities
	into a global scalar or vector.
	</P>
	<P>Only global quantities from a compute can be used for thermodynamic
	output. The user-defined ID of the compute is used along with an
	optional subscript as part of the <A HREF = "thermo_style.html">thermo_style</A>
	command. E.g. c_myTemp outputs the single scalar value generated by
	the compute; c_myTemp[2] would output the 2nd vector value.
	</P>
	<P><A HREF = "fix.html">Fixes</A> can also generate global scalar or vector values
	which can be output with thermodynamic output, e.g. the energy of an
	indenter's interaction with the simulation atoms. These values are
	accessed via the same format as a compute's values, as f_ID or
	f_ID[N]. See the doc pages for individual fix commands to see which
	ones generate global values that can be output with thermodynamic
	info. The <A HREF = "fix_ave_time.html">fix ave/time</A> command generates
	time-averaged global quantities which can be accessed for
	thermodynamic output.
	</P>
	<P>Input script variables of various kinds are defined by the
	<A HREF = "variable.html">variable</A> command. All kinds except the atom-style
	variable can be used for thermodynamic output. A variable with name
	"abc" is referenced in a thermo_style command as v_abc.
	</P>
	<P>The variable formula defined in the input script can contain math
	functions (add, exp, etc), atom values (x[N], fx[N]), groups
	quantities (mass(), vcm(), etc), references to thermodynamic
	quantities (e.g. temp, volume, etc), or references to other variables
	or <A HREF = "compute.html">computes</A> or <A HREF = "fix.html">fixes</A>. Thus a variable is
	the most general way to define some quantity you want calculated and
	output with thermodynamic info.
	</P>
	<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). There is also a <A HREF = "dump.html">dump
	custom</A> format where you specify what values are output with
	each atom. Pre-defined keywords can be specified (e.g. tag, type, x,
	etc). 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> provides the
	values to be output.
	</P>
	<P><A HREF = "compute.html">Computes</A> that generate per-atom values can be accessed
	by the dump custom command. These are computes that have the word
	"atom" in their style name, e.g. ke/atom, stress/atom, etc. The
	values are accessed as c_myKE for a scalar per-atom quantity or as
	c_myStress[2] for a component of a vector per-atom quantity. The
	<A HREF = "compute_variable_atom.html">compute variable/atom</A> command takes a
	user-defined atom-style <A HREF = "variable.html">variable</A> as input and
	calculates its value for each atom. Since this compute can be
	accessed by the dump custom command, this is a general way to define
	some quantity you want calculated and output in a dump file.
	</P>
	<P><A HREF = "fix.html">Fixes</A> can also generate per-atom values to output to dump
	files. For example, the <A HREF = "fix_ave_atom.html">fix ave/atom</A> command
	calculates time-averages of compute quantities. As indicated in the
	preceeding paragraph, a <A HREF = "compute.html">compute quantity</A> can be a
	calculated value such as <A HREF = "compute_epair_atom.html">energy</A> or
	<A HREF = "compute_stress_atom.html">stress</A> or it can be a value calculated by
	an atom-style <A HREF = "variable.html">variable</A>, or it can be an <A HREF = "compute_attribute_atom.html">atom
	attribute</A> such as velocity or force.
	These per-atom fix values are accessed by the <A HREF = "dump.html">dump custom</A>
	command as f_myKE for a scalar per-atom quantity or as f_myStress[2]
	for a component of a vector per-atom quantity.
	</P>
	<P>Three other fixes are of particular note for output: <A HREF = "fix_ave_time.html">fix
	ave/time</A>, <A HREF = "fix_ave_spatial.html">fix ave/spatial</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 time-averaging of
	global quantities like temperature or pressure. The global quantities
	are calculated by a <A HREF = "compute.html">compute</A> or a <A HREF = "fix.html">fix</A>. The
	compute or fix must generate global scalar or vector quantities. Note
	that this includes the "compute sum" command which computes global
	values by summing per-atom quantities. The time-averaged values
	generated by <A HREF = "fix_ave_time.html">fix ave/time</A> can be written directly
	to a file and/or accessed by any output command that uses fixes as a
	source of input, e.g. the <A HREF = "thermo_style.html">thermo_style custom</A>
	command. Fix ave/time options allow for running cummulative averages
	or moving time-windowed averages to be output.
	</P>
	<P>The <A HREF = "fix_ave_spatial.html">fix ave/spatial</A> command enables
	spatial-averaging of per-atom quantities like per-atom energy or
	stress. The per-atom quantities can be atom density (mass or number)
	or be calculated by a <A HREF = "compute.html">compute</A> or a <A HREF = "fix.html">fix</A>. The
	compute or fix must generate per-atom scalar or vector quantities.
	Note that if you use the <A HREF = "fix_ave_atom.html">fix ave/atom</A> command with
	fix ave/spatial, it means you are effectively calculating a time
	average of a spatial average of a time-averaged per-atom quantity.
	The time-averaged values generated by <A HREF = "fix_ave_spatial.html">fix
	ave/spatial</A> can be written directly to a file
	and/or accessed by any output command that uses fixes as a source of
	input, e.g. the <A HREF = "thermo_style.html">thermo_style custom</A> command. Fix
	ave/spatial options allow for running cummulative averages or moving
	time-windowed averages to be output.
	</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 periodically during a running
	simulation. Since the line can contain one or more
	<A HREF = "variable.html">variable</A> quantities, this command is a means to output
	desired calculated quantities that are not part of thermodynamic or
	dump file output.
	</P>
	<HR>

	<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 = "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 = "Jorgensen"></A>

	<P><B>(Jorgensen)</B> Jorgensen, Chandrasekhar, Madura, Impey, Klein, J Chem
	Phys, 79, 926 (1983).
	</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