forked from lijiext/lammps
2661 lines
129 KiB
HTML
2661 lines
129 KiB
HTML
<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>
|
|
6.23 <A HREF = "#howto_23">Using chunks to calculate system properties</A><BR>
|
|
6.24 <A HREF = "#howto_24">Setting parameters for the kspace_style pppm/disp command</A><BR>
|
|
6.25 <A HREF = "#howto_25">Adiabatic core/shell model</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 directly (e.g. 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 for velocity biases indirectly, via an optional extra
|
|
argument, 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 which has two effects. First, the current calculated
|
|
temperature, which is compared to the requested target temperature, is
|
|
caluclated with the velocity bias removed. Second, the thermostat
|
|
adjusts only the thermal temperature component of the particle's
|
|
velocities, which are the velocities with the bias removed. The
|
|
removed bias is then added back to the adjusted 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>. Of you could thermostat only
|
|
the thermal temperature of a streaming flow of particles without
|
|
affecting the streaming velocity, by using <A HREF = "compute_temp_profile.html">compute
|
|
temp/profile</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_set_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
|
|
"set_variable" function can set an existing string-style variable to a
|
|
new value, so that subsequent LAMMPS commands can access the 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>
|
|
|
|
<A NAME = "howto_23"></A><H4>6.23 Using chunks to calculate system properties
|
|
</H4>
|
|
<P>In LAMMS, "chunks" are collections of atoms, as defined by the
|
|
<A HREF = "compute_chunk_atom.html">compute chunk/atom</A> command, which assigns
|
|
each atom to a chunk ID (or to no chunk at all). The number of chunks
|
|
and the assignment of chunk IDs to atoms can be static or change over
|
|
time. Examples of "chunks" are molecules or spatial bins or atoms
|
|
with similar values (e.g. coordination number or potential energy).
|
|
</P>
|
|
<P>The per-atom chunk IDs can be used as input to two other kinds of
|
|
commands, to calculate various properties of a system:
|
|
</P>
|
|
<UL><LI><A HREF = "fix_ave_chunk.html">fix ave/chunk</A>
|
|
<LI>any of the <A HREF = "compute.html">compute */chunk</A> commands
|
|
</UL>
|
|
<P>Here, each of the 3 kinds of chunk-related commands is briefly
|
|
overviewed. Then some examples are given of how to compute different
|
|
properties with chunk commands.
|
|
</P>
|
|
<H5>Compute chunk/atom command:
|
|
</H5>
|
|
<P>This compute can assign atoms to chunks of various styles. Only atoms
|
|
in the specified group and optional specified region are assigned to a
|
|
chunk. Here are some possible chunk definitions:
|
|
</P>
|
|
<DIV ALIGN=center><TABLE BORDER=1 >
|
|
<TR><TD >atoms in same molecule </TD><TD > chunk ID = molecule ID </TD></TR>
|
|
<TR><TD >atoms of same atom type </TD><TD > chunk ID = atom type </TD></TR>
|
|
<TR><TD >all atoms with same atom property (charge, radius, etc) </TD><TD > chunk ID = output of compute property/atom </TD></TR>
|
|
<TR><TD >atoms in same cluster </TD><TD > chunk ID = output of <A HREF = "compute_cluster_atom.html">compute cluster/atom</A> command </TD></TR>
|
|
<TR><TD >atoms in same spatial bin </TD><TD > chunk ID = bin ID </TD></TR>
|
|
<TR><TD >atoms in same rigid body </TD><TD > chunk ID = molecule ID used to define rigid bodies </TD></TR>
|
|
<TR><TD >atoms with similar potential energy </TD><TD > chunk ID = output of <A HREF = "compute_pe_atom.html">compute pe/atom</A> </TD></TR>
|
|
<TR><TD >atoms with same local defect structure </TD><TD > chunk ID = output of <A HREF = "compute_centro_atom.html">compute centro/atom</A> or <A HREF = "compute_coord_atom.html">compute coord/atom</A> command
|
|
</TD></TR></TABLE></DIV>
|
|
|
|
<P>Note that chunk IDs are integer values, so for atom properties or
|
|
computes that produce a floating point value, they will be truncated
|
|
to an integer. You could also use the compute in a variable that
|
|
scales the floating point value to spread it across multiple intergers.
|
|
</P>
|
|
<P>Spatial bins can be of various kinds, e.g. 1d bins = slabs, 2d bins =
|
|
pencils, 3d bins = boxes, spherical bins, cylindrical bins.
|
|
</P>
|
|
<P>This compute also calculates the number of chunks <I>Nchunk</I>, which is
|
|
used by other commands to tally per-chunk data. <I>Nchunk</I> can be a
|
|
static value or change over time (e.g. the number of clusters). The
|
|
chunk ID for an individual atom can also be static (e.g. a molecule
|
|
ID), or dynamic (e.g. what spatial bin an atom is in as it moves).
|
|
</P>
|
|
<P>Note that this compute allows the per-atom output of other
|
|
<A HREF = "compute.html">computes</A>, <A HREF = "fix.html">fixes</A>, and
|
|
<A HREF = "variable.html">variables</A> to be used to define chunk IDs for each
|
|
atom. This means you can write your own compute or fix to output a
|
|
per-atom quantity to use as chunk ID. See
|
|
<A HREF = "Section_modify.html">Section_modify</A> of the documentation for how to
|
|
do this. You can also define a <A HREF = "variable.html">per-atom variable</A> in
|
|
the input script that uses a formula to generate a chunk ID for each
|
|
atom.
|
|
</P>
|
|
<H5>Fix ave/chunk command:
|
|
</H5>
|
|
<P>This fix takes the ID of a <A HREF = "compute_chunk_atom.html">compute
|
|
chunk/atom</A> command as input. For each chunk,
|
|
it then sums one or more specified per-atom values over the atoms in
|
|
each chunk. The per-atom values can be any atom property, such as
|
|
velocity, force, charge, potential energy, kinetic energy, stress,
|
|
etc. Additional keywords are defined for per-chunk properties like
|
|
density and temperature. More generally any per-atom value generated
|
|
by other <A HREF = "compute.html">computes</A>, <A HREF = "fix.html">fixes</A>, and <A HREF = "variable.html">per-atom
|
|
variables</A>, can be summed over atoms in each chunk.
|
|
</P>
|
|
<P>Similar to other averaging fixes, this fix allows the summed per-chunk
|
|
values to be time-averaged in various ways, and output to a file. The
|
|
fix produces a global array as output with one row of values per
|
|
chunk.
|
|
</P>
|
|
<H5>Compute */chunk commands:
|
|
</H5>
|
|
<P>Currently the following computes operate on chunks of atoms to produce
|
|
per-chunk values.
|
|
</P>
|
|
<UL><LI><A HREF = "compute_com_chunk.html">compute com/chunk</A>
|
|
<LI><A HREF = "compute_gyration_chunk.html">compute gyration/chunk</A>
|
|
<LI><A HREF = "compute_inertia_chunk.html">compute inertia/chunk</A>
|
|
<LI><A HREF = "compute_msd_chunk.html">compute msd/chunk</A>
|
|
<LI><A HREF = "compute_property_chunk.html">compute property/chunk</A>
|
|
<LI><A HREF = "compute_temp_chunk.html">compute temp/chunk</A>
|
|
<LI><A HREF = "compute_vcm_chunk.html">compute torque/chunk</A>
|
|
<LI><A HREF = "compute_vcm_chunk.html">compute vcm/chunk</A>
|
|
</UL>
|
|
<P>They each take the ID of a <A HREF = "compute_chunk_atom.html">compute
|
|
chunk/atom</A> command as input. As their names
|
|
indicate, they calculate the center-of-mass, radius of gyration,
|
|
moments of inertia, mean-squared displacement, temperature, torque,
|
|
and velocity of center-of-mass for each chunk of atoms. The <A HREF = "compute_property_chunk.html">compute
|
|
property/chunk</A> command can tally the
|
|
count of atoms in each chunk and extract other per-chunk properties.
|
|
</P>
|
|
<P>The reason these various calculations are not part of the <A HREF = "fix_ave_chunk.html">fix
|
|
ave/chunk command</A>, is that each requires a more
|
|
complicated operation than simply summing and averaging over per-atom
|
|
values in each chunk. For example, many of them require calculation
|
|
of a center of mass, which requires summing mass*position over the
|
|
atoms and then dividing by summed mass.
|
|
</P>
|
|
<P>All of these computes produce a global vector or global array as
|
|
output, wih one or more values per chunk. They can be used
|
|
in various ways:
|
|
</P>
|
|
<UL><LI>As input to the <A HREF = "fix_ave_time.html">fix ave/time</A> command, which can
|
|
write the values to a file and optionally time average them.
|
|
|
|
<LI>As input to the <A HREF = "fix_ave_histo.html">fix ave/histo</A> command to
|
|
histogram values across chunks. E.g. a histogram of cluster sizes or
|
|
molecule diffusion rates.
|
|
|
|
<LI>As input to special functions of <A HREF = "variable.html">equal-style
|
|
variables</A>, like sum() and max(). E.g. to find the
|
|
largest cluster or fastest diffusing molecule.
|
|
</UL>
|
|
<H5>Example calculations with chunks
|
|
</H5>
|
|
<P>Here are eaxmples using chunk commands to calculate various
|
|
properties:
|
|
</P>
|
|
<P>(1) Average velocity in each of 1000 2d spatial bins:
|
|
</P>
|
|
<PRE>compute cc1 all chunk/atom bin/2d x 0.0 0.1 y lower 0.01 units reduced
|
|
fix 1 all ave/chunk 100 10 1000 cc1 vx vy file tmp.out
|
|
</PRE>
|
|
<P>(2) Temperature in each spatial bin, after subtracting a flow
|
|
velocity:
|
|
</P>
|
|
<PRE>compute cc1 all chunk/atom bin/2d x 0.0 0.1 y lower 0.1 units reduced
|
|
compute vbias all temp/profile 1 0 0 y 10
|
|
fix 1 all ave/chunk 100 10 1000 cc1 temp bias vbias file tmp.out
|
|
</PRE>
|
|
<P>(3) Center of mass of each molecule:
|
|
</P>
|
|
<PRE>compute cc1 all chunk/atom molecule
|
|
compute myChunk all com/chunk cc1
|
|
fix 1 all ave/time 100 1 100 c_myChunk file tmp.out mode vector
|
|
</PRE>
|
|
<P>(4) Total force on each molecule and ave/max across all molecules:
|
|
</P>
|
|
<PRE>compute cc1 all chunk/atom molecule
|
|
fix 1 all ave/chunk 1000 1 1000 cc1 fx fy fz file tmp.out
|
|
variable xave equal ave(f_1<B>2</B>)
|
|
variable xmax equal max(f_1<B>2</B>)
|
|
thermo 1000
|
|
thermo_style custom step temp v_xave v_xmax
|
|
</PRE>
|
|
<P>(5) Histogram of cluster sizes:
|
|
</P>
|
|
<PRE>compute cluster all cluster/atom 1.0
|
|
compute cc1 all chunk/atom c_cluster compress yes
|
|
compute size all property/chunk cc1 count
|
|
fix 1 all ave/histo 100 1 100 0 20 20 c_size mode vector ave running beyond ignore file tmp.histo
|
|
</PRE>
|
|
<HR>
|
|
|
|
<A NAME = "howto_24"></A><H4>6.24 Setting parameters for the <A HREF = "kspace_style.html">kspace_style pppm/disp</A> command
|
|
</H4>
|
|
<P>The PPPM method computes interactions by splitting the pair potential
|
|
into two parts, one of which is computed in a normal pairwise fashion,
|
|
the so-called real-space part, and one of which is computed using the
|
|
Fourier transform, the so called reciprocal-space or kspace part. For
|
|
both parts, the potential is not computed exactly but is approximated.
|
|
Thus, there is an error in both parts of the computation, the
|
|
real-space and the kspace error. The just mentioned facts are true
|
|
both for the PPPM for Coulomb as well as dispersion interactions. The
|
|
deciding difference - and also the reason why the parameters for
|
|
pppm/disp have to be selected with more care - is the impact of the
|
|
errors on the results: The kspace error of the PPPM for Coulomb and
|
|
dispersion interaction and the real-space error of the PPPM for
|
|
Coulomb interaction have the character of noise. In contrast, the
|
|
real-space error of the PPPM for dispersion has a clear physical
|
|
interpretation: the underprediction of cohesion. As a consequence, the
|
|
real-space error has a much stronger effect than the kspace error on
|
|
simulation results for pppm/disp. Parameters must thus be chosen in a
|
|
way that this error is much smaller than the kspace error.
|
|
</P>
|
|
<P>When using pppm/disp and not making any specifications on the PPPM
|
|
parameters via the kspace modify command, parameters will be tuned
|
|
such that the real-space error and the kspace error are equal. This
|
|
will result in simulations that are either inaccurate or slow, both of
|
|
which is not desirable. For selecting parameters for the pppm/disp
|
|
that provide fast and accurate simulations, there are two approaches,
|
|
which both have their up- and downsides.
|
|
</P>
|
|
<P>The first approach is to set desired real-space an kspace accuracies
|
|
via the <I>kspace_modify force/disp/real</I> and <I>kspace_modify
|
|
force/disp/kspace</I> commands. Note that the accuracies have to be
|
|
specified in force units and are thus dependend on the chosen unit
|
|
settings. For real units, 0.0001 and 0.002 seem to provide reasonable
|
|
accurate and efficient computations for the real-space and kspace
|
|
accuracies. 0.002 and 0.05 work well for most systems using lj
|
|
units. PPPM parameters will be generated based on the desired
|
|
accuracies. The upside of this approach is that it usually provides a
|
|
good set of parameters and will work for both the <I>kspace_modify diff
|
|
ad</I> and <I>kspace_modify diff ik</I> options. The downside of the method
|
|
is that setting the PPPM parameters will take some time during the
|
|
initialization of the simulation.
|
|
</P>
|
|
<P>The second approach is to set the parameters for the pppm/disp
|
|
explicitly using the <I>kspace_modify mesh/disp</I>, <I>kspace_modify
|
|
order/disp</I>, and <I>kspace_modify gewald/disp</I> commands. This approach
|
|
requires a more experienced user who understands well the impact of
|
|
the choice of parameters on the simulation accuracy and
|
|
performance. This approach provides a fast initialization of the
|
|
simulation. However, it is sensitive to errors: A combination of
|
|
parameters that will perform well for one system might result in
|
|
far-from-optimal conditions for other simulations. For example,
|
|
parametes that provide accurate and fast computations for
|
|
all-atomistic force fields can provide insufficient accuracy or
|
|
united-atomistic force fields (which is related to that the latter
|
|
typically have larger dispersion coefficients).
|
|
</P>
|
|
<P>To avoid inaccurate or inefficient simulations, the pppm/disp stops
|
|
simulations with an error message if no action is taken to control the
|
|
PPPM parameters. If the automatic parameter generation is desired and
|
|
real-space and kspace accuracies are desired to be equal, this error
|
|
message can be suppressed using the <I>kspace_modify disp/auto yes</I>
|
|
command.
|
|
</P>
|
|
<P>A reasonable approach that combines the upsides of both methods is to
|
|
make the first run using the <I>kspace_modify force/disp/real</I> and
|
|
<I>kspace_modify force/disp/kspace</I> commands, write down the PPPM
|
|
parameters from the outut, and specify these parameters using the
|
|
second approach in subsequent runs (which have the same composition,
|
|
force field, and approximately the same volume).
|
|
</P>
|
|
<P>Concerning the performance of the pppm/disp there are two more things
|
|
to consider. The first is that when using the pppm/disp, the cutoff
|
|
parameter does no longer affect the accuracy of the simulation
|
|
(subject to that gewald/disp is adjusted when changing the cutoff).
|
|
The performance can thus be increased by examining different values
|
|
for the cutoff parameter. A lower bound for the cutoff is only set by
|
|
the truncation error of the repulsive term of pair potentials.
|
|
</P>
|
|
<P>The second is that the mixing rule of the pair style has an impact on
|
|
the computation time when using the pppm/disp. Fastest computations
|
|
are achieved when using the geometric mixing rule. Using the
|
|
arithmetic mixing rule substantially increases the computational cost.
|
|
The computational overhead can be reduced using the <I>kspace_modify
|
|
mix/disp geom</I> and <I>kspace_modify splittol</I> commands. The first
|
|
command simply enforces geometric mixing of the dispersion
|
|
coeffiecients in kspace computations. This introduces some error in
|
|
the computations but will also significantly speed-up the
|
|
simulations. The second keyword sets the accuracy with which the
|
|
dispersion coefficients are approximated using a matrix factorization
|
|
approach. This may result in better accuracy then using the first
|
|
command, but will usually also not provide an equally good increase of
|
|
efficiency.
|
|
</P>
|
|
<P>Finally, pppm/disp can also be used when no mixing rules apply.
|
|
This can be achieved using the <I>kspace_modify mix/disp none</I> command.
|
|
Note that the code does not check automatically whether any mixing
|
|
rule is fulfilled. If mixing rules do not apply, the user will have
|
|
to specify this command explicitly.
|
|
</P>
|
|
<HR>
|
|
|
|
<A NAME = "howto_25"></A><H4>6.25 Adiabatic core/shell model
|
|
</H4>
|
|
<P>The adiabatic core-shell model by <A HREF = "#MitchellFinchham">Mitchell and
|
|
Finchham</A> is a simple method for adding
|
|
polarizability to a system. In order to mimic the electron shell of
|
|
an ion, a ghost atom is attached to it. This way the ions are split
|
|
into a core and a shell where the latter is meant to react to the
|
|
electrostatic environment inducing polarizability.
|
|
</P>
|
|
<P>Technically, shells are attached to the cores by a spring force f =
|
|
k*r where k is a parametrized spring constant and r is the distance
|
|
between the core and the shell. The charges of the core and the shell
|
|
add up to the ion charge, thus q(ion) = q(core) + q(shell). In a
|
|
similar fashion the mass of the ion is distributed on the core and the
|
|
shell with the core having the larger mass.
|
|
</P>
|
|
<P>To run this model in LAMMPS, <A HREF = "atom_style.html">atom_style</A> <I>full</I> can
|
|
be used since atom charge and bonds are needed. Each kind of
|
|
core/shell pair requires two atom types and a bond type. The core and
|
|
shell of a core/shell pair should be bonded to each other with a
|
|
harmonic bond that provides the spring force. For example, a data file
|
|
for NaCl, as found in examples/coreshell, has this format:
|
|
</P>
|
|
<PRE>432 atoms # core and shell atoms
|
|
216 bonds # number of core/shell springs
|
|
</PRE>
|
|
<PRE>4 atom types # 2 cores and 2 shells for Na and Cl
|
|
2 bond types
|
|
</PRE>
|
|
<PRE>0.0 24.09597 xlo xhi
|
|
0.0 24.09597 ylo yhi
|
|
0.0 24.09597 zlo zhi
|
|
</PRE>
|
|
<PRE>Masses # core/shell mass ratio = 0.1
|
|
</PRE>
|
|
<PRE>1 20.690784 # Na core
|
|
2 31.90500 # Cl core
|
|
3 2.298976 # Na shell
|
|
4 3.54500 # Cl shell
|
|
</PRE>
|
|
<PRE>Atoms
|
|
</PRE>
|
|
<PRE>1 1 2 1.5005 0.00000000 0.00000000 0.00000000 # core of core/shell pair 1
|
|
2 1 4 -2.5005 0.00000000 0.00000000 0.00000000 # shell of core/shell pair 1
|
|
3 2 1 1.5056 4.01599500 4.01599500 4.01599500 # core of core/shell pair 2
|
|
4 2 3 -0.5056 4.01599500 4.01599500 4.01599500 # shell of core/shell pair 2
|
|
(...)
|
|
</PRE>
|
|
<PRE>Bonds # Bond topology for spring forces
|
|
</PRE>
|
|
<PRE>1 2 1 2 # spring for core/shell pair 1
|
|
2 2 3 4 # spring for core/shell pair 2
|
|
(...)
|
|
</PRE>
|
|
<P>Non-Coulombic (e.g. Lennard-Jones) pairwise interactions are only
|
|
defined between the shells. Coulombic interactions are defined
|
|
between all cores and shells. If desired, additional bonds can be
|
|
specified between cores.
|
|
</P>
|
|
<P>The <A HREF = "special_bonds.html">special_bonds</A> command should be used to
|
|
turn-off the Coulombic interaction within core/shell pairs, since that
|
|
interaction is set by the bond spring. This is done using the
|
|
<A HREF = "special_bonds.html">special_bonds</A> command with a 1-2 weight = 0.0,
|
|
which is the default value.
|
|
</P>
|
|
<P>Since the core/shell model permits distances of r = 0.0 between the
|
|
core and shell, a pair style with a "cs" suffix needs to be used to
|
|
implement a valid long-range Coulombic correction. Several such pair
|
|
styles are provided in the CORESHELL package. See <A HREF = "pair_cs.html">this doc
|
|
page</A> for details. All of the core/shell enabled pair
|
|
styles require the use of a long-range Coulombic solver, as specified
|
|
by the <A HREF = "kspace_style.html">kspace_style</A> command. Either the PPPM or
|
|
Ewald solvers can be used.
|
|
</P>
|
|
<P>For the NaCL example problem, these pair style and bond style settings
|
|
are used:
|
|
</P>
|
|
<PRE>pair_style born/coul/long/cs 20.0 20.0
|
|
pair_coeff * * 0.0 1.000 0.00 0.00 0.00
|
|
pair_coeff 3 3 487.0 0.23768 0.00 1.05 0.50 #Na-Na
|
|
pair_coeff 3 4 145134.0 0.23768 0.00 6.99 8.70 #Na-Cl
|
|
pair_coeff 4 4 405774.0 0.23768 0.00 72.40 145.40 #Cl-Cl
|
|
</PRE>
|
|
<PRE>bond_style harmonic
|
|
bond_coeff 1 63.014 0.0
|
|
bond_coeff 2 25.724 0.0
|
|
</PRE>
|
|
<P>When running dynamics with the adiabatic core/shell model, the
|
|
following issues should be considered. Since the relative motion of
|
|
the core and shell particles corresponds to the polarization, typical
|
|
thermostats can alter the polarization behaviour, meaining the shell
|
|
will not react freely to its electrostatic environment. Therefore
|
|
it's typically desirable to decouple the relative motion of the
|
|
core/shell pair, which is an imaginary degree of freedom, from the
|
|
real physical system. To do that, the <A HREF = "compute_temp_cs.html">compute
|
|
temp/cs</A> command can be used, in conjunction with
|
|
any of the thermostat fixes, such as <A HREF = "fix_nh.html">fix nvt</A> or <A HREF = "fix_langevin">fix
|
|
langevin</A>. This compute uses the center-of-mass velocity
|
|
of the core/shell pairs to calculate a temperature, and insures that
|
|
velocity is what is rescaled for thermostatting purposes. The
|
|
<A HREF = "compute_temp_cs.html">compute temp/cs</A> command requires input of two
|
|
groups, one for the core atoms, another for the shell atoms. These
|
|
can be defined using the <A HREF = "group.html">group <I>type</I></A> command. Note that
|
|
to perform thermostatting using this definition of temperature, the
|
|
<A HREF = "fix_modify.html">fix modify temp</A> command should be used to assign the
|
|
comptue to the thermostat fix. Likewise the <A HREF = "thermo_modify.html">thermo_modify
|
|
temp</A> command can be used to make this temperature
|
|
be output for the overall system.
|
|
</P>
|
|
<P>For the NaCl example, this can be done as follows:
|
|
</P>
|
|
<PRE>group cores type 1 2
|
|
group shells type 3 4
|
|
compute CSequ all temp/cs cores shells
|
|
fix thermoberendsen all temp/berendsen 1427 1427 0.4 # thermostat for the true physical system
|
|
fix thermostatequ all nve # integrator as needed for the berendsen thermostat
|
|
fix_modify thermoberendsen temp CSequ
|
|
thermo_modify temp CSequ # output of center-of-mass derived temperature
|
|
</PRE>
|
|
<P>When intializing the velocities of a system with core/shell pairs, it
|
|
is also desirable to not introduce energy into the relative motion of
|
|
the core/shell particles, but only assign a center-of-mass velocity to
|
|
the pairs. This can be done by using the <I>bias</I> keyword of the
|
|
<A HREF = "velocity.html">velocity create</A> command and assigning the <A HREF = "compute_temp_cs.html">compute
|
|
temp/cs</A> command to the <I>temp</I> keyword of the
|
|
<A HREF = "velocity.html">velocity</A> commmand, e.g.
|
|
</P>
|
|
<PRE>velocity all create 1427 134 bias yes temp CSequ
|
|
velocity all scale 1427 temp CSequ
|
|
</PRE>
|
|
<P>It is important to note that the polarizability of the core/shell
|
|
pairs is based on their relative motion. Therefore the choice of
|
|
spring force and mass ratio need to ensure much faster relative motion
|
|
of the 2 atoms within the core/shell pair than their center-of-mass
|
|
velocity. This allow the shells to effectively react instantaneously
|
|
to the electrostatic environment. This fast movement also limits the
|
|
timestep size that can be used.
|
|
</P>
|
|
<P>Additionally, the mass mismatch of the core and shell particles means
|
|
that only a small amount of energy is transfered to the decoupled
|
|
imaginary degrees of freedom. However, this transfer will typically
|
|
lead to a a small drift in total energy over time. This internal
|
|
energy can be monitored using the <A HREF = "compute_chunk_atom.html">compute
|
|
chunk/atom</A> and <A HREF = "compute_temp_chunk.html">compute
|
|
temp/chunk</A> commands. The internal kinetic
|
|
energies of each core/shell pair can then be summed using the sum()
|
|
special functino of the <A HREF = "variable.html">variable</A> command. Or they can
|
|
be time/averaged and output using the <A HREF = "fix_ave_time.html">fix ave/time</A>
|
|
command. To use these commands, each core/shell pair must be defined
|
|
as a "chunk". If each core/shell pair is defined as its own molecule,
|
|
the molecule ID can be used to define the chunks. If cores are bonded
|
|
to each other to form larger molecules, then another way to define the
|
|
chunks is to use the <A HREF = "fix_property_atom.html">fix property/atom</A> to
|
|
assign a core/shell ID to each atom via a special field in the data
|
|
file read by the <A HREF = "read_data.html">read_data</A> command. This field can
|
|
then be accessed by the <A HREF = "compute_property_atom.html">compute
|
|
property/atom</A> command, to use as input to
|
|
the <A HREF = "compute_chunk_atom.html">compute chunk/atom</A> command to define the
|
|
core/shell pairs as chunks.
|
|
</P>
|
|
<P>For example,
|
|
</P>
|
|
<PRE>fix csinfo all property/atom i_CSID # property/atom command
|
|
read_data NaCl_CS_x0.1_prop.data fix csinfo NULL CS-Info # atom property added in the data-file
|
|
compute prop all property/atom i_CSID
|
|
compute cs_chunk all chunk/atom c_prop
|
|
compute cstherm all temp/chunk cs_chunk temp internal com yes cdof 3.0 # note the chosen degrees of freedom for the core/shell pairs
|
|
fix ave_chunk all ave/time 10 1 10 c_cstherm file chunk.dump mode vector
|
|
</PRE>
|
|
<P>The additional section in the date file would be formatted like this:
|
|
</P>
|
|
<PRE>CS-Info # header of additional section
|
|
</PRE>
|
|
<PRE>1 1 # column 1 = atom ID, column 2 = core/shell ID
|
|
2 1
|
|
3 2
|
|
4 2
|
|
5 3
|
|
6 3
|
|
7 4
|
|
8 4
|
|
(...)
|
|
</PRE>
|
|
<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>
|
|
<A NAME = "MitchellFinchham"></A>
|
|
|
|
<P><B>(Mitchell and Finchham)</B> Mitchell, Finchham, J Phys Condensed Matter,
|
|
5, 1031-1038 (1993).
|
|
</P>
|
|
</HTML>
|